2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Discrete Fourier Transform Methods %
22 % Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization %
23 % dedicated to making software imaging solutions freely available. %
25 % You may not use this file except in compliance with the License. You may %
26 % obtain a copy of the License at %
28 % http://www.imagemagick.org/script/license.php %
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. %
36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 #include "MagickCore/studio.h"
46 #include "MagickCore/attribute.h"
47 #include "MagickCore/cache.h"
48 #include "MagickCore/image.h"
49 #include "MagickCore/image-private.h"
50 #include "MagickCore/list.h"
51 #include "MagickCore/fourier.h"
52 #include "MagickCore/log.h"
53 #include "MagickCore/memory_.h"
54 #include "MagickCore/monitor.h"
55 #include "MagickCore/pixel-accessor.h"
56 #include "MagickCore/property.h"
57 #include "MagickCore/quantum-private.h"
58 #include "MagickCore/thread-private.h"
59 #if defined(MAGICKCORE_FFTW_DELEGATE)
60 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
64 #if !defined(MAGICKCORE_HAVE_CABS)
65 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
67 #if !defined(MAGICKCORE_HAVE_CARG)
68 #define carg(z) (atan2(cimag(z),creal(z)))
70 #if !defined(MAGICKCORE_HAVE_CIMAG)
71 #define cimag(z) (z[1])
73 #if !defined(MAGICKCORE_HAVE_CREAL)
74 #define creal(z) (z[0])
81 typedef struct _FourierInfo
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
102 % F o r w a r d F o u r i e r T r a n s f o r m I m a g e %
106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 % ForwardFourierTransformImage() implements the discrete Fourier transform
109 % (DFT) of the image either as a magnitude / phase or real / imaginary image
112 % The format of the ForwadFourierTransformImage method is:
114 % Image *ForwardFourierTransformImage(const Image *image,
115 % const MagickBooleanType modulus,ExceptionInfo *exception)
117 % A description of each parameter follows:
119 % o image: the image.
121 % o modulus: if true, return as transform as a magnitude / phase pair
122 % otherwise a real / imaginary image pair.
124 % o exception: return any errors or warnings in this structure.
128 #if defined(MAGICKCORE_FFTW_DELEGATE)
130 static MagickBooleanType RollFourier(const size_t width,const size_t height,
131 const ssize_t x_offset,const ssize_t y_offset,double *fourier)
146 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
148 roll=(double *) AcquireQuantumMemory((size_t) height,width*sizeof(*roll));
149 if (roll == (double *) NULL)
152 for (y=0L; y < (ssize_t) height; y++)
155 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
157 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
159 for (x=0L; x < (ssize_t) width; x++)
162 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
164 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
166 roll[v*width+u]=fourier[i++];
169 (void) CopyMagickMemory(fourier,roll,height*width*sizeof(*roll));
170 roll=(double *) RelinquishMagickMemory(roll);
174 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
175 const size_t height,double *source,double *destination)
190 center=(ssize_t) floor((double) width/2L)+1L;
191 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
192 if (status == MagickFalse)
194 for (y=0L; y < (ssize_t) height; y++)
195 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
196 destination[width*y+x+width/2L]=source[center*y+x];
197 for (y=1; y < (ssize_t) height; y++)
198 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
199 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
200 for (x=0L; x < (ssize_t) (width/2L); x++)
201 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
205 static void CorrectPhaseLHS(const size_t width,const size_t height,
214 for (y=0L; y < (ssize_t) height; y++)
215 for (x=0L; x < (ssize_t) (width/2L); x++)
216 fourier[y*width+x]*=(-1.0);
219 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
220 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
247 magnitude_image=GetFirstImageInList(image);
248 phase_image=GetNextImageInList(image);
249 if (phase_image == (Image *) NULL)
251 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
252 "ImageSequenceRequired","`%s'",image->filename);
256 Create "Fourier Transform" image from constituent arrays.
258 magnitude_source=(double *) AcquireQuantumMemory((size_t)
259 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
260 if (magnitude_source == (double *) NULL)
262 (void) ResetMagickMemory(magnitude_source,0,fourier_info->height*
263 fourier_info->width*sizeof(*magnitude_source));
264 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
265 fourier_info->width*sizeof(*phase_source));
266 if (phase_source == (double *) NULL)
268 (void) ThrowMagickException(exception,GetMagickModule(),
269 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
270 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
273 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
274 magnitude,magnitude_source);
275 if (status != MagickFalse)
276 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
278 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
279 if (fourier_info->modulus != MagickFalse)
282 for (y=0L; y < (ssize_t) fourier_info->height; y++)
283 for (x=0L; x < (ssize_t) fourier_info->width; x++)
285 phase_source[i]/=(2.0*MagickPI);
286 phase_source[i]+=0.5;
290 magnitude_view=AcquireCacheView(magnitude_image);
291 phase_view=AcquireCacheView(phase_image);
293 for (y=0L; y < (ssize_t) fourier_info->height; y++)
295 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
297 if (q == (Quantum *) NULL)
299 for (x=0L; x < (ssize_t) fourier_info->width; x++)
301 switch (fourier_info->channel)
303 case RedPixelChannel:
306 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
307 magnitude_source[i]),q);
310 case GreenPixelChannel:
312 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
313 magnitude_source[i]),q);
316 case BluePixelChannel:
318 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
319 magnitude_source[i]),q);
322 case BlackPixelChannel:
324 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
325 magnitude_source[i]),q);
328 case AlphaPixelChannel:
330 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
331 magnitude_source[i]),q);
336 q+=GetPixelChannels(magnitude_image);
338 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
339 if (status == MagickFalse)
343 for (y=0L; y < (ssize_t) fourier_info->height; y++)
345 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
347 if (q == (Quantum *) NULL)
349 for (x=0L; x < (ssize_t) fourier_info->width; x++)
351 switch (fourier_info->channel)
353 case RedPixelChannel:
356 SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
360 case GreenPixelChannel:
362 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
366 case BluePixelChannel:
368 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
372 case BlackPixelChannel:
374 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
378 case AlphaPixelChannel:
380 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
386 q+=GetPixelChannels(phase_image);
388 status=SyncCacheViewAuthenticPixels(phase_view,exception);
389 if (status == MagickFalse)
392 phase_view=DestroyCacheView(phase_view);
393 magnitude_view=DestroyCacheView(magnitude_view);
394 phase_source=(double *) RelinquishMagickMemory(phase_source);
395 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
399 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
400 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
415 register const Quantum
426 Generate the forward Fourier transform.
428 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
429 fourier_info->width*sizeof(*source));
430 if (source == (double *) NULL)
432 (void) ThrowMagickException(exception,GetMagickModule(),
433 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
436 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
439 image_view=AcquireCacheView(image);
440 for (y=0L; y < (ssize_t) fourier_info->height; y++)
442 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
444 if (p == (const Quantum *) NULL)
446 for (x=0L; x < (ssize_t) fourier_info->width; x++)
448 switch (fourier_info->channel)
450 case RedPixelChannel:
453 source[i]=QuantumScale*GetPixelRed(image,p);
456 case GreenPixelChannel:
458 source[i]=QuantumScale*GetPixelGreen(image,p);
461 case BluePixelChannel:
463 source[i]=QuantumScale*GetPixelBlue(image,p);
466 case BlackPixelChannel:
468 source[i]=QuantumScale*GetPixelBlack(image,p);
471 case AlphaPixelChannel:
473 source[i]=QuantumScale*GetPixelAlpha(image,p);
478 p+=GetPixelChannels(image);
481 image_view=DestroyCacheView(image_view);
482 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
483 fourier_info->center*sizeof(*fourier));
484 if (fourier == (fftw_complex *) NULL)
486 (void) ThrowMagickException(exception,GetMagickModule(),
487 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
488 source=(double *) RelinquishMagickMemory(source);
491 #if defined(MAGICKCORE_OPENMP_SUPPORT)
492 #pragma omp critical (MagickCore_ForwardFourierTransform)
494 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
495 source,fourier,FFTW_ESTIMATE);
496 fftw_execute(fftw_r2c_plan);
497 fftw_destroy_plan(fftw_r2c_plan);
498 source=(double *) RelinquishMagickMemory(source);
500 Normalize Fourier transform.
502 n=(double) fourier_info->width*(double) fourier_info->width;
504 for (y=0L; y < (ssize_t) fourier_info->height; y++)
505 for (x=0L; x < (ssize_t) fourier_info->center; x++)
507 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
516 Generate magnitude and phase (or real and imaginary).
519 if (fourier_info->modulus != MagickFalse)
520 for (y=0L; y < (ssize_t) fourier_info->height; y++)
521 for (x=0L; x < (ssize_t) fourier_info->center; x++)
523 magnitude[i]=cabs(fourier[i]);
524 phase[i]=carg(fourier[i]);
528 for (y=0L; y < (ssize_t) fourier_info->height; y++)
529 for (x=0L; x < (ssize_t) fourier_info->center; x++)
531 magnitude[i]=creal(fourier[i]);
532 phase[i]=cimag(fourier[i]);
535 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
539 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
540 const PixelChannel channel,const MagickBooleanType modulus,
541 Image *fourier_image,ExceptionInfo *exception)
559 fourier_info.width=image->columns;
560 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
561 ((image->rows % 2) != 0))
563 extent=image->columns < image->rows ? image->rows : image->columns;
564 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
566 fourier_info.height=fourier_info.width;
567 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
568 fourier_info.channel=channel;
569 fourier_info.modulus=modulus;
570 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
571 fourier_info.center*sizeof(*magnitude));
572 if (magnitude == (double *) NULL)
574 (void) ThrowMagickException(exception,GetMagickModule(),
575 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
578 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
579 fourier_info.center*sizeof(*phase));
580 if (phase == (double *) NULL)
582 (void) ThrowMagickException(exception,GetMagickModule(),
583 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
584 magnitude=(double *) RelinquishMagickMemory(magnitude);
587 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
588 fourier_info.center*sizeof(*fourier));
589 if (fourier == (fftw_complex *) NULL)
591 (void) ThrowMagickException(exception,GetMagickModule(),
592 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
593 phase=(double *) RelinquishMagickMemory(phase);
594 magnitude=(double *) RelinquishMagickMemory(magnitude);
597 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
598 if (status != MagickFalse)
599 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
601 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
602 phase=(double *) RelinquishMagickMemory(phase);
603 magnitude=(double *) RelinquishMagickMemory(magnitude);
608 MagickExport Image *ForwardFourierTransformImage(const Image *image,
609 const MagickBooleanType modulus,ExceptionInfo *exception)
614 fourier_image=NewImageList();
615 #if !defined(MAGICKCORE_FFTW_DELEGATE)
617 (void) ThrowMagickException(exception,GetMagickModule(),
618 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
629 width=image->columns;
630 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
631 ((image->rows % 2) != 0))
633 extent=image->columns < image->rows ? image->rows : image->columns;
634 width=(extent & 0x01) == 1 ? extent+1UL : extent;
636 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
637 if (magnitude_image != (Image *) NULL)
642 magnitude_image->storage_class=DirectClass;
643 magnitude_image->depth=32UL;
644 phase_image=CloneImage(image,width,width,MagickFalse,exception);
645 if (phase_image == (Image *) NULL)
646 magnitude_image=DestroyImage(magnitude_image);
653 phase_image->storage_class=DirectClass;
654 phase_image->depth=32UL;
655 AppendImageToList(&fourier_image,magnitude_image);
656 AppendImageToList(&fourier_image,phase_image);
658 is_gray=IsImageGray(image,exception);
659 #if defined(MAGICKCORE_OPENMP_SUPPORT)
660 #pragma omp parallel sections
663 #if defined(MAGICKCORE_OPENMP_SUPPORT)
670 if (is_gray != MagickFalse)
671 thread_status=ForwardFourierTransformChannel(image,
672 GrayPixelChannel,modulus,fourier_image,exception);
674 thread_status=ForwardFourierTransformChannel(image,
675 RedPixelChannel,modulus,fourier_image,exception);
676 if (thread_status == MagickFalse)
677 status=thread_status;
679 #if defined(MAGICKCORE_OPENMP_SUPPORT)
686 thread_status=MagickTrue;
687 if (is_gray == MagickFalse)
688 thread_status=ForwardFourierTransformChannel(image,
689 GreenPixelChannel,modulus,fourier_image,exception);
690 if (thread_status == MagickFalse)
691 status=thread_status;
693 #if defined(MAGICKCORE_OPENMP_SUPPORT)
700 thread_status=MagickTrue;
701 if (is_gray == MagickFalse)
702 thread_status=ForwardFourierTransformChannel(image,
703 BluePixelChannel,modulus,fourier_image,exception);
704 if (thread_status == MagickFalse)
705 status=thread_status;
707 #if defined(MAGICKCORE_OPENMP_SUPPORT)
714 thread_status=MagickTrue;
715 if (image->colorspace == CMYKColorspace)
716 thread_status=ForwardFourierTransformChannel(image,
717 BlackPixelChannel,modulus,fourier_image,exception);
718 if (thread_status == MagickFalse)
719 status=thread_status;
721 #if defined(MAGICKCORE_OPENMP_SUPPORT)
728 thread_status=MagickTrue;
729 if (image->matte != MagickFalse)
730 thread_status=ForwardFourierTransformChannel(image,
731 AlphaPixelChannel,modulus,fourier_image,exception);
732 if (thread_status == MagickFalse)
733 status=thread_status;
736 if (status == MagickFalse)
737 fourier_image=DestroyImageList(fourier_image);
743 return(fourier_image);
747 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
751 % 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 %
755 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
757 % InverseFourierTransformImage() implements the inverse discrete Fourier
758 % transform (DFT) of the image either as a magnitude / phase or real /
759 % imaginary image pair.
761 % The format of the InverseFourierTransformImage method is:
763 % Image *InverseFourierTransformImage(const Image *magnitude_image,
764 % const Image *phase_image,const MagickBooleanType modulus,
765 % ExceptionInfo *exception)
767 % A description of each parameter follows:
769 % o magnitude_image: the magnitude or real image.
771 % o phase_image: the phase or imaginary image.
773 % o modulus: if true, return transform as a magnitude / phase pair
774 % otherwise a real / imaginary image pair.
776 % o exception: return any errors or warnings in this structure.
780 #if defined(MAGICKCORE_FFTW_DELEGATE)
781 static MagickBooleanType InverseQuadrantSwap(const size_t width,
782 const size_t height,const double *source,double *destination)
794 center=(ssize_t) floor((double) width/2.0)+1L;
795 for (y=1L; y < (ssize_t) height; y++)
796 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
797 destination[center*(height-y)-x+width/2L]=source[y*width+x];
798 for (y=0L; y < (ssize_t) height; y++)
799 destination[center*y]=source[y*width+width/2L];
800 for (x=0L; x < center; x++)
801 destination[x]=source[center-x-1L];
802 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
805 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
806 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
807 ExceptionInfo *exception)
822 register const Quantum
833 Inverse fourier - read image and break down into a double array.
835 magnitude_source=(double *) AcquireQuantumMemory((size_t)
836 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
837 if (magnitude_source == (double *) NULL)
839 (void) ThrowMagickException(exception,GetMagickModule(),
840 ResourceLimitError,"MemoryAllocationFailed","`%s'",
841 magnitude_image->filename);
844 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
845 fourier_info->width*sizeof(*phase_source));
846 if (phase_source == (double *) NULL)
848 (void) ThrowMagickException(exception,GetMagickModule(),
849 ResourceLimitError,"MemoryAllocationFailed","`%s'",
850 magnitude_image->filename);
851 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
855 magnitude_view=AcquireCacheView(magnitude_image);
856 for (y=0L; y < (ssize_t) fourier_info->height; y++)
858 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
860 if (p == (const Quantum *) NULL)
862 for (x=0L; x < (ssize_t) fourier_info->width; x++)
864 switch (fourier_info->channel)
866 case RedPixelChannel:
869 magnitude_source[i]=QuantumScale*GetPixelRed(magnitude_image,p);
872 case GreenPixelChannel:
874 magnitude_source[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
877 case BluePixelChannel:
879 magnitude_source[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
882 case BlackPixelChannel:
884 magnitude_source[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
887 case AlphaPixelChannel:
889 magnitude_source[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
894 p+=GetPixelChannels(magnitude_image);
898 phase_view=AcquireCacheView(phase_image);
899 for (y=0L; y < (ssize_t) fourier_info->height; y++)
901 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
903 if (p == (const Quantum *) NULL)
905 for (x=0L; x < (ssize_t) fourier_info->width; x++)
907 switch (fourier_info->channel)
909 case RedPixelChannel:
912 phase_source[i]=QuantumScale*GetPixelRed(phase_image,p);
915 case GreenPixelChannel:
917 phase_source[i]=QuantumScale*GetPixelGreen(phase_image,p);
920 case BluePixelChannel:
922 phase_source[i]=QuantumScale*GetPixelBlue(phase_image,p);
925 case BlackPixelChannel:
927 phase_source[i]=QuantumScale*GetPixelBlack(phase_image,p);
930 case AlphaPixelChannel:
932 phase_source[i]=QuantumScale*GetPixelAlpha(phase_image,p);
937 p+=GetPixelChannels(phase_image);
940 if (fourier_info->modulus != MagickFalse)
943 for (y=0L; y < (ssize_t) fourier_info->height; y++)
944 for (x=0L; x < (ssize_t) fourier_info->width; x++)
946 phase_source[i]-=0.5;
947 phase_source[i]*=(2.0*MagickPI);
951 magnitude_view=DestroyCacheView(magnitude_view);
952 phase_view=DestroyCacheView(phase_view);
953 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
954 fourier_info->center*sizeof(*magnitude));
955 if (magnitude == (double *) NULL)
957 (void) ThrowMagickException(exception,GetMagickModule(),
958 ResourceLimitError,"MemoryAllocationFailed","`%s'",
959 magnitude_image->filename);
960 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
961 phase_source=(double *) RelinquishMagickMemory(phase_source);
964 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
965 magnitude_source,magnitude);
966 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
967 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
968 fourier_info->width*sizeof(*phase));
969 if (phase == (double *) NULL)
971 (void) ThrowMagickException(exception,GetMagickModule(),
972 ResourceLimitError,"MemoryAllocationFailed","`%s'",
973 magnitude_image->filename);
974 phase_source=(double *) RelinquishMagickMemory(phase_source);
977 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
978 if (status != MagickFalse)
979 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
981 phase_source=(double *) RelinquishMagickMemory(phase_source);
986 if (fourier_info->modulus != MagickFalse)
987 for (y=0L; y < (ssize_t) fourier_info->height; y++)
988 for (x=0L; x < (ssize_t) fourier_info->center; x++)
990 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
991 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
993 fourier[i][0]=magnitude[i]*cos(phase[i]);
994 fourier[i][1]=magnitude[i]*sin(phase[i]);
999 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1000 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1002 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1003 fourier[i]=magnitude[i]+I*phase[i];
1005 fourier[i][0]=magnitude[i];
1006 fourier[i][1]=phase[i];
1010 phase=(double *) RelinquishMagickMemory(phase);
1011 magnitude=(double *) RelinquishMagickMemory(magnitude);
1015 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1016 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1037 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1038 fourier_info->width*sizeof(*source));
1039 if (source == (double *) NULL)
1041 (void) ThrowMagickException(exception,GetMagickModule(),
1042 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1043 return(MagickFalse);
1045 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1046 #pragma omp critical (MagickCore_InverseFourierTransform)
1049 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1050 fourier,source,FFTW_ESTIMATE);
1051 fftw_execute(fftw_c2r_plan);
1052 fftw_destroy_plan(fftw_c2r_plan);
1055 image_view=AcquireCacheView(image);
1056 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1058 if (y >= (ssize_t) image->rows)
1060 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1061 image->columns ? image->columns : fourier_info->width,1UL,exception);
1062 if (q == (Quantum *) NULL)
1064 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1066 switch (fourier_info->channel)
1068 case RedPixelChannel:
1071 SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
1074 case GreenPixelChannel:
1076 SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
1079 case BluePixelChannel:
1081 SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
1084 case BlackPixelChannel:
1086 SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
1089 case AlphaPixelChannel:
1091 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
1096 q+=GetPixelChannels(image);
1098 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1101 image_view=DestroyCacheView(image_view);
1102 source=(double *) RelinquishMagickMemory(source);
1106 static MagickBooleanType InverseFourierTransformChannel(
1107 const Image *magnitude_image,const Image *phase_image,
1108 const PixelChannel channel,const MagickBooleanType modulus,
1109 Image *fourier_image,ExceptionInfo *exception)
1127 fourier_info.width=magnitude_image->columns;
1128 if ((magnitude_image->columns != magnitude_image->rows) ||
1129 ((magnitude_image->columns % 2) != 0) ||
1130 ((magnitude_image->rows % 2) != 0))
1132 extent=magnitude_image->columns < magnitude_image->rows ?
1133 magnitude_image->rows : magnitude_image->columns;
1134 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1136 fourier_info.height=fourier_info.width;
1137 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
1138 fourier_info.channel=channel;
1139 fourier_info.modulus=modulus;
1140 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1141 fourier_info.center*sizeof(*magnitude));
1142 if (magnitude == (double *) NULL)
1144 (void) ThrowMagickException(exception,GetMagickModule(),
1145 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1146 magnitude_image->filename);
1147 return(MagickFalse);
1149 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1150 fourier_info.center*sizeof(*phase));
1151 if (phase == (double *) NULL)
1153 (void) ThrowMagickException(exception,GetMagickModule(),
1154 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1155 magnitude_image->filename);
1156 magnitude=(double *) RelinquishMagickMemory(magnitude);
1157 return(MagickFalse);
1159 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1160 fourier_info.center*sizeof(*fourier));
1161 if (fourier == (fftw_complex *) NULL)
1163 (void) ThrowMagickException(exception,GetMagickModule(),
1164 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1165 magnitude_image->filename);
1166 phase=(double *) RelinquishMagickMemory(phase);
1167 magnitude=(double *) RelinquishMagickMemory(magnitude);
1168 return(MagickFalse);
1170 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1172 if (status != MagickFalse)
1173 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1175 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1176 phase=(double *) RelinquishMagickMemory(phase);
1177 magnitude=(double *) RelinquishMagickMemory(magnitude);
1182 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1183 const Image *phase_image,const MagickBooleanType modulus,
1184 ExceptionInfo *exception)
1189 assert(magnitude_image != (Image *) NULL);
1190 assert(magnitude_image->signature == MagickSignature);
1191 if (magnitude_image->debug != MagickFalse)
1192 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1193 magnitude_image->filename);
1194 if (phase_image == (Image *) NULL)
1196 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1197 "ImageSequenceRequired","`%s'",magnitude_image->filename);
1198 return((Image *) NULL);
1200 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1201 fourier_image=(Image *) NULL;
1203 (void) ThrowMagickException(exception,GetMagickModule(),
1204 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1205 magnitude_image->filename);
1208 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1209 magnitude_image->rows,MagickFalse,exception);
1210 if (fourier_image != (Image *) NULL)
1217 is_gray=IsImageGray(magnitude_image,exception);
1218 if (is_gray != MagickFalse)
1219 is_gray=IsImageGray(phase_image,exception);
1220 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1221 #pragma omp parallel sections
1224 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1231 if (is_gray != MagickFalse)
1232 thread_status=InverseFourierTransformChannel(magnitude_image,
1233 phase_image,GrayPixelChannel,modulus,fourier_image,exception);
1235 thread_status=InverseFourierTransformChannel(magnitude_image,
1236 phase_image,RedPixelChannel,modulus,fourier_image,exception);
1237 if (thread_status == MagickFalse)
1238 status=thread_status;
1240 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1247 thread_status=MagickTrue;
1248 if (is_gray == MagickFalse)
1249 thread_status=InverseFourierTransformChannel(magnitude_image,
1250 phase_image,GreenPixelChannel,modulus,fourier_image,exception);
1251 if (thread_status == MagickFalse)
1252 status=thread_status;
1254 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1261 thread_status=MagickTrue;
1262 if (is_gray == MagickFalse)
1263 thread_status=InverseFourierTransformChannel(magnitude_image,
1264 phase_image,BluePixelChannel,modulus,fourier_image,exception);
1265 if (thread_status == MagickFalse)
1266 status=thread_status;
1268 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1275 thread_status=MagickTrue;
1276 if (magnitude_image->colorspace == CMYKColorspace)
1277 thread_status=InverseFourierTransformChannel(magnitude_image,
1278 phase_image,BlackPixelChannel,modulus,fourier_image,exception);
1279 if (thread_status == MagickFalse)
1280 status=thread_status;
1282 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1289 thread_status=MagickTrue;
1290 if (magnitude_image->matte != MagickFalse)
1291 thread_status=InverseFourierTransformChannel(magnitude_image,
1292 phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
1293 if (thread_status == MagickFalse)
1294 status=thread_status;
1297 if (status == MagickFalse)
1298 fourier_image=DestroyImage(fourier_image);
1303 return(fourier_image);