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