1 /**
2 * Copyright (c) 2011, The University of Southampton and the individual contributors.
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without modification,
6 * are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * * Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * * Neither the name of the University of Southampton nor the names of its
16 * contributors may be used to endorse or promote products derived from this
17 * software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
23 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
26 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30 package org.openimaj.image.processing.algorithm;
31
32 import org.openimaj.image.FImage;
33
34 import edu.emory.mathcs.jtransforms.fft.FloatFFT_2D;
35
36 /**
37 * Perform forward and inverse Fast Fourier Transforms on image data. This class
38 * computes the result of the transform in polar form (i.e. in terms of phase
39 * and magnitude). If you need the result of the transform in the raw complex
40 * form, then use the {@link FourierTransformComplex} class.
41 *
42 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
43 *
44 */
45 public class FourierTransform {
46 private FImage phase;
47 private FImage magnitude;
48 private boolean centre;
49
50 /**
51 * Construct Fourier Transform by performing a forward transform on the
52 * given image. If the centre option is set, the FFT will be re-ordered so
53 * that the DC component is in the centre.
54 *
55 * @param image
56 * the image to transform
57 * @param centre
58 * should the FFT be reordered so the centre is DC component
59 */
60 public FourierTransform(FImage image, boolean centre) {
61 this.centre = centre;
62
63 process(image);
64 }
65
66 /**
67 * Construct Fourier Transform object from the given magnitude and phase
68 * images in the frequency domain. The resultant object can then be used to
69 * construct the image using the {@link #inverse()} method.
70 *
71 * @param magnitude
72 * the magnitude image
73 * @param phase
74 * the phase image
75 * @param centre
76 * is the DC component in the image centre?
77 */
78 public FourierTransform(FImage magnitude, FImage phase, boolean centre) {
79 this.centre = centre;
80 this.magnitude = magnitude;
81 this.phase = phase;
82 }
83
84 /**
85 * Prepare data for a input to the FFT, padding if necessary.
86 *
87 * @param input
88 * input data
89 * @param rs
90 * desired number of rows
91 * @param cs
92 * desired number of columns
93 * @param centre
94 * if true, then the data will be prepared so that the DC
95 * component is centred.
96 * @return prepared data
97 */
98 public static float[][] prepareData(FImage input, int rs, int cs, boolean centre) {
99 return prepareData(input.pixels, rs, cs, centre);
100 }
101
102 /**
103 * Prepare data for a input to the FFT, padding if necessary. The data is
104 * prepared as a packed 1D array.
105 *
106 * @param input
107 * input data
108 * @param rs
109 * desired number of rows
110 * @param cs
111 * desired number of columns
112 * @param centre
113 * if true, then the data will be prepared so that the DC
114 * component is centred.
115 * @return prepared data
116 */
117 public static float[] prepareData1d(FImage input, int rs, int cs, boolean centre) {
118 return prepareData1d(input.pixels, rs, cs, centre);
119 }
120
121 /**
122 * Prepare data for a input to the FFT, padding if necessary.
123 *
124 * @param input
125 * input data
126 * @param rs
127 * desired number of rows
128 * @param cs
129 * desired number of columns
130 * @param centre
131 * if true, then the data will be prepared so that the DC
132 * component is centered.
133 * @return prepared data
134 */
135 public static float[][] prepareData(float[][] input, int rs, int cs, boolean centre) {
136 final float[][] prepared = new float[rs][cs * 2];
137
138 if (centre) {
139 for (int r = 0; r < Math.min(rs, input.length); r++) {
140 for (int c = 0; c < Math.min(cs, input[0].length); c++) {
141 prepared[r][c * 2] = input[r][c] * (1 - 2 * ((r + c) % 2));
142 }
143 }
144 } else {
145 for (int r = 0; r < Math.min(rs, input.length); r++) {
146 for (int c = 0; c < Math.min(cs, input[0].length); c++) {
147 prepared[r][c * 2] = input[r][c];
148 }
149 }
150 }
151
152 return prepared;
153 }
154
155 /**
156 * Prepare data for a input to the FFT, padding if necessary. The data is
157 * prepared as a packed 1D array.
158 *
159 * @param input
160 * input data
161 * @param rs
162 * desired number of rows
163 * @param cs
164 * desired number of columns
165 * @param centre
166 * if true, then the data will be prepared so that the DC
167 * component is centered.
168 * @return prepared data
169 */
170 public static float[] prepareData1d(float[][] input, int rs, int cs, boolean centre) {
171 final float[] prepared = new float[rs * cs * 2];
172
173 if (centre) {
174 for (int r = 0; r < Math.min(rs, input.length); r++) {
175 for (int c = 0; c < Math.min(cs, input[0].length); c++) {
176 prepared[r * 2 * cs + 2 * c] = input[r][c] * (1 - 2 * ((r + c) % 2));
177 }
178 }
179 } else {
180 for (int r = 0; r < Math.min(rs, input.length); r++) {
181 for (int c = 0; c < Math.min(cs, input[0].length); c++) {
182 prepared[r * 2 * cs + 2 * c] = input[r][c];
183 }
184 }
185 }
186
187 return prepared;
188 }
189
190 /**
191 * Extract the actual data from prepared data. The output image must have
192 * the same number of rows as the prepared data, and half the number of
193 * columns.
194 *
195 * @param prepared
196 * the prepared data
197 * @param output
198 * the output
199 * @param centre
200 * if true, then the data will be prepared so that the DC
201 * component is centered.
202 */
203 public static void unprepareData(float[][] prepared, FImage output, boolean centre) {
204 unprepareData(prepared, output.pixels, centre);
205 }
206
207 /**
208 * Extract the actual data from prepared data. The output image must have
209 * the same number of rows as the prepared data, and half the number of
210 * columns.
211 *
212 * @param prepared
213 * the prepared data
214 * @param output
215 * the output
216 * @param centre
217 * if true, then the data will be prepared so that the DC
218 * component is centered.
219 */
220 public static void unprepareData(float[] prepared, FImage output, boolean centre) {
221 unprepareData(prepared, output.pixels, centre);
222 }
223
224 /**
225 * Extract the actual data from prepared data. The output array must have
226 * the same number of rows as the prepared data, and half the number of
227 * columns.
228 *
229 * @param prepared
230 * the prepared data
231 * @param output
232 * the output
233 * @param centre
234 * if true, then the data will be prepared so that the DC
235 * component is centered.
236 */
237 public static void unprepareData(float[][] prepared, float[][] output, boolean centre) {
238 final int rs = output.length;
239 final int cs = output[0].length;
240
241 if (centre) {
242 for (int r = 0; r < rs; r++) {
243 for (int c = 0; c < cs; c++) {
244 output[r][c] = prepared[r][c * 2] * (1 - 2 * ((r + c) % 2));
245 }
246 }
247 } else {
248 for (int r = 0; r < rs; r++) {
249 for (int c = 0; c < cs; c++) {
250 output[r][c] = prepared[r][c * 2];
251 }
252 }
253 }
254 }
255
256 /**
257 * Extract the actual data from prepared data. The output array must have
258 * the same number of rows as the prepared data, and half the number of
259 * columns.
260 *
261 * @param prepared
262 * the prepared data
263 * @param output
264 * the output
265 * @param centre
266 * if true, then the data will be prepared so that the DC
267 * component is centered.
268 */
269 public static void unprepareData(float[] prepared, float[][] output, boolean centre) {
270 final int rs = output.length;
271 final int cs = output[0].length;
272
273 if (centre) {
274 for (int r = 0; r < rs; r++) {
275 for (int c = 0; c < cs; c++) {
276 output[r][c] = prepared[r * 2 * cs + 2 * c] * (1 - 2 * ((r + c) % 2));
277 }
278 }
279 } else {
280 for (int r = 0; r < rs; r++) {
281 for (int c = 0; c < cs; c++) {
282 output[r][c] = prepared[r * 2 * cs + 2 * c];
283 }
284 }
285 }
286 }
287
288 private void process(FImage image) {
289 final int cs = image.getCols();
290 final int rs = image.getRows();
291
292 phase = new FImage(cs, rs);
293 magnitude = new FImage(cs, rs);
294
295 final FloatFFT_2D fft = new FloatFFT_2D(rs, cs);
296 final float[][] prepared = prepareData(image.pixels, rs, cs, centre);
297
298 fft.complexForward(prepared);
299
300 for (int y = 0; y < rs; y++) {
301 for (int x = 0; x < cs; x++) {
302 final float re = prepared[y][x * 2];
303 final float im = prepared[y][1 + x * 2];
304
305 phase.pixels[y][x] = (float) Math.atan2(im, re);
306 magnitude.pixels[y][x] = (float) Math.sqrt(re * re + im * im);
307 }
308 }
309 }
310
311 /**
312 * Perform the inverse FFT using the underlying magnitude and phase images.
313 * The resultant reconstructed image may need normalisation.
314 *
315 * @return the reconstructed image
316 */
317 public FImage inverse() {
318 final int cs = magnitude.getCols();
319 final int rs = magnitude.getRows();
320
321 final FloatFFT_2D fft = new FloatFFT_2D(rs, cs);
322 final float[][] prepared = new float[rs][cs * 2];
323 for (int y = 0; y < rs; y++) {
324 for (int x = 0; x < cs; x++) {
325 final float p = phase.pixels[y][x];
326 final float m = magnitude.pixels[y][x];
327
328 final float re = (float) (m * Math.cos(p));
329 final float im = (float) (m * Math.sin(p));
330
331 prepared[y][x * 2] = re;
332 prepared[y][1 + x * 2] = im;
333 }
334 }
335
336 fft.complexInverse(prepared, true);
337
338 final FImage image = new FImage(cs, rs);
339 unprepareData(prepared, image, centre);
340
341 return image;
342 }
343
344 /**
345 * Get a log-normalised copy of the magnitude image suitable for displaying.
346 *
347 * @return a log-normalised copy of the magnitude image
348 */
349 public FImage getLogNormalisedMagnitude() {
350 final FImage im = magnitude.clone();
351
352 for (int y = 0; y < im.height; y++) {
353 for (int x = 0; x < im.width; x++) {
354 im.pixels[y][x] = (float) Math.log(im.pixels[y][x] + 1);
355 }
356 }
357
358 return im.normalise();
359 }
360
361 /**
362 * @return the phase image
363 */
364 public FImage getPhase() {
365 return phase;
366 }
367
368 /**
369 * @return the magnitude image
370 */
371 public FImage getMagnitude() {
372 return magnitude;
373 }
374
375 /**
376 * @return true if the DC component is in the centre; false otherwise
377 */
378 public boolean isCentre() {
379 return centre;
380 }
381 }