blob: b7fbf67679b6a470d7cc058988213d599096a5e0 [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*/
cristy4c08aed2011-07-01 19:47:50 +000045#include "MagickCore/studio.h"
46#include "MagickCore/attribute.h"
47#include "MagickCore/cache.h"
48#include "MagickCore/image.h"
49#include "MagickCore/image-private.h"
50#include "MagickCore/list.h"
51#include "MagickCore/fourier.h"
52#include "MagickCore/log.h"
53#include "MagickCore/memory_.h"
54#include "MagickCore/monitor.h"
55#include "MagickCore/pixel-accessor.h"
56#include "MagickCore/property.h"
57#include "MagickCore/quantum-private.h"
58#include "MagickCore/thread-private.h"
cristy3ed852e2009-09-05 21:47:34 +000059#if defined(MAGICKCORE_FFTW_DELEGATE)
cristy56ed31c2010-03-22 00:46:21 +000060#if defined(MAGICKCORE_HAVE_COMPLEX_H)
cristy3ed852e2009-09-05 21:47:34 +000061#include <complex.h>
cristy56ed31c2010-03-22 00:46:21 +000062#endif
cristy3ed852e2009-09-05 21:47:34 +000063#include <fftw3.h>
cristy47b022b2011-01-18 22:29:21 +000064#if !defined(MAGICKCORE_HAVE_CABS)
65#define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
66#endif
67#if !defined(MAGICKCORE_HAVE_CARG)
cristy4da3ba32011-02-07 16:58:11 +000068#define carg(z) (atan2(cimag(z),creal(z)))
cristy47b022b2011-01-18 22:29:21 +000069#endif
70#if !defined(MAGICKCORE_HAVE_CIMAG)
71#define cimag(z) (z[1])
72#endif
73#if !defined(MAGICKCORE_HAVE_CREAL)
74#define creal(z) (z[0])
75#endif
cristy3ed852e2009-09-05 21:47:34 +000076#endif
77
78/*
79 Typedef declarations.
80*/
81typedef struct _FourierInfo
82{
83 ChannelType
84 channel;
85
86 MagickBooleanType
87 modulus;
88
cristybb503372010-05-27 20:51:26 +000089 size_t
cristy3ed852e2009-09-05 21:47:34 +000090 width,
91 height;
92
cristybb503372010-05-27 20:51:26 +000093 ssize_t
cristy3ed852e2009-09-05 21:47:34 +000094 center;
95} FourierInfo;
96
97/*
98%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99% %
100% %
101% %
102% 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 %
103% %
104% %
105% %
106%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107%
108% ForwardFourierTransformImage() implements the discrete Fourier transform
109% (DFT) of the image either as a magnitude / phase or real / imaginary image
110% pair.
111%
112% The format of the ForwadFourierTransformImage method is:
113%
114% Image *ForwardFourierTransformImage(const Image *image,
115% const MagickBooleanType modulus,ExceptionInfo *exception)
116%
117% A description of each parameter follows:
118%
119% o image: the image.
120%
121% o modulus: if true, return as transform as a magnitude / phase pair
122% otherwise a real / imaginary image pair.
123%
124% o exception: return any errors or warnings in this structure.
125%
126*/
127
128#if defined(MAGICKCORE_FFTW_DELEGATE)
129
cristyc4ea4a42011-01-24 01:43:30 +0000130static MagickBooleanType RollFourier(const size_t width,const size_t height,
131 const ssize_t x_offset,const ssize_t y_offset,double *fourier)
cristy3ed852e2009-09-05 21:47:34 +0000132{
133 double
134 *roll;
135
cristyc4ea4a42011-01-24 01:43:30 +0000136 register ssize_t
137 i,
138 x;
139
cristybb503372010-05-27 20:51:26 +0000140 ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000141 u,
142 v,
143 y;
144
cristy3ed852e2009-09-05 21:47:34 +0000145 /*
cristy2da05352010-09-15 19:22:17 +0000146 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
cristy3ed852e2009-09-05 21:47:34 +0000147 */
cristyc4ea4a42011-01-24 01:43:30 +0000148 roll=(double *) AcquireQuantumMemory((size_t) height,width*sizeof(*roll));
cristy3ed852e2009-09-05 21:47:34 +0000149 if (roll == (double *) NULL)
150 return(MagickFalse);
151 i=0L;
cristybb503372010-05-27 20:51:26 +0000152 for (y=0L; y < (ssize_t) height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000153 {
154 if (y_offset < 0L)
cristybb503372010-05-27 20:51:26 +0000155 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
cristy3ed852e2009-09-05 21:47:34 +0000156 else
cristybb503372010-05-27 20:51:26 +0000157 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
cristy3ed852e2009-09-05 21:47:34 +0000158 y+y_offset;
cristybb503372010-05-27 20:51:26 +0000159 for (x=0L; x < (ssize_t) width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000160 {
161 if (x_offset < 0L)
cristybb503372010-05-27 20:51:26 +0000162 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
cristy3ed852e2009-09-05 21:47:34 +0000163 else
cristybb503372010-05-27 20:51:26 +0000164 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
cristy3ed852e2009-09-05 21:47:34 +0000165 x+x_offset;
166 roll[v*width+u]=fourier[i++];
cristyc4ea4a42011-01-24 01:43:30 +0000167 }
cristy3ed852e2009-09-05 21:47:34 +0000168 }
cristyc4ea4a42011-01-24 01:43:30 +0000169 (void) CopyMagickMemory(fourier,roll,height*width*sizeof(*roll));
cristy3ed852e2009-09-05 21:47:34 +0000170 roll=(double *) RelinquishMagickMemory(roll);
171 return(MagickTrue);
172}
173
cristybb503372010-05-27 20:51:26 +0000174static MagickBooleanType ForwardQuadrantSwap(const size_t width,
175 const size_t height,double *source,double *destination)
cristy3ed852e2009-09-05 21:47:34 +0000176{
cristy3ed852e2009-09-05 21:47:34 +0000177 MagickBooleanType
178 status;
179
cristybb503372010-05-27 20:51:26 +0000180 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000181 x;
182
cristyc4ea4a42011-01-24 01:43:30 +0000183 ssize_t
184 center,
185 y;
186
cristy3ed852e2009-09-05 21:47:34 +0000187 /*
188 Swap quadrants.
189 */
cristybb503372010-05-27 20:51:26 +0000190 center=(ssize_t) floor((double) width/2L)+1L;
191 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
cristy3ed852e2009-09-05 21:47:34 +0000192 if (status == MagickFalse)
193 return(MagickFalse);
cristybb503372010-05-27 20:51:26 +0000194 for (y=0L; y < (ssize_t) height; y++)
195 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000196 destination[width*y+x+width/2L]=source[center*y+x];
cristybb503372010-05-27 20:51:26 +0000197 for (y=1; y < (ssize_t) height; y++)
198 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000199 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
cristybb503372010-05-27 20:51:26 +0000200 for (x=0L; x < (ssize_t) (width/2L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000201 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
202 return(MagickTrue);
203}
204
cristyc4ea4a42011-01-24 01:43:30 +0000205static void CorrectPhaseLHS(const size_t width,const size_t height,
206 double *fourier)
cristy3ed852e2009-09-05 21:47:34 +0000207{
cristybb503372010-05-27 20:51:26 +0000208 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000209 x;
210
cristy9d314ff2011-03-09 01:30:28 +0000211 ssize_t
212 y;
213
cristybb503372010-05-27 20:51:26 +0000214 for (y=0L; y < (ssize_t) height; y++)
215 for (x=0L; x < (ssize_t) (width/2L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000216 fourier[y*width+x]*=(-1.0);
217}
218
219static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
220 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
221{
222 CacheView
223 *magnitude_view,
224 *phase_view;
225
226 double
227 *magnitude_source,
228 *phase_source;
229
230 Image
231 *magnitude_image,
232 *phase_image;
233
cristy3ed852e2009-09-05 21:47:34 +0000234 MagickBooleanType
235 status;
236
cristybb503372010-05-27 20:51:26 +0000237 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000238 x;
239
cristy4c08aed2011-07-01 19:47:50 +0000240 register Quantum
cristy3ed852e2009-09-05 21:47:34 +0000241 *q;
242
cristyc4ea4a42011-01-24 01:43:30 +0000243 ssize_t
244 i,
245 y;
246
cristy3ed852e2009-09-05 21:47:34 +0000247 magnitude_image=GetFirstImageInList(image);
248 phase_image=GetNextImageInList(image);
249 if (phase_image == (Image *) NULL)
250 {
251 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
252 "ImageSequenceRequired","`%s'",image->filename);
253 return(MagickFalse);
254 }
255 /*
256 Create "Fourier Transform" image from constituent arrays.
257 */
258 magnitude_source=(double *) AcquireQuantumMemory((size_t)
259 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
260 if (magnitude_source == (double *) NULL)
261 return(MagickFalse);
cristyc4ea4a42011-01-24 01:43:30 +0000262 (void) ResetMagickMemory(magnitude_source,0,fourier_info->height*
263 fourier_info->width*sizeof(*magnitude_source));
cristy3ed852e2009-09-05 21:47:34 +0000264 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
265 fourier_info->width*sizeof(*phase_source));
cristyc4ea4a42011-01-24 01:43:30 +0000266 if (phase_source == (double *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000267 {
268 (void) ThrowMagickException(exception,GetMagickModule(),
269 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
270 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
271 return(MagickFalse);
272 }
273 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
274 magnitude,magnitude_source);
275 if (status != MagickFalse)
276 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
277 phase_source);
278 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
279 if (fourier_info->modulus != MagickFalse)
280 {
281 i=0L;
cristybb503372010-05-27 20:51:26 +0000282 for (y=0L; y < (ssize_t) fourier_info->height; y++)
283 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000284 {
285 phase_source[i]/=(2.0*MagickPI);
286 phase_source[i]+=0.5;
287 i++;
288 }
289 }
290 magnitude_view=AcquireCacheView(magnitude_image);
291 phase_view=AcquireCacheView(phase_image);
292 i=0L;
cristybb503372010-05-27 20:51:26 +0000293 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000294 {
295 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
296 exception);
cristy4c08aed2011-07-01 19:47:50 +0000297 if (q == (const Quantum *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000298 break;
cristybb503372010-05-27 20:51:26 +0000299 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000300 {
301 switch (fourier_info->channel)
302 {
303 case RedChannel:
304 default:
305 {
cristy4c08aed2011-07-01 19:47:50 +0000306 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
307 magnitude_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000308 break;
309 }
310 case GreenChannel:
311 {
cristy4c08aed2011-07-01 19:47:50 +0000312 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
313 magnitude_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000314 break;
315 }
316 case BlueChannel:
317 {
cristy4c08aed2011-07-01 19:47:50 +0000318 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
319 magnitude_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000320 break;
321 }
cristy4c08aed2011-07-01 19:47:50 +0000322 case BlackChannel:
cristy3ed852e2009-09-05 21:47:34 +0000323 {
cristy4c08aed2011-07-01 19:47:50 +0000324 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
325 magnitude_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000326 break;
327 }
cristy4c08aed2011-07-01 19:47:50 +0000328 case AlphaChannel:
cristy3ed852e2009-09-05 21:47:34 +0000329 {
cristy4c08aed2011-07-01 19:47:50 +0000330 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
331 magnitude_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000332 break;
333 }
334 case GrayChannels:
335 {
cristy4c08aed2011-07-01 19:47:50 +0000336 SetPixelGray(magnitude_image,ClampToQuantum(QuantumRange*
337 magnitude_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000338 break;
339 }
340 }
341 i++;
cristyed231572011-07-14 02:18:59 +0000342 q+=GetPixelChannels(magnitude_image);
cristy3ed852e2009-09-05 21:47:34 +0000343 }
344 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
345 if (status == MagickFalse)
346 break;
347 }
348 i=0L;
cristybb503372010-05-27 20:51:26 +0000349 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000350 {
351 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
352 exception);
cristy4c08aed2011-07-01 19:47:50 +0000353 if (q == (const Quantum *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000354 break;
cristybb503372010-05-27 20:51:26 +0000355 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000356 {
357 switch (fourier_info->channel)
358 {
359 case RedChannel:
360 default:
361 {
cristy4c08aed2011-07-01 19:47:50 +0000362 SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
363 phase_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000364 break;
365 }
366 case GreenChannel:
367 {
cristy4c08aed2011-07-01 19:47:50 +0000368 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
369 phase_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000370 break;
371 }
372 case BlueChannel:
373 {
cristy4c08aed2011-07-01 19:47:50 +0000374 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
375 phase_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000376 break;
377 }
cristy4c08aed2011-07-01 19:47:50 +0000378 case BlackChannel:
cristy3ed852e2009-09-05 21:47:34 +0000379 {
cristy4c08aed2011-07-01 19:47:50 +0000380 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
381 phase_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000382 break;
383 }
cristy4c08aed2011-07-01 19:47:50 +0000384 case AlphaChannel:
cristy3ed852e2009-09-05 21:47:34 +0000385 {
cristy4c08aed2011-07-01 19:47:50 +0000386 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
387 phase_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000388 break;
389 }
390 case GrayChannels:
391 {
cristy4c08aed2011-07-01 19:47:50 +0000392 SetPixelGray(phase_image,ClampToQuantum(QuantumRange*
393 phase_source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +0000394 break;
395 }
396 }
397 i++;
cristyed231572011-07-14 02:18:59 +0000398 q+=GetPixelChannels(phase_image);
cristy3ed852e2009-09-05 21:47:34 +0000399 }
400 status=SyncCacheViewAuthenticPixels(phase_view,exception);
401 if (status == MagickFalse)
402 break;
403 }
404 phase_view=DestroyCacheView(phase_view);
405 magnitude_view=DestroyCacheView(magnitude_view);
406 phase_source=(double *) RelinquishMagickMemory(phase_source);
407 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
408 return(status);
409}
410
411static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
412 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
413{
414 CacheView
415 *image_view;
416
417 double
418 n,
419 *source;
420
421 fftw_complex
422 *fourier;
423
424 fftw_plan
425 fftw_r2c_plan;
426
cristy4c08aed2011-07-01 19:47:50 +0000427 register const Quantum
cristy3ed852e2009-09-05 21:47:34 +0000428 *p;
429
cristybb503372010-05-27 20:51:26 +0000430 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000431 i,
432 x;
433
cristyc4ea4a42011-01-24 01:43:30 +0000434 ssize_t
435 y;
436
cristy3ed852e2009-09-05 21:47:34 +0000437 /*
438 Generate the forward Fourier transform.
439 */
440 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
441 fourier_info->width*sizeof(*source));
442 if (source == (double *) NULL)
443 {
444 (void) ThrowMagickException(exception,GetMagickModule(),
445 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
446 return(MagickFalse);
447 }
cristyc4ea4a42011-01-24 01:43:30 +0000448 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
cristy3ed852e2009-09-05 21:47:34 +0000449 sizeof(*source));
450 i=0L;
451 image_view=AcquireCacheView(image);
cristybb503372010-05-27 20:51:26 +0000452 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000453 {
454 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
455 exception);
cristy4c08aed2011-07-01 19:47:50 +0000456 if (p == (const Quantum *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000457 break;
cristybb503372010-05-27 20:51:26 +0000458 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000459 {
460 switch (fourier_info->channel)
461 {
462 case RedChannel:
463 default:
464 {
cristy4c08aed2011-07-01 19:47:50 +0000465 source[i]=QuantumScale*GetPixelRed(image,p);
cristy3ed852e2009-09-05 21:47:34 +0000466 break;
467 }
468 case GreenChannel:
469 {
cristy4c08aed2011-07-01 19:47:50 +0000470 source[i]=QuantumScale*GetPixelGreen(image,p);
cristy3ed852e2009-09-05 21:47:34 +0000471 break;
472 }
473 case BlueChannel:
474 {
cristy4c08aed2011-07-01 19:47:50 +0000475 source[i]=QuantumScale*GetPixelBlue(image,p);
cristy3ed852e2009-09-05 21:47:34 +0000476 break;
477 }
cristy4c08aed2011-07-01 19:47:50 +0000478 case BlackChannel:
cristy3ed852e2009-09-05 21:47:34 +0000479 {
cristy4c08aed2011-07-01 19:47:50 +0000480 source[i]=QuantumScale*GetPixelBlack(image,p);
cristy3ed852e2009-09-05 21:47:34 +0000481 break;
482 }
cristy4c08aed2011-07-01 19:47:50 +0000483 case AlphaChannel:
cristy3ed852e2009-09-05 21:47:34 +0000484 {
cristy4c08aed2011-07-01 19:47:50 +0000485 source[i]=QuantumScale*GetPixelAlpha(image,p);
cristy3ed852e2009-09-05 21:47:34 +0000486 break;
487 }
488 case GrayChannels:
489 {
cristy4c08aed2011-07-01 19:47:50 +0000490 source[i]=QuantumScale*GetPixelGray(image,p);
cristy3ed852e2009-09-05 21:47:34 +0000491 break;
492 }
493 }
494 i++;
cristyed231572011-07-14 02:18:59 +0000495 p+=GetPixelChannels(image);
cristy3ed852e2009-09-05 21:47:34 +0000496 }
497 }
498 image_view=DestroyCacheView(image_view);
cristyb41ee102010-10-04 16:46:15 +0000499 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
cristy3ed852e2009-09-05 21:47:34 +0000500 fourier_info->center*sizeof(*fourier));
501 if (fourier == (fftw_complex *) NULL)
502 {
503 (void) ThrowMagickException(exception,GetMagickModule(),
504 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
505 source=(double *) RelinquishMagickMemory(source);
506 return(MagickFalse);
507 }
cristyb5d5f722009-11-04 03:03:49 +0000508#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristy3ed852e2009-09-05 21:47:34 +0000509 #pragma omp critical (MagickCore_ForwardFourierTransform)
510#endif
511 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
512 source,fourier,FFTW_ESTIMATE);
513 fftw_execute(fftw_r2c_plan);
514 fftw_destroy_plan(fftw_r2c_plan);
515 source=(double *) RelinquishMagickMemory(source);
516 /*
517 Normalize Fourier transform.
518 */
519 n=(double) fourier_info->width*(double) fourier_info->width;
520 i=0L;
cristybb503372010-05-27 20:51:26 +0000521 for (y=0L; y < (ssize_t) fourier_info->height; y++)
522 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy56ed31c2010-03-22 00:46:21 +0000523 {
524#if defined(MAGICKCORE_HAVE_COMPLEX_H)
525 fourier[i]/=n;
526#else
527 fourier[i][0]/=n;
528 fourier[i][1]/=n;
529#endif
530 i++;
531 }
cristy3ed852e2009-09-05 21:47:34 +0000532 /*
533 Generate magnitude and phase (or real and imaginary).
534 */
535 i=0L;
536 if (fourier_info->modulus != MagickFalse)
cristybb503372010-05-27 20:51:26 +0000537 for (y=0L; y < (ssize_t) fourier_info->height; y++)
538 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +0000539 {
540 magnitude[i]=cabs(fourier[i]);
541 phase[i]=carg(fourier[i]);
542 i++;
543 }
544 else
cristybb503372010-05-27 20:51:26 +0000545 for (y=0L; y < (ssize_t) fourier_info->height; y++)
546 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +0000547 {
548 magnitude[i]=creal(fourier[i]);
549 phase[i]=cimag(fourier[i]);
550 i++;
551 }
cristyb41ee102010-10-04 16:46:15 +0000552 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
cristy3ed852e2009-09-05 21:47:34 +0000553 return(MagickTrue);
554}
555
556static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
557 const ChannelType channel,const MagickBooleanType modulus,
558 Image *fourier_image,ExceptionInfo *exception)
559{
560 double
561 *magnitude,
562 *phase;
563
564 fftw_complex
565 *fourier;
566
cristy56ed31c2010-03-22 00:46:21 +0000567 FourierInfo
568 fourier_info;
569
cristyc4ea4a42011-01-24 01:43:30 +0000570 MagickBooleanType
571 status;
572
cristy3ed852e2009-09-05 21:47:34 +0000573 size_t
574 extent;
575
576 fourier_info.width=image->columns;
577 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
578 ((image->rows % 2) != 0))
579 {
580 extent=image->columns < image->rows ? image->rows : image->columns;
581 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
582 }
583 fourier_info.height=fourier_info.width;
cristybb503372010-05-27 20:51:26 +0000584 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
cristy3ed852e2009-09-05 21:47:34 +0000585 fourier_info.channel=channel;
586 fourier_info.modulus=modulus;
587 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
588 fourier_info.center*sizeof(*magnitude));
589 if (magnitude == (double *) NULL)
590 {
591 (void) ThrowMagickException(exception,GetMagickModule(),
592 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
593 return(MagickFalse);
594 }
595 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
596 fourier_info.center*sizeof(*phase));
597 if (phase == (double *) NULL)
598 {
599 (void) ThrowMagickException(exception,GetMagickModule(),
600 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
601 magnitude=(double *) RelinquishMagickMemory(magnitude);
602 return(MagickFalse);
603 }
cristyb41ee102010-10-04 16:46:15 +0000604 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
cristy3ed852e2009-09-05 21:47:34 +0000605 fourier_info.center*sizeof(*fourier));
606 if (fourier == (fftw_complex *) NULL)
607 {
608 (void) ThrowMagickException(exception,GetMagickModule(),
609 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
610 phase=(double *) RelinquishMagickMemory(phase);
611 magnitude=(double *) RelinquishMagickMemory(magnitude);
612 return(MagickFalse);
613 }
614 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
615 if (status != MagickFalse)
616 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
617 exception);
cristyb41ee102010-10-04 16:46:15 +0000618 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
cristy3ed852e2009-09-05 21:47:34 +0000619 phase=(double *) RelinquishMagickMemory(phase);
620 magnitude=(double *) RelinquishMagickMemory(magnitude);
621 return(status);
622}
623#endif
624
625MagickExport Image *ForwardFourierTransformImage(const Image *image,
626 const MagickBooleanType modulus,ExceptionInfo *exception)
627{
628 Image
629 *fourier_image;
630
631 fourier_image=NewImageList();
632#if !defined(MAGICKCORE_FFTW_DELEGATE)
633 (void) modulus;
634 (void) ThrowMagickException(exception,GetMagickModule(),
635 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
636 image->filename);
637#else
638 {
639 Image
640 *magnitude_image;
641
cristybb503372010-05-27 20:51:26 +0000642 size_t
cristy3ed852e2009-09-05 21:47:34 +0000643 extent,
644 width;
645
646 width=image->columns;
647 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
648 ((image->rows % 2) != 0))
649 {
650 extent=image->columns < image->rows ? image->rows : image->columns;
651 width=(extent & 0x01) == 1 ? extent+1UL : extent;
652 }
653 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
654 if (magnitude_image != (Image *) NULL)
655 {
656 Image
657 *phase_image;
658
659 magnitude_image->storage_class=DirectClass;
660 magnitude_image->depth=32UL;
661 phase_image=CloneImage(image,width,width,MagickFalse,exception);
662 if (phase_image == (Image *) NULL)
663 magnitude_image=DestroyImage(magnitude_image);
664 else
665 {
666 MagickBooleanType
667 is_gray,
668 status;
669
cristy3ed852e2009-09-05 21:47:34 +0000670 phase_image->storage_class=DirectClass;
671 phase_image->depth=32UL;
672 AppendImageToList(&fourier_image,magnitude_image);
673 AppendImageToList(&fourier_image,phase_image);
674 status=MagickTrue;
cristy4c08aed2011-07-01 19:47:50 +0000675 is_gray=IsImageGray(image,exception);
cristyb5d5f722009-11-04 03:03:49 +0000676#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristyb34ef052010-10-07 00:12:05 +0000677 #pragma omp parallel sections
cristy3ed852e2009-09-05 21:47:34 +0000678#endif
cristy3ed852e2009-09-05 21:47:34 +0000679 {
cristyb34ef052010-10-07 00:12:05 +0000680#if defined(MAGICKCORE_OPENMP_SUPPORT)
681 #pragma omp section
682#endif
cristy3ed852e2009-09-05 21:47:34 +0000683 {
cristyb34ef052010-10-07 00:12:05 +0000684 MagickBooleanType
685 thread_status;
686
687 if (is_gray != MagickFalse)
688 thread_status=ForwardFourierTransformChannel(image,
689 GrayChannels,modulus,fourier_image,exception);
690 else
691 thread_status=ForwardFourierTransformChannel(image,
692 RedChannel,modulus,fourier_image,exception);
693 if (thread_status == MagickFalse)
694 status=thread_status;
cristy3ed852e2009-09-05 21:47:34 +0000695 }
cristyb34ef052010-10-07 00:12:05 +0000696#if defined(MAGICKCORE_OPENMP_SUPPORT)
697 #pragma omp section
698#endif
699 {
700 MagickBooleanType
701 thread_status;
702
703 thread_status=MagickTrue;
704 if (is_gray == MagickFalse)
705 thread_status=ForwardFourierTransformChannel(image,
706 GreenChannel,modulus,fourier_image,exception);
707 if (thread_status == MagickFalse)
708 status=thread_status;
709 }
710#if defined(MAGICKCORE_OPENMP_SUPPORT)
711 #pragma omp section
712#endif
713 {
714 MagickBooleanType
715 thread_status;
716
717 thread_status=MagickTrue;
718 if (is_gray == MagickFalse)
719 thread_status=ForwardFourierTransformChannel(image,
720 BlueChannel,modulus,fourier_image,exception);
721 if (thread_status == MagickFalse)
722 status=thread_status;
723 }
724#if defined(MAGICKCORE_OPENMP_SUPPORT)
725 #pragma omp section
726#endif
727 {
728 MagickBooleanType
729 thread_status;
730
731 thread_status=MagickTrue;
cristy4c08aed2011-07-01 19:47:50 +0000732 if (image->colorspace == CMYKColorspace)
cristyb34ef052010-10-07 00:12:05 +0000733 thread_status=ForwardFourierTransformChannel(image,
cristy4c08aed2011-07-01 19:47:50 +0000734 BlackChannel,modulus,fourier_image,exception);
cristyb34ef052010-10-07 00:12:05 +0000735 if (thread_status == MagickFalse)
736 status=thread_status;
737 }
738#if defined(MAGICKCORE_OPENMP_SUPPORT)
739 #pragma omp section
740#endif
741 {
742 MagickBooleanType
743 thread_status;
744
745 thread_status=MagickTrue;
cristy4c08aed2011-07-01 19:47:50 +0000746 if (image->matte != MagickFalse)
cristyb34ef052010-10-07 00:12:05 +0000747 thread_status=ForwardFourierTransformChannel(image,
cristy4c08aed2011-07-01 19:47:50 +0000748 AlphaChannel,modulus,fourier_image,exception);
cristyb34ef052010-10-07 00:12:05 +0000749 if (thread_status == MagickFalse)
750 status=thread_status;
751 }
cristy3ed852e2009-09-05 21:47:34 +0000752 }
753 if (status == MagickFalse)
754 fourier_image=DestroyImageList(fourier_image);
755 fftw_cleanup();
756 }
757 }
758 }
759#endif
760 return(fourier_image);
761}
762
763/*
764%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
765% %
766% %
767% %
768% 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 %
769% %
770% %
771% %
772%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
773%
774% InverseFourierTransformImage() implements the inverse discrete Fourier
775% transform (DFT) of the image either as a magnitude / phase or real /
776% imaginary image pair.
777%
778% The format of the InverseFourierTransformImage method is:
779%
cristyc9550792009-11-13 20:05:42 +0000780% Image *InverseFourierTransformImage(const Image *magnitude_image,
781% const Image *phase_image,const MagickBooleanType modulus,
782% ExceptionInfo *exception)
cristy3ed852e2009-09-05 21:47:34 +0000783%
784% A description of each parameter follows:
785%
cristyc9550792009-11-13 20:05:42 +0000786% o magnitude_image: the magnitude or real image.
787%
788% o phase_image: the phase or imaginary image.
cristy3ed852e2009-09-05 21:47:34 +0000789%
790% o modulus: if true, return transform as a magnitude / phase pair
791% otherwise a real / imaginary image pair.
792%
793% o exception: return any errors or warnings in this structure.
794%
795*/
796
797#if defined(MAGICKCORE_FFTW_DELEGATE)
cristybb503372010-05-27 20:51:26 +0000798static MagickBooleanType InverseQuadrantSwap(const size_t width,
799 const size_t height,const double *source,double *destination)
cristy3ed852e2009-09-05 21:47:34 +0000800{
cristyc4ea4a42011-01-24 01:43:30 +0000801 register ssize_t
802 x;
803
cristybb503372010-05-27 20:51:26 +0000804 ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000805 center,
806 y;
807
cristy3ed852e2009-09-05 21:47:34 +0000808 /*
809 Swap quadrants.
810 */
cristybb503372010-05-27 20:51:26 +0000811 center=(ssize_t) floor((double) width/2.0)+1L;
812 for (y=1L; y < (ssize_t) height; y++)
813 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
cristy3ed852e2009-09-05 21:47:34 +0000814 destination[center*(height-y)-x+width/2L]=source[y*width+x];
cristybb503372010-05-27 20:51:26 +0000815 for (y=0L; y < (ssize_t) height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000816 destination[center*y]=source[y*width+width/2L];
817 for (x=0L; x < center; x++)
818 destination[x]=source[center-x-1L];
cristybb503372010-05-27 20:51:26 +0000819 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
cristy3ed852e2009-09-05 21:47:34 +0000820}
821
822static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
cristyc9550792009-11-13 20:05:42 +0000823 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
824 ExceptionInfo *exception)
cristy3ed852e2009-09-05 21:47:34 +0000825{
826 CacheView
827 *magnitude_view,
828 *phase_view;
829
830 double
831 *magnitude,
832 *phase,
833 *magnitude_source,
834 *phase_source;
835
cristy3ed852e2009-09-05 21:47:34 +0000836 MagickBooleanType
837 status;
838
cristy4c08aed2011-07-01 19:47:50 +0000839 register const Quantum
cristy3ed852e2009-09-05 21:47:34 +0000840 *p;
841
cristybb503372010-05-27 20:51:26 +0000842 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +0000843 i,
844 x;
845
cristyc4ea4a42011-01-24 01:43:30 +0000846 ssize_t
847 y;
848
cristy3ed852e2009-09-05 21:47:34 +0000849 /*
850 Inverse fourier - read image and break down into a double array.
851 */
cristy3ed852e2009-09-05 21:47:34 +0000852 magnitude_source=(double *) AcquireQuantumMemory((size_t)
853 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
854 if (magnitude_source == (double *) NULL)
855 {
856 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000857 ResourceLimitError,"MemoryAllocationFailed","`%s'",
858 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +0000859 return(MagickFalse);
860 }
861 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
cristyc4ea4a42011-01-24 01:43:30 +0000862 fourier_info->width*sizeof(*phase_source));
cristy3ed852e2009-09-05 21:47:34 +0000863 if (phase_source == (double *) NULL)
864 {
865 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000866 ResourceLimitError,"MemoryAllocationFailed","`%s'",
867 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +0000868 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
869 return(MagickFalse);
870 }
871 i=0L;
872 magnitude_view=AcquireCacheView(magnitude_image);
cristybb503372010-05-27 20:51:26 +0000873 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +0000874 {
875 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
876 exception);
cristy4c08aed2011-07-01 19:47:50 +0000877 if (p == (const Quantum *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000878 break;
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 {
cristy4c08aed2011-07-01 19:47:50 +0000886 magnitude_source[i]=QuantumScale*GetPixelRed(magnitude_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000887 break;
888 }
889 case GreenChannel:
890 {
cristy4c08aed2011-07-01 19:47:50 +0000891 magnitude_source[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000892 break;
893 }
894 case BlueChannel:
895 {
cristy4c08aed2011-07-01 19:47:50 +0000896 magnitude_source[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000897 break;
898 }
cristy4c08aed2011-07-01 19:47:50 +0000899 case BlackChannel:
cristy3ed852e2009-09-05 21:47:34 +0000900 {
cristy4c08aed2011-07-01 19:47:50 +0000901 magnitude_source[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000902 break;
903 }
cristy4c08aed2011-07-01 19:47:50 +0000904 case AlphaChannel:
cristy3ed852e2009-09-05 21:47:34 +0000905 {
cristy4c08aed2011-07-01 19:47:50 +0000906 magnitude_source[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000907 break;
908 }
909 case GrayChannels:
910 {
cristy4c08aed2011-07-01 19:47:50 +0000911 magnitude_source[i]=QuantumScale*GetPixelGray(magnitude_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000912 break;
913 }
914 }
915 i++;
cristyed231572011-07-14 02:18:59 +0000916 p+=GetPixelChannels(magnitude_image);
cristy3ed852e2009-09-05 21:47:34 +0000917 }
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);
cristy4c08aed2011-07-01 19:47:50 +0000925 if (p == (const Quantum *) NULL)
cristy3ed852e2009-09-05 21:47:34 +0000926 break;
cristybb503372010-05-27 20:51:26 +0000927 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000928 {
929 switch (fourier_info->channel)
930 {
931 case RedChannel:
932 default:
933 {
cristy4c08aed2011-07-01 19:47:50 +0000934 phase_source[i]=QuantumScale*GetPixelRed(phase_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000935 break;
936 }
937 case GreenChannel:
938 {
cristy4c08aed2011-07-01 19:47:50 +0000939 phase_source[i]=QuantumScale*GetPixelGreen(phase_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000940 break;
941 }
942 case BlueChannel:
943 {
cristy4c08aed2011-07-01 19:47:50 +0000944 phase_source[i]=QuantumScale*GetPixelBlue(phase_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000945 break;
946 }
cristy4c08aed2011-07-01 19:47:50 +0000947 case BlackChannel:
cristy3ed852e2009-09-05 21:47:34 +0000948 {
cristy4c08aed2011-07-01 19:47:50 +0000949 phase_source[i]=QuantumScale*GetPixelBlack(phase_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000950 break;
951 }
cristy4c08aed2011-07-01 19:47:50 +0000952 case AlphaChannel:
cristy3ed852e2009-09-05 21:47:34 +0000953 {
cristy4c08aed2011-07-01 19:47:50 +0000954 phase_source[i]=QuantumScale*GetPixelAlpha(phase_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000955 break;
956 }
957 case GrayChannels:
958 {
cristy4c08aed2011-07-01 19:47:50 +0000959 phase_source[i]=QuantumScale*GetPixelGray(phase_image,p);
cristy3ed852e2009-09-05 21:47:34 +0000960 break;
961 }
962 }
963 i++;
cristyed231572011-07-14 02:18:59 +0000964 p+=GetPixelChannels(phase_image);
cristy3ed852e2009-09-05 21:47:34 +0000965 }
966 }
967 if (fourier_info->modulus != MagickFalse)
968 {
969 i=0L;
cristybb503372010-05-27 20:51:26 +0000970 for (y=0L; y < (ssize_t) fourier_info->height; y++)
971 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +0000972 {
973 phase_source[i]-=0.5;
974 phase_source[i]*=(2.0*MagickPI);
975 i++;
976 }
977 }
978 magnitude_view=DestroyCacheView(magnitude_view);
979 phase_view=DestroyCacheView(phase_view);
980 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
981 fourier_info->center*sizeof(*magnitude));
982 if (magnitude == (double *) NULL)
983 {
984 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000985 ResourceLimitError,"MemoryAllocationFailed","`%s'",
986 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +0000987 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
988 phase_source=(double *) RelinquishMagickMemory(phase_source);
989 return(MagickFalse);
990 }
991 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
992 magnitude_source,magnitude);
993 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
cristyc4ea4a42011-01-24 01:43:30 +0000994 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
995 fourier_info->width*sizeof(*phase));
cristy3ed852e2009-09-05 21:47:34 +0000996 if (phase == (double *) NULL)
997 {
998 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +0000999 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1000 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001001 phase_source=(double *) RelinquishMagickMemory(phase_source);
1002 return(MagickFalse);
1003 }
1004 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1005 if (status != MagickFalse)
1006 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1007 phase_source,phase);
1008 phase_source=(double *) RelinquishMagickMemory(phase_source);
1009 /*
1010 Merge two sets.
1011 */
1012 i=0L;
1013 if (fourier_info->modulus != MagickFalse)
cristybb503372010-05-27 20:51:26 +00001014 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1015 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +00001016 {
cristy56ed31c2010-03-22 00:46:21 +00001017#if defined(MAGICKCORE_HAVE_COMPLEX_H)
cristy3ed852e2009-09-05 21:47:34 +00001018 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
cristy56ed31c2010-03-22 00:46:21 +00001019#else
1020 fourier[i][0]=magnitude[i]*cos(phase[i]);
1021 fourier[i][1]=magnitude[i]*sin(phase[i]);
1022#endif
cristy3ed852e2009-09-05 21:47:34 +00001023 i++;
1024 }
1025 else
cristybb503372010-05-27 20:51:26 +00001026 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1027 for (x=0L; x < (ssize_t) fourier_info->center; x++)
cristy3ed852e2009-09-05 21:47:34 +00001028 {
cristy56ed31c2010-03-22 00:46:21 +00001029#if defined(MAGICKCORE_HAVE_COMPLEX_H)
cristy3ed852e2009-09-05 21:47:34 +00001030 fourier[i]=magnitude[i]+I*phase[i];
cristy56ed31c2010-03-22 00:46:21 +00001031#else
1032 fourier[i][0]=magnitude[i];
1033 fourier[i][1]=phase[i];
1034#endif
cristy3ed852e2009-09-05 21:47:34 +00001035 i++;
1036 }
1037 phase=(double *) RelinquishMagickMemory(phase);
1038 magnitude=(double *) RelinquishMagickMemory(magnitude);
1039 return(status);
1040}
1041
1042static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1043 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1044{
1045 CacheView
1046 *image_view;
1047
1048 double
1049 *source;
1050
1051 fftw_plan
1052 fftw_c2r_plan;
1053
cristy4c08aed2011-07-01 19:47:50 +00001054 register Quantum
cristyc4ea4a42011-01-24 01:43:30 +00001055 *q;
1056
cristybb503372010-05-27 20:51:26 +00001057 register ssize_t
cristy3ed852e2009-09-05 21:47:34 +00001058 i,
1059 x;
1060
cristyc4ea4a42011-01-24 01:43:30 +00001061 ssize_t
1062 y;
cristy3ed852e2009-09-05 21:47:34 +00001063
cristyc4ea4a42011-01-24 01:43:30 +00001064 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1065 fourier_info->width*sizeof(*source));
cristy3ed852e2009-09-05 21:47:34 +00001066 if (source == (double *) NULL)
1067 {
1068 (void) ThrowMagickException(exception,GetMagickModule(),
1069 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1070 return(MagickFalse);
1071 }
cristyb5d5f722009-11-04 03:03:49 +00001072#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristy3ed852e2009-09-05 21:47:34 +00001073 #pragma omp critical (MagickCore_InverseFourierTransform)
1074#endif
cristydf079ac2010-09-10 01:45:44 +00001075 {
1076 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1077 fourier,source,FFTW_ESTIMATE);
1078 fftw_execute(fftw_c2r_plan);
1079 fftw_destroy_plan(fftw_c2r_plan);
1080 }
cristy3ed852e2009-09-05 21:47:34 +00001081 i=0L;
1082 image_view=AcquireCacheView(image);
cristybb503372010-05-27 20:51:26 +00001083 for (y=0L; y < (ssize_t) fourier_info->height; y++)
cristy3ed852e2009-09-05 21:47:34 +00001084 {
cristy85812052010-09-14 17:56:15 +00001085 if (y >= (ssize_t) image->rows)
1086 break;
1087 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1088 image->columns ? image->columns : fourier_info->width,1UL,exception);
cristy4c08aed2011-07-01 19:47:50 +00001089 if (q == (const Quantum *) NULL)
cristy3ed852e2009-09-05 21:47:34 +00001090 break;
cristybb503372010-05-27 20:51:26 +00001091 for (x=0L; x < (ssize_t) fourier_info->width; x++)
cristy3ed852e2009-09-05 21:47:34 +00001092 {
1093 switch (fourier_info->channel)
1094 {
1095 case RedChannel:
1096 default:
1097 {
cristy4c08aed2011-07-01 19:47:50 +00001098 SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +00001099 break;
1100 }
1101 case GreenChannel:
1102 {
cristy4c08aed2011-07-01 19:47:50 +00001103 SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +00001104 break;
1105 }
1106 case BlueChannel:
1107 {
cristy4c08aed2011-07-01 19:47:50 +00001108 SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +00001109 break;
1110 }
cristy4c08aed2011-07-01 19:47:50 +00001111 case BlackChannel:
cristy3ed852e2009-09-05 21:47:34 +00001112 {
cristy4c08aed2011-07-01 19:47:50 +00001113 SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +00001114 break;
1115 }
cristy4c08aed2011-07-01 19:47:50 +00001116 case AlphaChannel:
cristy3ed852e2009-09-05 21:47:34 +00001117 {
cristy4c08aed2011-07-01 19:47:50 +00001118 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +00001119 break;
1120 }
1121 case GrayChannels:
1122 {
cristy4c08aed2011-07-01 19:47:50 +00001123 SetPixelGray(image,ClampToQuantum(QuantumRange*source[i]),q);
cristy3ed852e2009-09-05 21:47:34 +00001124 break;
1125 }
1126 }
1127 i++;
cristyed231572011-07-14 02:18:59 +00001128 q+=GetPixelChannels(image);
cristy3ed852e2009-09-05 21:47:34 +00001129 }
1130 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1131 break;
1132 }
1133 image_view=DestroyCacheView(image_view);
1134 source=(double *) RelinquishMagickMemory(source);
1135 return(MagickTrue);
1136}
1137
cristyc9550792009-11-13 20:05:42 +00001138static MagickBooleanType InverseFourierTransformChannel(
1139 const Image *magnitude_image,const Image *phase_image,
cristy3ed852e2009-09-05 21:47:34 +00001140 const ChannelType channel,const MagickBooleanType modulus,
1141 Image *fourier_image,ExceptionInfo *exception)
1142{
1143 double
1144 *magnitude,
1145 *phase;
1146
1147 fftw_complex
1148 *fourier;
1149
1150 FourierInfo
1151 fourier_info;
1152
1153 MagickBooleanType
1154 status;
1155
1156 size_t
1157 extent;
1158
cristyc9550792009-11-13 20:05:42 +00001159 fourier_info.width=magnitude_image->columns;
1160 if ((magnitude_image->columns != magnitude_image->rows) ||
1161 ((magnitude_image->columns % 2) != 0) ||
1162 ((magnitude_image->rows % 2) != 0))
cristy3ed852e2009-09-05 21:47:34 +00001163 {
cristyc9550792009-11-13 20:05:42 +00001164 extent=magnitude_image->columns < magnitude_image->rows ?
1165 magnitude_image->rows : magnitude_image->columns;
cristy3ed852e2009-09-05 21:47:34 +00001166 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1167 }
1168 fourier_info.height=fourier_info.width;
cristybb503372010-05-27 20:51:26 +00001169 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
cristy3ed852e2009-09-05 21:47:34 +00001170 fourier_info.channel=channel;
1171 fourier_info.modulus=modulus;
1172 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
cristyc4ea4a42011-01-24 01:43:30 +00001173 fourier_info.center*sizeof(*magnitude));
cristy3ed852e2009-09-05 21:47:34 +00001174 if (magnitude == (double *) NULL)
1175 {
1176 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001177 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1178 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001179 return(MagickFalse);
1180 }
1181 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
cristyc4ea4a42011-01-24 01:43:30 +00001182 fourier_info.center*sizeof(*phase));
cristy3ed852e2009-09-05 21:47:34 +00001183 if (phase == (double *) NULL)
1184 {
1185 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001186 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1187 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001188 magnitude=(double *) RelinquishMagickMemory(magnitude);
1189 return(MagickFalse);
1190 }
cristyb41ee102010-10-04 16:46:15 +00001191 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
cristy3ed852e2009-09-05 21:47:34 +00001192 fourier_info.center*sizeof(*fourier));
1193 if (fourier == (fftw_complex *) NULL)
1194 {
1195 (void) ThrowMagickException(exception,GetMagickModule(),
cristyc9550792009-11-13 20:05:42 +00001196 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1197 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001198 phase=(double *) RelinquishMagickMemory(phase);
1199 magnitude=(double *) RelinquishMagickMemory(magnitude);
1200 return(MagickFalse);
1201 }
cristyc9550792009-11-13 20:05:42 +00001202 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1203 exception);
cristy3ed852e2009-09-05 21:47:34 +00001204 if (status != MagickFalse)
1205 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1206 exception);
cristyb41ee102010-10-04 16:46:15 +00001207 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
cristy3ed852e2009-09-05 21:47:34 +00001208 phase=(double *) RelinquishMagickMemory(phase);
1209 magnitude=(double *) RelinquishMagickMemory(magnitude);
1210 return(status);
1211}
1212#endif
1213
cristyc9550792009-11-13 20:05:42 +00001214MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1215 const Image *phase_image,const MagickBooleanType modulus,
1216 ExceptionInfo *exception)
cristy3ed852e2009-09-05 21:47:34 +00001217{
1218 Image
1219 *fourier_image;
1220
cristyc9550792009-11-13 20:05:42 +00001221 assert(magnitude_image != (Image *) NULL);
1222 assert(magnitude_image->signature == MagickSignature);
1223 if (magnitude_image->debug != MagickFalse)
1224 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1225 magnitude_image->filename);
1226 if (phase_image == (Image *) NULL)
1227 {
1228 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1229 "ImageSequenceRequired","`%s'",magnitude_image->filename);
cristy9372a152009-11-18 01:42:16 +00001230 return((Image *) NULL);
cristyc9550792009-11-13 20:05:42 +00001231 }
cristy3ed852e2009-09-05 21:47:34 +00001232#if !defined(MAGICKCORE_FFTW_DELEGATE)
1233 fourier_image=(Image *) NULL;
1234 (void) modulus;
1235 (void) ThrowMagickException(exception,GetMagickModule(),
1236 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
cristyc9550792009-11-13 20:05:42 +00001237 magnitude_image->filename);
cristy3ed852e2009-09-05 21:47:34 +00001238#else
1239 {
cristyc9550792009-11-13 20:05:42 +00001240 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1241 magnitude_image->rows,MagickFalse,exception);
cristy3ed852e2009-09-05 21:47:34 +00001242 if (fourier_image != (Image *) NULL)
1243 {
1244 MagickBooleanType
1245 is_gray,
1246 status;
1247
cristy3ed852e2009-09-05 21:47:34 +00001248 status=MagickTrue;
cristy4c08aed2011-07-01 19:47:50 +00001249 is_gray=IsImageGray(magnitude_image,exception);
cristyc9550792009-11-13 20:05:42 +00001250 if (is_gray != MagickFalse)
cristy4c08aed2011-07-01 19:47:50 +00001251 is_gray=IsImageGray(phase_image,exception);
cristyb5d5f722009-11-04 03:03:49 +00001252#if defined(MAGICKCORE_OPENMP_SUPPORT)
cristyb34ef052010-10-07 00:12:05 +00001253 #pragma omp parallel sections
cristy3ed852e2009-09-05 21:47:34 +00001254#endif
cristy3ed852e2009-09-05 21:47:34 +00001255 {
cristyb34ef052010-10-07 00:12:05 +00001256#if defined(MAGICKCORE_OPENMP_SUPPORT)
1257 #pragma omp section
1258#endif
cristy3ed852e2009-09-05 21:47:34 +00001259 {
cristyb34ef052010-10-07 00:12:05 +00001260 MagickBooleanType
1261 thread_status;
1262
1263 if (is_gray != MagickFalse)
1264 thread_status=InverseFourierTransformChannel(magnitude_image,
1265 phase_image,GrayChannels,modulus,fourier_image,exception);
1266 else
cristyc9550792009-11-13 20:05:42 +00001267 thread_status=InverseFourierTransformChannel(magnitude_image,
1268 phase_image,RedChannel,modulus,fourier_image,exception);
cristyb34ef052010-10-07 00:12:05 +00001269 if (thread_status == MagickFalse)
1270 status=thread_status;
cristy3ed852e2009-09-05 21:47:34 +00001271 }
cristyb34ef052010-10-07 00:12:05 +00001272#if defined(MAGICKCORE_OPENMP_SUPPORT)
1273 #pragma omp section
1274#endif
1275 {
1276 MagickBooleanType
1277 thread_status;
1278
1279 thread_status=MagickTrue;
1280 if (is_gray == MagickFalse)
1281 thread_status=InverseFourierTransformChannel(magnitude_image,
1282 phase_image,GreenChannel,modulus,fourier_image,exception);
1283 if (thread_status == MagickFalse)
1284 status=thread_status;
1285 }
1286#if defined(MAGICKCORE_OPENMP_SUPPORT)
1287 #pragma omp section
1288#endif
1289 {
1290 MagickBooleanType
1291 thread_status;
1292
1293 thread_status=MagickTrue;
1294 if (is_gray == MagickFalse)
1295 thread_status=InverseFourierTransformChannel(magnitude_image,
1296 phase_image,BlueChannel,modulus,fourier_image,exception);
1297 if (thread_status == MagickFalse)
1298 status=thread_status;
1299 }
1300#if defined(MAGICKCORE_OPENMP_SUPPORT)
1301 #pragma omp section
1302#endif
1303 {
1304 MagickBooleanType
1305 thread_status;
1306
1307 thread_status=MagickTrue;
cristy4c08aed2011-07-01 19:47:50 +00001308 if (magnitude_image->colorspace == CMYKColorspace)
cristyb34ef052010-10-07 00:12:05 +00001309 thread_status=InverseFourierTransformChannel(magnitude_image,
cristy4c08aed2011-07-01 19:47:50 +00001310 phase_image,BlackChannel,modulus,fourier_image,exception);
cristyb34ef052010-10-07 00:12:05 +00001311 if (thread_status == MagickFalse)
1312 status=thread_status;
1313 }
1314#if defined(MAGICKCORE_OPENMP_SUPPORT)
1315 #pragma omp section
1316#endif
1317 {
1318 MagickBooleanType
1319 thread_status;
1320
1321 thread_status=MagickTrue;
cristy4c08aed2011-07-01 19:47:50 +00001322 if (magnitude_image->matte != MagickFalse)
cristyb34ef052010-10-07 00:12:05 +00001323 thread_status=InverseFourierTransformChannel(magnitude_image,
cristy4c08aed2011-07-01 19:47:50 +00001324 phase_image,AlphaChannel,modulus,fourier_image,exception);
cristyb34ef052010-10-07 00:12:05 +00001325 if (thread_status == MagickFalse)
1326 status=thread_status;
1327 }
cristy3ed852e2009-09-05 21:47:34 +00001328 }
1329 if (status == MagickFalse)
1330 fourier_image=DestroyImage(fourier_image);
1331 }
cristyb34ef052010-10-07 00:12:05 +00001332 fftw_cleanup();
cristy3ed852e2009-09-05 21:47:34 +00001333 }
1334#endif
1335 return(fourier_image);
1336}