blob: 806777bb184163f69a5b77eef8f0c3c60fbf76dc [file] [log] [blame]
cristy3ed852e2009-09-05 21:47:34 +00001/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
6% FFFFF OOO U U RRRR IIIII EEEEE RRRR %
7% F O O U U R R I E R R %
8% FFF O O U U RRRR I EEE RRRR %
9% F O O U U R R I E R R %
10% F OOO UUU R R IIIII EEEEE R R %
11% %
12% %
13% MagickCore Discrete Fourier Transform Methods %
14% %
15% Software Design %
16% Sean Burke %
17% Fred Weinhaus %
18% John Cristy %
19% July 2009 %
20% %
21% %
cristy7e41fe82010-12-04 23:12:08 +000022% Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization %
cristy3ed852e2009-09-05 21:47:34 +000023% dedicated to making software imaging solutions freely available. %
24% %
25% You may not use this file except in compliance with the License. You may %
26% obtain a copy of the License at %
27% %
28% http://www.imagemagick.org/script/license.php %
29% %
30% Unless required by applicable law or agreed to in writing, software %
31% distributed under the License is distributed on an "AS IS" BASIS, %
32% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
33% See the License for the specific language governing permissions and %
34% limitations under the License. %
35% %
36%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37%
38%
39%
40*/
41
42/*
43 Include declarations.
44*/
45#include "magick/studio.h"
46#include "magick/cache.h"
47#include "magick/image.h"
48#include "magick/image-private.h"
49#include "magick/list.h"
50#include "magick/fourier.h"
51#include "magick/log.h"
52#include "magick/memory_.h"
53#include "magick/monitor.h"
54#include "magick/property.h"
cristye3038982010-05-17 02:20:52 +000055#include "magick/quantum-private.h"
cristy3ed852e2009-09-05 21:47:34 +000056#include "magick/thread-private.h"
57#if defined(MAGICKCORE_FFTW_DELEGATE)
cristy56ed31c2010-03-22 00:46:21 +000058#if defined(MAGICKCORE_HAVE_COMPLEX_H)
cristy3ed852e2009-09-05 21:47:34 +000059#include <complex.h>
cristy56ed31c2010-03-22 00:46:21 +000060#endif
cristy3ed852e2009-09-05 21:47:34 +000061#include <fftw3.h>
cristy47b022b2011-01-18 22:29:21 +000062#if !defined(MAGICKCORE_HAVE_CABS)
63#define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
64#endif
65#if !defined(MAGICKCORE_HAVE_CARG)
cristy4da3ba32011-02-07 16:58:11 +000066#define carg(z) (atan2(cimag(z),creal(z)))
cristy47b022b2011-01-18 22:29:21 +000067#endif
68#if !defined(MAGICKCORE_HAVE_CIMAG)
69#define cimag(z) (z[1])
70#endif
71#if !defined(MAGICKCORE_HAVE_CREAL)
72#define creal(z) (z[0])
73#endif
cristy3ed852e2009-09-05 21:47:34 +000074#endif
75
76/*
77 Typedef declarations.
78*/
79typedef struct _FourierInfo
80{
81 ChannelType
82 channel;
83
84 MagickBooleanType
85 modulus;
86
cristybb503372010-05-27 20:51:26 +000087 size_t
cristy3ed852e2009-09-05 21:47:34 +000088 width,
89 height;
90
cristybb503372010-05-27 20:51:26 +000091 ssize_t
cristy3ed852e2009-09-05 21:47:34 +000092 center;
93} FourierInfo;
94
95/*
96%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
97% %
98% %
99% %
100% F o r w a r d F o u r i e r T r a n s f o r m I m a g e %
101% %
102% %
103% %
104%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105%
106% ForwardFourierTransformImage() implements the discrete Fourier transform
107% (DFT) of the image either as a magnitude / phase or real / imaginary image
108% pair.
109%
110% The format of the ForwadFourierTransformImage method is:
111%
112% Image *ForwardFourierTransformImage(const Image *image,
113% const MagickBooleanType modulus,ExceptionInfo *exception)
114%
115% A description of each parameter follows:
116%
117% o image: the image.
118%
119% o modulus: if true, return as transform as a magnitude / phase pair
120% otherwise a real / imaginary image pair.
121%
122% o exception: return any errors or warnings in this structure.
123%
124*/
125
126#if defined(MAGICKCORE_FFTW_DELEGATE)
127
cristyc4ea4a42011-01-24 01:43:30 +0000128static MagickBooleanType RollFourier(const size_t width,const size_t height,
129 const ssize_t x_offset,const ssize_t y_offset,double *fourier)
cristy3ed852e2009-09-05 21:47:34 +0000130{
131 double
132 *roll;
133
cristyc4ea4a42011-01-24 01:43:30 +0000134 register ssize_t
135 i,
136 x;
137
cristybb503372010-05-27 20:51:26 +0000138 ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000139 u,
140 v,
141 y;
142
cristy3ed852e2009-09-05 21:47:34 +0000143 /*
cristy2da05352010-09-15 19:22:17 +0000144 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
cristy3ed852e2009-09-05 21:47:34 +0000145 */
cristyc4ea4a42011-01-24 01:43:30 +0000146 roll=(double *) AcquireQuantumMemory((size_t) height,width*sizeof(*roll));
cristy3ed852e2009-09-05 21:47:34 +0000147 if (roll == (double *) NULL)
148 return(MagickFalse);
149 i=0L;
cristybb503372010-05-27 20:51:26 +0000150 for (y=0L; y < (ssize_t) height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000151 {
152 if (y_offset < 0L)
cristybb503372010-05-27 20:51:26 +0000153 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
cristy3ed852e2009-09-05 21:47:34 +0000154 else
cristybb503372010-05-27 20:51:26 +0000155 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
cristy3ed852e2009-09-05 21:47:34 +0000156 y+y_offset;
cristybb503372010-05-27 20:51:26 +0000157 for (x=0L; x < (ssize_t) width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000158 {
159 if (x_offset < 0L)
cristybb503372010-05-27 20:51:26 +0000160 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
cristy3ed852e2009-09-05 21:47:34 +0000161 else
cristybb503372010-05-27 20:51:26 +0000162 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
cristy3ed852e2009-09-05 21:47:34 +0000163 x+x_offset;
164 roll[v*width+u]=fourier[i++];
cristyc4ea4a42011-01-24 01:43:30 +0000165 }
cristy3ed852e2009-09-05 21:47:34 +0000166 }
cristyc4ea4a42011-01-24 01:43:30 +0000167 (void) CopyMagickMemory(fourier,roll,height*width*sizeof(*roll));
cristy3ed852e2009-09-05 21:47:34 +0000168 roll=(double *) RelinquishMagickMemory(roll);
169 return(MagickTrue);
170}
171
cristybb503372010-05-27 20:51:26 +0000172static MagickBooleanType ForwardQuadrantSwap(const size_t width,
173 const size_t height,double *source,double *destination)
cristy3ed852e2009-09-05 21:47:34 +0000174{
cristy3ed852e2009-09-05 21:47:34 +0000175 MagickBooleanType
176 status;
177
cristybb503372010-05-27 20:51:26 +0000178 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000179 x;
180
cristyc4ea4a42011-01-24 01:43:30 +0000181 ssize_t
182 center,
183 y;
184
cristy3ed852e2009-09-05 21:47:34 +0000185 /*
186 Swap quadrants.
187 */
cristybb503372010-05-27 20:51:26 +0000188 center=(ssize_t) floor((double) width/2L)+1L;
189 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
cristy3ed852e2009-09-05 21:47:34 +0000190 if (status == MagickFalse)
191 return(MagickFalse);
cristybb503372010-05-27 20:51:26 +0000192 for (y=0L; y < (ssize_t) height; y++)
193 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000194 destination[width*y+x+width/2L]=source[center*y+x];
cristybb503372010-05-27 20:51:26 +0000195 for (y=1; y < (ssize_t) height; y++)
196 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000197 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
cristybb503372010-05-27 20:51:26 +0000198 for (x=0L; x < (ssize_t) (width/2L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000199 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
200 return(MagickTrue);
201}
202
cristyc4ea4a42011-01-24 01:43:30 +0000203static void CorrectPhaseLHS(const size_t width,const size_t height,
204 double *fourier)
cristy3ed852e2009-09-05 21:47:34 +0000205{
cristybb503372010-05-27 20:51:26 +0000206 ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000207 y;
208
cristybb503372010-05-27 20:51:26 +0000209 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000210 x;
211
cristybb503372010-05-27 20:51:26 +0000212 for (y=0L; y < (ssize_t) height; y++)
213 for (x=0L; x < (ssize_t) (width/2L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000214 fourier[y*width+x]*=(-1.0);
215}
216
217static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
218 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
219{
220 CacheView
221 *magnitude_view,
222 *phase_view;
223
224 double
225 *magnitude_source,
226 *phase_source;
227
228 Image
229 *magnitude_image,
230 *phase_image;
231
cristy3ed852e2009-09-05 21:47:34 +0000232 MagickBooleanType
233 status;
234
235 register IndexPacket
236 *indexes;
237
cristybb503372010-05-27 20:51:26 +0000238 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000239 x;
240
241 register PixelPacket
242 *q;
243
cristyc4ea4a42011-01-24 01:43:30 +0000244 ssize_t
245 i,
246 y;
247
cristy3ed852e2009-09-05 21:47:34 +0000248 magnitude_image=GetFirstImageInList(image);
249 phase_image=GetNextImageInList(image);
250 if (phase_image == (Image *) NULL)
251 {
252 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
253 "ImageSequenceRequired","`%s'",image->filename);
254 return(MagickFalse);
255 }
256 /*
257 Create "Fourier Transform" image from constituent arrays.
258 */
259 magnitude_source=(double *) AcquireQuantumMemory((size_t)
260 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
261 if (magnitude_source == (double *) NULL)
262 return(MagickFalse);
cristyc4ea4a42011-01-24 01:43:30 +0000263 (void) ResetMagickMemory(magnitude_source,0,fourier_info->height*
264 fourier_info->width*sizeof(*magnitude_source));
cristy3ed852e2009-09-05 21:47:34 +0000265 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
266 fourier_info->width*sizeof(*phase_source));
cristyc4ea4a42011-01-24 01:43:30 +0000267 if (phase_source == (double *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000268 {
269 (void) ThrowMagickException(exception,GetMagickModule(),
270 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
271 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
272 return(MagickFalse);
273 }
274 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
275 magnitude,magnitude_source);
276 if (status != MagickFalse)
277 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
278 phase_source);
279 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
280 if (fourier_info->modulus != MagickFalse)
281 {
282 i=0L;
cristybb503372010-05-27 20:51:26 +0000283 for (y=0L; y < (ssize_t) fourier_info->height; y++)
284 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000285 {
286 phase_source[i]/=(2.0*MagickPI);
287 phase_source[i]+=0.5;
288 i++;
289 }
290 }
291 magnitude_view=AcquireCacheView(magnitude_image);
292 phase_view=AcquireCacheView(phase_image);
293 i=0L;
cristybb503372010-05-27 20:51:26 +0000294 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000295 {
296 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
297 exception);
298 if (q == (PixelPacket *) NULL)
299 break;
300 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
cristybb503372010-05-27 20:51:26 +0000301 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000302 {
303 switch (fourier_info->channel)
304 {
305 case RedChannel:
306 default:
307 {
cristye85007d2010-06-06 22:51:36 +0000308 q->red=ClampToQuantum(QuantumRange*magnitude_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000309 break;
310 }
311 case GreenChannel:
312 {
cristye85007d2010-06-06 22:51:36 +0000313 q->green=ClampToQuantum(QuantumRange*magnitude_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000314 break;
315 }
316 case BlueChannel:
317 {
cristye85007d2010-06-06 22:51:36 +0000318 q->blue=ClampToQuantum(QuantumRange*magnitude_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000319 break;
320 }
321 case OpacityChannel:
322 {
cristye85007d2010-06-06 22:51:36 +0000323 q->opacity=ClampToQuantum(QuantumRange*magnitude_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000324 break;
325 }
326 case IndexChannel:
327 {
cristye85007d2010-06-06 22:51:36 +0000328 indexes[x]=ClampToQuantum(QuantumRange*magnitude_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000329 break;
330 }
331 case GrayChannels:
332 {
cristyc4ea4a42011-01-24 01:43:30 +0000333 SetGrayPixelComponent(q,ClampToQuantum(QuantumRange*
334 magnitude_source[i]));
cristy3ed852e2009-09-05 21:47:34 +0000335 break;
336 }
337 }
338 i++;
339 q++;
340 }
341 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
342 if (status == MagickFalse)
343 break;
344 }
345 i=0L;
cristybb503372010-05-27 20:51:26 +0000346 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000347 {
348 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
349 exception);
350 if (q == (PixelPacket *) NULL)
351 break;
352 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
cristybb503372010-05-27 20:51:26 +0000353 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000354 {
355 switch (fourier_info->channel)
356 {
357 case RedChannel:
358 default:
359 {
cristye85007d2010-06-06 22:51:36 +0000360 q->red=ClampToQuantum(QuantumRange*phase_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000361 break;
362 }
363 case GreenChannel:
364 {
cristye85007d2010-06-06 22:51:36 +0000365 q->green=ClampToQuantum(QuantumRange*phase_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000366 break;
367 }
368 case BlueChannel:
369 {
cristye85007d2010-06-06 22:51:36 +0000370 q->blue=ClampToQuantum(QuantumRange*phase_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000371 break;
372 }
373 case OpacityChannel:
374 {
cristye85007d2010-06-06 22:51:36 +0000375 q->opacity=ClampToQuantum(QuantumRange*phase_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000376 break;
377 }
378 case IndexChannel:
379 {
cristye85007d2010-06-06 22:51:36 +0000380 indexes[x]=ClampToQuantum(QuantumRange*phase_source[i]);
cristy3ed852e2009-09-05 21:47:34 +0000381 break;
382 }
383 case GrayChannels:
384 {
cristyc4ea4a42011-01-24 01:43:30 +0000385 SetGrayPixelComponent(q,ClampToQuantum(QuantumRange*phase_source[i]));
cristy3ed852e2009-09-05 21:47:34 +0000386 break;
387 }
388 }
389 i++;
390 q++;
391 }
392 status=SyncCacheViewAuthenticPixels(phase_view,exception);
393 if (status == MagickFalse)
394 break;
395 }
396 phase_view=DestroyCacheView(phase_view);
397 magnitude_view=DestroyCacheView(magnitude_view);
398 phase_source=(double *) RelinquishMagickMemory(phase_source);
399 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
400 return(status);
401}
402
403static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
404 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
405{
406 CacheView
407 *image_view;
408
409 double
410 n,
411 *source;
412
413 fftw_complex
414 *fourier;
415
416 fftw_plan
417 fftw_r2c_plan;
418
cristy3ed852e2009-09-05 21:47:34 +0000419 register const IndexPacket
420 *indexes;
421
422 register const PixelPacket
423 *p;
424
cristybb503372010-05-27 20:51:26 +0000425 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000426 i,
427 x;
428
cristyc4ea4a42011-01-24 01:43:30 +0000429 ssize_t
430 y;
431
cristy3ed852e2009-09-05 21:47:34 +0000432 /*
433 Generate the forward Fourier transform.
434 */
435 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
436 fourier_info->width*sizeof(*source));
437 if (source == (double *) NULL)
438 {
439 (void) ThrowMagickException(exception,GetMagickModule(),
440 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
441 return(MagickFalse);
442 }
cristyc4ea4a42011-01-24 01:43:30 +0000443 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
cristy3ed852e2009-09-05 21:47:34 +0000444 sizeof(*source));
445 i=0L;
446 image_view=AcquireCacheView(image);
cristybb503372010-05-27 20:51:26 +0000447 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000448 {
449 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
450 exception);
451 if (p == (const PixelPacket *) NULL)
452 break;
453 indexes=GetCacheViewVirtualIndexQueue(image_view);
cristybb503372010-05-27 20:51:26 +0000454 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000455 {
456 switch (fourier_info->channel)
457 {
458 case RedChannel:
459 default:
460 {
cristyce70c172010-01-07 17:15:30 +0000461 source[i]=QuantumScale*GetRedPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000462 break;
463 }
464 case GreenChannel:
465 {
cristyce70c172010-01-07 17:15:30 +0000466 source[i]=QuantumScale*GetGreenPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000467 break;
468 }
469 case BlueChannel:
470 {
cristyce70c172010-01-07 17:15:30 +0000471 source[i]=QuantumScale*GetBluePixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000472 break;
473 }
474 case OpacityChannel:
475 {
cristyce70c172010-01-07 17:15:30 +0000476 source[i]=QuantumScale*GetOpacityPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000477 break;
478 }
479 case IndexChannel:
480 {
481 source[i]=QuantumScale*indexes[x];
482 break;
483 }
484 case GrayChannels:
485 {
cristy04ec04d2011-01-24 14:18:55 +0000486 source[i]=QuantumScale*GetGrayPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000487 break;
488 }
489 }
490 i++;
491 p++;
492 }
493 }
494 image_view=DestroyCacheView(image_view);
cristyb41ee102010-10-04 16:46:15 +0000495 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
cristy3ed852e2009-09-05 21:47:34 +0000496 fourier_info->center*sizeof(*fourier));
497 if (fourier == (fftw_complex *) NULL)
498 {
499 (void) ThrowMagickException(exception,GetMagickModule(),
500 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
501 source=(double *) RelinquishMagickMemory(source);
502 return(MagickFalse);
503 }
cristyb5d5f722009-11-04 03:03:49 +0000504#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristy3ed852e2009-09-05 21:47:34 +0000505 #pragma omp critical (MagickCore_ForwardFourierTransform)
506#endif
507 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
508 source,fourier,FFTW_ESTIMATE);
509 fftw_execute(fftw_r2c_plan);
510 fftw_destroy_plan(fftw_r2c_plan);
511 source=(double *) RelinquishMagickMemory(source);
512 /*
513 Normalize Fourier transform.
514 */
515 n=(double) fourier_info->width*(double) fourier_info->width;
516 i=0L;
cristybb503372010-05-27 20:51:26 +0000517 for (y=0L; y < (ssize_t) fourier_info->height; y++)
518 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy56ed31c2010-03-22 00:46:21 +0000519 {
520#if defined(MAGICKCORE_HAVE_COMPLEX_H)
521 fourier[i]/=n;
522#else
523 fourier[i][0]/=n;
524 fourier[i][1]/=n;
525#endif
526 i++;
527 }
cristy3ed852e2009-09-05 21:47:34 +0000528 /*
529 Generate magnitude and phase (or real and imaginary).
530 */
531 i=0L;
532 if (fourier_info->modulus != MagickFalse)
cristybb503372010-05-27 20:51:26 +0000533 for (y=0L; y < (ssize_t) fourier_info->height; y++)
534 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +0000535 {
536 magnitude[i]=cabs(fourier[i]);
537 phase[i]=carg(fourier[i]);
538 i++;
539 }
540 else
cristybb503372010-05-27 20:51:26 +0000541 for (y=0L; y < (ssize_t) fourier_info->height; y++)
542 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +0000543 {
544 magnitude[i]=creal(fourier[i]);
545 phase[i]=cimag(fourier[i]);
546 i++;
547 }
cristyb41ee102010-10-04 16:46:15 +0000548 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
cristy3ed852e2009-09-05 21:47:34 +0000549 return(MagickTrue);
550}
551
552static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
553 const ChannelType channel,const MagickBooleanType modulus,
554 Image *fourier_image,ExceptionInfo *exception)
555{
556 double
557 *magnitude,
558 *phase;
559
560 fftw_complex
561 *fourier;
562
cristy56ed31c2010-03-22 00:46:21 +0000563 FourierInfo
564 fourier_info;
565
cristyc4ea4a42011-01-24 01:43:30 +0000566 MagickBooleanType
567 status;
568
cristy3ed852e2009-09-05 21:47:34 +0000569 size_t
570 extent;
571
572 fourier_info.width=image->columns;
573 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
574 ((image->rows % 2) != 0))
575 {
576 extent=image->columns < image->rows ? image->rows : image->columns;
577 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
578 }
579 fourier_info.height=fourier_info.width;
cristybb503372010-05-27 20:51:26 +0000580 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
cristy3ed852e2009-09-05 21:47:34 +0000581 fourier_info.channel=channel;
582 fourier_info.modulus=modulus;
583 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
584 fourier_info.center*sizeof(*magnitude));
585 if (magnitude == (double *) NULL)
586 {
587 (void) ThrowMagickException(exception,GetMagickModule(),
588 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
589 return(MagickFalse);
590 }
591 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
592 fourier_info.center*sizeof(*phase));
593 if (phase == (double *) NULL)
594 {
595 (void) ThrowMagickException(exception,GetMagickModule(),
596 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
597 magnitude=(double *) RelinquishMagickMemory(magnitude);
598 return(MagickFalse);
599 }
cristyb41ee102010-10-04 16:46:15 +0000600 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
cristy3ed852e2009-09-05 21:47:34 +0000601 fourier_info.center*sizeof(*fourier));
602 if (fourier == (fftw_complex *) NULL)
603 {
604 (void) ThrowMagickException(exception,GetMagickModule(),
605 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
606 phase=(double *) RelinquishMagickMemory(phase);
607 magnitude=(double *) RelinquishMagickMemory(magnitude);
608 return(MagickFalse);
609 }
610 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
611 if (status != MagickFalse)
612 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
613 exception);
cristyb41ee102010-10-04 16:46:15 +0000614 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
cristy3ed852e2009-09-05 21:47:34 +0000615 phase=(double *) RelinquishMagickMemory(phase);
616 magnitude=(double *) RelinquishMagickMemory(magnitude);
617 return(status);
618}
619#endif
620
621MagickExport Image *ForwardFourierTransformImage(const Image *image,
622 const MagickBooleanType modulus,ExceptionInfo *exception)
623{
624 Image
625 *fourier_image;
626
627 fourier_image=NewImageList();
628#if !defined(MAGICKCORE_FFTW_DELEGATE)
629 (void) modulus;
630 (void) ThrowMagickException(exception,GetMagickModule(),
631 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
632 image->filename);
633#else
634 {
635 Image
636 *magnitude_image;
637
cristybb503372010-05-27 20:51:26 +0000638 size_t
cristy3ed852e2009-09-05 21:47:34 +0000639 extent,
640 width;
641
642 width=image->columns;
643 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
644 ((image->rows % 2) != 0))
645 {
646 extent=image->columns < image->rows ? image->rows : image->columns;
647 width=(extent & 0x01) == 1 ? extent+1UL : extent;
648 }
649 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
650 if (magnitude_image != (Image *) NULL)
651 {
652 Image
653 *phase_image;
654
655 magnitude_image->storage_class=DirectClass;
656 magnitude_image->depth=32UL;
657 phase_image=CloneImage(image,width,width,MagickFalse,exception);
658 if (phase_image == (Image *) NULL)
659 magnitude_image=DestroyImage(magnitude_image);
660 else
661 {
662 MagickBooleanType
663 is_gray,
664 status;
665
cristy3ed852e2009-09-05 21:47:34 +0000666 phase_image->storage_class=DirectClass;
667 phase_image->depth=32UL;
668 AppendImageToList(&fourier_image,magnitude_image);
669 AppendImageToList(&fourier_image,phase_image);
670 status=MagickTrue;
671 is_gray=IsGrayImage(image,exception);
cristyb5d5f722009-11-04 03:03:49 +0000672#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristyb34ef052010-10-07 00:12:05 +0000673 #pragma omp parallel sections
cristy3ed852e2009-09-05 21:47:34 +0000674#endif
cristy3ed852e2009-09-05 21:47:34 +0000675 {
cristyb34ef052010-10-07 00:12:05 +0000676#if defined(MAGICKCORE_OPENMP_SUPPORT)
677 #pragma omp section
678#endif
cristy3ed852e2009-09-05 21:47:34 +0000679 {
cristyb34ef052010-10-07 00:12:05 +0000680 MagickBooleanType
681 thread_status;
682
683 if (is_gray != MagickFalse)
684 thread_status=ForwardFourierTransformChannel(image,
685 GrayChannels,modulus,fourier_image,exception);
686 else
687 thread_status=ForwardFourierTransformChannel(image,
688 RedChannel,modulus,fourier_image,exception);
689 if (thread_status == MagickFalse)
690 status=thread_status;
cristy3ed852e2009-09-05 21:47:34 +0000691 }
cristyb34ef052010-10-07 00:12:05 +0000692#if defined(MAGICKCORE_OPENMP_SUPPORT)
693 #pragma omp section
694#endif
695 {
696 MagickBooleanType
697 thread_status;
698
699 thread_status=MagickTrue;
700 if (is_gray == MagickFalse)
701 thread_status=ForwardFourierTransformChannel(image,
702 GreenChannel,modulus,fourier_image,exception);
703 if (thread_status == MagickFalse)
704 status=thread_status;
705 }
706#if defined(MAGICKCORE_OPENMP_SUPPORT)
707 #pragma omp section
708#endif
709 {
710 MagickBooleanType
711 thread_status;
712
713 thread_status=MagickTrue;
714 if (is_gray == MagickFalse)
715 thread_status=ForwardFourierTransformChannel(image,
716 BlueChannel,modulus,fourier_image,exception);
717 if (thread_status == MagickFalse)
718 status=thread_status;
719 }
720#if defined(MAGICKCORE_OPENMP_SUPPORT)
721 #pragma omp section
722#endif
723 {
724 MagickBooleanType
725 thread_status;
726
727 thread_status=MagickTrue;
728 if (image->matte != MagickFalse)
729 thread_status=ForwardFourierTransformChannel(image,
730 OpacityChannel,modulus,fourier_image,exception);
731 if (thread_status == MagickFalse)
732 status=thread_status;
733 }
734#if defined(MAGICKCORE_OPENMP_SUPPORT)
735 #pragma omp section
736#endif
737 {
738 MagickBooleanType
739 thread_status;
740
741 thread_status=MagickTrue;
742 if (image->colorspace == CMYKColorspace)
743 thread_status=ForwardFourierTransformChannel(image,
744 IndexChannel,modulus,fourier_image,exception);
745 if (thread_status == MagickFalse)
746 status=thread_status;
747 }
cristy3ed852e2009-09-05 21:47:34 +0000748 }
749 if (status == MagickFalse)
750 fourier_image=DestroyImageList(fourier_image);
751 fftw_cleanup();
752 }
753 }
754 }
755#endif
756 return(fourier_image);
757}
758
759/*
760%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
761% %
762% %
763% %
764% I n v e r s e F o u r i e r T r a n s f o r m I m a g e %
765% %
766% %
767% %
768%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
769%
770% InverseFourierTransformImage() implements the inverse discrete Fourier
771% transform (DFT) of the image either as a magnitude / phase or real /
772% imaginary image pair.
773%
774% The format of the InverseFourierTransformImage method is:
775%
cristyc9550792009-11-13 20:05:42 +0000776% Image *InverseFourierTransformImage(const Image *magnitude_image,
777% const Image *phase_image,const MagickBooleanType modulus,
778% ExceptionInfo *exception)
cristy3ed852e2009-09-05 21:47:34 +0000779%
780% A description of each parameter follows:
781%
cristyc9550792009-11-13 20:05:42 +0000782% o magnitude_image: the magnitude or real image.
783%
784% o phase_image: the phase or imaginary image.
cristy3ed852e2009-09-05 21:47:34 +0000785%
786% o modulus: if true, return transform as a magnitude / phase pair
787% otherwise a real / imaginary image pair.
788%
789% o exception: return any errors or warnings in this structure.
790%
791*/
792
793#if defined(MAGICKCORE_FFTW_DELEGATE)
cristybb503372010-05-27 20:51:26 +0000794static MagickBooleanType InverseQuadrantSwap(const size_t width,
795 const size_t height,const double *source,double *destination)
cristy3ed852e2009-09-05 21:47:34 +0000796{
cristyc4ea4a42011-01-24 01:43:30 +0000797 register ssize_t
798 x;
799
cristybb503372010-05-27 20:51:26 +0000800 ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000801 center,
802 y;
803
cristy3ed852e2009-09-05 21:47:34 +0000804 /*
805 Swap quadrants.
806 */
cristybb503372010-05-27 20:51:26 +0000807 center=(ssize_t) floor((double) width/2.0)+1L;
808 for (y=1L; y < (ssize_t) height; y++)
809 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000810 destination[center*(height-y)-x+width/2L]=source[y*width+x];
cristybb503372010-05-27 20:51:26 +0000811 for (y=0L; y < (ssize_t) height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000812 destination[center*y]=source[y*width+width/2L];
813 for (x=0L; x < center; x++)
814 destination[x]=source[center-x-1L];
cristybb503372010-05-27 20:51:26 +0000815 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
cristy3ed852e2009-09-05 21:47:34 +0000816}
817
818static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
cristyc9550792009-11-13 20:05:42 +0000819 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
820 ExceptionInfo *exception)
cristy3ed852e2009-09-05 21:47:34 +0000821{
822 CacheView
823 *magnitude_view,
824 *phase_view;
825
826 double
827 *magnitude,
828 *phase,
829 *magnitude_source,
830 *phase_source;
831
cristy3ed852e2009-09-05 21:47:34 +0000832 MagickBooleanType
833 status;
834
835 register const IndexPacket
836 *indexes;
837
838 register const PixelPacket
839 *p;
840
cristybb503372010-05-27 20:51:26 +0000841 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000842 i,
843 x;
844
cristyc4ea4a42011-01-24 01:43:30 +0000845 ssize_t
846 y;
847
cristy3ed852e2009-09-05 21:47:34 +0000848 /*
849 Inverse fourier - read image and break down into a double array.
850 */
cristy3ed852e2009-09-05 21:47:34 +0000851 magnitude_source=(double *) AcquireQuantumMemory((size_t)
852 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
853 if (magnitude_source == (double *) NULL)
854 {
855 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000856 ResourceLimitError,"MemoryAllocationFailed","`%s'",
857 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +0000858 return(MagickFalse);
859 }
860 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
cristyc4ea4a42011-01-24 01:43:30 +0000861 fourier_info->width*sizeof(*phase_source));
cristy3ed852e2009-09-05 21:47:34 +0000862 if (phase_source == (double *) NULL)
863 {
864 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000865 ResourceLimitError,"MemoryAllocationFailed","`%s'",
866 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +0000867 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
868 return(MagickFalse);
869 }
870 i=0L;
871 magnitude_view=AcquireCacheView(magnitude_image);
cristybb503372010-05-27 20:51:26 +0000872 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000873 {
874 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
875 exception);
876 if (p == (const PixelPacket *) NULL)
877 break;
878 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
cristybb503372010-05-27 20:51:26 +0000879 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000880 {
881 switch (fourier_info->channel)
882 {
883 case RedChannel:
884 default:
885 {
cristyce70c172010-01-07 17:15:30 +0000886 magnitude_source[i]=QuantumScale*GetRedPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000887 break;
888 }
889 case GreenChannel:
890 {
cristyce70c172010-01-07 17:15:30 +0000891 magnitude_source[i]=QuantumScale*GetGreenPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000892 break;
893 }
894 case BlueChannel:
895 {
cristyce70c172010-01-07 17:15:30 +0000896 magnitude_source[i]=QuantumScale*GetBluePixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000897 break;
898 }
899 case OpacityChannel:
900 {
cristyce70c172010-01-07 17:15:30 +0000901 magnitude_source[i]=QuantumScale*GetOpacityPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000902 break;
903 }
904 case IndexChannel:
905 {
906 magnitude_source[i]=QuantumScale*indexes[x];
907 break;
908 }
909 case GrayChannels:
910 {
cristyc4ea4a42011-01-24 01:43:30 +0000911 magnitude_source[i]=QuantumScale*GetGrayPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000912 break;
913 }
914 }
915 i++;
916 p++;
917 }
918 }
919 i=0L;
920 phase_view=AcquireCacheView(phase_image);
cristybb503372010-05-27 20:51:26 +0000921 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000922 {
923 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
924 exception);
925 if (p == (const PixelPacket *) NULL)
926 break;
927 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
cristybb503372010-05-27 20:51:26 +0000928 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000929 {
930 switch (fourier_info->channel)
931 {
932 case RedChannel:
933 default:
934 {
cristyce70c172010-01-07 17:15:30 +0000935 phase_source[i]=QuantumScale*GetRedPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000936 break;
937 }
938 case GreenChannel:
939 {
cristyce70c172010-01-07 17:15:30 +0000940 phase_source[i]=QuantumScale*GetGreenPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000941 break;
942 }
943 case BlueChannel:
944 {
cristyce70c172010-01-07 17:15:30 +0000945 phase_source[i]=QuantumScale*GetBluePixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000946 break;
947 }
948 case OpacityChannel:
949 {
cristyce70c172010-01-07 17:15:30 +0000950 phase_source[i]=QuantumScale*GetOpacityPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000951 break;
952 }
953 case IndexChannel:
954 {
955 phase_source[i]=QuantumScale*indexes[x];
956 break;
957 }
958 case GrayChannels:
959 {
cristyc4ea4a42011-01-24 01:43:30 +0000960 phase_source[i]=QuantumScale*GetGrayPixelComponent(p);
cristy3ed852e2009-09-05 21:47:34 +0000961 break;
962 }
963 }
964 i++;
965 p++;
966 }
967 }
968 if (fourier_info->modulus != MagickFalse)
969 {
970 i=0L;
cristybb503372010-05-27 20:51:26 +0000971 for (y=0L; y < (ssize_t) fourier_info->height; y++)
972 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000973 {
974 phase_source[i]-=0.5;
975 phase_source[i]*=(2.0*MagickPI);
976 i++;
977 }
978 }
979 magnitude_view=DestroyCacheView(magnitude_view);
980 phase_view=DestroyCacheView(phase_view);
981 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
982 fourier_info->center*sizeof(*magnitude));
983 if (magnitude == (double *) NULL)
984 {
985 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000986 ResourceLimitError,"MemoryAllocationFailed","`%s'",
987 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +0000988 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
989 phase_source=(double *) RelinquishMagickMemory(phase_source);
990 return(MagickFalse);
991 }
992 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
993 magnitude_source,magnitude);
994 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
cristyc4ea4a42011-01-24 01:43:30 +0000995 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
996 fourier_info->width*sizeof(*phase));
cristy3ed852e2009-09-05 21:47:34 +0000997 if (phase == (double *) NULL)
998 {
999 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001000 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1001 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001002 phase_source=(double *) RelinquishMagickMemory(phase_source);
1003 return(MagickFalse);
1004 }
1005 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1006 if (status != MagickFalse)
1007 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1008 phase_source,phase);
1009 phase_source=(double *) RelinquishMagickMemory(phase_source);
1010 /*
1011 Merge two sets.
1012 */
1013 i=0L;
1014 if (fourier_info->modulus != MagickFalse)
cristybb503372010-05-27 20:51:26 +00001015 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1016 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +00001017 {
cristy56ed31c2010-03-22 00:46:21 +00001018#if defined(MAGICKCORE_HAVE_COMPLEX_H)
cristy3ed852e2009-09-05 21:47:34 +00001019 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
cristy56ed31c2010-03-22 00:46:21 +00001020#else
1021 fourier[i][0]=magnitude[i]*cos(phase[i]);
1022 fourier[i][1]=magnitude[i]*sin(phase[i]);
1023#endif
cristy3ed852e2009-09-05 21:47:34 +00001024 i++;
1025 }
1026 else
cristybb503372010-05-27 20:51:26 +00001027 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1028 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +00001029 {
cristy56ed31c2010-03-22 00:46:21 +00001030#if defined(MAGICKCORE_HAVE_COMPLEX_H)
cristy3ed852e2009-09-05 21:47:34 +00001031 fourier[i]=magnitude[i]+I*phase[i];
cristy56ed31c2010-03-22 00:46:21 +00001032#else
1033 fourier[i][0]=magnitude[i];
1034 fourier[i][1]=phase[i];
1035#endif
cristy3ed852e2009-09-05 21:47:34 +00001036 i++;
1037 }
1038 phase=(double *) RelinquishMagickMemory(phase);
1039 magnitude=(double *) RelinquishMagickMemory(magnitude);
1040 return(status);
1041}
1042
1043static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1044 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1045{
1046 CacheView
1047 *image_view;
1048
1049 double
1050 *source;
1051
1052 fftw_plan
1053 fftw_c2r_plan;
1054
cristy3ed852e2009-09-05 21:47:34 +00001055 register IndexPacket
1056 *indexes;
1057
cristyc4ea4a42011-01-24 01:43:30 +00001058 register PixelPacket
1059 *q;
1060
cristybb503372010-05-27 20:51:26 +00001061 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +00001062 i,
1063 x;
1064
cristyc4ea4a42011-01-24 01:43:30 +00001065 ssize_t
1066 y;
cristy3ed852e2009-09-05 21:47:34 +00001067
cristyc4ea4a42011-01-24 01:43:30 +00001068 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1069 fourier_info->width*sizeof(*source));
cristy3ed852e2009-09-05 21:47:34 +00001070 if (source == (double *) NULL)
1071 {
1072 (void) ThrowMagickException(exception,GetMagickModule(),
1073 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1074 return(MagickFalse);
1075 }
cristyb5d5f722009-11-04 03:03:49 +00001076#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristy3ed852e2009-09-05 21:47:34 +00001077 #pragma omp critical (MagickCore_InverseFourierTransform)
1078#endif
cristydf079ac2010-09-10 01:45:44 +00001079 {
1080 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1081 fourier,source,FFTW_ESTIMATE);
1082 fftw_execute(fftw_c2r_plan);
1083 fftw_destroy_plan(fftw_c2r_plan);
1084 }
cristy3ed852e2009-09-05 21:47:34 +00001085 i=0L;
1086 image_view=AcquireCacheView(image);
cristybb503372010-05-27 20:51:26 +00001087 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +00001088 {
cristy85812052010-09-14 17:56:15 +00001089 if (y >= (ssize_t) image->rows)
1090 break;
1091 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1092 image->columns ? image->columns : fourier_info->width,1UL,exception);
cristy3ed852e2009-09-05 21:47:34 +00001093 if (q == (PixelPacket *) NULL)
1094 break;
1095 indexes=GetCacheViewAuthenticIndexQueue(image_view);
cristybb503372010-05-27 20:51:26 +00001096 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +00001097 {
1098 switch (fourier_info->channel)
1099 {
1100 case RedChannel:
1101 default:
1102 {
cristye85007d2010-06-06 22:51:36 +00001103 q->red=ClampToQuantum(QuantumRange*source[i]);
cristy3ed852e2009-09-05 21:47:34 +00001104 break;
1105 }
1106 case GreenChannel:
1107 {
cristye85007d2010-06-06 22:51:36 +00001108 q->green=ClampToQuantum(QuantumRange*source[i]);
cristy3ed852e2009-09-05 21:47:34 +00001109 break;
1110 }
1111 case BlueChannel:
1112 {
cristye85007d2010-06-06 22:51:36 +00001113 q->blue=ClampToQuantum(QuantumRange*source[i]);
cristy3ed852e2009-09-05 21:47:34 +00001114 break;
1115 }
1116 case OpacityChannel:
1117 {
cristye85007d2010-06-06 22:51:36 +00001118 q->opacity=ClampToQuantum(QuantumRange*source[i]);
cristy3ed852e2009-09-05 21:47:34 +00001119 break;
1120 }
1121 case IndexChannel:
1122 {
cristye85007d2010-06-06 22:51:36 +00001123 indexes[x]=ClampToQuantum(QuantumRange*source[i]);
cristy3ed852e2009-09-05 21:47:34 +00001124 break;
1125 }
1126 case GrayChannels:
1127 {
cristyc4ea4a42011-01-24 01:43:30 +00001128 SetGrayPixelComponent(q,ClampToQuantum(QuantumRange*source[i]));
cristy3ed852e2009-09-05 21:47:34 +00001129 break;
1130 }
1131 }
1132 i++;
1133 q++;
1134 }
1135 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1136 break;
1137 }
1138 image_view=DestroyCacheView(image_view);
1139 source=(double *) RelinquishMagickMemory(source);
1140 return(MagickTrue);
1141}
1142
cristyc9550792009-11-13 20:05:42 +00001143static MagickBooleanType InverseFourierTransformChannel(
1144 const Image *magnitude_image,const Image *phase_image,
cristy3ed852e2009-09-05 21:47:34 +00001145 const ChannelType channel,const MagickBooleanType modulus,
1146 Image *fourier_image,ExceptionInfo *exception)
1147{
1148 double
1149 *magnitude,
1150 *phase;
1151
1152 fftw_complex
1153 *fourier;
1154
1155 FourierInfo
1156 fourier_info;
1157
1158 MagickBooleanType
1159 status;
1160
1161 size_t
1162 extent;
1163
cristyc9550792009-11-13 20:05:42 +00001164 fourier_info.width=magnitude_image->columns;
1165 if ((magnitude_image->columns != magnitude_image->rows) ||
1166 ((magnitude_image->columns % 2) != 0) ||
1167 ((magnitude_image->rows % 2) != 0))
cristy3ed852e2009-09-05 21:47:34 +00001168 {
cristyc9550792009-11-13 20:05:42 +00001169 extent=magnitude_image->columns < magnitude_image->rows ?
1170 magnitude_image->rows : magnitude_image->columns;
cristy3ed852e2009-09-05 21:47:34 +00001171 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1172 }
1173 fourier_info.height=fourier_info.width;
cristybb503372010-05-27 20:51:26 +00001174 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
cristy3ed852e2009-09-05 21:47:34 +00001175 fourier_info.channel=channel;
1176 fourier_info.modulus=modulus;
1177 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
cristyc4ea4a42011-01-24 01:43:30 +00001178 fourier_info.center*sizeof(*magnitude));
cristy3ed852e2009-09-05 21:47:34 +00001179 if (magnitude == (double *) NULL)
1180 {
1181 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001182 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1183 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001184 return(MagickFalse);
1185 }
1186 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
cristyc4ea4a42011-01-24 01:43:30 +00001187 fourier_info.center*sizeof(*phase));
cristy3ed852e2009-09-05 21:47:34 +00001188 if (phase == (double *) NULL)
1189 {
1190 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001191 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1192 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001193 magnitude=(double *) RelinquishMagickMemory(magnitude);
1194 return(MagickFalse);
1195 }
cristyb41ee102010-10-04 16:46:15 +00001196 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
cristy3ed852e2009-09-05 21:47:34 +00001197 fourier_info.center*sizeof(*fourier));
1198 if (fourier == (fftw_complex *) NULL)
1199 {
1200 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001201 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1202 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001203 phase=(double *) RelinquishMagickMemory(phase);
1204 magnitude=(double *) RelinquishMagickMemory(magnitude);
1205 return(MagickFalse);
1206 }
cristyc9550792009-11-13 20:05:42 +00001207 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1208 exception);
cristy3ed852e2009-09-05 21:47:34 +00001209 if (status != MagickFalse)
1210 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1211 exception);
cristyb41ee102010-10-04 16:46:15 +00001212 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
cristy3ed852e2009-09-05 21:47:34 +00001213 phase=(double *) RelinquishMagickMemory(phase);
1214 magnitude=(double *) RelinquishMagickMemory(magnitude);
1215 return(status);
1216}
1217#endif
1218
cristyc9550792009-11-13 20:05:42 +00001219MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1220 const Image *phase_image,const MagickBooleanType modulus,
1221 ExceptionInfo *exception)
cristy3ed852e2009-09-05 21:47:34 +00001222{
1223 Image
1224 *fourier_image;
1225
cristyc9550792009-11-13 20:05:42 +00001226 assert(magnitude_image != (Image *) NULL);
1227 assert(magnitude_image->signature == MagickSignature);
1228 if (magnitude_image->debug != MagickFalse)
1229 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1230 magnitude_image->filename);
1231 if (phase_image == (Image *) NULL)
1232 {
1233 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1234 "ImageSequenceRequired","`%s'",magnitude_image->filename);
cristy9372a152009-11-18 01:42:16 +00001235 return((Image *) NULL);
cristyc9550792009-11-13 20:05:42 +00001236 }
cristy3ed852e2009-09-05 21:47:34 +00001237#if !defined(MAGICKCORE_FFTW_DELEGATE)
1238 fourier_image=(Image *) NULL;
1239 (void) modulus;
1240 (void) ThrowMagickException(exception,GetMagickModule(),
1241 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
cristyc9550792009-11-13 20:05:42 +00001242 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001243#else
1244 {
cristyc9550792009-11-13 20:05:42 +00001245 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1246 magnitude_image->rows,MagickFalse,exception);
cristy3ed852e2009-09-05 21:47:34 +00001247 if (fourier_image != (Image *) NULL)
1248 {
1249 MagickBooleanType
1250 is_gray,
1251 status;
1252
cristy3ed852e2009-09-05 21:47:34 +00001253 status=MagickTrue;
cristyc9550792009-11-13 20:05:42 +00001254 is_gray=IsGrayImage(magnitude_image,exception);
1255 if (is_gray != MagickFalse)
1256 is_gray=IsGrayImage(phase_image,exception);
cristyb5d5f722009-11-04 03:03:49 +00001257#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristyb34ef052010-10-07 00:12:05 +00001258 #pragma omp parallel sections
cristy3ed852e2009-09-05 21:47:34 +00001259#endif
cristy3ed852e2009-09-05 21:47:34 +00001260 {
cristyb34ef052010-10-07 00:12:05 +00001261#if defined(MAGICKCORE_OPENMP_SUPPORT)
1262 #pragma omp section
1263#endif
cristy3ed852e2009-09-05 21:47:34 +00001264 {
cristyb34ef052010-10-07 00:12:05 +00001265 MagickBooleanType
1266 thread_status;
1267
1268 if (is_gray != MagickFalse)
1269 thread_status=InverseFourierTransformChannel(magnitude_image,
1270 phase_image,GrayChannels,modulus,fourier_image,exception);
1271 else
cristyc9550792009-11-13 20:05:42 +00001272 thread_status=InverseFourierTransformChannel(magnitude_image,
1273 phase_image,RedChannel,modulus,fourier_image,exception);
cristyb34ef052010-10-07 00:12:05 +00001274 if (thread_status == MagickFalse)
1275 status=thread_status;
cristy3ed852e2009-09-05 21:47:34 +00001276 }
cristyb34ef052010-10-07 00:12:05 +00001277#if defined(MAGICKCORE_OPENMP_SUPPORT)
1278 #pragma omp section
1279#endif
1280 {
1281 MagickBooleanType
1282 thread_status;
1283
1284 thread_status=MagickTrue;
1285 if (is_gray == MagickFalse)
1286 thread_status=InverseFourierTransformChannel(magnitude_image,
1287 phase_image,GreenChannel,modulus,fourier_image,exception);
1288 if (thread_status == MagickFalse)
1289 status=thread_status;
1290 }
1291#if defined(MAGICKCORE_OPENMP_SUPPORT)
1292 #pragma omp section
1293#endif
1294 {
1295 MagickBooleanType
1296 thread_status;
1297
1298 thread_status=MagickTrue;
1299 if (is_gray == MagickFalse)
1300 thread_status=InverseFourierTransformChannel(magnitude_image,
1301 phase_image,BlueChannel,modulus,fourier_image,exception);
1302 if (thread_status == MagickFalse)
1303 status=thread_status;
1304 }
1305#if defined(MAGICKCORE_OPENMP_SUPPORT)
1306 #pragma omp section
1307#endif
1308 {
1309 MagickBooleanType
1310 thread_status;
1311
1312 thread_status=MagickTrue;
1313 if (magnitude_image->matte != MagickFalse)
1314 thread_status=InverseFourierTransformChannel(magnitude_image,
1315 phase_image,OpacityChannel,modulus,fourier_image,exception);
1316 if (thread_status == MagickFalse)
1317 status=thread_status;
1318 }
1319#if defined(MAGICKCORE_OPENMP_SUPPORT)
1320 #pragma omp section
1321#endif
1322 {
1323 MagickBooleanType
1324 thread_status;
1325
1326 thread_status=MagickTrue;
1327 if (magnitude_image->colorspace == CMYKColorspace)
1328 thread_status=InverseFourierTransformChannel(magnitude_image,
1329 phase_image,IndexChannel,modulus,fourier_image,exception);
1330 if (thread_status == MagickFalse)
1331 status=thread_status;
1332 }
cristy3ed852e2009-09-05 21:47:34 +00001333 }
1334 if (status == MagickFalse)
1335 fourier_image=DestroyImage(fourier_image);
1336 }
cristyb34ef052010-10-07 00:12:05 +00001337 fftw_cleanup();
cristy3ed852e2009-09-05 21:47:34 +00001338 }
1339#endif
1340 return(fourier_image);
1341}