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 "magick/studio.h"
46 #include "magick/cache.h"
47 #include "magick/image.h"
48 #include "magick/image-private.h"
49 #include "magick/list.h"
50 #include "magick/fourier.h"
51 #include "magick/log.h"
52 #include "magick/memory_.h"
53 #include "magick/monitor.h"
54 #include "magick/property.h"
55 #include "magick/quantum-private.h"
56 #include "magick/thread-private.h"
57 #if defined(MAGICKCORE_FFTW_DELEGATE)
58 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
62 #if !defined(MAGICKCORE_HAVE_CABS)
63 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
65 #if !defined(MAGICKCORE_HAVE_CARG)
66 #define carg(z) (atan2(cimag(z[1]),creal(z[0])))
68 #if !defined(MAGICKCORE_HAVE_CIMAG)
69 #define cimag(z) (z[1])
71 #if !defined(MAGICKCORE_HAVE_CREAL)
72 #define creal(z) (z[0])
79 typedef struct _FourierInfo
96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100 % 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 %
104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 % ForwardFourierTransformImage() implements the discrete Fourier transform
107 % (DFT) of the image either as a magnitude / phase or real / imaginary image
110 % The format of the ForwadFourierTransformImage method is:
112 % Image *ForwardFourierTransformImage(const Image *image,
113 % const MagickBooleanType modulus,ExceptionInfo *exception)
115 % A description of each parameter follows:
117 % o image: the image.
119 % o modulus: if true, return as transform as a magnitude / phase pair
120 % otherwise a real / imaginary image pair.
122 % o exception: return any errors or warnings in this structure.
126 #if defined(MAGICKCORE_FFTW_DELEGATE)
128 static MagickBooleanType RollFourier(const size_t width,const size_t height,
129 const ssize_t x_offset,const ssize_t y_offset,double *fourier)
144 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
146 roll=(double *) AcquireQuantumMemory((size_t) height,width*sizeof(*roll));
147 if (roll == (double *) NULL)
150 for (y=0L; y < (ssize_t) height; y++)
153 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
155 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
157 for (x=0L; x < (ssize_t) width; x++)
160 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
162 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
164 roll[v*width+u]=fourier[i++];
167 (void) CopyMagickMemory(fourier,roll,height*width*sizeof(*roll));
168 roll=(double *) RelinquishMagickMemory(roll);
172 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
173 const size_t height,double *source,double *destination)
188 center=(ssize_t) floor((double) width/2L)+1L;
189 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
190 if (status == MagickFalse)
192 for (y=0L; y < (ssize_t) height; y++)
193 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
194 destination[width*y+x+width/2L]=source[center*y+x];
195 for (y=1; y < (ssize_t) height; y++)
196 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
197 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
198 for (x=0L; x < (ssize_t) (width/2L); x++)
199 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
203 static void CorrectPhaseLHS(const size_t width,const size_t height,
212 for (y=0L; y < (ssize_t) height; y++)
213 for (x=0L; x < (ssize_t) (width/2L); x++)
214 fourier[y*width+x]*=(-1.0);
217 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
218 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
248 magnitude_image=GetFirstImageInList(image);
249 phase_image=GetNextImageInList(image);
250 if (phase_image == (Image *) NULL)
252 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
253 "ImageSequenceRequired","`%s'",image->filename);
257 Create "Fourier Transform" image from constituent arrays.
259 magnitude_source=(double *) AcquireQuantumMemory((size_t)
260 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
261 if (magnitude_source == (double *) NULL)
263 (void) ResetMagickMemory(magnitude_source,0,fourier_info->height*
264 fourier_info->width*sizeof(*magnitude_source));
265 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
266 fourier_info->width*sizeof(*phase_source));
267 if (phase_source == (double *) NULL)
269 (void) ThrowMagickException(exception,GetMagickModule(),
270 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
271 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
274 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
275 magnitude,magnitude_source);
276 if (status != MagickFalse)
277 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
279 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
280 if (fourier_info->modulus != MagickFalse)
283 for (y=0L; y < (ssize_t) fourier_info->height; y++)
284 for (x=0L; x < (ssize_t) fourier_info->width; x++)
286 phase_source[i]/=(2.0*MagickPI);
287 phase_source[i]+=0.5;
291 magnitude_view=AcquireCacheView(magnitude_image);
292 phase_view=AcquireCacheView(phase_image);
294 for (y=0L; y < (ssize_t) fourier_info->height; y++)
296 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
298 if (q == (PixelPacket *) NULL)
300 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
301 for (x=0L; x < (ssize_t) fourier_info->width; x++)
303 switch (fourier_info->channel)
308 q->red=ClampToQuantum(QuantumRange*magnitude_source[i]);
313 q->green=ClampToQuantum(QuantumRange*magnitude_source[i]);
318 q->blue=ClampToQuantum(QuantumRange*magnitude_source[i]);
323 q->opacity=ClampToQuantum(QuantumRange*magnitude_source[i]);
328 indexes[x]=ClampToQuantum(QuantumRange*magnitude_source[i]);
333 SetGrayPixelComponent(q,ClampToQuantum(QuantumRange*
334 magnitude_source[i]));
341 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
342 if (status == MagickFalse)
346 for (y=0L; y < (ssize_t) fourier_info->height; y++)
348 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
350 if (q == (PixelPacket *) NULL)
352 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
353 for (x=0L; x < (ssize_t) fourier_info->width; x++)
355 switch (fourier_info->channel)
360 q->red=ClampToQuantum(QuantumRange*phase_source[i]);
365 q->green=ClampToQuantum(QuantumRange*phase_source[i]);
370 q->blue=ClampToQuantum(QuantumRange*phase_source[i]);
375 q->opacity=ClampToQuantum(QuantumRange*phase_source[i]);
380 indexes[x]=ClampToQuantum(QuantumRange*phase_source[i]);
385 SetGrayPixelComponent(q,ClampToQuantum(QuantumRange*phase_source[i]));
392 status=SyncCacheViewAuthenticPixels(phase_view,exception);
393 if (status == MagickFalse)
396 phase_view=DestroyCacheView(phase_view);
397 magnitude_view=DestroyCacheView(magnitude_view);
398 phase_source=(double *) RelinquishMagickMemory(phase_source);
399 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
403 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
404 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
419 register const IndexPacket
422 register const PixelPacket
433 Generate the forward Fourier transform.
435 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
436 fourier_info->width*sizeof(*source));
437 if (source == (double *) NULL)
439 (void) ThrowMagickException(exception,GetMagickModule(),
440 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
443 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
446 image_view=AcquireCacheView(image);
447 for (y=0L; y < (ssize_t) fourier_info->height; y++)
449 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
451 if (p == (const PixelPacket *) NULL)
453 indexes=GetCacheViewVirtualIndexQueue(image_view);
454 for (x=0L; x < (ssize_t) fourier_info->width; x++)
456 switch (fourier_info->channel)
461 source[i]=QuantumScale*GetRedPixelComponent(p);
466 source[i]=QuantumScale*GetGreenPixelComponent(p);
471 source[i]=QuantumScale*GetBluePixelComponent(p);
476 source[i]=QuantumScale*GetOpacityPixelComponent(p);
481 source[i]=QuantumScale*indexes[x];
486 source[i]=QuantumScale*GetGraGrayyPixelComponent(p);
494 image_view=DestroyCacheView(image_view);
495 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
496 fourier_info->center*sizeof(*fourier));
497 if (fourier == (fftw_complex *) NULL)
499 (void) ThrowMagickException(exception,GetMagickModule(),
500 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
501 source=(double *) RelinquishMagickMemory(source);
504 #if defined(MAGICKCORE_OPENMP_SUPPORT)
505 #pragma omp critical (MagickCore_ForwardFourierTransform)
507 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
508 source,fourier,FFTW_ESTIMATE);
509 fftw_execute(fftw_r2c_plan);
510 fftw_destroy_plan(fftw_r2c_plan);
511 source=(double *) RelinquishMagickMemory(source);
513 Normalize Fourier transform.
515 n=(double) fourier_info->width*(double) fourier_info->width;
517 for (y=0L; y < (ssize_t) fourier_info->height; y++)
518 for (x=0L; x < (ssize_t) fourier_info->center; x++)
520 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
529 Generate magnitude and phase (or real and imaginary).
532 if (fourier_info->modulus != MagickFalse)
533 for (y=0L; y < (ssize_t) fourier_info->height; y++)
534 for (x=0L; x < (ssize_t) fourier_info->center; x++)
536 magnitude[i]=cabs(fourier[i]);
537 phase[i]=carg(fourier[i]);
541 for (y=0L; y < (ssize_t) fourier_info->height; y++)
542 for (x=0L; x < (ssize_t) fourier_info->center; x++)
544 magnitude[i]=creal(fourier[i]);
545 phase[i]=cimag(fourier[i]);
548 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
552 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
553 const ChannelType channel,const MagickBooleanType modulus,
554 Image *fourier_image,ExceptionInfo *exception)
572 fourier_info.width=image->columns;
573 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
574 ((image->rows % 2) != 0))
576 extent=image->columns < image->rows ? image->rows : image->columns;
577 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
579 fourier_info.height=fourier_info.width;
580 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
581 fourier_info.channel=channel;
582 fourier_info.modulus=modulus;
583 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
584 fourier_info.center*sizeof(*magnitude));
585 if (magnitude == (double *) NULL)
587 (void) ThrowMagickException(exception,GetMagickModule(),
588 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
591 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
592 fourier_info.center*sizeof(*phase));
593 if (phase == (double *) NULL)
595 (void) ThrowMagickException(exception,GetMagickModule(),
596 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
597 magnitude=(double *) RelinquishMagickMemory(magnitude);
600 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
601 fourier_info.center*sizeof(*fourier));
602 if (fourier == (fftw_complex *) NULL)
604 (void) ThrowMagickException(exception,GetMagickModule(),
605 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
606 phase=(double *) RelinquishMagickMemory(phase);
607 magnitude=(double *) RelinquishMagickMemory(magnitude);
610 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
611 if (status != MagickFalse)
612 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
614 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
615 phase=(double *) RelinquishMagickMemory(phase);
616 magnitude=(double *) RelinquishMagickMemory(magnitude);
621 MagickExport Image *ForwardFourierTransformImage(const Image *image,
622 const MagickBooleanType modulus,ExceptionInfo *exception)
627 fourier_image=NewImageList();
628 #if !defined(MAGICKCORE_FFTW_DELEGATE)
630 (void) ThrowMagickException(exception,GetMagickModule(),
631 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
642 width=image->columns;
643 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
644 ((image->rows % 2) != 0))
646 extent=image->columns < image->rows ? image->rows : image->columns;
647 width=(extent & 0x01) == 1 ? extent+1UL : extent;
649 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
650 if (magnitude_image != (Image *) NULL)
655 magnitude_image->storage_class=DirectClass;
656 magnitude_image->depth=32UL;
657 phase_image=CloneImage(image,width,width,MagickFalse,exception);
658 if (phase_image == (Image *) NULL)
659 magnitude_image=DestroyImage(magnitude_image);
666 phase_image->storage_class=DirectClass;
667 phase_image->depth=32UL;
668 AppendImageToList(&fourier_image,magnitude_image);
669 AppendImageToList(&fourier_image,phase_image);
671 is_gray=IsGrayImage(image,exception);
672 #if defined(MAGICKCORE_OPENMP_SUPPORT)
673 #pragma omp parallel sections
676 #if defined(MAGICKCORE_OPENMP_SUPPORT)
683 if (is_gray != MagickFalse)
684 thread_status=ForwardFourierTransformChannel(image,
685 GrayChannels,modulus,fourier_image,exception);
687 thread_status=ForwardFourierTransformChannel(image,
688 RedChannel,modulus,fourier_image,exception);
689 if (thread_status == MagickFalse)
690 status=thread_status;
692 #if defined(MAGICKCORE_OPENMP_SUPPORT)
699 thread_status=MagickTrue;
700 if (is_gray == MagickFalse)
701 thread_status=ForwardFourierTransformChannel(image,
702 GreenChannel,modulus,fourier_image,exception);
703 if (thread_status == MagickFalse)
704 status=thread_status;
706 #if defined(MAGICKCORE_OPENMP_SUPPORT)
713 thread_status=MagickTrue;
714 if (is_gray == MagickFalse)
715 thread_status=ForwardFourierTransformChannel(image,
716 BlueChannel,modulus,fourier_image,exception);
717 if (thread_status == MagickFalse)
718 status=thread_status;
720 #if defined(MAGICKCORE_OPENMP_SUPPORT)
727 thread_status=MagickTrue;
728 if (image->matte != MagickFalse)
729 thread_status=ForwardFourierTransformChannel(image,
730 OpacityChannel,modulus,fourier_image,exception);
731 if (thread_status == MagickFalse)
732 status=thread_status;
734 #if defined(MAGICKCORE_OPENMP_SUPPORT)
741 thread_status=MagickTrue;
742 if (image->colorspace == CMYKColorspace)
743 thread_status=ForwardFourierTransformChannel(image,
744 IndexChannel,modulus,fourier_image,exception);
745 if (thread_status == MagickFalse)
746 status=thread_status;
749 if (status == MagickFalse)
750 fourier_image=DestroyImageList(fourier_image);
756 return(fourier_image);
760 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
764 % 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 %
768 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
770 % InverseFourierTransformImage() implements the inverse discrete Fourier
771 % transform (DFT) of the image either as a magnitude / phase or real /
772 % imaginary image pair.
774 % The format of the InverseFourierTransformImage method is:
776 % Image *InverseFourierTransformImage(const Image *magnitude_image,
777 % const Image *phase_image,const MagickBooleanType modulus,
778 % ExceptionInfo *exception)
780 % A description of each parameter follows:
782 % o magnitude_image: the magnitude or real image.
784 % o phase_image: the phase or imaginary image.
786 % o modulus: if true, return transform as a magnitude / phase pair
787 % otherwise a real / imaginary image pair.
789 % o exception: return any errors or warnings in this structure.
793 #if defined(MAGICKCORE_FFTW_DELEGATE)
794 static MagickBooleanType InverseQuadrantSwap(const size_t width,
795 const size_t height,const double *source,double *destination)
807 center=(ssize_t) floor((double) width/2.0)+1L;
808 for (y=1L; y < (ssize_t) height; y++)
809 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
810 destination[center*(height-y)-x+width/2L]=source[y*width+x];
811 for (y=0L; y < (ssize_t) height; y++)
812 destination[center*y]=source[y*width+width/2L];
813 for (x=0L; x < center; x++)
814 destination[x]=source[center-x-1L];
815 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
818 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
819 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
820 ExceptionInfo *exception)
835 register const IndexPacket
838 register const PixelPacket
849 Inverse fourier - read image and break down into a double array.
851 magnitude_source=(double *) AcquireQuantumMemory((size_t)
852 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
853 if (magnitude_source == (double *) NULL)
855 (void) ThrowMagickException(exception,GetMagickModule(),
856 ResourceLimitError,"MemoryAllocationFailed","`%s'",
857 magnitude_image->filename);
860 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
861 fourier_info->width*sizeof(*phase_source));
862 if (phase_source == (double *) NULL)
864 (void) ThrowMagickException(exception,GetMagickModule(),
865 ResourceLimitError,"MemoryAllocationFailed","`%s'",
866 magnitude_image->filename);
867 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
871 magnitude_view=AcquireCacheView(magnitude_image);
872 for (y=0L; y < (ssize_t) fourier_info->height; y++)
874 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
876 if (p == (const PixelPacket *) NULL)
878 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
879 for (x=0L; x < (ssize_t) fourier_info->width; x++)
881 switch (fourier_info->channel)
886 magnitude_source[i]=QuantumScale*GetRedPixelComponent(p);
891 magnitude_source[i]=QuantumScale*GetGreenPixelComponent(p);
896 magnitude_source[i]=QuantumScale*GetBluePixelComponent(p);
901 magnitude_source[i]=QuantumScale*GetOpacityPixelComponent(p);
906 magnitude_source[i]=QuantumScale*indexes[x];
911 magnitude_source[i]=QuantumScale*GetGrayPixelComponent(p);
920 phase_view=AcquireCacheView(phase_image);
921 for (y=0L; y < (ssize_t) fourier_info->height; y++)
923 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
925 if (p == (const PixelPacket *) NULL)
927 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
928 for (x=0L; x < (ssize_t) fourier_info->width; x++)
930 switch (fourier_info->channel)
935 phase_source[i]=QuantumScale*GetRedPixelComponent(p);
940 phase_source[i]=QuantumScale*GetGreenPixelComponent(p);
945 phase_source[i]=QuantumScale*GetBluePixelComponent(p);
950 phase_source[i]=QuantumScale*GetOpacityPixelComponent(p);
955 phase_source[i]=QuantumScale*indexes[x];
960 phase_source[i]=QuantumScale*GetGrayPixelComponent(p);
968 if (fourier_info->modulus != MagickFalse)
971 for (y=0L; y < (ssize_t) fourier_info->height; y++)
972 for (x=0L; x < (ssize_t) fourier_info->width; x++)
974 phase_source[i]-=0.5;
975 phase_source[i]*=(2.0*MagickPI);
979 magnitude_view=DestroyCacheView(magnitude_view);
980 phase_view=DestroyCacheView(phase_view);
981 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
982 fourier_info->center*sizeof(*magnitude));
983 if (magnitude == (double *) NULL)
985 (void) ThrowMagickException(exception,GetMagickModule(),
986 ResourceLimitError,"MemoryAllocationFailed","`%s'",
987 magnitude_image->filename);
988 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
989 phase_source=(double *) RelinquishMagickMemory(phase_source);
992 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
993 magnitude_source,magnitude);
994 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
995 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
996 fourier_info->width*sizeof(*phase));
997 if (phase == (double *) NULL)
999 (void) ThrowMagickException(exception,GetMagickModule(),
1000 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1001 magnitude_image->filename);
1002 phase_source=(double *) RelinquishMagickMemory(phase_source);
1003 return(MagickFalse);
1005 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1006 if (status != MagickFalse)
1007 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1008 phase_source,phase);
1009 phase_source=(double *) RelinquishMagickMemory(phase_source);
1014 if (fourier_info->modulus != MagickFalse)
1015 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1016 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1018 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1019 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
1021 fourier[i][0]=magnitude[i]*cos(phase[i]);
1022 fourier[i][1]=magnitude[i]*sin(phase[i]);
1027 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1028 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1030 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1031 fourier[i]=magnitude[i]+I*phase[i];
1033 fourier[i][0]=magnitude[i];
1034 fourier[i][1]=phase[i];
1038 phase=(double *) RelinquishMagickMemory(phase);
1039 magnitude=(double *) RelinquishMagickMemory(magnitude);
1043 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1044 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1055 register IndexPacket
1058 register PixelPacket
1068 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1069 fourier_info->width*sizeof(*source));
1070 if (source == (double *) NULL)
1072 (void) ThrowMagickException(exception,GetMagickModule(),
1073 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1074 return(MagickFalse);
1076 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1077 #pragma omp critical (MagickCore_InverseFourierTransform)
1080 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1081 fourier,source,FFTW_ESTIMATE);
1082 fftw_execute(fftw_c2r_plan);
1083 fftw_destroy_plan(fftw_c2r_plan);
1086 image_view=AcquireCacheView(image);
1087 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1089 if (y >= (ssize_t) image->rows)
1091 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1092 image->columns ? image->columns : fourier_info->width,1UL,exception);
1093 if (q == (PixelPacket *) NULL)
1095 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1096 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1098 switch (fourier_info->channel)
1103 q->red=ClampToQuantum(QuantumRange*source[i]);
1108 q->green=ClampToQuantum(QuantumRange*source[i]);
1113 q->blue=ClampToQuantum(QuantumRange*source[i]);
1116 case OpacityChannel:
1118 q->opacity=ClampToQuantum(QuantumRange*source[i]);
1123 indexes[x]=ClampToQuantum(QuantumRange*source[i]);
1128 SetGrayPixelComponent(q,ClampToQuantum(QuantumRange*source[i]));
1135 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1138 image_view=DestroyCacheView(image_view);
1139 source=(double *) RelinquishMagickMemory(source);
1143 static MagickBooleanType InverseFourierTransformChannel(
1144 const Image *magnitude_image,const Image *phase_image,
1145 const ChannelType channel,const MagickBooleanType modulus,
1146 Image *fourier_image,ExceptionInfo *exception)
1164 fourier_info.width=magnitude_image->columns;
1165 if ((magnitude_image->columns != magnitude_image->rows) ||
1166 ((magnitude_image->columns % 2) != 0) ||
1167 ((magnitude_image->rows % 2) != 0))
1169 extent=magnitude_image->columns < magnitude_image->rows ?
1170 magnitude_image->rows : magnitude_image->columns;
1171 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1173 fourier_info.height=fourier_info.width;
1174 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
1175 fourier_info.channel=channel;
1176 fourier_info.modulus=modulus;
1177 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1178 fourier_info.center*sizeof(*magnitude));
1179 if (magnitude == (double *) NULL)
1181 (void) ThrowMagickException(exception,GetMagickModule(),
1182 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1183 magnitude_image->filename);
1184 return(MagickFalse);
1186 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1187 fourier_info.center*sizeof(*phase));
1188 if (phase == (double *) NULL)
1190 (void) ThrowMagickException(exception,GetMagickModule(),
1191 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1192 magnitude_image->filename);
1193 magnitude=(double *) RelinquishMagickMemory(magnitude);
1194 return(MagickFalse);
1196 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1197 fourier_info.center*sizeof(*fourier));
1198 if (fourier == (fftw_complex *) NULL)
1200 (void) ThrowMagickException(exception,GetMagickModule(),
1201 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1202 magnitude_image->filename);
1203 phase=(double *) RelinquishMagickMemory(phase);
1204 magnitude=(double *) RelinquishMagickMemory(magnitude);
1205 return(MagickFalse);
1207 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1209 if (status != MagickFalse)
1210 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1212 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1213 phase=(double *) RelinquishMagickMemory(phase);
1214 magnitude=(double *) RelinquishMagickMemory(magnitude);
1219 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1220 const Image *phase_image,const MagickBooleanType modulus,
1221 ExceptionInfo *exception)
1226 assert(magnitude_image != (Image *) NULL);
1227 assert(magnitude_image->signature == MagickSignature);
1228 if (magnitude_image->debug != MagickFalse)
1229 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1230 magnitude_image->filename);
1231 if (phase_image == (Image *) NULL)
1233 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1234 "ImageSequenceRequired","`%s'",magnitude_image->filename);
1235 return((Image *) NULL);
1237 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1238 fourier_image=(Image *) NULL;
1240 (void) ThrowMagickException(exception,GetMagickModule(),
1241 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1242 magnitude_image->filename);
1245 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1246 magnitude_image->rows,MagickFalse,exception);
1247 if (fourier_image != (Image *) NULL)
1254 is_gray=IsGrayImage(magnitude_image,exception);
1255 if (is_gray != MagickFalse)
1256 is_gray=IsGrayImage(phase_image,exception);
1257 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1258 #pragma omp parallel sections
1261 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1268 if (is_gray != MagickFalse)
1269 thread_status=InverseFourierTransformChannel(magnitude_image,
1270 phase_image,GrayChannels,modulus,fourier_image,exception);
1272 thread_status=InverseFourierTransformChannel(magnitude_image,
1273 phase_image,RedChannel,modulus,fourier_image,exception);
1274 if (thread_status == MagickFalse)
1275 status=thread_status;
1277 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1284 thread_status=MagickTrue;
1285 if (is_gray == MagickFalse)
1286 thread_status=InverseFourierTransformChannel(magnitude_image,
1287 phase_image,GreenChannel,modulus,fourier_image,exception);
1288 if (thread_status == MagickFalse)
1289 status=thread_status;
1291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1298 thread_status=MagickTrue;
1299 if (is_gray == MagickFalse)
1300 thread_status=InverseFourierTransformChannel(magnitude_image,
1301 phase_image,BlueChannel,modulus,fourier_image,exception);
1302 if (thread_status == MagickFalse)
1303 status=thread_status;
1305 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1312 thread_status=MagickTrue;
1313 if (magnitude_image->matte != MagickFalse)
1314 thread_status=InverseFourierTransformChannel(magnitude_image,
1315 phase_image,OpacityChannel,modulus,fourier_image,exception);
1316 if (thread_status == MagickFalse)
1317 status=thread_status;
1319 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1326 thread_status=MagickTrue;
1327 if (magnitude_image->colorspace == CMYKColorspace)
1328 thread_status=InverseFourierTransformChannel(magnitude_image,
1329 phase_image,IndexChannel,modulus,fourier_image,exception);
1330 if (thread_status == MagickFalse)
1331 status=thread_status;
1334 if (status == MagickFalse)
1335 fourier_image=DestroyImage(fourier_image);
1340 return(fourier_image);