MagickCore  6.9.12-94
Convert, Edit, Or Compose Bitmap Images
 All Data Structures
fourier.c
1 /*
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 % Cristy %
19 % July 2009 %
20 % %
21 % %
22 % Copyright 1999 ImageMagick Studio LLC, a non-profit organization %
23 % 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 % https://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/artifact.h"
47 #include "magick/attribute.h"
48 #include "magick/blob.h"
49 #include "magick/cache.h"
50 #include "magick/image.h"
51 #include "magick/image-private.h"
52 #include "magick/list.h"
53 #include "magick/fourier.h"
54 #include "magick/log.h"
55 #include "magick/memory_.h"
56 #include "magick/monitor.h"
57 #include "magick/monitor-private.h"
58 #include "magick/pixel-accessor.h"
59 #include "magick/pixel-private.h"
60 #include "magick/property.h"
61 #include "magick/quantum-private.h"
62 #include "magick/resource_.h"
63 #include "magick/string-private.h"
64 #include "magick/thread-private.h"
65 #if defined(MAGICKCORE_FFTW_DELEGATE) && !defined(__cplusplus) && !defined(c_plusplus)
66 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
67 #include <complex.h>
68 #endif
69 #include <fftw3.h>
70 #if !defined(MAGICKCORE_HAVE_CABS)
71 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
72 #endif
73 #if !defined(MAGICKCORE_HAVE_CARG)
74 #define carg(z) (atan2(cimag(z),creal(z)))
75 #endif
76 #if !defined(MAGICKCORE_HAVE_CIMAG)
77 #define cimag(z) (z[1])
78 #endif
79 #if !defined(MAGICKCORE_HAVE_CREAL)
80 #define creal(z) (z[0])
81 #endif
82 #endif
83 
84 /*
85  Typedef declarations.
86 */
87 typedef struct _FourierInfo
88 {
89  ChannelType
90  channel;
91 
92  MagickBooleanType
93  modulus;
94 
95  size_t
96  width,
97  height;
98 
99  ssize_t
100  center;
101 } FourierInfo;
102 
103 /*
104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 % %
106 % %
107 % %
108 % C o m p l e x I m a g e s %
109 % %
110 % %
111 % %
112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113 %
114 % ComplexImages() performs complex mathematics on an image sequence.
115 %
116 % The format of the ComplexImages method is:
117 %
118 % MagickBooleanType ComplexImages(Image *images,const ComplexOperator op,
119 % ExceptionInfo *exception)
120 %
121 % A description of each parameter follows:
122 %
123 % o image: the image.
124 %
125 % o op: A complex operator.
126 %
127 % o exception: return any errors or warnings in this structure.
128 %
129 */
130 MagickExport Image *ComplexImages(const Image *images,const ComplexOperator op,
131  ExceptionInfo *exception)
132 {
133 #define ComplexImageTag "Complex/Image"
134 
135  CacheView
136  *Ai_view,
137  *Ar_view,
138  *Bi_view,
139  *Br_view,
140  *Ci_view,
141  *Cr_view;
142 
143  const char
144  *artifact;
145 
146  const Image
147  *Ai_image,
148  *Ar_image,
149  *Bi_image,
150  *Br_image;
151 
152  double
153  snr;
154 
155  Image
156  *Ci_image,
157  *complex_images,
158  *Cr_image,
159  *image;
160 
161  MagickBooleanType
162  status;
163 
164  MagickOffsetType
165  progress;
166 
167  size_t
168  columns,
169  rows;
170 
171  ssize_t
172  y;
173 
174  assert(images != (Image *) NULL);
175  assert(images->signature == MagickCoreSignature);
176  assert(exception != (ExceptionInfo *) NULL);
177  assert(exception->signature == MagickCoreSignature);
178  if (IsEventLogging() != MagickFalse)
179  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
180  if (images->next == (Image *) NULL)
181  {
182  (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
183  "ImageSequenceRequired","`%s'",images->filename);
184  return((Image *) NULL);
185  }
186  image=CloneImage(images,0,0,MagickTrue,exception);
187  if (image == (Image *) NULL)
188  return((Image *) NULL);
189  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
190  {
191  image=DestroyImageList(image);
192  return(image);
193  }
194  image->depth=32UL;
195  complex_images=NewImageList();
196  AppendImageToList(&complex_images,image);
197  image=CloneImage(images->next,0,0,MagickTrue,exception);
198  if (image == (Image *) NULL)
199  {
200  complex_images=DestroyImageList(complex_images);
201  return(complex_images);
202  }
203  if (SetImageStorageClass(image,DirectClass) == MagickFalse)
204  {
205  image=DestroyImageList(image);
206  return(image);
207  }
208  image->depth=32UL;
209  AppendImageToList(&complex_images,image);
210  /*
211  Apply complex mathematics to image pixels.
212  */
213  artifact=GetImageArtifact(image,"complex:snr");
214  snr=0.0;
215  if (artifact != (const char *) NULL)
216  snr=StringToDouble(artifact,(char **) NULL);
217  Ar_image=images;
218  Ai_image=images->next;
219  Br_image=images;
220  Bi_image=images->next;
221  if ((images->next->next != (Image *) NULL) &&
222  (images->next->next->next != (Image *) NULL))
223  {
224  Br_image=images->next->next;
225  Bi_image=images->next->next->next;
226  }
227  Cr_image=complex_images;
228  Ci_image=complex_images->next;
229  Ar_view=AcquireVirtualCacheView(Ar_image,exception);
230  Ai_view=AcquireVirtualCacheView(Ai_image,exception);
231  Br_view=AcquireVirtualCacheView(Br_image,exception);
232  Bi_view=AcquireVirtualCacheView(Bi_image,exception);
233  Cr_view=AcquireAuthenticCacheView(Cr_image,exception);
234  Ci_view=AcquireAuthenticCacheView(Ci_image,exception);
235  status=MagickTrue;
236  progress=0;
237  columns=MagickMin(Cr_image->columns,Ci_image->columns);
238  rows=MagickMin(Cr_image->rows,Ci_image->rows);
239 #if defined(MAGICKCORE_OPENMP_SUPPORT)
240  #pragma omp parallel for schedule(static) shared(progress,status) \
241  magick_number_threads(Cr_image,complex_images,rows,1L)
242 #endif
243  for (y=0; y < (ssize_t) rows; y++)
244  {
245  const PixelPacket
246  *magick_restrict Ai,
247  *magick_restrict Ar,
248  *magick_restrict Bi,
249  *magick_restrict Br;
250 
252  *magick_restrict Ci,
253  *magick_restrict Cr;
254 
255  ssize_t
256  x;
257 
258  if (status == MagickFalse)
259  continue;
260  Ar=GetCacheViewVirtualPixels(Ar_view,0,y,columns,1,exception);
261  Ai=GetCacheViewVirtualPixels(Ai_view,0,y,columns,1,exception);
262  Br=GetCacheViewVirtualPixels(Br_view,0,y,columns,1,exception);
263  Bi=GetCacheViewVirtualPixels(Bi_view,0,y,columns,1,exception);
264  Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,columns,1,exception);
265  Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,columns,1,exception);
266  if ((Ar == (const PixelPacket *) NULL) ||
267  (Ai == (const PixelPacket *) NULL) ||
268  (Br == (const PixelPacket *) NULL) ||
269  (Bi == (const PixelPacket *) NULL) ||
270  (Cr == (PixelPacket *) NULL) || (Ci == (PixelPacket *) NULL))
271  {
272  status=MagickFalse;
273  continue;
274  }
275  for (x=0; x < (ssize_t) columns; x++)
276  {
278  ai = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Ai)),
279  (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Ai)),
280  (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Ai)),
281  (Quantum) (image->matte != MagickFalse ? QuantumScale*
282  (MagickRealType) GetPixelOpacity(Ai) : (MagickRealType)
283  OpaqueOpacity), 0 },
284  ar = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Ar)),
285  (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Ar)),
286  (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Ar)),
287  (Quantum) (image->matte != MagickFalse ? QuantumScale*
288  (MagickRealType) GetPixelOpacity(Ar) : (MagickRealType)
289  OpaqueOpacity), 0 },
290  bi = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Bi)),
291  (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Bi)),
292  (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Bi)),
293  (Quantum) (image->matte != MagickFalse ? QuantumScale*
294  (MagickRealType) GetPixelOpacity(Bi) : (MagickRealType)
295  OpaqueOpacity), 0 },
296  br = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Br)),
297  (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Br)),
298  (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Br)),
299  (Quantum) (image->matte != MagickFalse ? QuantumScale*
300  (MagickRealType) GetPixelOpacity(Br) : (MagickRealType)
301  OpaqueOpacity), 0 },
302  ci,
303  cr;
304 
305  switch (op)
306  {
307  case AddComplexOperator:
308  {
309  cr.red=ar.red+br.red;
310  ci.red=ai.red+bi.red;
311  cr.green=ar.green+br.green;
312  ci.green=ai.green+bi.green;
313  cr.blue=ar.blue+br.blue;
314  ci.blue=ai.blue+bi.blue;
315  cr.opacity=ar.opacity+br.opacity;
316  ci.opacity=ai.opacity+bi.opacity;
317  break;
318  }
319  case ConjugateComplexOperator:
320  default:
321  {
322  cr.red=ar.red;
323  ci.red=(-ai.red);
324  cr.green=ar.green;
325  ci.green=(-ai.green);
326  cr.blue=ar.blue;
327  ci.blue=(-ai.blue);
328  cr.opacity=ar.opacity;
329  ci.opacity=(-ai.opacity);
330  break;
331  }
332  case DivideComplexOperator:
333  {
334  double
335  gamma;
336 
337  gamma=PerceptibleReciprocal((double) br.red*(double) br.red+
338  (double) bi.red*(double) bi.red+snr);
339  cr.red=gamma*((double) ar.red*(double) br.red+(double) ai.red*
340  (double) bi.red);
341  ci.red=gamma*((double) ai.red*(double) br.red-(double) ar.red*
342  (double) bi.red);
343  gamma=PerceptibleReciprocal((double) br.green*(double) br.green+
344  (double) bi.green*(double) bi.green+snr);
345  cr.green=gamma*((double) ar.green*(double) br.green+
346  (double) ai.green*(double) bi.green);
347  ci.green=gamma*((double) ai.green*(double) br.green-
348  (double) ar.green*(double) bi.green);
349  gamma=PerceptibleReciprocal((double) br.blue*(double) br.blue+
350  (double) bi.blue*(double) bi.blue+snr);
351  cr.blue=gamma*((double) ar.blue*(double) br.blue+
352  (double) ai.blue*(double) bi.blue);
353  ci.blue=gamma*((double) ai.blue*(double) br.blue-
354  (double) ar.blue*(double) bi.blue);
355  gamma=PerceptibleReciprocal((double) br.opacity*(double) br.opacity+
356  (double) bi.opacity*(double) bi.opacity+snr);
357  cr.opacity=gamma*((double) ar.opacity*(double) br.opacity+
358  (double) ai.opacity*(double) bi.opacity);
359  ci.opacity=gamma*((double) ai.opacity*(double) br.opacity-
360  (double) ar.opacity*(double) bi.opacity);
361  break;
362  }
363  case MagnitudePhaseComplexOperator:
364  {
365  cr.red=sqrt((double) ar.red*(double) ar.red+(double) ai.red*
366  (double) ai.red);
367  ci.red=atan2((double) ai.red,(double) ar.red)/(2.0*MagickPI)+0.5;
368  cr.green=sqrt((double) ar.green*(double) ar.green+(double) ai.green*
369  (double) ai.green);
370  ci.green=atan2((double) ai.green,(double) ar.green)/(2.0*MagickPI)+
371  0.5;
372  cr.blue=sqrt((double) ar.blue*(double) ar.blue+(double) ai.blue*
373  (double) ai.blue);
374  ci.blue=atan2((double) ai.blue,(double) ar.blue)/(2.0*MagickPI)+0.5;
375  cr.opacity=sqrt((double) ar.opacity*(double) ar.opacity+
376  (double) ai.opacity*(double) ai.opacity);
377  ci.opacity=atan2((double) ai.opacity,(double) ar.opacity)/
378  (2.0*MagickPI)+0.5;
379  break;
380  }
381  case MultiplyComplexOperator:
382  {
383  cr.red=((double) ar.red*(double) br.red-(double) ai.red*
384  (double) bi.red);
385  ci.red=((double) ai.red*(double) br.red+(double) ar.red*
386  (double) bi.red);
387  cr.green=((double) ar.green*(double) br.green-(double) ai.green*
388  (double) bi.green);
389  ci.green=((double) ai.green*(double) br.green+(double) ar.green*
390  (double) bi.green);
391  cr.blue=((double) ar.blue*(double) br.blue-(double) ai.blue*
392  (double) bi.blue);
393  ci.blue=((double) ai.blue*(double) br.blue+(double) ar.blue*
394  (double) bi.blue);
395  cr.opacity=((double) ar.opacity*(double) br.opacity-
396  (double) ai.opacity*(double) bi.opacity);
397  ci.opacity=((double) ai.opacity*(double) br.opacity+
398  (double) ar.opacity*(double) bi.opacity);
399  break;
400  }
401  case RealImaginaryComplexOperator:
402  {
403  cr.red=(double) ar.red*cos(2.0*MagickPI*((double) ai.red-0.5));
404  ci.red=(double) ar.red*sin(2.0*MagickPI*((double) ai.red-0.5));
405  cr.green=(double) ar.green*cos(2.0*MagickPI*((double) ai.green-0.5));
406  ci.green=(double) ar.green*sin(2.0*MagickPI*((double) ai.green-0.5));
407  cr.blue=(double) ar.blue*cos(2.0*MagickPI*((double) ai.blue-0.5));
408  ci.blue=(double) ar.blue*sin(2.0*MagickPI*((double) ai.blue-0.5));
409  cr.opacity=(double) ar.opacity*cos(2.0*MagickPI*((double) ai.opacity-
410  0.5));
411  ci.opacity=(double) ar.opacity*sin(2.0*MagickPI*((double) ai.opacity-
412  0.5));
413  break;
414  }
415  case SubtractComplexOperator:
416  {
417  cr.red=ar.red-br.red;
418  ci.red=ai.red-bi.red;
419  cr.green=ar.green-br.green;
420  ci.green=ai.green-bi.green;
421  cr.blue=ar.blue-br.blue;
422  ci.blue=ai.blue-bi.blue;
423  cr.opacity=ar.opacity-br.opacity;
424  ci.opacity=ai.opacity-bi.opacity;
425  break;
426  }
427  }
428  Cr->red=(MagickRealType) QuantumRange*(MagickRealType) cr.red;
429  Ci->red=(MagickRealType) QuantumRange*(MagickRealType) ci.red;
430  Cr->green=(MagickRealType) QuantumRange*(MagickRealType) cr.green;
431  Ci->green=(MagickRealType) QuantumRange*(MagickRealType) ci.green;
432  Cr->blue=(MagickRealType) QuantumRange*(MagickRealType) cr.blue;
433  Ci->blue=(MagickRealType) QuantumRange*(MagickRealType) ci.blue;
434  if (images->matte != MagickFalse)
435  {
436  Cr->opacity=(MagickRealType) QuantumRange*(MagickRealType) cr.opacity;
437  Ci->opacity=(MagickRealType) QuantumRange*(MagickRealType) ci.opacity;
438  }
439  Ar++;
440  Ai++;
441  Br++;
442  Bi++;
443  Cr++;
444  Ci++;
445  }
446  if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
447  status=MagickFalse;
448  if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
449  status=MagickFalse;
450  if (images->progress_monitor != (MagickProgressMonitor) NULL)
451  {
452  MagickBooleanType
453  proceed;
454 
455 #if defined(MAGICKCORE_OPENMP_SUPPORT)
456  #pragma omp atomic
457 #endif
458  progress++;
459  proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
460  if (proceed == MagickFalse)
461  status=MagickFalse;
462  }
463  }
464  Cr_view=DestroyCacheView(Cr_view);
465  Ci_view=DestroyCacheView(Ci_view);
466  Br_view=DestroyCacheView(Br_view);
467  Bi_view=DestroyCacheView(Bi_view);
468  Ar_view=DestroyCacheView(Ar_view);
469  Ai_view=DestroyCacheView(Ai_view);
470  if (status == MagickFalse)
471  complex_images=DestroyImageList(complex_images);
472  return(complex_images);
473 }
474 
475 /*
476 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
477 % %
478 % %
479 % %
480 % 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 %
481 % %
482 % %
483 % %
484 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
485 %
486 % ForwardFourierTransformImage() implements the discrete Fourier transform
487 % (DFT) of the image either as a magnitude / phase or real / imaginary image
488 % pair.
489 %
490 % The format of the ForwardFourierTransformImage method is:
491 %
492 % Image *ForwardFourierTransformImage(const Image *image,
493 % const MagickBooleanType modulus,ExceptionInfo *exception)
494 %
495 % A description of each parameter follows:
496 %
497 % o image: the image.
498 %
499 % o modulus: if true, return as transform as a magnitude / phase pair
500 % otherwise a real / imaginary image pair.
501 %
502 % o exception: return any errors or warnings in this structure.
503 %
504 */
505 
506 #if defined(MAGICKCORE_FFTW_DELEGATE) && !defined(__cplusplus) && !defined(c_plusplus)
507 
508 static MagickBooleanType RollFourier(const size_t width,const size_t height,
509  const ssize_t x_offset,const ssize_t y_offset,double *roll_pixels)
510 {
511  double
512  *source_pixels;
513 
514  MemoryInfo
515  *source_info;
516 
517  ssize_t
518  i,
519  u,
520  v,
521  x,
522  y;
523 
524  /*
525  Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
526  */
527  source_info=AcquireVirtualMemory(width,height*sizeof(*source_pixels));
528  if (source_info == (MemoryInfo *) NULL)
529  return(MagickFalse);
530  source_pixels=(double *) GetVirtualMemoryBlob(source_info);
531  i=0L;
532  for (y=0L; y < (ssize_t) height; y++)
533  {
534  if (y_offset < 0L)
535  v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
536  else
537  v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
538  y+y_offset;
539  for (x=0L; x < (ssize_t) width; x++)
540  {
541  if (x_offset < 0L)
542  u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
543  else
544  u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
545  x+x_offset;
546  source_pixels[v*width+u]=roll_pixels[i++];
547  }
548  }
549  (void) memcpy(roll_pixels,source_pixels,height*width*sizeof(*source_pixels));
550  source_info=RelinquishVirtualMemory(source_info);
551  return(MagickTrue);
552 }
553 
554 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
555  const size_t height,double *source_pixels,double *forward_pixels)
556 {
557  MagickBooleanType
558  status;
559 
560  ssize_t
561  center,
562  x,
563  y;
564 
565  /*
566  Swap quadrants.
567  */
568  center=(ssize_t) (width/2L)+1L;
569  status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,
570  source_pixels);
571  if (status == MagickFalse)
572  return(MagickFalse);
573  for (y=0L; y < (ssize_t) height; y++)
574  for (x=0L; x < (ssize_t) (width/2L); x++)
575  forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
576  for (y=1; y < (ssize_t) height; y++)
577  for (x=0L; x < (ssize_t) (width/2L); x++)
578  forward_pixels[(height-y)*width+width/2L-x-1L]=
579  source_pixels[y*center+x+1L];
580  for (x=0L; x < (ssize_t) (width/2L); x++)
581  forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
582  return(MagickTrue);
583 }
584 
585 static void CorrectPhaseLHS(const size_t width,const size_t height,
586  double *fourier_pixels)
587 {
588  ssize_t
589  x,
590  y;
591 
592  for (y=0L; y < (ssize_t) height; y++)
593  for (x=0L; x < (ssize_t) (width/2L); x++)
594  fourier_pixels[y*width+x]*=(-1.0);
595 }
596 
597 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
598  Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
599 {
600  CacheView
601  *magnitude_view,
602  *phase_view;
603 
604  double
605  *magnitude_pixels,
606  *phase_pixels;
607 
608  Image
609  *magnitude_image,
610  *phase_image;
611 
612  MagickBooleanType
613  status;
614 
615  MemoryInfo
616  *magnitude_info,
617  *phase_info;
618 
619  IndexPacket
620  *indexes;
621 
623  *q;
624 
625  ssize_t
626  i,
627  x,
628  y;
629 
630  magnitude_image=GetFirstImageInList(image);
631  phase_image=GetNextImageInList(image);
632  if (phase_image == (Image *) NULL)
633  {
634  (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
635  "ImageSequenceRequired","`%s'",image->filename);
636  return(MagickFalse);
637  }
638  /*
639  Create "Fourier Transform" image from constituent arrays.
640  */
641  magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
642  sizeof(*magnitude_pixels));
643  phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
644  sizeof(*phase_pixels));
645  if ((magnitude_info == (MemoryInfo *) NULL) ||
646  (phase_info == (MemoryInfo *) NULL))
647  {
648  if (phase_info != (MemoryInfo *) NULL)
649  phase_info=RelinquishVirtualMemory(phase_info);
650  if (magnitude_info != (MemoryInfo *) NULL)
651  magnitude_info=RelinquishVirtualMemory(magnitude_info);
652  (void) ThrowMagickException(exception,GetMagickModule(),
653  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
654  return(MagickFalse);
655  }
656  magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
657  (void) memset(magnitude_pixels,0,fourier_info->width*
658  fourier_info->height*sizeof(*magnitude_pixels));
659  phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
660  (void) memset(phase_pixels,0,fourier_info->width*
661  fourier_info->height*sizeof(*phase_pixels));
662  status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
663  magnitude,magnitude_pixels);
664  if (status != MagickFalse)
665  status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
666  phase_pixels);
667  CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
668  if (fourier_info->modulus != MagickFalse)
669  {
670  i=0L;
671  for (y=0L; y < (ssize_t) fourier_info->height; y++)
672  for (x=0L; x < (ssize_t) fourier_info->width; x++)
673  {
674  phase_pixels[i]/=(2.0*MagickPI);
675  phase_pixels[i]+=0.5;
676  i++;
677  }
678  }
679  magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
680  i=0L;
681  for (y=0L; y < (ssize_t) fourier_info->height; y++)
682  {
683  q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
684  exception);
685  if (q == (PixelPacket *) NULL)
686  break;
687  indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
688  for (x=0L; x < (ssize_t) fourier_info->width; x++)
689  {
690  switch (fourier_info->channel)
691  {
692  case RedChannel:
693  default:
694  {
695  SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
696  break;
697  }
698  case GreenChannel:
699  {
700  SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
701  break;
702  }
703  case BlueChannel:
704  {
705  SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
706  break;
707  }
708  case OpacityChannel:
709  {
710  SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
711  break;
712  }
713  case IndexChannel:
714  {
715  SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType) QuantumRange*
716  magnitude_pixels[i]));
717  break;
718  }
719  case GrayChannels:
720  {
721  SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
722  break;
723  }
724  }
725  i++;
726  q++;
727  }
728  status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
729  if (status == MagickFalse)
730  break;
731  }
732  magnitude_view=DestroyCacheView(magnitude_view);
733  i=0L;
734  phase_view=AcquireAuthenticCacheView(phase_image,exception);
735  for (y=0L; y < (ssize_t) fourier_info->height; y++)
736  {
737  q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
738  exception);
739  if (q == (PixelPacket *) NULL)
740  break;
741  indexes=GetCacheViewAuthenticIndexQueue(phase_view);
742  for (x=0L; x < (ssize_t) fourier_info->width; x++)
743  {
744  switch (fourier_info->channel)
745  {
746  case RedChannel:
747  default:
748  {
749  SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
750  break;
751  }
752  case GreenChannel:
753  {
754  SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
755  break;
756  }
757  case BlueChannel:
758  {
759  SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
760  break;
761  }
762  case OpacityChannel:
763  {
764  SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
765  break;
766  }
767  case IndexChannel:
768  {
769  SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
770  break;
771  }
772  case GrayChannels:
773  {
774  SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
775  break;
776  }
777  }
778  i++;
779  q++;
780  }
781  status=SyncCacheViewAuthenticPixels(phase_view,exception);
782  if (status == MagickFalse)
783  break;
784  }
785  phase_view=DestroyCacheView(phase_view);
786  phase_info=RelinquishVirtualMemory(phase_info);
787  magnitude_info=RelinquishVirtualMemory(magnitude_info);
788  return(status);
789 }
790 
791 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
792  const Image *image,double *magnitude_pixels,double *phase_pixels,
793  ExceptionInfo *exception)
794 {
795  CacheView
796  *image_view;
797 
798  const char
799  *value;
800 
801  double
802  *source_pixels;
803 
804  fftw_complex
805  *forward_pixels;
806 
807 
808  fftw_plan
809  fftw_r2c_plan;
810 
811  MemoryInfo
812  *forward_info,
813  *source_info;
814 
815  const IndexPacket
816  *indexes;
817 
818  const PixelPacket
819  *p;
820 
821  ssize_t
822  i,
823  x,
824  y;
825 
826  /*
827  Generate the forward Fourier transform.
828  */
829  if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
830  {
831  (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
832  "WidthOrHeightExceedsLimit","`%s'",image->filename);
833  return(MagickFalse);
834  }
835  source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
836  sizeof(*source_pixels));
837  if (source_info == (MemoryInfo *) NULL)
838  {
839  (void) ThrowMagickException(exception,GetMagickModule(),
840  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
841  return(MagickFalse);
842  }
843  source_pixels=(double *) GetVirtualMemoryBlob(source_info);
844  memset(source_pixels,0,fourier_info->width*fourier_info->height*
845  sizeof(*source_pixels));
846  i=0L;
847  image_view=AcquireVirtualCacheView(image,exception);
848  for (y=0L; y < (ssize_t) fourier_info->height; y++)
849  {
850  p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
851  exception);
852  if (p == (const PixelPacket *) NULL)
853  break;
854  indexes=GetCacheViewVirtualIndexQueue(image_view);
855  for (x=0L; x < (ssize_t) fourier_info->width; x++)
856  {
857  switch (fourier_info->channel)
858  {
859  case RedChannel:
860  default:
861  {
862  source_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
863  break;
864  }
865  case GreenChannel:
866  {
867  source_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
868  break;
869  }
870  case BlueChannel:
871  {
872  source_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
873  break;
874  }
875  case OpacityChannel:
876  {
877  source_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
878  break;
879  }
880  case IndexChannel:
881  {
882  source_pixels[i]=QuantumScale*(MagickRealType)
883  GetPixelIndex(indexes+x);
884  break;
885  }
886  case GrayChannels:
887  {
888  source_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
889  break;
890  }
891  }
892  i++;
893  p++;
894  }
895  }
896  image_view=DestroyCacheView(image_view);
897  forward_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/2+
898  1)*sizeof(*forward_pixels));
899  if (forward_info == (MemoryInfo *) NULL)
900  {
901  (void) ThrowMagickException(exception,GetMagickModule(),
902  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
903  source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
904  return(MagickFalse);
905  }
906  forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
907 #if defined(MAGICKCORE_OPENMP_SUPPORT)
908  #pragma omp critical (MagickCore_ForwardFourierTransform)
909 #endif
910  fftw_r2c_plan=fftw_plan_dft_r2c_2d((int) fourier_info->width,
911  (int) fourier_info->height,source_pixels,forward_pixels,FFTW_ESTIMATE);
912  fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
913  fftw_destroy_plan(fftw_r2c_plan);
914  source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
915  value=GetImageArtifact(image,"fourier:normalize");
916  if ((value == (const char *) NULL) || (LocaleCompare(value,"forward") == 0))
917  {
918  double
919  gamma;
920 
921  /*
922  Normalize forward transform.
923  */
924  i=0L;
925  gamma=PerceptibleReciprocal((double) fourier_info->width*
926  fourier_info->height);
927  for (y=0L; y < (ssize_t) fourier_info->height; y++)
928  for (x=0L; x < (ssize_t) fourier_info->center; x++)
929  {
930 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
931  forward_pixels[i]*=gamma;
932 #else
933  forward_pixels[i][0]*=gamma;
934  forward_pixels[i][1]*=gamma;
935 #endif
936  i++;
937  }
938  }
939  /*
940  Generate magnitude and phase (or real and imaginary).
941  */
942  i=0L;
943  if (fourier_info->modulus != MagickFalse)
944  for (y=0L; y < (ssize_t) fourier_info->height; y++)
945  for (x=0L; x < (ssize_t) fourier_info->center; x++)
946  {
947  magnitude_pixels[i]=cabs(forward_pixels[i]);
948  phase_pixels[i]=carg(forward_pixels[i]);
949  i++;
950  }
951  else
952  for (y=0L; y < (ssize_t) fourier_info->height; y++)
953  for (x=0L; x < (ssize_t) fourier_info->center; x++)
954  {
955  magnitude_pixels[i]=creal(forward_pixels[i]);
956  phase_pixels[i]=cimag(forward_pixels[i]);
957  i++;
958  }
959  forward_info=(MemoryInfo *) RelinquishVirtualMemory(forward_info);
960  return(MagickTrue);
961 }
962 
963 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
964  const ChannelType channel,const MagickBooleanType modulus,
965  Image *fourier_image,ExceptionInfo *exception)
966 {
967  double
968  *magnitude_pixels,
969  *phase_pixels;
970 
972  fourier_info;
973 
974  MagickBooleanType
975  status;
976 
977  MemoryInfo
978  *magnitude_info,
979  *phase_info;
980 
981  fourier_info.width=image->columns;
982  fourier_info.height=image->rows;
983  if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
984  ((image->rows % 2) != 0))
985  {
986  size_t extent=image->columns < image->rows ? image->rows : image->columns;
987  fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
988  }
989  fourier_info.height=fourier_info.width;
990  fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
991  fourier_info.channel=channel;
992  fourier_info.modulus=modulus;
993  magnitude_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
994  1)*sizeof(*magnitude_pixels));
995  phase_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+1)*
996  sizeof(*phase_pixels));
997  if ((magnitude_info == (MemoryInfo *) NULL) ||
998  (phase_info == (MemoryInfo *) NULL))
999  {
1000  if (phase_info != (MemoryInfo *) NULL)
1001  phase_info=RelinquishVirtualMemory(phase_info);
1002  if (magnitude_info == (MemoryInfo *) NULL)
1003  magnitude_info=RelinquishVirtualMemory(magnitude_info);
1004  (void) ThrowMagickException(exception,GetMagickModule(),
1005  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1006  return(MagickFalse);
1007  }
1008  magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
1009  phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
1010  status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
1011  phase_pixels,exception);
1012  if (status != MagickFalse)
1013  status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
1014  phase_pixels,exception);
1015  phase_info=RelinquishVirtualMemory(phase_info);
1016  magnitude_info=RelinquishVirtualMemory(magnitude_info);
1017  return(status);
1018 }
1019 #endif
1020 
1021 MagickExport Image *ForwardFourierTransformImage(const Image *image,
1022  const MagickBooleanType modulus,ExceptionInfo *exception)
1023 {
1024  Image
1025  *fourier_image;
1026 
1027  fourier_image=NewImageList();
1028 #if !defined(MAGICKCORE_FFTW_DELEGATE) || defined(__cplusplus) || defined(c_plusplus)
1029  (void) modulus;
1030  (void) ThrowMagickException(exception,GetMagickModule(),
1031  MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1032  image->filename);
1033 #else
1034  {
1035  Image
1036  *magnitude_image;
1037 
1038  size_t
1039  height,
1040  width;
1041 
1042  width=image->columns;
1043  height=image->rows;
1044  if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
1045  ((image->rows % 2) != 0))
1046  {
1047  size_t extent=image->columns < image->rows ? image->rows :
1048  image->columns;
1049  width=(extent & 0x01) == 1 ? extent+1UL : extent;
1050  }
1051  height=width;
1052  magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
1053  if (magnitude_image != (Image *) NULL)
1054  {
1055  Image
1056  *phase_image;
1057 
1058  magnitude_image->storage_class=DirectClass;
1059  magnitude_image->depth=32UL;
1060  phase_image=CloneImage(image,width,height,MagickTrue,exception);
1061  if (phase_image == (Image *) NULL)
1062  magnitude_image=DestroyImage(magnitude_image);
1063  else
1064  {
1065  MagickBooleanType
1066  is_gray,
1067  status;
1068 
1069  phase_image->storage_class=DirectClass;
1070  phase_image->depth=32UL;
1071  AppendImageToList(&fourier_image,magnitude_image);
1072  AppendImageToList(&fourier_image,phase_image);
1073  status=MagickTrue;
1074  is_gray=IsGrayImage(image,exception);
1075 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1076  #pragma omp parallel sections
1077 #endif
1078  {
1079 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1080  #pragma omp section
1081 #endif
1082  {
1083  MagickBooleanType
1084  thread_status;
1085 
1086  if (is_gray != MagickFalse)
1087  thread_status=ForwardFourierTransformChannel(image,
1088  GrayChannels,modulus,fourier_image,exception);
1089  else
1090  thread_status=ForwardFourierTransformChannel(image,RedChannel,
1091  modulus,fourier_image,exception);
1092  if (thread_status == MagickFalse)
1093  status=thread_status;
1094  }
1095 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1096  #pragma omp section
1097 #endif
1098  {
1099  MagickBooleanType
1100  thread_status;
1101 
1102  thread_status=MagickTrue;
1103  if (is_gray == MagickFalse)
1104  thread_status=ForwardFourierTransformChannel(image,
1105  GreenChannel,modulus,fourier_image,exception);
1106  if (thread_status == MagickFalse)
1107  status=thread_status;
1108  }
1109 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1110  #pragma omp section
1111 #endif
1112  {
1113  MagickBooleanType
1114  thread_status;
1115 
1116  thread_status=MagickTrue;
1117  if (is_gray == MagickFalse)
1118  thread_status=ForwardFourierTransformChannel(image,
1119  BlueChannel,modulus,fourier_image,exception);
1120  if (thread_status == MagickFalse)
1121  status=thread_status;
1122  }
1123 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1124  #pragma omp section
1125 #endif
1126  {
1127  MagickBooleanType
1128  thread_status;
1129 
1130  thread_status=MagickTrue;
1131  if (image->matte != MagickFalse)
1132  thread_status=ForwardFourierTransformChannel(image,
1133  OpacityChannel,modulus,fourier_image,exception);
1134  if (thread_status == MagickFalse)
1135  status=thread_status;
1136  }
1137 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1138  #pragma omp section
1139 #endif
1140  {
1141  MagickBooleanType
1142  thread_status;
1143 
1144  thread_status=MagickTrue;
1145  if (image->colorspace == CMYKColorspace)
1146  thread_status=ForwardFourierTransformChannel(image,
1147  IndexChannel,modulus,fourier_image,exception);
1148  if (thread_status == MagickFalse)
1149  status=thread_status;
1150  }
1151  }
1152  if (status == MagickFalse)
1153  fourier_image=DestroyImageList(fourier_image);
1154  fftw_cleanup();
1155  }
1156  }
1157  }
1158 #endif
1159  return(fourier_image);
1160 }
1161 
1162 /*
1163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1164 % %
1165 % %
1166 % %
1167 % 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 %
1168 % %
1169 % %
1170 % %
1171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1172 %
1173 % InverseFourierTransformImage() implements the inverse discrete Fourier
1174 % transform (DFT) of the image either as a magnitude / phase or real /
1175 % imaginary image pair.
1176 %
1177 % The format of the InverseFourierTransformImage method is:
1178 %
1179 % Image *InverseFourierTransformImage(const Image *magnitude_image,
1180 % const Image *phase_image,const MagickBooleanType modulus,
1181 % ExceptionInfo *exception)
1182 %
1183 % A description of each parameter follows:
1184 %
1185 % o magnitude_image: the magnitude or real image.
1186 %
1187 % o phase_image: the phase or imaginary image.
1188 %
1189 % o modulus: if true, return transform as a magnitude / phase pair
1190 % otherwise a real / imaginary image pair.
1191 %
1192 % o exception: return any errors or warnings in this structure.
1193 %
1194 */
1195 
1196 #if defined(MAGICKCORE_FFTW_DELEGATE) && !defined(__cplusplus) && !defined(c_plusplus)
1197 static MagickBooleanType InverseQuadrantSwap(const size_t width,
1198  const size_t height,const double *source,double *destination)
1199 {
1200  ssize_t
1201  center,
1202  x,
1203  y;
1204 
1205  /*
1206  Swap quadrants.
1207  */
1208  center=(ssize_t) (width/2L)+1L;
1209  for (y=1L; y < (ssize_t) height; y++)
1210  for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1211  destination[(height-y)*center-x+width/2L]=source[y*width+x];
1212  for (y=0L; y < (ssize_t) height; y++)
1213  destination[y*center]=source[y*width+width/2L];
1214  for (x=0L; x < center; x++)
1215  destination[x]=source[center-x-1L];
1216  return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1217 }
1218 
1219 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
1220  const Image *magnitude_image,const Image *phase_image,
1221  fftw_complex *fourier_pixels,ExceptionInfo *exception)
1222 {
1223  CacheView
1224  *magnitude_view,
1225  *phase_view;
1226 
1227  double
1228  *inverse_pixels,
1229  *magnitude_pixels,
1230  *phase_pixels;
1231 
1232  MagickBooleanType
1233  status;
1234 
1235  MemoryInfo
1236  *inverse_info,
1237  *magnitude_info,
1238  *phase_info;
1239 
1240  const IndexPacket
1241  *indexes;
1242 
1243  const PixelPacket
1244  *p;
1245 
1246  ssize_t
1247  i,
1248  x,
1249  y;
1250 
1251  /*
1252  Inverse Fourier - read image and break down into a double array.
1253  */
1254  magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1255  sizeof(*magnitude_pixels));
1256  phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1257  sizeof(*phase_pixels));
1258  inverse_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/
1259  2+1)*sizeof(*inverse_pixels));
1260  if ((magnitude_info == (MemoryInfo *) NULL) ||
1261  (phase_info == (MemoryInfo *) NULL) ||
1262  (inverse_info == (MemoryInfo *) NULL))
1263  {
1264  if (magnitude_info != (MemoryInfo *) NULL)
1265  magnitude_info=RelinquishVirtualMemory(magnitude_info);
1266  if (phase_info != (MemoryInfo *) NULL)
1267  phase_info=RelinquishVirtualMemory(phase_info);
1268  if (inverse_info != (MemoryInfo *) NULL)
1269  inverse_info=RelinquishVirtualMemory(inverse_info);
1270  (void) ThrowMagickException(exception,GetMagickModule(),
1271  ResourceLimitError,"MemoryAllocationFailed","`%s'",
1272  magnitude_image->filename);
1273  return(MagickFalse);
1274  }
1275  magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
1276  phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
1277  inverse_pixels=(double *) GetVirtualMemoryBlob(inverse_info);
1278  i=0L;
1279  magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1280  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1281  {
1282  p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1283  exception);
1284  if (p == (const PixelPacket *) NULL)
1285  break;
1286  indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
1287  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1288  {
1289  switch (fourier_info->channel)
1290  {
1291  case RedChannel:
1292  default:
1293  {
1294  magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
1295  break;
1296  }
1297  case GreenChannel:
1298  {
1299  magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
1300  break;
1301  }
1302  case BlueChannel:
1303  {
1304  magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
1305  break;
1306  }
1307  case OpacityChannel:
1308  {
1309  magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
1310  break;
1311  }
1312  case IndexChannel:
1313  {
1314  magnitude_pixels[i]=QuantumScale*(MagickRealType)
1315  GetPixelIndex(indexes+x);
1316  break;
1317  }
1318  case GrayChannels:
1319  {
1320  magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
1321  break;
1322  }
1323  }
1324  i++;
1325  p++;
1326  }
1327  }
1328  magnitude_view=DestroyCacheView(magnitude_view);
1329  status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1330  magnitude_pixels,inverse_pixels);
1331  (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1332  fourier_info->center*sizeof(*magnitude_pixels));
1333  i=0L;
1334  phase_view=AcquireVirtualCacheView(phase_image,exception);
1335  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1336  {
1337  p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1338  exception);
1339  if (p == (const PixelPacket *) NULL)
1340  break;
1341  indexes=GetCacheViewAuthenticIndexQueue(phase_view);
1342  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1343  {
1344  switch (fourier_info->channel)
1345  {
1346  case RedChannel:
1347  default:
1348  {
1349  phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
1350  break;
1351  }
1352  case GreenChannel:
1353  {
1354  phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
1355  break;
1356  }
1357  case BlueChannel:
1358  {
1359  phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
1360  break;
1361  }
1362  case OpacityChannel:
1363  {
1364  phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
1365  break;
1366  }
1367  case IndexChannel:
1368  {
1369  phase_pixels[i]=QuantumScale*(MagickRealType)
1370  GetPixelIndex(indexes+x);
1371  break;
1372  }
1373  case GrayChannels:
1374  {
1375  phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
1376  break;
1377  }
1378  }
1379  i++;
1380  p++;
1381  }
1382  }
1383  if (fourier_info->modulus != MagickFalse)
1384  {
1385  i=0L;
1386  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1387  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1388  {
1389  phase_pixels[i]-=0.5;
1390  phase_pixels[i]*=(2.0*MagickPI);
1391  i++;
1392  }
1393  }
1394  phase_view=DestroyCacheView(phase_view);
1395  CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1396  if (status != MagickFalse)
1397  status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1398  phase_pixels,inverse_pixels);
1399  (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1400  fourier_info->center*sizeof(*phase_pixels));
1401  inverse_info=RelinquishVirtualMemory(inverse_info);
1402  /*
1403  Merge two sets.
1404  */
1405  i=0L;
1406  if (fourier_info->modulus != MagickFalse)
1407  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1408  for (x=0L; x < (ssize_t) fourier_info->center; x++)
1409  {
1410 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1411  fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+
1412  (MagickRealType) I*magnitude_pixels[i]*sin(phase_pixels[i]);
1413 #else
1414  fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1415  fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1416 #endif
1417  i++;
1418  }
1419  else
1420  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1421  for (x=0L; x < (ssize_t) fourier_info->center; x++)
1422  {
1423 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1424  fourier_pixels[i]=magnitude_pixels[i]+(MagickRealType) I*
1425  phase_pixels[i];
1426 #else
1427  fourier_pixels[i][0]=magnitude_pixels[i];
1428  fourier_pixels[i][1]=phase_pixels[i];
1429 #endif
1430  i++;
1431  }
1432  magnitude_info=RelinquishVirtualMemory(magnitude_info);
1433  phase_info=RelinquishVirtualMemory(phase_info);
1434  return(status);
1435 }
1436 
1437 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1438  fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception)
1439 {
1440  CacheView
1441  *image_view;
1442 
1443  double
1444  *source_pixels;
1445 
1446  const char
1447  *value;
1448 
1449  fftw_plan
1450  fftw_c2r_plan;
1451 
1452  MemoryInfo
1453  *source_info;
1454 
1455  IndexPacket
1456  *indexes;
1457 
1458  PixelPacket
1459  *q;
1460 
1461  ssize_t
1462  i,
1463  x,
1464  y;
1465 
1466  /*
1467  Generate the inverse Fourier transform.
1468  */
1469  if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
1470  {
1471  (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1472  "WidthOrHeightExceedsLimit","`%s'",image->filename);
1473  return(MagickFalse);
1474  }
1475  source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1476  sizeof(*source_pixels));
1477  if (source_info == (MemoryInfo *) NULL)
1478  {
1479  (void) ThrowMagickException(exception,GetMagickModule(),
1480  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1481  return(MagickFalse);
1482  }
1483  source_pixels=(double *) GetVirtualMemoryBlob(source_info);
1484  value=GetImageArtifact(image,"fourier:normalize");
1485  if (LocaleCompare(value,"inverse") == 0)
1486  {
1487  double
1488  gamma;
1489 
1490  /*
1491  Normalize inverse transform.
1492  */
1493  i=0L;
1494  gamma=PerceptibleReciprocal((double) fourier_info->width*
1495  fourier_info->height);
1496  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1497  for (x=0L; x < (ssize_t) fourier_info->center; x++)
1498  {
1499 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1500  fourier_pixels[i]*=gamma;
1501 #else
1502  fourier_pixels[i][0]*=gamma;
1503  fourier_pixels[i][1]*=gamma;
1504 #endif
1505  i++;
1506  }
1507  }
1508 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1509  #pragma omp critical (MagickCore_InverseFourierTransform)
1510 #endif
1511  fftw_c2r_plan=fftw_plan_dft_c2r_2d((int) fourier_info->width,
1512  (int) fourier_info->height,fourier_pixels,source_pixels,FFTW_ESTIMATE);
1513  fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1514  fftw_destroy_plan(fftw_c2r_plan);
1515  i=0L;
1516  image_view=AcquireAuthenticCacheView(image,exception);
1517  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1518  {
1519  if (y >= (ssize_t) image->rows)
1520  break;
1521  q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1522  image->columns ? image->columns : fourier_info->width,1UL,exception);
1523  if (q == (PixelPacket *) NULL)
1524  break;
1525  indexes=GetCacheViewAuthenticIndexQueue(image_view);
1526  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1527  {
1528  if (x < (ssize_t) image->columns)
1529  switch (fourier_info->channel)
1530  {
1531  case RedChannel:
1532  default:
1533  {
1534  SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
1535  source_pixels[i]));
1536  break;
1537  }
1538  case GreenChannel:
1539  {
1540  SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
1541  source_pixels[i]));
1542  break;
1543  }
1544  case BlueChannel:
1545  {
1546  SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
1547  source_pixels[i]));
1548  break;
1549  }
1550  case OpacityChannel:
1551  {
1552  SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*
1553  source_pixels[i]));
1554  break;
1555  }
1556  case IndexChannel:
1557  {
1558  SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType)
1559  QuantumRange*source_pixels[i]));
1560  break;
1561  }
1562  case GrayChannels:
1563  {
1564  SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*
1565  source_pixels[i]));
1566  break;
1567  }
1568  }
1569  i++;
1570  q++;
1571  }
1572  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1573  break;
1574  }
1575  image_view=DestroyCacheView(image_view);
1576  source_info=RelinquishVirtualMemory(source_info);
1577  return(MagickTrue);
1578 }
1579 
1580 static MagickBooleanType InverseFourierTransformChannel(
1581  const Image *magnitude_image,const Image *phase_image,
1582  const ChannelType channel,const MagickBooleanType modulus,
1583  Image *fourier_image,ExceptionInfo *exception)
1584 {
1585  fftw_complex
1586  *inverse_pixels;
1587 
1588  FourierInfo
1589  fourier_info;
1590 
1591  MagickBooleanType
1592  status;
1593 
1594  MemoryInfo
1595  *inverse_info;
1596 
1597  fourier_info.width=magnitude_image->columns;
1598  fourier_info.height=magnitude_image->rows;
1599  if ((magnitude_image->columns != magnitude_image->rows) ||
1600  ((magnitude_image->columns % 2) != 0) ||
1601  ((magnitude_image->rows % 2) != 0))
1602  {
1603  size_t extent=magnitude_image->columns < magnitude_image->rows ?
1604  magnitude_image->rows : magnitude_image->columns;
1605  fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1606  }
1607  fourier_info.height=fourier_info.width;
1608  fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1609  fourier_info.channel=channel;
1610  fourier_info.modulus=modulus;
1611  inverse_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1612  1)*sizeof(*inverse_pixels));
1613  if (inverse_info == (MemoryInfo *) NULL)
1614  {
1615  (void) ThrowMagickException(exception,GetMagickModule(),
1616  ResourceLimitError,"MemoryAllocationFailed","`%s'",
1617  magnitude_image->filename);
1618  return(MagickFalse);
1619  }
1620  inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1621  status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1622  inverse_pixels,exception);
1623  if (status != MagickFalse)
1624  status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1625  exception);
1626  inverse_info=RelinquishVirtualMemory(inverse_info);
1627  return(status);
1628 }
1629 #endif
1630 
1631 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1632  const Image *phase_image,const MagickBooleanType modulus,
1633  ExceptionInfo *exception)
1634 {
1635  Image
1636  *fourier_image;
1637 
1638  assert(magnitude_image != (Image *) NULL);
1639  assert(magnitude_image->signature == MagickCoreSignature);
1640  if (IsEventLogging() != MagickFalse)
1641  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1642  magnitude_image->filename);
1643  if (phase_image == (Image *) NULL)
1644  {
1645  (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1646  "ImageSequenceRequired","`%s'",magnitude_image->filename);
1647  return((Image *) NULL);
1648  }
1649 #if !defined(MAGICKCORE_FFTW_DELEGATE) || defined(__cplusplus) || defined(c_plusplus)
1650  fourier_image=(Image *) NULL;
1651  (void) modulus;
1652  (void) ThrowMagickException(exception,GetMagickModule(),
1653  MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1654  magnitude_image->filename);
1655 #else
1656  {
1657  fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1658  magnitude_image->rows,MagickTrue,exception);
1659  if (fourier_image != (Image *) NULL)
1660  {
1661  MagickBooleanType
1662  is_gray,
1663  status;
1664 
1665  status=MagickTrue;
1666  is_gray=IsGrayImage(magnitude_image,exception);
1667  if (is_gray != MagickFalse)
1668  is_gray=IsGrayImage(phase_image,exception);
1669 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1670  #pragma omp parallel sections
1671 #endif
1672  {
1673 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1674  #pragma omp section
1675 #endif
1676  {
1677  MagickBooleanType
1678  thread_status;
1679 
1680  if (is_gray != MagickFalse)
1681  thread_status=InverseFourierTransformChannel(magnitude_image,
1682  phase_image,GrayChannels,modulus,fourier_image,exception);
1683  else
1684  thread_status=InverseFourierTransformChannel(magnitude_image,
1685  phase_image,RedChannel,modulus,fourier_image,exception);
1686  if (thread_status == MagickFalse)
1687  status=thread_status;
1688  }
1689 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1690  #pragma omp section
1691 #endif
1692  {
1693  MagickBooleanType
1694  thread_status;
1695 
1696  thread_status=MagickTrue;
1697  if (is_gray == MagickFalse)
1698  thread_status=InverseFourierTransformChannel(magnitude_image,
1699  phase_image,GreenChannel,modulus,fourier_image,exception);
1700  if (thread_status == MagickFalse)
1701  status=thread_status;
1702  }
1703 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1704  #pragma omp section
1705 #endif
1706  {
1707  MagickBooleanType
1708  thread_status;
1709 
1710  thread_status=MagickTrue;
1711  if (is_gray == MagickFalse)
1712  thread_status=InverseFourierTransformChannel(magnitude_image,
1713  phase_image,BlueChannel,modulus,fourier_image,exception);
1714  if (thread_status == MagickFalse)
1715  status=thread_status;
1716  }
1717 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1718  #pragma omp section
1719 #endif
1720  {
1721  MagickBooleanType
1722  thread_status;
1723 
1724  thread_status=MagickTrue;
1725  if (magnitude_image->matte != MagickFalse)
1726  thread_status=InverseFourierTransformChannel(magnitude_image,
1727  phase_image,OpacityChannel,modulus,fourier_image,exception);
1728  if (thread_status == MagickFalse)
1729  status=thread_status;
1730  }
1731 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1732  #pragma omp section
1733 #endif
1734  {
1735  MagickBooleanType
1736  thread_status;
1737 
1738  thread_status=MagickTrue;
1739  if (magnitude_image->colorspace == CMYKColorspace)
1740  thread_status=InverseFourierTransformChannel(magnitude_image,
1741  phase_image,IndexChannel,modulus,fourier_image,exception);
1742  if (thread_status == MagickFalse)
1743  status=thread_status;
1744  }
1745  }
1746  if (status == MagickFalse)
1747  fourier_image=DestroyImage(fourier_image);
1748  }
1749  fftw_cleanup();
1750  }
1751 #endif
1752  return(fourier_image);
1753 }
Definition: image.h:152