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-2012 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/resource_.h"
59 #include "MagickCore/thread-private.h"
60 #if defined(MAGICKCORE_FFTW_DELEGATE)
61 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
65 #if !defined(MAGICKCORE_HAVE_CABS)
66 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
68 #if !defined(MAGICKCORE_HAVE_CARG)
69 #define carg(z) (atan2(cimag(z),creal(z)))
71 #if !defined(MAGICKCORE_HAVE_CIMAG)
72 #define cimag(z) (z[1])
74 #if !defined(MAGICKCORE_HAVE_CREAL)
75 #define creal(z) (z[0])
82 typedef struct _FourierInfo
99 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
103 % 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 %
107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
109 % ForwardFourierTransformImage() implements the discrete Fourier transform
110 % (DFT) of the image either as a magnitude / phase or real / imaginary image
113 % The format of the ForwadFourierTransformImage method is:
115 % Image *ForwardFourierTransformImage(const Image *image,
116 % const MagickBooleanType modulus,ExceptionInfo *exception)
118 % A description of each parameter follows:
120 % o image: the image.
122 % o modulus: if true, return as transform as a magnitude / phase pair
123 % otherwise a real / imaginary image pair.
125 % o exception: return any errors or warnings in this structure.
129 #if defined(MAGICKCORE_FFTW_DELEGATE)
131 static MagickBooleanType RollFourier(const size_t width,const size_t height,
132 const ssize_t x_offset,const ssize_t y_offset,double *fourier)
147 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
149 roll=(double *) AcquireQuantumMemory((size_t) height,width*sizeof(*roll));
150 if (roll == (double *) NULL)
153 for (y=0L; y < (ssize_t) height; y++)
156 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
158 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
160 for (x=0L; x < (ssize_t) width; x++)
163 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
165 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
167 roll[v*width+u]=fourier[i++];
170 (void) CopyMagickMemory(fourier,roll,height*width*sizeof(*roll));
171 roll=(double *) RelinquishMagickMemory(roll);
175 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
176 const size_t height,double *source,double *destination)
191 center=(ssize_t) floor((double) width/2L)+1L;
192 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
193 if (status == MagickFalse)
195 for (y=0L; y < (ssize_t) height; y++)
196 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
197 destination[width*y+x+width/2L]=source[center*y+x];
198 for (y=1; y < (ssize_t) height; y++)
199 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
200 destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
201 for (x=0L; x < (ssize_t) (width/2L); x++)
202 destination[-x+width/2L-1L]=destination[x+width/2L+1L];
206 static void CorrectPhaseLHS(const size_t width,const size_t height,
215 for (y=0L; y < (ssize_t) height; y++)
216 for (x=0L; x < (ssize_t) (width/2L); x++)
217 fourier[y*width+x]*=(-1.0);
220 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
221 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
248 magnitude_image=GetFirstImageInList(image);
249 phase_image=GetNextImageInList(image);
250 if (phase_image == (Image *) NULL)
252 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
253 "TwoOrMoreImagesRequired","'%s'",image->filename);
257 Create "Fourier Transform" image from constituent arrays.
259 magnitude_source=(double *) AcquireQuantumMemory((size_t)
260 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
261 if (magnitude_source == (double *) NULL)
263 (void) ResetMagickMemory(magnitude_source,0,fourier_info->height*
264 fourier_info->width*sizeof(*magnitude_source));
265 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
266 fourier_info->width*sizeof(*phase_source));
267 if (phase_source == (double *) NULL)
269 (void) ThrowMagickException(exception,GetMagickModule(),
270 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
271 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
274 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
275 magnitude,magnitude_source);
276 if (status != MagickFalse)
277 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
279 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
280 if (fourier_info->modulus != MagickFalse)
283 for (y=0L; y < (ssize_t) fourier_info->height; y++)
284 for (x=0L; x < (ssize_t) fourier_info->width; x++)
286 phase_source[i]/=(2.0*MagickPI);
287 phase_source[i]+=0.5;
291 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
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)
303 case RedPixelChannel:
306 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
307 magnitude_source[i]),q);
310 case GreenPixelChannel:
312 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
313 magnitude_source[i]),q);
316 case BluePixelChannel:
318 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
319 magnitude_source[i]),q);
322 case BlackPixelChannel:
324 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
325 magnitude_source[i]),q);
328 case AlphaPixelChannel:
330 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
331 magnitude_source[i]),q);
336 q+=GetPixelChannels(magnitude_image);
338 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
339 if (status == MagickFalse)
342 magnitude_view=DestroyCacheView(magnitude_view);
344 phase_view=AcquireAuthenticCacheView(phase_image,exception);
345 for (y=0L; y < (ssize_t) fourier_info->height; y++)
347 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
349 if (q == (Quantum *) NULL)
351 for (x=0L; x < (ssize_t) fourier_info->width; x++)
353 switch (fourier_info->channel)
355 case RedPixelChannel:
358 SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
362 case GreenPixelChannel:
364 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
368 case BluePixelChannel:
370 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
374 case BlackPixelChannel:
376 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
380 case AlphaPixelChannel:
382 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
388 q+=GetPixelChannels(phase_image);
390 status=SyncCacheViewAuthenticPixels(phase_view,exception);
391 if (status == MagickFalse)
394 phase_view=DestroyCacheView(phase_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)
416 register const Quantum
427 Generate the forward Fourier transform.
429 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
430 fourier_info->width*sizeof(*source));
431 if (source == (double *) NULL)
433 (void) ThrowMagickException(exception,GetMagickModule(),
434 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
437 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
440 image_view=AcquireVirtualCacheView(image,exception);
441 for (y=0L; y < (ssize_t) fourier_info->height; y++)
443 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
445 if (p == (const Quantum *) NULL)
447 for (x=0L; x < (ssize_t) fourier_info->width; x++)
449 switch (fourier_info->channel)
451 case RedPixelChannel:
454 source[i]=QuantumScale*GetPixelRed(image,p);
457 case GreenPixelChannel:
459 source[i]=QuantumScale*GetPixelGreen(image,p);
462 case BluePixelChannel:
464 source[i]=QuantumScale*GetPixelBlue(image,p);
467 case BlackPixelChannel:
469 source[i]=QuantumScale*GetPixelBlack(image,p);
472 case AlphaPixelChannel:
474 source[i]=QuantumScale*GetPixelAlpha(image,p);
479 p+=GetPixelChannels(image);
482 image_view=DestroyCacheView(image_view);
483 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
484 fourier_info->center*sizeof(*fourier));
485 if (fourier == (fftw_complex *) NULL)
487 (void) ThrowMagickException(exception,GetMagickModule(),
488 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
489 source=(double *) RelinquishMagickMemory(source);
492 #if defined(MAGICKCORE_OPENMP_SUPPORT)
493 #pragma omp critical (MagickCore_ForwardFourierTransform)
495 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
496 source,fourier,FFTW_ESTIMATE);
497 fftw_execute(fftw_r2c_plan);
498 fftw_destroy_plan(fftw_r2c_plan);
499 source=(double *) RelinquishMagickMemory(source);
501 Normalize Fourier transform.
503 n=(double) fourier_info->width*(double) fourier_info->width;
505 for (y=0L; y < (ssize_t) fourier_info->height; y++)
506 for (x=0L; x < (ssize_t) fourier_info->center; x++)
508 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
517 Generate magnitude and phase (or real and imaginary).
520 if (fourier_info->modulus != MagickFalse)
521 for (y=0L; y < (ssize_t) fourier_info->height; y++)
522 for (x=0L; x < (ssize_t) fourier_info->center; x++)
524 magnitude[i]=cabs(fourier[i]);
525 phase[i]=carg(fourier[i]);
529 for (y=0L; y < (ssize_t) fourier_info->height; y++)
530 for (x=0L; x < (ssize_t) fourier_info->center; x++)
532 magnitude[i]=creal(fourier[i]);
533 phase[i]=cimag(fourier[i]);
536 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
540 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
541 const PixelChannel channel,const MagickBooleanType modulus,
542 Image *fourier_image,ExceptionInfo *exception)
560 fourier_info.width=image->columns;
561 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
562 ((image->rows % 2) != 0))
564 extent=image->columns < image->rows ? image->rows : image->columns;
565 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
567 fourier_info.height=fourier_info.width;
568 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
569 fourier_info.channel=channel;
570 fourier_info.modulus=modulus;
571 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
572 fourier_info.center*sizeof(*magnitude));
573 if (magnitude == (double *) NULL)
575 (void) ThrowMagickException(exception,GetMagickModule(),
576 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
579 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
580 fourier_info.center*sizeof(*phase));
581 if (phase == (double *) NULL)
583 (void) ThrowMagickException(exception,GetMagickModule(),
584 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
585 magnitude=(double *) RelinquishMagickMemory(magnitude);
588 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
589 fourier_info.center*sizeof(*fourier));
590 if (fourier == (fftw_complex *) NULL)
592 (void) ThrowMagickException(exception,GetMagickModule(),
593 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
594 phase=(double *) RelinquishMagickMemory(phase);
595 magnitude=(double *) RelinquishMagickMemory(magnitude);
598 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
599 if (status != MagickFalse)
600 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
602 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
603 phase=(double *) RelinquishMagickMemory(phase);
604 magnitude=(double *) RelinquishMagickMemory(magnitude);
609 MagickExport Image *ForwardFourierTransformImage(const Image *image,
610 const MagickBooleanType modulus,ExceptionInfo *exception)
615 fourier_image=NewImageList();
616 #if !defined(MAGICKCORE_FFTW_DELEGATE)
618 (void) ThrowMagickException(exception,GetMagickModule(),
619 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
630 width=image->columns;
631 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
632 ((image->rows % 2) != 0))
634 extent=image->columns < image->rows ? image->rows : image->columns;
635 width=(extent & 0x01) == 1 ? extent+1UL : extent;
637 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
638 if (magnitude_image != (Image *) NULL)
643 magnitude_image->storage_class=DirectClass;
644 magnitude_image->depth=32UL;
645 phase_image=CloneImage(image,width,width,MagickFalse,exception);
646 if (phase_image == (Image *) NULL)
647 magnitude_image=DestroyImage(magnitude_image);
654 phase_image->storage_class=DirectClass;
655 phase_image->depth=32UL;
656 AppendImageToList(&fourier_image,magnitude_image);
657 AppendImageToList(&fourier_image,phase_image);
659 is_gray=IsImageGray(image,exception);
660 #if defined(MAGICKCORE_OPENMP_SUPPORT)
661 #pragma omp parallel sections
664 #if defined(MAGICKCORE_OPENMP_SUPPORT)
671 if (is_gray != MagickFalse)
672 thread_status=ForwardFourierTransformChannel(image,
673 GrayPixelChannel,modulus,fourier_image,exception);
675 thread_status=ForwardFourierTransformChannel(image,
676 RedPixelChannel,modulus,fourier_image,exception);
677 if (thread_status == MagickFalse)
678 status=thread_status;
680 #if defined(MAGICKCORE_OPENMP_SUPPORT)
687 thread_status=MagickTrue;
688 if (is_gray == MagickFalse)
689 thread_status=ForwardFourierTransformChannel(image,
690 GreenPixelChannel,modulus,fourier_image,exception);
691 if (thread_status == MagickFalse)
692 status=thread_status;
694 #if defined(MAGICKCORE_OPENMP_SUPPORT)
701 thread_status=MagickTrue;
702 if (is_gray == MagickFalse)
703 thread_status=ForwardFourierTransformChannel(image,
704 BluePixelChannel,modulus,fourier_image,exception);
705 if (thread_status == MagickFalse)
706 status=thread_status;
708 #if defined(MAGICKCORE_OPENMP_SUPPORT)
715 thread_status=MagickTrue;
716 if (image->colorspace == CMYKColorspace)
717 thread_status=ForwardFourierTransformChannel(image,
718 BlackPixelChannel,modulus,fourier_image,exception);
719 if (thread_status == MagickFalse)
720 status=thread_status;
722 #if defined(MAGICKCORE_OPENMP_SUPPORT)
729 thread_status=MagickTrue;
730 if (image->matte != MagickFalse)
731 thread_status=ForwardFourierTransformChannel(image,
732 AlphaPixelChannel,modulus,fourier_image,exception);
733 if (thread_status == MagickFalse)
734 status=thread_status;
737 if (status == MagickFalse)
738 fourier_image=DestroyImageList(fourier_image);
744 return(fourier_image);
748 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
752 % 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 %
756 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
758 % InverseFourierTransformImage() implements the inverse discrete Fourier
759 % transform (DFT) of the image either as a magnitude / phase or real /
760 % imaginary image pair.
762 % The format of the InverseFourierTransformImage method is:
764 % Image *InverseFourierTransformImage(const Image *magnitude_image,
765 % const Image *phase_image,const MagickBooleanType modulus,
766 % ExceptionInfo *exception)
768 % A description of each parameter follows:
770 % o magnitude_image: the magnitude or real image.
772 % o phase_image: the phase or imaginary image.
774 % o modulus: if true, return transform as a magnitude / phase pair
775 % otherwise a real / imaginary image pair.
777 % o exception: return any errors or warnings in this structure.
781 #if defined(MAGICKCORE_FFTW_DELEGATE)
782 static MagickBooleanType InverseQuadrantSwap(const size_t width,
783 const size_t height,const double *source,double *destination)
795 center=(ssize_t) floor((double) width/2L)+1L;
796 for (y=1L; y < (ssize_t) height; y++)
797 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
798 destination[center*(height-y)-x+width/2L]=source[y*width+x];
799 for (y=0L; y < (ssize_t) height; y++)
800 destination[center*y]=source[y*width+width/2L];
801 for (x=0L; x < center; x++)
802 destination[x]=source[center-x-1L];
803 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
806 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
807 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
808 ExceptionInfo *exception)
823 register const Quantum
834 Inverse fourier - read image and break down into a double array.
836 magnitude_source=(double *) AcquireQuantumMemory((size_t)
837 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
838 if (magnitude_source == (double *) NULL)
840 (void) ThrowMagickException(exception,GetMagickModule(),
841 ResourceLimitError,"MemoryAllocationFailed","'%s'",
842 magnitude_image->filename);
845 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
846 fourier_info->width*sizeof(*phase_source));
847 if (phase_source == (double *) NULL)
849 (void) ThrowMagickException(exception,GetMagickModule(),
850 ResourceLimitError,"MemoryAllocationFailed","'%s'",
851 magnitude_image->filename);
852 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
856 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
857 for (y=0L; y < (ssize_t) fourier_info->height; y++)
859 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
861 if (p == (const Quantum *) NULL)
863 for (x=0L; x < (ssize_t) fourier_info->width; x++)
865 switch (fourier_info->channel)
867 case RedPixelChannel:
870 magnitude_source[i]=QuantumScale*GetPixelRed(magnitude_image,p);
873 case GreenPixelChannel:
875 magnitude_source[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
878 case BluePixelChannel:
880 magnitude_source[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
883 case BlackPixelChannel:
885 magnitude_source[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
888 case AlphaPixelChannel:
890 magnitude_source[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
895 p+=GetPixelChannels(magnitude_image);
899 phase_view=AcquireVirtualCacheView(phase_image,exception);
900 for (y=0L; y < (ssize_t) fourier_info->height; y++)
902 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
904 if (p == (const Quantum *) NULL)
906 for (x=0L; x < (ssize_t) fourier_info->width; x++)
908 switch (fourier_info->channel)
910 case RedPixelChannel:
913 phase_source[i]=QuantumScale*GetPixelRed(phase_image,p);
916 case GreenPixelChannel:
918 phase_source[i]=QuantumScale*GetPixelGreen(phase_image,p);
921 case BluePixelChannel:
923 phase_source[i]=QuantumScale*GetPixelBlue(phase_image,p);
926 case BlackPixelChannel:
928 phase_source[i]=QuantumScale*GetPixelBlack(phase_image,p);
931 case AlphaPixelChannel:
933 phase_source[i]=QuantumScale*GetPixelAlpha(phase_image,p);
938 p+=GetPixelChannels(phase_image);
941 if (fourier_info->modulus != MagickFalse)
944 for (y=0L; y < (ssize_t) fourier_info->height; y++)
945 for (x=0L; x < (ssize_t) fourier_info->width; x++)
947 phase_source[i]-=0.5;
948 phase_source[i]*=(2.0*MagickPI);
952 magnitude_view=DestroyCacheView(magnitude_view);
953 phase_view=DestroyCacheView(phase_view);
954 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
955 fourier_info->center*sizeof(*magnitude));
956 if (magnitude == (double *) NULL)
958 (void) ThrowMagickException(exception,GetMagickModule(),
959 ResourceLimitError,"MemoryAllocationFailed","'%s'",
960 magnitude_image->filename);
961 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
962 phase_source=(double *) RelinquishMagickMemory(phase_source);
965 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
966 magnitude_source,magnitude);
967 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
968 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
969 fourier_info->width*sizeof(*phase));
970 if (phase == (double *) NULL)
972 (void) ThrowMagickException(exception,GetMagickModule(),
973 ResourceLimitError,"MemoryAllocationFailed","'%s'",
974 magnitude_image->filename);
975 phase_source=(double *) RelinquishMagickMemory(phase_source);
978 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
979 if (status != MagickFalse)
980 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
982 phase_source=(double *) RelinquishMagickMemory(phase_source);
987 if (fourier_info->modulus != MagickFalse)
988 for (y=0L; y < (ssize_t) fourier_info->height; y++)
989 for (x=0L; x < (ssize_t) fourier_info->center; x++)
991 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
992 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
994 fourier[i][0]=magnitude[i]*cos(phase[i]);
995 fourier[i][1]=magnitude[i]*sin(phase[i]);
1000 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1001 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1003 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1004 fourier[i]=magnitude[i]+I*phase[i];
1006 fourier[i][0]=magnitude[i];
1007 fourier[i][1]=phase[i];
1011 phase=(double *) RelinquishMagickMemory(phase);
1012 magnitude=(double *) RelinquishMagickMemory(magnitude);
1016 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1017 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1038 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1039 fourier_info->width*sizeof(*source));
1040 if (source == (double *) NULL)
1042 (void) ThrowMagickException(exception,GetMagickModule(),
1043 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
1044 return(MagickFalse);
1046 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1047 #pragma omp critical (MagickCore_InverseFourierTransform)
1050 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1051 fourier,source,FFTW_ESTIMATE);
1052 fftw_execute(fftw_c2r_plan);
1053 fftw_destroy_plan(fftw_c2r_plan);
1056 image_view=AcquireAuthenticCacheView(image,exception);
1057 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1059 if (y >= (ssize_t) image->rows)
1061 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1062 image->columns ? image->columns : fourier_info->width,1UL,exception);
1063 if (q == (Quantum *) NULL)
1065 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1067 if (x < (ssize_t) image->columns)
1068 switch (fourier_info->channel)
1070 case RedPixelChannel:
1073 SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
1076 case GreenPixelChannel:
1078 SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
1081 case BluePixelChannel:
1083 SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
1086 case BlackPixelChannel:
1088 SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
1091 case AlphaPixelChannel:
1093 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
1098 q+=GetPixelChannels(image);
1100 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1103 image_view=DestroyCacheView(image_view);
1104 source=(double *) RelinquishMagickMemory(source);
1108 static MagickBooleanType InverseFourierTransformChannel(
1109 const Image *magnitude_image,const Image *phase_image,
1110 const PixelChannel channel,const MagickBooleanType modulus,
1111 Image *fourier_image,ExceptionInfo *exception)
1129 fourier_info.width=magnitude_image->columns;
1130 if ((magnitude_image->columns != magnitude_image->rows) ||
1131 ((magnitude_image->columns % 2) != 0) ||
1132 ((magnitude_image->rows % 2) != 0))
1134 extent=magnitude_image->columns < magnitude_image->rows ?
1135 magnitude_image->rows : magnitude_image->columns;
1136 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1138 fourier_info.height=fourier_info.width;
1139 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
1140 fourier_info.channel=channel;
1141 fourier_info.modulus=modulus;
1142 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1143 fourier_info.center*sizeof(*magnitude));
1144 if (magnitude == (double *) NULL)
1146 (void) ThrowMagickException(exception,GetMagickModule(),
1147 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1148 magnitude_image->filename);
1149 return(MagickFalse);
1151 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1152 fourier_info.center*sizeof(*phase));
1153 if (phase == (double *) NULL)
1155 (void) ThrowMagickException(exception,GetMagickModule(),
1156 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1157 magnitude_image->filename);
1158 magnitude=(double *) RelinquishMagickMemory(magnitude);
1159 return(MagickFalse);
1161 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1162 fourier_info.center*sizeof(*fourier));
1163 if (fourier == (fftw_complex *) NULL)
1165 (void) ThrowMagickException(exception,GetMagickModule(),
1166 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1167 magnitude_image->filename);
1168 phase=(double *) RelinquishMagickMemory(phase);
1169 magnitude=(double *) RelinquishMagickMemory(magnitude);
1170 return(MagickFalse);
1172 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1174 if (status != MagickFalse)
1175 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1177 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1178 phase=(double *) RelinquishMagickMemory(phase);
1179 magnitude=(double *) RelinquishMagickMemory(magnitude);
1184 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1185 const Image *phase_image,const MagickBooleanType modulus,
1186 ExceptionInfo *exception)
1191 assert(magnitude_image != (Image *) NULL);
1192 assert(magnitude_image->signature == MagickSignature);
1193 if (magnitude_image->debug != MagickFalse)
1194 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1195 magnitude_image->filename);
1196 if (phase_image == (Image *) NULL)
1198 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1199 "TwoOrMoreImagesRequired","'%s'",magnitude_image->filename);
1200 return((Image *) NULL);
1202 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1203 fourier_image=(Image *) NULL;
1205 (void) ThrowMagickException(exception,GetMagickModule(),
1206 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
1207 magnitude_image->filename);
1210 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1211 magnitude_image->rows,MagickFalse,exception);
1212 if (fourier_image != (Image *) NULL)
1219 is_gray=IsImageGray(magnitude_image,exception);
1220 if (is_gray != MagickFalse)
1221 is_gray=IsImageGray(phase_image,exception);
1222 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1223 #pragma omp parallel sections
1226 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1233 if (is_gray != MagickFalse)
1234 thread_status=InverseFourierTransformChannel(magnitude_image,
1235 phase_image,GrayPixelChannel,modulus,fourier_image,exception);
1237 thread_status=InverseFourierTransformChannel(magnitude_image,
1238 phase_image,RedPixelChannel,modulus,fourier_image,exception);
1239 if (thread_status == MagickFalse)
1240 status=thread_status;
1242 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1249 thread_status=MagickTrue;
1250 if (is_gray == MagickFalse)
1251 thread_status=InverseFourierTransformChannel(magnitude_image,
1252 phase_image,GreenPixelChannel,modulus,fourier_image,exception);
1253 if (thread_status == MagickFalse)
1254 status=thread_status;
1256 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1263 thread_status=MagickTrue;
1264 if (is_gray == MagickFalse)
1265 thread_status=InverseFourierTransformChannel(magnitude_image,
1266 phase_image,BluePixelChannel,modulus,fourier_image,exception);
1267 if (thread_status == MagickFalse)
1268 status=thread_status;
1270 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1277 thread_status=MagickTrue;
1278 if (magnitude_image->colorspace == CMYKColorspace)
1279 thread_status=InverseFourierTransformChannel(magnitude_image,
1280 phase_image,BlackPixelChannel,modulus,fourier_image,exception);
1281 if (thread_status == MagickFalse)
1282 status=thread_status;
1284 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1291 thread_status=MagickTrue;
1292 if (magnitude_image->matte != MagickFalse)
1293 thread_status=InverseFourierTransformChannel(magnitude_image,
1294 phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
1295 if (thread_status == MagickFalse)
1296 status=thread_status;
1299 if (status == MagickFalse)
1300 fourier_image=DestroyImage(fourier_image);
1305 return(fourier_image);