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