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/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 *magnitude_image,
739 % const Image *phase_image,const MagickBooleanType modulus,
740 % ExceptionInfo *exception)
742 % A description of each parameter follows:
744 % o magnitude_image: the magnitude or real image.
746 % o phase_image: the phase or imaginary image.
748 % o modulus: if true, return transform as a magnitude / phase pair
749 % otherwise a real / imaginary image pair.
751 % o exception: return any errors or warnings in this structure.
755 #if defined(MAGICKCORE_FFTW_DELEGATE)
756 static MagickBooleanType InverseQuadrantSwap(const unsigned long width,
757 const unsigned long height,const double *source,double *destination)
769 center=(long) floor((double) width/2.0)+1L;
770 for (y=1L; y < (long) height; y++)
771 for (x=0L; x < (long) (width/2L+1L); x++)
772 destination[center*(height-y)-x+width/2L]=source[y*width+x];
773 for (y=0L; y < (long) height; y++)
774 destination[center*y]=source[y*width+width/2L];
775 for (x=0L; x < center; x++)
776 destination[x]=source[center-x-1L];
777 return(RollFourier(center,height,0L,(long) height/-2L,destination));
780 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
781 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
782 ExceptionInfo *exception)
800 register const IndexPacket
803 register const PixelPacket
811 Inverse fourier - read image and break down into a double array.
813 magnitude_source=(double *) AcquireQuantumMemory((size_t)
814 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
815 if (magnitude_source == (double *) NULL)
817 (void) ThrowMagickException(exception,GetMagickModule(),
818 ResourceLimitError,"MemoryAllocationFailed","`%s'",
819 magnitude_image->filename);
822 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
823 fourier_info->height*sizeof(*phase_source));
824 if (phase_source == (double *) NULL)
826 (void) ThrowMagickException(exception,GetMagickModule(),
827 ResourceLimitError,"MemoryAllocationFailed","`%s'",
828 magnitude_image->filename);
829 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
833 magnitude_view=AcquireCacheView(magnitude_image);
834 for (y=0L; y < (long) fourier_info->height; y++)
836 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
838 if (p == (const PixelPacket *) NULL)
840 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
841 for (x=0L; x < (long) fourier_info->width; x++)
843 switch (fourier_info->channel)
848 magnitude_source[i]=QuantumScale*p->red;
853 magnitude_source[i]=QuantumScale*p->green;
858 magnitude_source[i]=QuantumScale*p->blue;
863 magnitude_source[i]=QuantumScale*p->opacity;
868 magnitude_source[i]=QuantumScale*indexes[x];
873 magnitude_source[i]=QuantumScale*p->red;
882 phase_view=AcquireCacheView(phase_image);
883 for (y=0L; y < (long) fourier_info->height; y++)
885 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
887 if (p == (const PixelPacket *) NULL)
889 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
890 for (x=0L; x < (long) fourier_info->width; x++)
892 switch (fourier_info->channel)
897 phase_source[i]=QuantumScale*p->red;
902 phase_source[i]=QuantumScale*p->green;
907 phase_source[i]=QuantumScale*p->blue;
912 phase_source[i]=QuantumScale*p->opacity;
917 phase_source[i]=QuantumScale*indexes[x];
922 phase_source[i]=QuantumScale*p->red;
930 if (fourier_info->modulus != MagickFalse)
933 for (y=0L; y < (long) fourier_info->height; y++)
934 for (x=0L; x < (long) fourier_info->width; x++)
936 phase_source[i]-=0.5;
937 phase_source[i]*=(2.0*MagickPI);
941 magnitude_view=DestroyCacheView(magnitude_view);
942 phase_view=DestroyCacheView(phase_view);
943 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
944 fourier_info->center*sizeof(*magnitude));
945 if (magnitude == (double *) NULL)
947 (void) ThrowMagickException(exception,GetMagickModule(),
948 ResourceLimitError,"MemoryAllocationFailed","`%s'",
949 magnitude_image->filename);
950 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
951 phase_source=(double *) RelinquishMagickMemory(phase_source);
954 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
955 magnitude_source,magnitude);
956 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
957 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
958 fourier_info->height*sizeof(*phase));
959 if (phase == (double *) NULL)
961 (void) ThrowMagickException(exception,GetMagickModule(),
962 ResourceLimitError,"MemoryAllocationFailed","`%s'",
963 magnitude_image->filename);
964 phase_source=(double *) RelinquishMagickMemory(phase_source);
967 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
968 if (status != MagickFalse)
969 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
971 phase_source=(double *) RelinquishMagickMemory(phase_source);
976 if (fourier_info->modulus != MagickFalse)
977 for (y=0L; y < (long) fourier_info->height; y++)
978 for (x=0L; x < (long) fourier_info->center; x++)
980 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
984 for (y=0L; y < (long) fourier_info->height; y++)
985 for (x=0L; x < (long) fourier_info->center; x++)
987 fourier[i]=magnitude[i]+I*phase[i];
990 phase=(double *) RelinquishMagickMemory(phase);
991 magnitude=(double *) RelinquishMagickMemory(magnitude);
995 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
996 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1010 register IndexPacket
1017 register PixelPacket
1020 source=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
1021 fourier_info->height*sizeof(double));
1022 if (source == (double *) NULL)
1024 (void) ThrowMagickException(exception,GetMagickModule(),
1025 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1026 return(MagickFalse);
1028 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1029 #pragma omp critical (MagickCore_InverseFourierTransform)
1031 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1032 fourier,source,FFTW_ESTIMATE);
1033 fftw_execute(fftw_c2r_plan);
1034 fftw_destroy_plan(fftw_c2r_plan);
1036 image_view=AcquireCacheView(image);
1037 for (y=0L; y < (long) fourier_info->height; y++)
1039 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width,1UL,
1041 if (q == (PixelPacket *) NULL)
1043 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1044 for (x=0L; x < (long) fourier_info->width; x++)
1046 switch (fourier_info->channel)
1051 q->red=RoundToQuantum(QuantumRange*source[i]);
1056 q->green=RoundToQuantum(QuantumRange*source[i]);
1061 q->blue=RoundToQuantum(QuantumRange*source[i]);
1064 case OpacityChannel:
1066 q->opacity=RoundToQuantum(QuantumRange*source[i]);
1071 indexes[x]=RoundToQuantum(QuantumRange*source[i]);
1076 q->red=RoundToQuantum(QuantumRange*source[i]);
1085 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1088 image_view=DestroyCacheView(image_view);
1089 source=(double *) RelinquishMagickMemory(source);
1093 static MagickBooleanType InverseFourierTransformChannel(
1094 const Image *magnitude_image,const Image *phase_image,
1095 const ChannelType channel,const MagickBooleanType modulus,
1096 Image *fourier_image,ExceptionInfo *exception)
1114 fourier_info.width=magnitude_image->columns;
1115 if ((magnitude_image->columns != magnitude_image->rows) ||
1116 ((magnitude_image->columns % 2) != 0) ||
1117 ((magnitude_image->rows % 2) != 0))
1119 extent=magnitude_image->columns < magnitude_image->rows ?
1120 magnitude_image->rows : magnitude_image->columns;
1121 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1123 fourier_info.height=fourier_info.width;
1124 fourier_info.center=(long) floor((double) fourier_info.width/2.0)+1L;
1125 fourier_info.channel=channel;
1126 fourier_info.modulus=modulus;
1127 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1128 fourier_info.center*sizeof(double));
1129 if (magnitude == (double *) NULL)
1131 (void) ThrowMagickException(exception,GetMagickModule(),
1132 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1133 magnitude_image->filename);
1134 return(MagickFalse);
1136 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1137 fourier_info.center*sizeof(double));
1138 if (phase == (double *) NULL)
1140 (void) ThrowMagickException(exception,GetMagickModule(),
1141 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1142 magnitude_image->filename);
1143 magnitude=(double *) RelinquishMagickMemory(magnitude);
1144 return(MagickFalse);
1146 fourier=(fftw_complex *) AcquireAlignedMemory((size_t) fourier_info.height,
1147 fourier_info.center*sizeof(*fourier));
1148 if (fourier == (fftw_complex *) NULL)
1150 (void) ThrowMagickException(exception,GetMagickModule(),
1151 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1152 magnitude_image->filename);
1153 phase=(double *) RelinquishMagickMemory(phase);
1154 magnitude=(double *) RelinquishMagickMemory(magnitude);
1155 return(MagickFalse);
1157 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1159 if (status != MagickFalse)
1160 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1162 fourier=(fftw_complex *) RelinquishAlignedMemory(fourier);
1163 phase=(double *) RelinquishMagickMemory(phase);
1164 magnitude=(double *) RelinquishMagickMemory(magnitude);
1169 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1170 const Image *phase_image,const MagickBooleanType modulus,
1171 ExceptionInfo *exception)
1176 assert(magnitude_image != (Image *) NULL);
1177 assert(magnitude_image->signature == MagickSignature);
1178 if (magnitude_image->debug != MagickFalse)
1179 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1180 magnitude_image->filename);
1181 if (phase_image == (Image *) NULL)
1183 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1184 "ImageSequenceRequired","`%s'",magnitude_image->filename);
1185 return((Image *) NULL);
1187 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1188 fourier_image=(Image *) NULL;
1190 (void) ThrowMagickException(exception,GetMagickModule(),
1191 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1192 magnitude_image->filename);
1195 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1196 magnitude_image->rows,MagickFalse,exception);
1197 if (fourier_image != (Image *) NULL)
1207 is_gray=IsGrayImage(magnitude_image,exception);
1208 if (is_gray != MagickFalse)
1209 is_gray=IsGrayImage(phase_image,exception);
1210 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1211 #pragma omp parallel for schedule(dynamic,4) shared(status)
1213 for (i=0L; i < 5L; i++)
1218 thread_status=MagickTrue;
1223 if (is_gray != MagickFalse)
1225 thread_status=InverseFourierTransformChannel(magnitude_image,
1226 phase_image,GrayChannels,modulus,fourier_image,exception);
1229 thread_status=InverseFourierTransformChannel(magnitude_image,
1230 phase_image,RedChannel,modulus,fourier_image,exception);
1235 if (is_gray == MagickFalse)
1236 thread_status=InverseFourierTransformChannel(magnitude_image,
1237 phase_image,GreenChannel,modulus,fourier_image,exception);
1242 if (is_gray == MagickFalse)
1243 thread_status=InverseFourierTransformChannel(magnitude_image,
1244 phase_image,BlueChannel,modulus,fourier_image,exception);
1249 if (magnitude_image->matte != MagickFalse)
1250 thread_status=InverseFourierTransformChannel(magnitude_image,
1251 phase_image,OpacityChannel,modulus,fourier_image,exception);
1256 if (magnitude_image->colorspace == CMYKColorspace)
1257 thread_status=InverseFourierTransformChannel(magnitude_image,
1258 phase_image,IndexChannel,modulus,fourier_image,exception);
1262 if (thread_status == MagickFalse)
1263 status=thread_status;
1265 if (status == MagickFalse)
1266 fourier_image=DestroyImage(fourier_image);
1271 return(fourier_image);