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