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)
306 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
307 magnitude_source[i]),q);
312 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
313 magnitude_source[i]),q);
318 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
319 magnitude_source[i]),q);
324 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
325 magnitude_source[i]),q);
330 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
331 magnitude_source[i]),q);
336 SetPixelGray(magnitude_image,ClampToQuantum(QuantumRange*
337 magnitude_source[i]),q);
342 q+=GetPixelChannels(magnitude_image);
344 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
345 if (status == MagickFalse)
349 for (y=0L; y < (ssize_t) fourier_info->height; y++)
351 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
353 if (q == (Quantum *) NULL)
355 for (x=0L; x < (ssize_t) fourier_info->width; x++)
357 switch (fourier_info->channel)
362 SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
368 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
374 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
380 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
386 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
392 SetPixelGray(phase_image,ClampToQuantum(QuantumRange*
398 q+=GetPixelChannels(phase_image);
400 status=SyncCacheViewAuthenticPixels(phase_view,exception);
401 if (status == MagickFalse)
404 phase_view=DestroyCacheView(phase_view);
405 magnitude_view=DestroyCacheView(magnitude_view);
406 phase_source=(double *) RelinquishMagickMemory(phase_source);
407 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
411 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
412 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
427 register const Quantum
438 Generate the forward Fourier transform.
440 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
441 fourier_info->width*sizeof(*source));
442 if (source == (double *) NULL)
444 (void) ThrowMagickException(exception,GetMagickModule(),
445 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
448 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
451 image_view=AcquireCacheView(image);
452 for (y=0L; y < (ssize_t) fourier_info->height; y++)
454 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
456 if (p == (const Quantum *) NULL)
458 for (x=0L; x < (ssize_t) fourier_info->width; x++)
460 switch (fourier_info->channel)
465 source[i]=QuantumScale*GetPixelRed(image,p);
470 source[i]=QuantumScale*GetPixelGreen(image,p);
475 source[i]=QuantumScale*GetPixelBlue(image,p);
480 source[i]=QuantumScale*GetPixelBlack(image,p);
485 source[i]=QuantumScale*GetPixelAlpha(image,p);
490 source[i]=QuantumScale*GetPixelGray(image,p);
495 p+=GetPixelChannels(image);
498 image_view=DestroyCacheView(image_view);
499 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
500 fourier_info->center*sizeof(*fourier));
501 if (fourier == (fftw_complex *) NULL)
503 (void) ThrowMagickException(exception,GetMagickModule(),
504 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
505 source=(double *) RelinquishMagickMemory(source);
508 #if defined(MAGICKCORE_OPENMP_SUPPORT)
509 #pragma omp critical (MagickCore_ForwardFourierTransform)
511 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
512 source,fourier,FFTW_ESTIMATE);
513 fftw_execute(fftw_r2c_plan);
514 fftw_destroy_plan(fftw_r2c_plan);
515 source=(double *) RelinquishMagickMemory(source);
517 Normalize Fourier transform.
519 n=(double) fourier_info->width*(double) fourier_info->width;
521 for (y=0L; y < (ssize_t) fourier_info->height; y++)
522 for (x=0L; x < (ssize_t) fourier_info->center; x++)
524 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
533 Generate magnitude and phase (or real and imaginary).
536 if (fourier_info->modulus != MagickFalse)
537 for (y=0L; y < (ssize_t) fourier_info->height; y++)
538 for (x=0L; x < (ssize_t) fourier_info->center; x++)
540 magnitude[i]=cabs(fourier[i]);
541 phase[i]=carg(fourier[i]);
545 for (y=0L; y < (ssize_t) fourier_info->height; y++)
546 for (x=0L; x < (ssize_t) fourier_info->center; x++)
548 magnitude[i]=creal(fourier[i]);
549 phase[i]=cimag(fourier[i]);
552 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
556 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
557 const ChannelType channel,const MagickBooleanType modulus,
558 Image *fourier_image,ExceptionInfo *exception)
576 fourier_info.width=image->columns;
577 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
578 ((image->rows % 2) != 0))
580 extent=image->columns < image->rows ? image->rows : image->columns;
581 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
583 fourier_info.height=fourier_info.width;
584 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
585 fourier_info.channel=channel;
586 fourier_info.modulus=modulus;
587 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
588 fourier_info.center*sizeof(*magnitude));
589 if (magnitude == (double *) NULL)
591 (void) ThrowMagickException(exception,GetMagickModule(),
592 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
595 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
596 fourier_info.center*sizeof(*phase));
597 if (phase == (double *) NULL)
599 (void) ThrowMagickException(exception,GetMagickModule(),
600 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
601 magnitude=(double *) RelinquishMagickMemory(magnitude);
604 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
605 fourier_info.center*sizeof(*fourier));
606 if (fourier == (fftw_complex *) NULL)
608 (void) ThrowMagickException(exception,GetMagickModule(),
609 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
610 phase=(double *) RelinquishMagickMemory(phase);
611 magnitude=(double *) RelinquishMagickMemory(magnitude);
614 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
615 if (status != MagickFalse)
616 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
618 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
619 phase=(double *) RelinquishMagickMemory(phase);
620 magnitude=(double *) RelinquishMagickMemory(magnitude);
625 MagickExport Image *ForwardFourierTransformImage(const Image *image,
626 const MagickBooleanType modulus,ExceptionInfo *exception)
631 fourier_image=NewImageList();
632 #if !defined(MAGICKCORE_FFTW_DELEGATE)
634 (void) ThrowMagickException(exception,GetMagickModule(),
635 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
646 width=image->columns;
647 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
648 ((image->rows % 2) != 0))
650 extent=image->columns < image->rows ? image->rows : image->columns;
651 width=(extent & 0x01) == 1 ? extent+1UL : extent;
653 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
654 if (magnitude_image != (Image *) NULL)
659 magnitude_image->storage_class=DirectClass;
660 magnitude_image->depth=32UL;
661 phase_image=CloneImage(image,width,width,MagickFalse,exception);
662 if (phase_image == (Image *) NULL)
663 magnitude_image=DestroyImage(magnitude_image);
670 phase_image->storage_class=DirectClass;
671 phase_image->depth=32UL;
672 AppendImageToList(&fourier_image,magnitude_image);
673 AppendImageToList(&fourier_image,phase_image);
675 is_gray=IsImageGray(image,exception);
676 #if defined(MAGICKCORE_OPENMP_SUPPORT)
677 #pragma omp parallel sections
680 #if defined(MAGICKCORE_OPENMP_SUPPORT)
687 if (is_gray != MagickFalse)
688 thread_status=ForwardFourierTransformChannel(image,
689 GrayChannels,modulus,fourier_image,exception);
691 thread_status=ForwardFourierTransformChannel(image,
692 RedChannel,modulus,fourier_image,exception);
693 if (thread_status == MagickFalse)
694 status=thread_status;
696 #if defined(MAGICKCORE_OPENMP_SUPPORT)
703 thread_status=MagickTrue;
704 if (is_gray == MagickFalse)
705 thread_status=ForwardFourierTransformChannel(image,
706 GreenChannel,modulus,fourier_image,exception);
707 if (thread_status == MagickFalse)
708 status=thread_status;
710 #if defined(MAGICKCORE_OPENMP_SUPPORT)
717 thread_status=MagickTrue;
718 if (is_gray == MagickFalse)
719 thread_status=ForwardFourierTransformChannel(image,
720 BlueChannel,modulus,fourier_image,exception);
721 if (thread_status == MagickFalse)
722 status=thread_status;
724 #if defined(MAGICKCORE_OPENMP_SUPPORT)
731 thread_status=MagickTrue;
732 if (image->colorspace == CMYKColorspace)
733 thread_status=ForwardFourierTransformChannel(image,
734 BlackChannel,modulus,fourier_image,exception);
735 if (thread_status == MagickFalse)
736 status=thread_status;
738 #if defined(MAGICKCORE_OPENMP_SUPPORT)
745 thread_status=MagickTrue;
746 if (image->matte != MagickFalse)
747 thread_status=ForwardFourierTransformChannel(image,
748 AlphaChannel,modulus,fourier_image,exception);
749 if (thread_status == MagickFalse)
750 status=thread_status;
753 if (status == MagickFalse)
754 fourier_image=DestroyImageList(fourier_image);
760 return(fourier_image);
764 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
768 % 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 %
772 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
774 % InverseFourierTransformImage() implements the inverse discrete Fourier
775 % transform (DFT) of the image either as a magnitude / phase or real /
776 % imaginary image pair.
778 % The format of the InverseFourierTransformImage method is:
780 % Image *InverseFourierTransformImage(const Image *magnitude_image,
781 % const Image *phase_image,const MagickBooleanType modulus,
782 % ExceptionInfo *exception)
784 % A description of each parameter follows:
786 % o magnitude_image: the magnitude or real image.
788 % o phase_image: the phase or imaginary image.
790 % o modulus: if true, return transform as a magnitude / phase pair
791 % otherwise a real / imaginary image pair.
793 % o exception: return any errors or warnings in this structure.
797 #if defined(MAGICKCORE_FFTW_DELEGATE)
798 static MagickBooleanType InverseQuadrantSwap(const size_t width,
799 const size_t height,const double *source,double *destination)
811 center=(ssize_t) floor((double) width/2.0)+1L;
812 for (y=1L; y < (ssize_t) height; y++)
813 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
814 destination[center*(height-y)-x+width/2L]=source[y*width+x];
815 for (y=0L; y < (ssize_t) height; y++)
816 destination[center*y]=source[y*width+width/2L];
817 for (x=0L; x < center; x++)
818 destination[x]=source[center-x-1L];
819 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
822 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
823 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
824 ExceptionInfo *exception)
839 register const Quantum
850 Inverse fourier - read image and break down into a double array.
852 magnitude_source=(double *) AcquireQuantumMemory((size_t)
853 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
854 if (magnitude_source == (double *) NULL)
856 (void) ThrowMagickException(exception,GetMagickModule(),
857 ResourceLimitError,"MemoryAllocationFailed","`%s'",
858 magnitude_image->filename);
861 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
862 fourier_info->width*sizeof(*phase_source));
863 if (phase_source == (double *) NULL)
865 (void) ThrowMagickException(exception,GetMagickModule(),
866 ResourceLimitError,"MemoryAllocationFailed","`%s'",
867 magnitude_image->filename);
868 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
872 magnitude_view=AcquireCacheView(magnitude_image);
873 for (y=0L; y < (ssize_t) fourier_info->height; y++)
875 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
877 if (p == (const Quantum *) NULL)
879 for (x=0L; x < (ssize_t) fourier_info->width; x++)
881 switch (fourier_info->channel)
886 magnitude_source[i]=QuantumScale*GetPixelRed(magnitude_image,p);
891 magnitude_source[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
896 magnitude_source[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
901 magnitude_source[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
906 magnitude_source[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
911 magnitude_source[i]=QuantumScale*GetPixelGray(magnitude_image,p);
916 p+=GetPixelChannels(magnitude_image);
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 Quantum *) NULL)
927 for (x=0L; x < (ssize_t) fourier_info->width; x++)
929 switch (fourier_info->channel)
934 phase_source[i]=QuantumScale*GetPixelRed(phase_image,p);
939 phase_source[i]=QuantumScale*GetPixelGreen(phase_image,p);
944 phase_source[i]=QuantumScale*GetPixelBlue(phase_image,p);
949 phase_source[i]=QuantumScale*GetPixelBlack(phase_image,p);
954 phase_source[i]=QuantumScale*GetPixelAlpha(phase_image,p);
959 phase_source[i]=QuantumScale*GetPixelGray(phase_image,p);
964 p+=GetPixelChannels(phase_image);
967 if (fourier_info->modulus != MagickFalse)
970 for (y=0L; y < (ssize_t) fourier_info->height; y++)
971 for (x=0L; x < (ssize_t) fourier_info->width; x++)
973 phase_source[i]-=0.5;
974 phase_source[i]*=(2.0*MagickPI);
978 magnitude_view=DestroyCacheView(magnitude_view);
979 phase_view=DestroyCacheView(phase_view);
980 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
981 fourier_info->center*sizeof(*magnitude));
982 if (magnitude == (double *) NULL)
984 (void) ThrowMagickException(exception,GetMagickModule(),
985 ResourceLimitError,"MemoryAllocationFailed","`%s'",
986 magnitude_image->filename);
987 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
988 phase_source=(double *) RelinquishMagickMemory(phase_source);
991 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
992 magnitude_source,magnitude);
993 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
994 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
995 fourier_info->width*sizeof(*phase));
996 if (phase == (double *) NULL)
998 (void) ThrowMagickException(exception,GetMagickModule(),
999 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1000 magnitude_image->filename);
1001 phase_source=(double *) RelinquishMagickMemory(phase_source);
1002 return(MagickFalse);
1004 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1005 if (status != MagickFalse)
1006 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1007 phase_source,phase);
1008 phase_source=(double *) RelinquishMagickMemory(phase_source);
1013 if (fourier_info->modulus != MagickFalse)
1014 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1015 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1017 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1018 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
1020 fourier[i][0]=magnitude[i]*cos(phase[i]);
1021 fourier[i][1]=magnitude[i]*sin(phase[i]);
1026 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1027 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1029 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1030 fourier[i]=magnitude[i]+I*phase[i];
1032 fourier[i][0]=magnitude[i];
1033 fourier[i][1]=phase[i];
1037 phase=(double *) RelinquishMagickMemory(phase);
1038 magnitude=(double *) RelinquishMagickMemory(magnitude);
1042 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1043 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1064 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1065 fourier_info->width*sizeof(*source));
1066 if (source == (double *) NULL)
1068 (void) ThrowMagickException(exception,GetMagickModule(),
1069 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1070 return(MagickFalse);
1072 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1073 #pragma omp critical (MagickCore_InverseFourierTransform)
1076 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1077 fourier,source,FFTW_ESTIMATE);
1078 fftw_execute(fftw_c2r_plan);
1079 fftw_destroy_plan(fftw_c2r_plan);
1082 image_view=AcquireCacheView(image);
1083 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1085 if (y >= (ssize_t) image->rows)
1087 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1088 image->columns ? image->columns : fourier_info->width,1UL,exception);
1089 if (q == (Quantum *) NULL)
1091 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1093 switch (fourier_info->channel)
1098 SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
1103 SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
1108 SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
1113 SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
1118 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
1123 SetPixelGray(image,ClampToQuantum(QuantumRange*source[i]),q);
1128 q+=GetPixelChannels(image);
1130 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1133 image_view=DestroyCacheView(image_view);
1134 source=(double *) RelinquishMagickMemory(source);
1138 static MagickBooleanType InverseFourierTransformChannel(
1139 const Image *magnitude_image,const Image *phase_image,
1140 const ChannelType channel,const MagickBooleanType modulus,
1141 Image *fourier_image,ExceptionInfo *exception)
1159 fourier_info.width=magnitude_image->columns;
1160 if ((magnitude_image->columns != magnitude_image->rows) ||
1161 ((magnitude_image->columns % 2) != 0) ||
1162 ((magnitude_image->rows % 2) != 0))
1164 extent=magnitude_image->columns < magnitude_image->rows ?
1165 magnitude_image->rows : magnitude_image->columns;
1166 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1168 fourier_info.height=fourier_info.width;
1169 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
1170 fourier_info.channel=channel;
1171 fourier_info.modulus=modulus;
1172 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1173 fourier_info.center*sizeof(*magnitude));
1174 if (magnitude == (double *) NULL)
1176 (void) ThrowMagickException(exception,GetMagickModule(),
1177 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1178 magnitude_image->filename);
1179 return(MagickFalse);
1181 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1182 fourier_info.center*sizeof(*phase));
1183 if (phase == (double *) NULL)
1185 (void) ThrowMagickException(exception,GetMagickModule(),
1186 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1187 magnitude_image->filename);
1188 magnitude=(double *) RelinquishMagickMemory(magnitude);
1189 return(MagickFalse);
1191 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1192 fourier_info.center*sizeof(*fourier));
1193 if (fourier == (fftw_complex *) NULL)
1195 (void) ThrowMagickException(exception,GetMagickModule(),
1196 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1197 magnitude_image->filename);
1198 phase=(double *) RelinquishMagickMemory(phase);
1199 magnitude=(double *) RelinquishMagickMemory(magnitude);
1200 return(MagickFalse);
1202 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1204 if (status != MagickFalse)
1205 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1207 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1208 phase=(double *) RelinquishMagickMemory(phase);
1209 magnitude=(double *) RelinquishMagickMemory(magnitude);
1214 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1215 const Image *phase_image,const MagickBooleanType modulus,
1216 ExceptionInfo *exception)
1221 assert(magnitude_image != (Image *) NULL);
1222 assert(magnitude_image->signature == MagickSignature);
1223 if (magnitude_image->debug != MagickFalse)
1224 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1225 magnitude_image->filename);
1226 if (phase_image == (Image *) NULL)
1228 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1229 "ImageSequenceRequired","`%s'",magnitude_image->filename);
1230 return((Image *) NULL);
1232 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1233 fourier_image=(Image *) NULL;
1235 (void) ThrowMagickException(exception,GetMagickModule(),
1236 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1237 magnitude_image->filename);
1240 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1241 magnitude_image->rows,MagickFalse,exception);
1242 if (fourier_image != (Image *) NULL)
1249 is_gray=IsImageGray(magnitude_image,exception);
1250 if (is_gray != MagickFalse)
1251 is_gray=IsImageGray(phase_image,exception);
1252 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1253 #pragma omp parallel sections
1256 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1263 if (is_gray != MagickFalse)
1264 thread_status=InverseFourierTransformChannel(magnitude_image,
1265 phase_image,GrayChannels,modulus,fourier_image,exception);
1267 thread_status=InverseFourierTransformChannel(magnitude_image,
1268 phase_image,RedChannel,modulus,fourier_image,exception);
1269 if (thread_status == MagickFalse)
1270 status=thread_status;
1272 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1279 thread_status=MagickTrue;
1280 if (is_gray == MagickFalse)
1281 thread_status=InverseFourierTransformChannel(magnitude_image,
1282 phase_image,GreenChannel,modulus,fourier_image,exception);
1283 if (thread_status == MagickFalse)
1284 status=thread_status;
1286 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1293 thread_status=MagickTrue;
1294 if (is_gray == MagickFalse)
1295 thread_status=InverseFourierTransformChannel(magnitude_image,
1296 phase_image,BlueChannel,modulus,fourier_image,exception);
1297 if (thread_status == MagickFalse)
1298 status=thread_status;
1300 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1307 thread_status=MagickTrue;
1308 if (magnitude_image->colorspace == CMYKColorspace)
1309 thread_status=InverseFourierTransformChannel(magnitude_image,
1310 phase_image,BlackChannel,modulus,fourier_image,exception);
1311 if (thread_status == MagickFalse)
1312 status=thread_status;
1314 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1321 thread_status=MagickTrue;
1322 if (magnitude_image->matte != MagickFalse)
1323 thread_status=InverseFourierTransformChannel(magnitude_image,
1324 phase_image,AlphaChannel,modulus,fourier_image,exception);
1325 if (thread_status == MagickFalse)
1326 status=thread_status;
1329 if (status == MagickFalse)
1330 fourier_image=DestroyImage(fourier_image);
1335 return(fourier_image);