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-2009 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/thread-private.h"
56 #if defined(MAGICKCORE_FFTW_DELEGATE)
64 typedef struct _FourierInfo
81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 % 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 %
89 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 % ForwardFourierTransformImage() implements the discrete Fourier transform
92 % (DFT) of the image either as a magnitude / phase or real / imaginary image
95 % The format of the ForwadFourierTransformImage method is:
97 % Image *ForwardFourierTransformImage(const Image *image,
98 % const MagickBooleanType modulus,ExceptionInfo *exception)
100 % A description of each parameter follows:
102 % o image: the image.
104 % o modulus: if true, return as transform as a magnitude / phase pair
105 % otherwise a real / imaginary image pair.
107 % o exception: return any errors or warnings in this structure.
111 #if defined(MAGICKCORE_FFTW_DELEGATE)
113 static MagickBooleanType RollFourier(const unsigned long width,
114 const unsigned long height,const long x_offset,const long y_offset,
130 Move the zero frequency (DC) from (0,0) to (width/2,height/2).
132 roll=(double *) AcquireQuantumMemory((size_t) width,height*sizeof(*roll));
133 if (roll == (double *) NULL)
136 for (y=0L; y < (long) height; y++)
139 v=((y+y_offset) < 0L) ? y+y_offset+(long) height : y+y_offset;
141 v=((y+y_offset) > ((long) height-1L)) ? y+y_offset-(long) height :
143 for (x=0L; x < (long) width; x++)
146 u=((x+x_offset) < 0L) ? x+x_offset+(long) width : x+x_offset;
148 u=((x+x_offset) > ((long) width-1L)) ? x+x_offset-(long) width :
150 roll[v*width+u]=fourier[i++];
153 (void) CopyMagickMemory(fourier,roll,width*height*sizeof(*roll));
154 roll=(double *) RelinquishMagickMemory(roll);
158 static MagickBooleanType ForwardQuadrantSwap(const unsigned long width,
159 const unsigned long height,double *source,double *destination)
174 center=(long) floor((double) width/2L)+1L;
175 status=RollFourier((unsigned long) center,height,0L,(long) height/2L,source);
176 if (status == MagickFalse)
178 for (y=0L; y < (long) height; y++)
179 for (x=0L; x < (long) (width/2L-1L); x++)
180 destination[width*y+x+width/2L]=source[center*y+x];
181 for (y=1; y < (long) height; y++)
182 for (x=0L; x < (long) (width/2L-1L); x++)
183 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
184 for (x=0L; x < (long) (width/2L); x++)
185 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
189 static void CorrectPhaseLHS(const unsigned long width,
190 const unsigned long height,double *fourier)
198 for (y=0L; y < (long) height; y++)
199 for (x=0L; x < (long) (width/2L); x++)
200 fourier[y*width+x]*=(-1.0);
203 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
204 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
234 magnitude_image=GetFirstImageInList(image);
235 phase_image=GetNextImageInList(image);
236 if (phase_image == (Image *) NULL)
238 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
239 "ImageSequenceRequired","`%s'",image->filename);
243 Create "Fourier Transform" image from constituent arrays.
245 magnitude_source=(double *) AcquireQuantumMemory((size_t)
246 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
247 if (magnitude_source == (double *) NULL)
249 (void) ResetMagickMemory(magnitude_source,0,fourier_info->width*
250 fourier_info->height*sizeof(*magnitude_source));
251 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
252 fourier_info->width*sizeof(*phase_source));
253 if (magnitude_source == (double *) NULL)
255 (void) ThrowMagickException(exception,GetMagickModule(),
256 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
257 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
260 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
261 magnitude,magnitude_source);
262 if (status != MagickFalse)
263 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
265 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
266 if (fourier_info->modulus != MagickFalse)
269 for (y=0L; y < (long) fourier_info->height; y++)
270 for (x=0L; x < (long) fourier_info->width; x++)
272 phase_source[i]/=(2.0*MagickPI);
273 phase_source[i]+=0.5;
277 magnitude_view=AcquireCacheView(magnitude_image);
278 phase_view=AcquireCacheView(phase_image);
280 for (y=0L; y < (long) fourier_info->height; y++)
282 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
284 if (q == (PixelPacket *) NULL)
286 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
287 for (x=0L; x < (long) fourier_info->width; x++)
289 switch (fourier_info->channel)
294 q->red=RoundToQuantum(QuantumRange*magnitude_source[i]);
299 q->green=RoundToQuantum(QuantumRange*magnitude_source[i]);
304 q->blue=RoundToQuantum(QuantumRange*magnitude_source[i]);
309 q->opacity=RoundToQuantum(QuantumRange*magnitude_source[i]);
314 indexes[x]=RoundToQuantum(QuantumRange*magnitude_source[i]);
319 q->red=RoundToQuantum(QuantumRange*magnitude_source[i]);
328 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
329 if (status == MagickFalse)
333 for (y=0L; y < (long) fourier_info->height; y++)
335 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
337 if (q == (PixelPacket *) NULL)
339 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
340 for (x=0L; x < (long) fourier_info->width; x++)
342 switch (fourier_info->channel)
347 q->red=RoundToQuantum(QuantumRange*phase_source[i]);
352 q->green=RoundToQuantum(QuantumRange*phase_source[i]);
357 q->blue=RoundToQuantum(QuantumRange*phase_source[i]);
362 q->opacity=RoundToQuantum(QuantumRange*phase_source[i]);
367 indexes[x]=RoundToQuantum(QuantumRange*phase_source[i]);
372 q->red=RoundToQuantum(QuantumRange*phase_source[i]);
381 status=SyncCacheViewAuthenticPixels(phase_view,exception);
382 if (status == MagickFalse)
385 phase_view=DestroyCacheView(phase_view);
386 magnitude_view=DestroyCacheView(magnitude_view);
387 phase_source=(double *) RelinquishMagickMemory(phase_source);
388 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
392 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
393 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
411 register const IndexPacket
414 register const PixelPacket
422 Generate the forward Fourier transform.
424 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
425 fourier_info->width*sizeof(*source));
426 if (source == (double *) NULL)
428 (void) ThrowMagickException(exception,GetMagickModule(),
429 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
432 ResetMagickMemory(source,0,fourier_info->width*fourier_info->height*
435 image_view=AcquireCacheView(image);
436 for (y=0L; y < (long) fourier_info->height; y++)
438 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
440 if (p == (const PixelPacket *) NULL)
442 indexes=GetCacheViewVirtualIndexQueue(image_view);
443 for (x=0L; x < (long) fourier_info->width; x++)
445 switch (fourier_info->channel)
450 source[i]=QuantumScale*p->red;
455 source[i]=QuantumScale*p->green;
460 source[i]=QuantumScale*p->blue;
465 source[i]=QuantumScale*p->opacity;
470 source[i]=QuantumScale*indexes[x];
475 source[i]=QuantumScale*p->red;
483 image_view=DestroyCacheView(image_view);
484 fourier=(fftw_complex *) AcquireAlignedMemory((size_t) fourier_info->height,
485 fourier_info->center*sizeof(*fourier));
486 if (fourier == (fftw_complex *) NULL)
488 (void) ThrowMagickException(exception,GetMagickModule(),
489 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
490 source=(double *) RelinquishMagickMemory(source);
493 #if defined(MAGICKCORE_OPENMP_SUPPORT)
494 #pragma omp critical (MagickCore_ForwardFourierTransform)
496 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
497 source,fourier,FFTW_ESTIMATE);
498 fftw_execute(fftw_r2c_plan);
499 fftw_destroy_plan(fftw_r2c_plan);
500 source=(double *) RelinquishMagickMemory(source);
502 Normalize Fourier transform.
504 n=(double) fourier_info->width*(double) fourier_info->width;
506 for (y=0L; y < (long) fourier_info->height; y++)
507 for (x=0L; x < (long) fourier_info->center; x++)
510 Generate magnitude and phase (or real and imaginary).
513 if (fourier_info->modulus != MagickFalse)
514 for (y=0L; y < (long) fourier_info->height; y++)
515 for (x=0L; x < (long) fourier_info->center; x++)
517 magnitude[i]=cabs(fourier[i]);
518 phase[i]=carg(fourier[i]);
522 for (y=0L; y < (long) fourier_info->height; y++)
523 for (x=0L; x < (long) fourier_info->center; x++)
525 magnitude[i]=creal(fourier[i]);
526 phase[i]=cimag(fourier[i]);
529 fourier=(fftw_complex *) RelinquishAlignedMemory(fourier);
533 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
534 const ChannelType channel,const MagickBooleanType modulus,
535 Image *fourier_image,ExceptionInfo *exception)
553 fourier_info.width=image->columns;
554 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
555 ((image->rows % 2) != 0))
557 extent=image->columns < image->rows ? image->rows : image->columns;
558 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
560 fourier_info.height=fourier_info.width;
561 fourier_info.center=(long) floor((double) fourier_info.width/2.0)+1L;
562 fourier_info.channel=channel;
563 fourier_info.modulus=modulus;
564 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
565 fourier_info.center*sizeof(*magnitude));
566 if (magnitude == (double *) NULL)
568 (void) ThrowMagickException(exception,GetMagickModule(),
569 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
572 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
573 fourier_info.center*sizeof(*phase));
574 if (phase == (double *) NULL)
576 (void) ThrowMagickException(exception,GetMagickModule(),
577 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
578 magnitude=(double *) RelinquishMagickMemory(magnitude);
581 fourier=(fftw_complex *) AcquireAlignedMemory((size_t) fourier_info.height,
582 fourier_info.center*sizeof(*fourier));
583 if (fourier == (fftw_complex *) NULL)
585 (void) ThrowMagickException(exception,GetMagickModule(),
586 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
587 phase=(double *) RelinquishMagickMemory(phase);
588 magnitude=(double *) RelinquishMagickMemory(magnitude);
591 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
592 if (status != MagickFalse)
593 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
595 fourier=(fftw_complex *) RelinquishAlignedMemory(fourier);
596 phase=(double *) RelinquishMagickMemory(phase);
597 magnitude=(double *) RelinquishMagickMemory(magnitude);
602 MagickExport Image *ForwardFourierTransformImage(const Image *image,
603 const MagickBooleanType modulus,ExceptionInfo *exception)
608 fourier_image=NewImageList();
609 #if !defined(MAGICKCORE_FFTW_DELEGATE)
611 (void) ThrowMagickException(exception,GetMagickModule(),
612 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
623 width=image->columns;
624 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
625 ((image->rows % 2) != 0))
627 extent=image->columns < image->rows ? image->rows : image->columns;
628 width=(extent & 0x01) == 1 ? extent+1UL : extent;
630 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
631 if (magnitude_image != (Image *) NULL)
636 magnitude_image->storage_class=DirectClass;
637 magnitude_image->depth=32UL;
638 phase_image=CloneImage(image,width,width,MagickFalse,exception);
639 if (phase_image == (Image *) NULL)
640 magnitude_image=DestroyImage(magnitude_image);
650 phase_image->storage_class=DirectClass;
651 phase_image->depth=32UL;
652 AppendImageToList(&fourier_image,magnitude_image);
653 AppendImageToList(&fourier_image,phase_image);
655 is_gray=IsGrayImage(image,exception);
656 #if defined(MAGICKCORE_OPENMP_SUPPORT)
657 #pragma omp parallel for schedule(dynamic,4) shared(status)
659 for (i=0L; i < 5L; i++)
664 thread_status=MagickTrue;
669 if (is_gray != MagickFalse)
671 thread_status=ForwardFourierTransformChannel(image,
672 GrayChannels,modulus,fourier_image,exception);
675 thread_status=ForwardFourierTransformChannel(image,RedChannel,
676 modulus,fourier_image,exception);
681 if (is_gray == MagickFalse)
682 thread_status=ForwardFourierTransformChannel(image,
683 GreenChannel,modulus,fourier_image,exception);
688 if (is_gray == MagickFalse)
689 thread_status=ForwardFourierTransformChannel(image,
690 BlueChannel,modulus,fourier_image,exception);
695 if (image->matte != MagickFalse)
696 thread_status=ForwardFourierTransformChannel(image,
697 OpacityChannel,modulus,fourier_image,exception);
702 if (image->colorspace == CMYKColorspace)
703 thread_status=ForwardFourierTransformChannel(image,
704 IndexChannel,modulus,fourier_image,exception);
708 if (thread_status == MagickFalse)
709 status=thread_status;
711 if (status == MagickFalse)
712 fourier_image=DestroyImageList(fourier_image);
718 return(fourier_image);
722 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
726 % 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 %
730 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
732 % InverseFourierTransformImage() implements the inverse discrete Fourier
733 % transform (DFT) of the image either as a magnitude / phase or real /
734 % imaginary image pair.
736 % The format of the InverseFourierTransformImage method is:
738 % Image *InverseFourierTransformImage(const Image *images,
739 % const MagickBooleanType modulus,ExceptionInfo *exception)
741 % A description of each parameter follows:
743 % o images: the image sequence.
745 % o modulus: if true, return transform as a magnitude / phase pair
746 % otherwise a real / imaginary image pair.
748 % o exception: return any errors or warnings in this structure.
752 #if defined(MAGICKCORE_FFTW_DELEGATE)
753 static MagickBooleanType InverseQuadrantSwap(const unsigned long width,
754 const unsigned long height,const double *source,double *destination)
766 center=(long) floor((double) width/2.0)+1L;
767 for (y=1L; y < (long) height; y++)
768 for (x=0L; x < (long) (width/2L+1L); x++)
769 destination[center*(height-y)-x+width/2L]=source[y*width+x];
770 for (y=0L; y < (long) height; y++)
771 destination[center*y]=source[y*width+width/2L];
772 for (x=0L; x < center; x++)
773 destination[x]=source[center-x-1L];
774 return(RollFourier(center,height,0L,(long) height/-2L,destination));
777 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
778 const Image *images,fftw_complex *fourier,ExceptionInfo *exception)
800 register const IndexPacket
803 register const PixelPacket
811 Inverse fourier - read image and break down into a double array.
813 assert(images != (Image *) NULL);
814 assert(images->signature == MagickSignature);
815 if (images->debug != MagickFalse)
816 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
817 magnitude_image=GetFirstImageInList(images),
818 phase_image=GetNextImageInList(images);
819 if (phase_image == (Image *) NULL)
821 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
822 "ImageSequenceRequired","`%s'",images->filename);
825 magnitude_source=(double *) AcquireQuantumMemory((size_t)
826 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
827 if (magnitude_source == (double *) NULL)
829 (void) ThrowMagickException(exception,GetMagickModule(),
830 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
833 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
834 fourier_info->height*sizeof(*phase_source));
835 if (phase_source == (double *) NULL)
837 (void) ThrowMagickException(exception,GetMagickModule(),
838 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
839 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
843 magnitude_view=AcquireCacheView(magnitude_image);
844 for (y=0L; y < (long) fourier_info->height; y++)
846 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
848 if (p == (const PixelPacket *) NULL)
850 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
851 for (x=0L; x < (long) fourier_info->width; x++)
853 switch (fourier_info->channel)
858 magnitude_source[i]=QuantumScale*p->red;
863 magnitude_source[i]=QuantumScale*p->green;
868 magnitude_source[i]=QuantumScale*p->blue;
873 magnitude_source[i]=QuantumScale*p->opacity;
878 magnitude_source[i]=QuantumScale*indexes[x];
883 magnitude_source[i]=QuantumScale*p->red;
892 phase_view=AcquireCacheView(phase_image);
893 for (y=0L; y < (long) fourier_info->height; y++)
895 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
897 if (p == (const PixelPacket *) NULL)
899 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
900 for (x=0L; x < (long) fourier_info->width; x++)
902 switch (fourier_info->channel)
907 phase_source[i]=QuantumScale*p->red;
912 phase_source[i]=QuantumScale*p->green;
917 phase_source[i]=QuantumScale*p->blue;
922 phase_source[i]=QuantumScale*p->opacity;
927 phase_source[i]=QuantumScale*indexes[x];
932 phase_source[i]=QuantumScale*p->red;
940 if (fourier_info->modulus != MagickFalse)
943 for (y=0L; y < (long) fourier_info->height; y++)
944 for (x=0L; x < (long) 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'",images->filename);
959 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
960 phase_source=(double *) RelinquishMagickMemory(phase_source);
963 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
964 magnitude_source,magnitude);
965 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
966 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
967 fourier_info->height*sizeof(*phase));
968 if (phase == (double *) NULL)
970 (void) ThrowMagickException(exception,GetMagickModule(),
971 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
972 phase_source=(double *) RelinquishMagickMemory(phase_source);
975 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
976 if (status != MagickFalse)
977 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
979 phase_source=(double *) RelinquishMagickMemory(phase_source);
984 if (fourier_info->modulus != MagickFalse)
985 for (y=0L; y < (long) fourier_info->height; y++)
986 for (x=0L; x < (long) fourier_info->center; x++)
988 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
992 for (y=0L; y < (long) fourier_info->height; y++)
993 for (x=0L; x < (long) fourier_info->center; x++)
995 fourier[i]=magnitude[i]+I*phase[i];
998 phase=(double *) RelinquishMagickMemory(phase);
999 magnitude=(double *) RelinquishMagickMemory(magnitude);
1003 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1004 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1018 register IndexPacket
1025 register PixelPacket
1028 source=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
1029 fourier_info->height*sizeof(double));
1030 if (source == (double *) NULL)
1032 (void) ThrowMagickException(exception,GetMagickModule(),
1033 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1034 return(MagickFalse);
1036 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1037 #pragma omp critical (MagickCore_InverseFourierTransform)
1039 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1040 fourier,source,FFTW_ESTIMATE);
1041 fftw_execute(fftw_c2r_plan);
1042 fftw_destroy_plan(fftw_c2r_plan);
1044 image_view=AcquireCacheView(image);
1045 for (y=0L; y < (long) fourier_info->height; y++)
1047 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width,1UL,
1049 if (q == (PixelPacket *) NULL)
1051 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1052 for (x=0L; x < (long) fourier_info->width; x++)
1054 switch (fourier_info->channel)
1059 q->red=RoundToQuantum(QuantumRange*source[i]);
1064 q->green=RoundToQuantum(QuantumRange*source[i]);
1069 q->blue=RoundToQuantum(QuantumRange*source[i]);
1072 case OpacityChannel:
1074 q->opacity=RoundToQuantum(QuantumRange*source[i]);
1079 indexes[x]=RoundToQuantum(QuantumRange*source[i]);
1084 q->red=RoundToQuantum(QuantumRange*source[i]);
1093 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1096 image_view=DestroyCacheView(image_view);
1097 source=(double *) RelinquishMagickMemory(source);
1101 static MagickBooleanType InverseFourierTransformChannel(const Image *images,
1102 const ChannelType channel,const MagickBooleanType modulus,
1103 Image *fourier_image,ExceptionInfo *exception)
1121 fourier_info.width=images->columns;
1122 if ((images->columns != images->rows) || ((images->columns % 2) != 0) ||
1123 ((images->rows % 2) != 0))
1125 extent=images->columns < images->rows ? images->rows : images->columns;
1126 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1128 fourier_info.height=fourier_info.width;
1129 fourier_info.center=(long) floor((double) fourier_info.width/2.0)+1L;
1130 fourier_info.channel=channel;
1131 fourier_info.modulus=modulus;
1132 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1133 fourier_info.center*sizeof(double));
1134 if (magnitude == (double *) NULL)
1136 (void) ThrowMagickException(exception,GetMagickModule(),
1137 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
1138 return(MagickFalse);
1140 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1141 fourier_info.center*sizeof(double));
1142 if (phase == (double *) NULL)
1144 (void) ThrowMagickException(exception,GetMagickModule(),
1145 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
1146 magnitude=(double *) RelinquishMagickMemory(magnitude);
1147 return(MagickFalse);
1149 fourier=(fftw_complex *) AcquireAlignedMemory((size_t) fourier_info.height,
1150 fourier_info.center*sizeof(*fourier));
1151 if (fourier == (fftw_complex *) NULL)
1153 (void) ThrowMagickException(exception,GetMagickModule(),
1154 ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
1155 phase=(double *) RelinquishMagickMemory(phase);
1156 magnitude=(double *) RelinquishMagickMemory(magnitude);
1157 return(MagickFalse);
1159 status=InverseFourier(&fourier_info,images,fourier,exception);
1160 if (status != MagickFalse)
1161 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1163 fourier=(fftw_complex *) RelinquishAlignedMemory(fourier);
1164 phase=(double *) RelinquishMagickMemory(phase);
1165 magnitude=(double *) RelinquishMagickMemory(magnitude);
1170 MagickExport Image *InverseFourierTransformImage(const Image *images,
1171 const MagickBooleanType modulus,ExceptionInfo *exception)
1176 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1177 fourier_image=(Image *) NULL;
1179 (void) ThrowMagickException(exception,GetMagickModule(),
1180 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1184 fourier_image=CloneImage(images,images->columns,images->rows,MagickFalse,
1186 if (fourier_image != (Image *) NULL)
1196 is_gray=IsGrayImage(images,exception);
1197 if ((is_gray != MagickFalse) && (images->next != (Image *) NULL))
1198 is_gray=IsGrayImage(images->next,exception);
1199 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1200 #pragma omp parallel for schedule(dynamic,4) shared(status)
1202 for (i=0L; i < 5L; i++)
1207 thread_status=MagickTrue;
1212 if (is_gray != MagickFalse)
1214 thread_status=InverseFourierTransformChannel(images,
1215 GrayChannels,modulus,fourier_image,exception);
1218 thread_status=InverseFourierTransformChannel(images,RedChannel,
1219 modulus,fourier_image,exception);
1224 if (is_gray == MagickFalse)
1225 thread_status=InverseFourierTransformChannel(images,
1226 GreenChannel,modulus,fourier_image,exception);
1231 if (is_gray == MagickFalse)
1232 thread_status=InverseFourierTransformChannel(images,BlueChannel,
1233 modulus,fourier_image,exception);
1238 if (images->matte != MagickFalse)
1239 thread_status=InverseFourierTransformChannel(images,
1240 OpacityChannel,modulus,fourier_image,exception);
1245 if (images->colorspace == CMYKColorspace)
1246 thread_status=InverseFourierTransformChannel(images,
1247 IndexChannel,modulus,fourier_image,exception);
1251 if (thread_status == MagickFalse)
1252 status=thread_status;
1254 if (status == MagickFalse)
1255 fourier_image=DestroyImage(fourier_image);
1260 return(fourier_image);