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-2010 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)
61 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
62 #define carg(z) (atan2(z[1],z[0]))
63 #define creal(z) (z[0])
64 #define cimag(z) (z[1])
72 typedef struct _FourierInfo
89 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93 % 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 %
97 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 % ForwardFourierTransformImage() implements the discrete Fourier transform
100 % (DFT) of the image either as a magnitude / phase or real / imaginary image
103 % The format of the ForwadFourierTransformImage method is:
105 % Image *ForwardFourierTransformImage(const Image *image,
106 % const MagickBooleanType modulus,ExceptionInfo *exception)
108 % A description of each parameter follows:
110 % o image: the image.
112 % o modulus: if true, return as transform as a magnitude / phase pair
113 % otherwise a real / imaginary image pair.
115 % o exception: return any errors or warnings in this structure.
119 #if defined(MAGICKCORE_FFTW_DELEGATE)
121 static MagickBooleanType RollFourier(const size_t width,
122 const size_t height,const ssize_t x_offset,const ssize_t y_offset,
138 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
140 roll=(double *) AcquireQuantumMemory((size_t) width,height*sizeof(*roll));
141 if (roll == (double *) NULL)
144 for (y=0L; y < (ssize_t) height; y++)
147 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
149 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
151 for (x=0L; x < (ssize_t) width; x++)
154 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
156 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
158 roll[v*width+u]=fourier[i++];
161 (void) CopyMagickMemory(fourier,roll,width*height*sizeof(*roll));
162 roll=(double *) RelinquishMagickMemory(roll);
166 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
167 const size_t height,double *source,double *destination)
182 center=(ssize_t) floor((double) width/2L)+1L;
183 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
184 if (status == MagickFalse)
186 for (y=0L; y < (ssize_t) height; y++)
187 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
188 destination[width*y+x+width/2L]=source[center*y+x];
189 for (y=1; y < (ssize_t) height; y++)
190 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
191 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
192 for (x=0L; x < (ssize_t) (width/2L); x++)
193 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
197 static void CorrectPhaseLHS(const size_t width,
198 const size_t height,double *fourier)
206 for (y=0L; y < (ssize_t) height; y++)
207 for (x=0L; x < (ssize_t) (width/2L); x++)
208 fourier[y*width+x]*=(-1.0);
211 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
212 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
242 magnitude_image=GetFirstImageInList(image);
243 phase_image=GetNextImageInList(image);
244 if (phase_image == (Image *) NULL)
246 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
247 "ImageSequenceRequired","`%s'",image->filename);
251 Create "Fourier Transform" image from constituent arrays.
253 magnitude_source=(double *) AcquireQuantumMemory((size_t)
254 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
255 if (magnitude_source == (double *) NULL)
257 (void) ResetMagickMemory(magnitude_source,0,fourier_info->width*
258 fourier_info->height*sizeof(*magnitude_source));
259 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
260 fourier_info->width*sizeof(*phase_source));
261 if (magnitude_source == (double *) NULL)
263 (void) ThrowMagickException(exception,GetMagickModule(),
264 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
265 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
268 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
269 magnitude,magnitude_source);
270 if (status != MagickFalse)
271 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
273 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
274 if (fourier_info->modulus != MagickFalse)
277 for (y=0L; y < (ssize_t) fourier_info->height; y++)
278 for (x=0L; x < (ssize_t) fourier_info->width; x++)
280 phase_source[i]/=(2.0*MagickPI);
281 phase_source[i]+=0.5;
285 magnitude_view=AcquireCacheView(magnitude_image);
286 phase_view=AcquireCacheView(phase_image);
288 for (y=0L; y < (ssize_t) fourier_info->height; y++)
290 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
292 if (q == (PixelPacket *) NULL)
294 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
295 for (x=0L; x < (ssize_t) fourier_info->width; x++)
297 switch (fourier_info->channel)
302 q->red=ClampToQuantum(QuantumRange*magnitude_source[i]);
307 q->green=ClampToQuantum(QuantumRange*magnitude_source[i]);
312 q->blue=ClampToQuantum(QuantumRange*magnitude_source[i]);
317 q->opacity=ClampToQuantum(QuantumRange*magnitude_source[i]);
322 indexes[x]=ClampToQuantum(QuantumRange*magnitude_source[i]);
327 q->red=ClampToQuantum(QuantumRange*magnitude_source[i]);
336 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
337 if (status == MagickFalse)
341 for (y=0L; y < (ssize_t) fourier_info->height; y++)
343 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
345 if (q == (PixelPacket *) NULL)
347 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
348 for (x=0L; x < (ssize_t) fourier_info->width; x++)
350 switch (fourier_info->channel)
355 q->red=ClampToQuantum(QuantumRange*phase_source[i]);
360 q->green=ClampToQuantum(QuantumRange*phase_source[i]);
365 q->blue=ClampToQuantum(QuantumRange*phase_source[i]);
370 q->opacity=ClampToQuantum(QuantumRange*phase_source[i]);
375 indexes[x]=ClampToQuantum(QuantumRange*phase_source[i]);
380 q->red=ClampToQuantum(QuantumRange*phase_source[i]);
389 status=SyncCacheViewAuthenticPixels(phase_view,exception);
390 if (status == MagickFalse)
393 phase_view=DestroyCacheView(phase_view);
394 magnitude_view=DestroyCacheView(magnitude_view);
395 phase_source=(double *) RelinquishMagickMemory(phase_source);
396 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
400 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
401 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
419 register const IndexPacket
422 register const PixelPacket
430 Generate the forward Fourier transform.
432 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
433 fourier_info->width*sizeof(*source));
434 if (source == (double *) NULL)
436 (void) ThrowMagickException(exception,GetMagickModule(),
437 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
440 ResetMagickMemory(source,0,fourier_info->width*fourier_info->height*
443 image_view=AcquireCacheView(image);
444 for (y=0L; y < (ssize_t) fourier_info->height; y++)
446 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
448 if (p == (const PixelPacket *) NULL)
450 indexes=GetCacheViewVirtualIndexQueue(image_view);
451 for (x=0L; x < (ssize_t) fourier_info->width; x++)
453 switch (fourier_info->channel)
458 source[i]=QuantumScale*GetRedPixelComponent(p);
463 source[i]=QuantumScale*GetGreenPixelComponent(p);
468 source[i]=QuantumScale*GetBluePixelComponent(p);
473 source[i]=QuantumScale*GetOpacityPixelComponent(p);
478 source[i]=QuantumScale*indexes[x];
483 source[i]=QuantumScale*GetRedPixelComponent(p);
491 image_view=DestroyCacheView(image_view);
492 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
493 fourier_info->center*sizeof(*fourier));
494 if (fourier == (fftw_complex *) NULL)
496 (void) ThrowMagickException(exception,GetMagickModule(),
497 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
498 source=(double *) RelinquishMagickMemory(source);
501 #if defined(MAGICKCORE_OPENMP_SUPPORT)
502 #pragma omp critical (MagickCore_ForwardFourierTransform)
504 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
505 source,fourier,FFTW_ESTIMATE);
506 fftw_execute(fftw_r2c_plan);
507 fftw_destroy_plan(fftw_r2c_plan);
508 source=(double *) RelinquishMagickMemory(source);
510 Normalize Fourier transform.
512 n=(double) fourier_info->width*(double) fourier_info->width;
514 for (y=0L; y < (ssize_t) fourier_info->height; y++)
515 for (x=0L; x < (ssize_t) fourier_info->center; x++)
517 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
526 Generate magnitude and phase (or real and imaginary).
529 if (fourier_info->modulus != MagickFalse)
530 for (y=0L; y < (ssize_t) fourier_info->height; y++)
531 for (x=0L; x < (ssize_t) fourier_info->center; x++)
533 magnitude[i]=cabs(fourier[i]);
534 phase[i]=carg(fourier[i]);
538 for (y=0L; y < (ssize_t) fourier_info->height; y++)
539 for (x=0L; x < (ssize_t) fourier_info->center; x++)
541 magnitude[i]=creal(fourier[i]);
542 phase[i]=cimag(fourier[i]);
545 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
549 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
550 const ChannelType channel,const MagickBooleanType modulus,
551 Image *fourier_image,ExceptionInfo *exception)
569 fourier_info.width=image->columns;
570 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
571 ((image->rows % 2) != 0))
573 extent=image->columns < image->rows ? image->rows : image->columns;
574 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
576 fourier_info.height=fourier_info.width;
577 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
578 fourier_info.channel=channel;
579 fourier_info.modulus=modulus;
580 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
581 fourier_info.center*sizeof(*magnitude));
582 if (magnitude == (double *) NULL)
584 (void) ThrowMagickException(exception,GetMagickModule(),
585 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
588 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
589 fourier_info.center*sizeof(*phase));
590 if (phase == (double *) NULL)
592 (void) ThrowMagickException(exception,GetMagickModule(),
593 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
594 magnitude=(double *) RelinquishMagickMemory(magnitude);
597 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
598 fourier_info.center*sizeof(*fourier));
599 if (fourier == (fftw_complex *) NULL)
601 (void) ThrowMagickException(exception,GetMagickModule(),
602 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
603 phase=(double *) RelinquishMagickMemory(phase);
604 magnitude=(double *) RelinquishMagickMemory(magnitude);
607 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
608 if (status != MagickFalse)
609 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
611 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
612 phase=(double *) RelinquishMagickMemory(phase);
613 magnitude=(double *) RelinquishMagickMemory(magnitude);
618 MagickExport Image *ForwardFourierTransformImage(const Image *image,
619 const MagickBooleanType modulus,ExceptionInfo *exception)
624 fourier_image=NewImageList();
625 #if !defined(MAGICKCORE_FFTW_DELEGATE)
627 (void) ThrowMagickException(exception,GetMagickModule(),
628 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
639 width=image->columns;
640 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
641 ((image->rows % 2) != 0))
643 extent=image->columns < image->rows ? image->rows : image->columns;
644 width=(extent & 0x01) == 1 ? extent+1UL : extent;
646 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
647 if (magnitude_image != (Image *) NULL)
652 magnitude_image->storage_class=DirectClass;
653 magnitude_image->depth=32UL;
654 phase_image=CloneImage(image,width,width,MagickFalse,exception);
655 if (phase_image == (Image *) NULL)
656 magnitude_image=DestroyImage(magnitude_image);
663 phase_image->storage_class=DirectClass;
664 phase_image->depth=32UL;
665 AppendImageToList(&fourier_image,magnitude_image);
666 AppendImageToList(&fourier_image,phase_image);
668 is_gray=IsGrayImage(image,exception);
669 #if defined(MAGICKCORE_OPENMP_SUPPORT)
670 #pragma omp parallel sections
673 #if defined(MAGICKCORE_OPENMP_SUPPORT)
680 if (is_gray != MagickFalse)
681 thread_status=ForwardFourierTransformChannel(image,
682 GrayChannels,modulus,fourier_image,exception);
684 thread_status=ForwardFourierTransformChannel(image,
685 RedChannel,modulus,fourier_image,exception);
686 if (thread_status == MagickFalse)
687 status=thread_status;
689 #if defined(MAGICKCORE_OPENMP_SUPPORT)
696 thread_status=MagickTrue;
697 if (is_gray == MagickFalse)
698 thread_status=ForwardFourierTransformChannel(image,
699 GreenChannel,modulus,fourier_image,exception);
700 if (thread_status == MagickFalse)
701 status=thread_status;
703 #if defined(MAGICKCORE_OPENMP_SUPPORT)
710 thread_status=MagickTrue;
711 if (is_gray == MagickFalse)
712 thread_status=ForwardFourierTransformChannel(image,
713 BlueChannel,modulus,fourier_image,exception);
714 if (thread_status == MagickFalse)
715 status=thread_status;
717 #if defined(MAGICKCORE_OPENMP_SUPPORT)
724 thread_status=MagickTrue;
725 if (image->matte != MagickFalse)
726 thread_status=ForwardFourierTransformChannel(image,
727 OpacityChannel,modulus,fourier_image,exception);
728 if (thread_status == MagickFalse)
729 status=thread_status;
731 #if defined(MAGICKCORE_OPENMP_SUPPORT)
738 thread_status=MagickTrue;
739 if (image->colorspace == CMYKColorspace)
740 thread_status=ForwardFourierTransformChannel(image,
741 IndexChannel,modulus,fourier_image,exception);
742 if (thread_status == MagickFalse)
743 status=thread_status;
746 if (status == MagickFalse)
747 fourier_image=DestroyImageList(fourier_image);
753 return(fourier_image);
757 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
761 % 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 %
765 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
767 % InverseFourierTransformImage() implements the inverse discrete Fourier
768 % transform (DFT) of the image either as a magnitude / phase or real /
769 % imaginary image pair.
771 % The format of the InverseFourierTransformImage method is:
773 % Image *InverseFourierTransformImage(const Image *magnitude_image,
774 % const Image *phase_image,const MagickBooleanType modulus,
775 % ExceptionInfo *exception)
777 % A description of each parameter follows:
779 % o magnitude_image: the magnitude or real image.
781 % o phase_image: the phase or imaginary image.
783 % o modulus: if true, return transform as a magnitude / phase pair
784 % otherwise a real / imaginary image pair.
786 % o exception: return any errors or warnings in this structure.
790 #if defined(MAGICKCORE_FFTW_DELEGATE)
791 static MagickBooleanType InverseQuadrantSwap(const size_t width,
792 const size_t height,const double *source,double *destination)
804 center=(ssize_t) floor((double) width/2.0)+1L;
805 for (y=1L; y < (ssize_t) height; y++)
806 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
807 destination[center*(height-y)-x+width/2L]=source[y*width+x];
808 for (y=0L; y < (ssize_t) height; y++)
809 destination[center*y]=source[y*width+width/2L];
810 for (x=0L; x < center; x++)
811 destination[x]=source[center-x-1L];
812 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
815 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
816 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
817 ExceptionInfo *exception)
835 register const IndexPacket
838 register const PixelPacket
846 Inverse fourier - read image and break down into a double array.
848 magnitude_source=(double *) AcquireQuantumMemory((size_t)
849 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
850 if (magnitude_source == (double *) NULL)
852 (void) ThrowMagickException(exception,GetMagickModule(),
853 ResourceLimitError,"MemoryAllocationFailed","`%s'",
854 magnitude_image->filename);
857 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
858 fourier_info->height*sizeof(*phase_source));
859 if (phase_source == (double *) NULL)
861 (void) ThrowMagickException(exception,GetMagickModule(),
862 ResourceLimitError,"MemoryAllocationFailed","`%s'",
863 magnitude_image->filename);
864 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
868 magnitude_view=AcquireCacheView(magnitude_image);
869 for (y=0L; y < (ssize_t) fourier_info->height; y++)
871 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
873 if (p == (const PixelPacket *) NULL)
875 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
876 for (x=0L; x < (ssize_t) fourier_info->width; x++)
878 switch (fourier_info->channel)
883 magnitude_source[i]=QuantumScale*GetRedPixelComponent(p);
888 magnitude_source[i]=QuantumScale*GetGreenPixelComponent(p);
893 magnitude_source[i]=QuantumScale*GetBluePixelComponent(p);
898 magnitude_source[i]=QuantumScale*GetOpacityPixelComponent(p);
903 magnitude_source[i]=QuantumScale*indexes[x];
908 magnitude_source[i]=QuantumScale*GetRedPixelComponent(p);
917 phase_view=AcquireCacheView(phase_image);
918 for (y=0L; y < (ssize_t) fourier_info->height; y++)
920 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
922 if (p == (const PixelPacket *) NULL)
924 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
925 for (x=0L; x < (ssize_t) fourier_info->width; x++)
927 switch (fourier_info->channel)
932 phase_source[i]=QuantumScale*GetRedPixelComponent(p);
937 phase_source[i]=QuantumScale*GetGreenPixelComponent(p);
942 phase_source[i]=QuantumScale*GetBluePixelComponent(p);
947 phase_source[i]=QuantumScale*GetOpacityPixelComponent(p);
952 phase_source[i]=QuantumScale*indexes[x];
957 phase_source[i]=QuantumScale*GetRedPixelComponent(p);
965 if (fourier_info->modulus != MagickFalse)
968 for (y=0L; y < (ssize_t) fourier_info->height; y++)
969 for (x=0L; x < (ssize_t) fourier_info->width; x++)
971 phase_source[i]-=0.5;
972 phase_source[i]*=(2.0*MagickPI);
976 magnitude_view=DestroyCacheView(magnitude_view);
977 phase_view=DestroyCacheView(phase_view);
978 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
979 fourier_info->center*sizeof(*magnitude));
980 if (magnitude == (double *) NULL)
982 (void) ThrowMagickException(exception,GetMagickModule(),
983 ResourceLimitError,"MemoryAllocationFailed","`%s'",
984 magnitude_image->filename);
985 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
986 phase_source=(double *) RelinquishMagickMemory(phase_source);
989 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
990 magnitude_source,magnitude);
991 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
992 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
993 fourier_info->height*sizeof(*phase));
994 if (phase == (double *) NULL)
996 (void) ThrowMagickException(exception,GetMagickModule(),
997 ResourceLimitError,"MemoryAllocationFailed","`%s'",
998 magnitude_image->filename);
999 phase_source=(double *) RelinquishMagickMemory(phase_source);
1000 return(MagickFalse);
1002 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1003 if (status != MagickFalse)
1004 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1005 phase_source,phase);
1006 phase_source=(double *) RelinquishMagickMemory(phase_source);
1011 if (fourier_info->modulus != MagickFalse)
1012 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1013 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1015 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1016 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
1018 fourier[i][0]=magnitude[i]*cos(phase[i]);
1019 fourier[i][1]=magnitude[i]*sin(phase[i]);
1024 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1025 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1027 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1028 fourier[i]=magnitude[i]+I*phase[i];
1030 fourier[i][0]=magnitude[i];
1031 fourier[i][1]=phase[i];
1035 phase=(double *) RelinquishMagickMemory(phase);
1036 magnitude=(double *) RelinquishMagickMemory(magnitude);
1040 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1041 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1055 register IndexPacket
1062 register PixelPacket
1065 source=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
1066 fourier_info->height*sizeof(double));
1067 if (source == (double *) NULL)
1069 (void) ThrowMagickException(exception,GetMagickModule(),
1070 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1071 return(MagickFalse);
1073 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1074 #pragma omp critical (MagickCore_InverseFourierTransform)
1077 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1078 fourier,source,FFTW_ESTIMATE);
1079 fftw_execute(fftw_c2r_plan);
1080 fftw_destroy_plan(fftw_c2r_plan);
1083 image_view=AcquireCacheView(image);
1084 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1086 if (y >= (ssize_t) image->rows)
1088 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1089 image->columns ? image->columns : fourier_info->width,1UL,exception);
1090 if (q == (PixelPacket *) NULL)
1092 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1093 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1095 switch (fourier_info->channel)
1100 q->red=ClampToQuantum(QuantumRange*source[i]);
1105 q->green=ClampToQuantum(QuantumRange*source[i]);
1110 q->blue=ClampToQuantum(QuantumRange*source[i]);
1113 case OpacityChannel:
1115 q->opacity=ClampToQuantum(QuantumRange*source[i]);
1120 indexes[x]=ClampToQuantum(QuantumRange*source[i]);
1125 q->red=ClampToQuantum(QuantumRange*source[i]);
1134 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1137 image_view=DestroyCacheView(image_view);
1138 source=(double *) RelinquishMagickMemory(source);
1142 static MagickBooleanType InverseFourierTransformChannel(
1143 const Image *magnitude_image,const Image *phase_image,
1144 const ChannelType channel,const MagickBooleanType modulus,
1145 Image *fourier_image,ExceptionInfo *exception)
1163 fourier_info.width=magnitude_image->columns;
1164 if ((magnitude_image->columns != magnitude_image->rows) ||
1165 ((magnitude_image->columns % 2) != 0) ||
1166 ((magnitude_image->rows % 2) != 0))
1168 extent=magnitude_image->columns < magnitude_image->rows ?
1169 magnitude_image->rows : magnitude_image->columns;
1170 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1172 fourier_info.height=fourier_info.width;
1173 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
1174 fourier_info.channel=channel;
1175 fourier_info.modulus=modulus;
1176 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1177 fourier_info.center*sizeof(double));
1178 if (magnitude == (double *) NULL)
1180 (void) ThrowMagickException(exception,GetMagickModule(),
1181 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1182 magnitude_image->filename);
1183 return(MagickFalse);
1185 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1186 fourier_info.center*sizeof(double));
1187 if (phase == (double *) NULL)
1189 (void) ThrowMagickException(exception,GetMagickModule(),
1190 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1191 magnitude_image->filename);
1192 magnitude=(double *) RelinquishMagickMemory(magnitude);
1193 return(MagickFalse);
1195 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1196 fourier_info.center*sizeof(*fourier));
1197 if (fourier == (fftw_complex *) NULL)
1199 (void) ThrowMagickException(exception,GetMagickModule(),
1200 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1201 magnitude_image->filename);
1202 phase=(double *) RelinquishMagickMemory(phase);
1203 magnitude=(double *) RelinquishMagickMemory(magnitude);
1204 return(MagickFalse);
1206 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1208 if (status != MagickFalse)
1209 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1211 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1212 phase=(double *) RelinquishMagickMemory(phase);
1213 magnitude=(double *) RelinquishMagickMemory(magnitude);
1218 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1219 const Image *phase_image,const MagickBooleanType modulus,
1220 ExceptionInfo *exception)
1225 assert(magnitude_image != (Image *) NULL);
1226 assert(magnitude_image->signature == MagickSignature);
1227 if (magnitude_image->debug != MagickFalse)
1228 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1229 magnitude_image->filename);
1230 if (phase_image == (Image *) NULL)
1232 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1233 "ImageSequenceRequired","`%s'",magnitude_image->filename);
1234 return((Image *) NULL);
1236 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1237 fourier_image=(Image *) NULL;
1239 (void) ThrowMagickException(exception,GetMagickModule(),
1240 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1241 magnitude_image->filename);
1244 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1245 magnitude_image->rows,MagickFalse,exception);
1246 if (fourier_image != (Image *) NULL)
1253 is_gray=IsGrayImage(magnitude_image,exception);
1254 if (is_gray != MagickFalse)
1255 is_gray=IsGrayImage(phase_image,exception);
1256 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1257 #pragma omp parallel sections
1260 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1267 if (is_gray != MagickFalse)
1268 thread_status=InverseFourierTransformChannel(magnitude_image,
1269 phase_image,GrayChannels,modulus,fourier_image,exception);
1271 thread_status=InverseFourierTransformChannel(magnitude_image,
1272 phase_image,RedChannel,modulus,fourier_image,exception);
1273 if (thread_status == MagickFalse)
1274 status=thread_status;
1276 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1283 thread_status=MagickTrue;
1284 if (is_gray == MagickFalse)
1285 thread_status=InverseFourierTransformChannel(magnitude_image,
1286 phase_image,GreenChannel,modulus,fourier_image,exception);
1287 if (thread_status == MagickFalse)
1288 status=thread_status;
1290 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1297 thread_status=MagickTrue;
1298 if (is_gray == MagickFalse)
1299 thread_status=InverseFourierTransformChannel(magnitude_image,
1300 phase_image,BlueChannel,modulus,fourier_image,exception);
1301 if (thread_status == MagickFalse)
1302 status=thread_status;
1304 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1311 thread_status=MagickTrue;
1312 if (magnitude_image->matte != MagickFalse)
1313 thread_status=InverseFourierTransformChannel(magnitude_image,
1314 phase_image,OpacityChannel,modulus,fourier_image,exception);
1315 if (thread_status == MagickFalse)
1316 status=thread_status;
1318 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1325 thread_status=MagickTrue;
1326 if (magnitude_image->colorspace == CMYKColorspace)
1327 thread_status=InverseFourierTransformChannel(magnitude_image,
1328 phase_image,IndexChannel,modulus,fourier_image,exception);
1329 if (thread_status == MagickFalse)
1330 status=thread_status;
1333 if (status == MagickFalse)
1334 fourier_image=DestroyImage(fourier_image);
1339 return(fourier_image);