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 (void) ResetMagickMemory(phase_source,0,fourier_info->height*
275 fourier_info->width*sizeof(*phase_source));
276 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
277 magnitude,magnitude_source);
278 if (status != MagickFalse)
279 status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
281 CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
282 if (fourier_info->modulus != MagickFalse)
285 for (y=0L; y < (ssize_t) fourier_info->height; y++)
286 for (x=0L; x < (ssize_t) fourier_info->width; x++)
288 phase_source[i]/=(2.0*MagickPI);
289 phase_source[i]+=0.5;
293 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
295 for (y=0L; y < (ssize_t) fourier_info->height; y++)
297 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
299 if (q == (Quantum *) NULL)
301 for (x=0L; x < (ssize_t) fourier_info->width; x++)
303 switch (fourier_info->channel)
305 case RedPixelChannel:
308 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
309 magnitude_source[i]),q);
312 case GreenPixelChannel:
314 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
315 magnitude_source[i]),q);
318 case BluePixelChannel:
320 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
321 magnitude_source[i]),q);
324 case BlackPixelChannel:
326 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
327 magnitude_source[i]),q);
330 case AlphaPixelChannel:
332 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
333 magnitude_source[i]),q);
338 q+=GetPixelChannels(magnitude_image);
340 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
341 if (status == MagickFalse)
344 magnitude_view=DestroyCacheView(magnitude_view);
346 phase_view=AcquireAuthenticCacheView(phase_image,exception);
347 for (y=0L; y < (ssize_t) fourier_info->height; y++)
349 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
351 if (q == (Quantum *) NULL)
353 for (x=0L; x < (ssize_t) fourier_info->width; x++)
355 switch (fourier_info->channel)
357 case RedPixelChannel:
360 SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
364 case GreenPixelChannel:
366 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
370 case BluePixelChannel:
372 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
376 case BlackPixelChannel:
378 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
382 case AlphaPixelChannel:
384 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
390 q+=GetPixelChannels(phase_image);
392 status=SyncCacheViewAuthenticPixels(phase_view,exception);
393 if (status == MagickFalse)
396 phase_view=DestroyCacheView(phase_view);
397 phase_source=(double *) RelinquishMagickMemory(phase_source);
398 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
402 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
403 const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
418 register const Quantum
429 Generate the forward Fourier transform.
431 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
432 fourier_info->width*sizeof(*source));
433 if (source == (double *) NULL)
435 (void) ThrowMagickException(exception,GetMagickModule(),
436 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
439 ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
442 image_view=AcquireVirtualCacheView(image,exception);
443 for (y=0L; y < (ssize_t) fourier_info->height; y++)
445 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
447 if (p == (const Quantum *) NULL)
449 for (x=0L; x < (ssize_t) fourier_info->width; x++)
451 switch (fourier_info->channel)
453 case RedPixelChannel:
456 source[i]=QuantumScale*GetPixelRed(image,p);
459 case GreenPixelChannel:
461 source[i]=QuantumScale*GetPixelGreen(image,p);
464 case BluePixelChannel:
466 source[i]=QuantumScale*GetPixelBlue(image,p);
469 case BlackPixelChannel:
471 source[i]=QuantumScale*GetPixelBlack(image,p);
474 case AlphaPixelChannel:
476 source[i]=QuantumScale*GetPixelAlpha(image,p);
481 p+=GetPixelChannels(image);
484 image_view=DestroyCacheView(image_view);
485 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
486 fourier_info->center*sizeof(*fourier));
487 if (fourier == (fftw_complex *) NULL)
489 (void) ThrowMagickException(exception,GetMagickModule(),
490 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
491 source=(double *) RelinquishMagickMemory(source);
494 #if defined(MAGICKCORE_OPENMP_SUPPORT)
495 #pragma omp critical (MagickCore_ForwardFourierTransform)
497 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->height,
498 source,fourier,FFTW_ESTIMATE);
499 fftw_execute(fftw_r2c_plan);
500 fftw_destroy_plan(fftw_r2c_plan);
501 source=(double *) RelinquishMagickMemory(source);
503 Normalize Fourier transform.
505 n=(double) fourier_info->width*(double) fourier_info->width;
507 for (y=0L; y < (ssize_t) fourier_info->height; y++)
508 for (x=0L; x < (ssize_t) fourier_info->center; x++)
510 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
519 Generate magnitude and phase (or real and imaginary).
522 if (fourier_info->modulus != MagickFalse)
523 for (y=0L; y < (ssize_t) fourier_info->height; y++)
524 for (x=0L; x < (ssize_t) fourier_info->center; x++)
526 magnitude[i]=cabs(fourier[i]);
527 phase[i]=carg(fourier[i]);
531 for (y=0L; y < (ssize_t) fourier_info->height; y++)
532 for (x=0L; x < (ssize_t) fourier_info->center; x++)
534 magnitude[i]=creal(fourier[i]);
535 phase[i]=cimag(fourier[i]);
538 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
542 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
543 const PixelChannel channel,const MagickBooleanType modulus,
544 Image *fourier_image,ExceptionInfo *exception)
562 fourier_info.width=image->columns;
563 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
564 ((image->rows % 2) != 0))
566 extent=image->columns < image->rows ? image->rows : image->columns;
567 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
569 fourier_info.height=fourier_info.width;
570 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
571 fourier_info.channel=channel;
572 fourier_info.modulus=modulus;
573 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
574 fourier_info.center*sizeof(*magnitude));
575 if (magnitude == (double *) NULL)
577 (void) ThrowMagickException(exception,GetMagickModule(),
578 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
581 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
582 fourier_info.center*sizeof(*phase));
583 if (phase == (double *) NULL)
585 (void) ThrowMagickException(exception,GetMagickModule(),
586 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
587 magnitude=(double *) RelinquishMagickMemory(magnitude);
590 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
591 fourier_info.center*sizeof(*fourier));
592 if (fourier == (fftw_complex *) NULL)
594 (void) ThrowMagickException(exception,GetMagickModule(),
595 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
596 phase=(double *) RelinquishMagickMemory(phase);
597 magnitude=(double *) RelinquishMagickMemory(magnitude);
600 status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
601 if (status != MagickFalse)
602 status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
604 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
605 phase=(double *) RelinquishMagickMemory(phase);
606 magnitude=(double *) RelinquishMagickMemory(magnitude);
611 MagickExport Image *ForwardFourierTransformImage(const Image *image,
612 const MagickBooleanType modulus,ExceptionInfo *exception)
617 fourier_image=NewImageList();
618 #if !defined(MAGICKCORE_FFTW_DELEGATE)
620 (void) ThrowMagickException(exception,GetMagickModule(),
621 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
632 width=image->columns;
633 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
634 ((image->rows % 2) != 0))
636 extent=image->columns < image->rows ? image->rows : image->columns;
637 width=(extent & 0x01) == 1 ? extent+1UL : extent;
639 magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
640 if (magnitude_image != (Image *) NULL)
645 magnitude_image->storage_class=DirectClass;
646 magnitude_image->depth=32UL;
647 phase_image=CloneImage(image,width,width,MagickFalse,exception);
648 if (phase_image == (Image *) NULL)
649 magnitude_image=DestroyImage(magnitude_image);
656 phase_image->storage_class=DirectClass;
657 phase_image->depth=32UL;
658 AppendImageToList(&fourier_image,magnitude_image);
659 AppendImageToList(&fourier_image,phase_image);
661 is_gray=IsImageGray(image,exception);
662 #if defined(MAGICKCORE_OPENMP_SUPPORT)
663 #pragma omp parallel sections
666 #if defined(MAGICKCORE_OPENMP_SUPPORT)
673 if (is_gray != MagickFalse)
674 thread_status=ForwardFourierTransformChannel(image,
675 GrayPixelChannel,modulus,fourier_image,exception);
677 thread_status=ForwardFourierTransformChannel(image,
678 RedPixelChannel,modulus,fourier_image,exception);
679 if (thread_status == MagickFalse)
680 status=thread_status;
682 #if defined(MAGICKCORE_OPENMP_SUPPORT)
689 thread_status=MagickTrue;
690 if (is_gray == MagickFalse)
691 thread_status=ForwardFourierTransformChannel(image,
692 GreenPixelChannel,modulus,fourier_image,exception);
693 if (thread_status == MagickFalse)
694 status=thread_status;
696 #if defined(MAGICKCORE_OPENMP_SUPPORT)
703 thread_status=MagickTrue;
704 if (is_gray == MagickFalse)
705 thread_status=ForwardFourierTransformChannel(image,
706 BluePixelChannel,modulus,fourier_image,exception);
707 if (thread_status == MagickFalse)
708 status=thread_status;
710 #if defined(MAGICKCORE_OPENMP_SUPPORT)
717 thread_status=MagickTrue;
718 if (image->colorspace == CMYKColorspace)
719 thread_status=ForwardFourierTransformChannel(image,
720 BlackPixelChannel,modulus,fourier_image,exception);
721 if (thread_status == MagickFalse)
722 status=thread_status;
724 #if defined(MAGICKCORE_OPENMP_SUPPORT)
731 thread_status=MagickTrue;
732 if (image->alpha_trait == BlendPixelTrait)
733 thread_status=ForwardFourierTransformChannel(image,
734 AlphaPixelChannel,modulus,fourier_image,exception);
735 if (thread_status == MagickFalse)
736 status=thread_status;
739 if (status == MagickFalse)
740 fourier_image=DestroyImageList(fourier_image);
746 return(fourier_image);
750 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
754 % 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 %
758 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
760 % InverseFourierTransformImage() implements the inverse discrete Fourier
761 % transform (DFT) of the image either as a magnitude / phase or real /
762 % imaginary image pair.
764 % The format of the InverseFourierTransformImage method is:
766 % Image *InverseFourierTransformImage(const Image *magnitude_image,
767 % const Image *phase_image,const MagickBooleanType modulus,
768 % ExceptionInfo *exception)
770 % A description of each parameter follows:
772 % o magnitude_image: the magnitude or real image.
774 % o phase_image: the phase or imaginary image.
776 % o modulus: if true, return transform as a magnitude / phase pair
777 % otherwise a real / imaginary image pair.
779 % o exception: return any errors or warnings in this structure.
783 #if defined(MAGICKCORE_FFTW_DELEGATE)
784 static MagickBooleanType InverseQuadrantSwap(const size_t width,
785 const size_t height,const double *source,double *destination)
797 center=(ssize_t) floor((double) width/2L)+1L;
798 for (y=1L; y < (ssize_t) height; y++)
799 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
800 destination[center*(height-y)-x+width/2L]=source[y*width+x];
801 for (y=0L; y < (ssize_t) height; y++)
802 destination[center*y]=source[y*width+width/2L];
803 for (x=0L; x < center; x++)
804 destination[x]=source[center-x-1L];
805 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
808 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
809 const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
810 ExceptionInfo *exception)
825 register const Quantum
836 Inverse fourier - read image and break down into a double array.
838 magnitude_source=(double *) AcquireQuantumMemory((size_t)
839 fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
840 if (magnitude_source == (double *) NULL)
842 (void) ThrowMagickException(exception,GetMagickModule(),
843 ResourceLimitError,"MemoryAllocationFailed","'%s'",
844 magnitude_image->filename);
847 phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
848 fourier_info->width*sizeof(*phase_source));
849 if (phase_source == (double *) NULL)
851 (void) ThrowMagickException(exception,GetMagickModule(),
852 ResourceLimitError,"MemoryAllocationFailed","'%s'",
853 magnitude_image->filename);
854 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
858 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
859 for (y=0L; y < (ssize_t) fourier_info->height; y++)
861 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
863 if (p == (const Quantum *) NULL)
865 for (x=0L; x < (ssize_t) fourier_info->width; x++)
867 switch (fourier_info->channel)
869 case RedPixelChannel:
872 magnitude_source[i]=QuantumScale*GetPixelRed(magnitude_image,p);
875 case GreenPixelChannel:
877 magnitude_source[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
880 case BluePixelChannel:
882 magnitude_source[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
885 case BlackPixelChannel:
887 magnitude_source[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
890 case AlphaPixelChannel:
892 magnitude_source[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
897 p+=GetPixelChannels(magnitude_image);
901 phase_view=AcquireVirtualCacheView(phase_image,exception);
902 for (y=0L; y < (ssize_t) fourier_info->height; y++)
904 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
906 if (p == (const Quantum *) NULL)
908 for (x=0L; x < (ssize_t) fourier_info->width; x++)
910 switch (fourier_info->channel)
912 case RedPixelChannel:
915 phase_source[i]=QuantumScale*GetPixelRed(phase_image,p);
918 case GreenPixelChannel:
920 phase_source[i]=QuantumScale*GetPixelGreen(phase_image,p);
923 case BluePixelChannel:
925 phase_source[i]=QuantumScale*GetPixelBlue(phase_image,p);
928 case BlackPixelChannel:
930 phase_source[i]=QuantumScale*GetPixelBlack(phase_image,p);
933 case AlphaPixelChannel:
935 phase_source[i]=QuantumScale*GetPixelAlpha(phase_image,p);
940 p+=GetPixelChannels(phase_image);
943 if (fourier_info->modulus != MagickFalse)
946 for (y=0L; y < (ssize_t) fourier_info->height; y++)
947 for (x=0L; x < (ssize_t) fourier_info->width; x++)
949 phase_source[i]-=0.5;
950 phase_source[i]*=(2.0*MagickPI);
954 magnitude_view=DestroyCacheView(magnitude_view);
955 phase_view=DestroyCacheView(phase_view);
956 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
957 fourier_info->center*sizeof(*magnitude));
958 if (magnitude == (double *) NULL)
960 (void) ThrowMagickException(exception,GetMagickModule(),
961 ResourceLimitError,"MemoryAllocationFailed","'%s'",
962 magnitude_image->filename);
963 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
964 phase_source=(double *) RelinquishMagickMemory(phase_source);
967 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
968 magnitude_source,magnitude);
969 magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
970 phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
971 fourier_info->width*sizeof(*phase));
972 if (phase == (double *) NULL)
974 (void) ThrowMagickException(exception,GetMagickModule(),
975 ResourceLimitError,"MemoryAllocationFailed","'%s'",
976 magnitude_image->filename);
977 phase_source=(double *) RelinquishMagickMemory(phase_source);
980 CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
981 if (status != MagickFalse)
982 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
984 phase_source=(double *) RelinquishMagickMemory(phase_source);
989 if (fourier_info->modulus != MagickFalse)
990 for (y=0L; y < (ssize_t) fourier_info->height; y++)
991 for (x=0L; x < (ssize_t) fourier_info->center; x++)
993 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
994 fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
996 fourier[i][0]=magnitude[i]*cos(phase[i]);
997 fourier[i][1]=magnitude[i]*sin(phase[i]);
1002 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1003 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1005 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1006 fourier[i]=magnitude[i]+I*phase[i];
1008 fourier[i][0]=magnitude[i];
1009 fourier[i][1]=phase[i];
1013 phase=(double *) RelinquishMagickMemory(phase);
1014 magnitude=(double *) RelinquishMagickMemory(magnitude);
1018 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1019 fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1040 source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1041 fourier_info->width*sizeof(*source));
1042 if (source == (double *) NULL)
1044 (void) ThrowMagickException(exception,GetMagickModule(),
1045 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
1046 return(MagickFalse);
1048 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1049 #pragma omp critical (MagickCore_InverseFourierTransform)
1052 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1053 fourier,source,FFTW_ESTIMATE);
1054 fftw_execute(fftw_c2r_plan);
1055 fftw_destroy_plan(fftw_c2r_plan);
1058 image_view=AcquireAuthenticCacheView(image,exception);
1059 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1061 if (y >= (ssize_t) image->rows)
1063 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1064 image->columns ? image->columns : fourier_info->width,1UL,exception);
1065 if (q == (Quantum *) NULL)
1067 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1069 if (x < (ssize_t) image->columns)
1070 switch (fourier_info->channel)
1072 case RedPixelChannel:
1075 SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
1078 case GreenPixelChannel:
1080 SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
1083 case BluePixelChannel:
1085 SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
1088 case BlackPixelChannel:
1090 SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
1093 case AlphaPixelChannel:
1095 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
1100 q+=GetPixelChannels(image);
1102 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1105 image_view=DestroyCacheView(image_view);
1106 source=(double *) RelinquishMagickMemory(source);
1110 static MagickBooleanType InverseFourierTransformChannel(
1111 const Image *magnitude_image,const Image *phase_image,
1112 const PixelChannel channel,const MagickBooleanType modulus,
1113 Image *fourier_image,ExceptionInfo *exception)
1131 fourier_info.width=magnitude_image->columns;
1132 if ((magnitude_image->columns != magnitude_image->rows) ||
1133 ((magnitude_image->columns % 2) != 0) ||
1134 ((magnitude_image->rows % 2) != 0))
1136 extent=magnitude_image->columns < magnitude_image->rows ?
1137 magnitude_image->rows : magnitude_image->columns;
1138 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1140 fourier_info.height=fourier_info.width;
1141 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
1142 fourier_info.channel=channel;
1143 fourier_info.modulus=modulus;
1144 magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1145 fourier_info.center*sizeof(*magnitude));
1146 if (magnitude == (double *) NULL)
1148 (void) ThrowMagickException(exception,GetMagickModule(),
1149 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1150 magnitude_image->filename);
1151 return(MagickFalse);
1153 phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1154 fourier_info.center*sizeof(*phase));
1155 if (phase == (double *) NULL)
1157 (void) ThrowMagickException(exception,GetMagickModule(),
1158 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1159 magnitude_image->filename);
1160 magnitude=(double *) RelinquishMagickMemory(magnitude);
1161 return(MagickFalse);
1163 fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1164 fourier_info.center*sizeof(*fourier));
1165 if (fourier == (fftw_complex *) NULL)
1167 (void) ThrowMagickException(exception,GetMagickModule(),
1168 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1169 magnitude_image->filename);
1170 phase=(double *) RelinquishMagickMemory(phase);
1171 magnitude=(double *) RelinquishMagickMemory(magnitude);
1172 return(MagickFalse);
1174 status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1176 if (status != MagickFalse)
1177 status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1179 fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1180 phase=(double *) RelinquishMagickMemory(phase);
1181 magnitude=(double *) RelinquishMagickMemory(magnitude);
1186 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1187 const Image *phase_image,const MagickBooleanType modulus,
1188 ExceptionInfo *exception)
1193 assert(magnitude_image != (Image *) NULL);
1194 assert(magnitude_image->signature == MagickSignature);
1195 if (magnitude_image->debug != MagickFalse)
1196 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1197 magnitude_image->filename);
1198 if (phase_image == (Image *) NULL)
1200 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1201 "TwoOrMoreImagesRequired","'%s'",magnitude_image->filename);
1202 return((Image *) NULL);
1204 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1205 fourier_image=(Image *) NULL;
1207 (void) ThrowMagickException(exception,GetMagickModule(),
1208 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
1209 magnitude_image->filename);
1212 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1213 magnitude_image->rows,MagickFalse,exception);
1214 if (fourier_image != (Image *) NULL)
1221 is_gray=IsImageGray(magnitude_image,exception);
1222 if (is_gray != MagickFalse)
1223 is_gray=IsImageGray(phase_image,exception);
1224 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1225 #pragma omp parallel sections
1228 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1235 if (is_gray != MagickFalse)
1236 thread_status=InverseFourierTransformChannel(magnitude_image,
1237 phase_image,GrayPixelChannel,modulus,fourier_image,exception);
1239 thread_status=InverseFourierTransformChannel(magnitude_image,
1240 phase_image,RedPixelChannel,modulus,fourier_image,exception);
1241 if (thread_status == MagickFalse)
1242 status=thread_status;
1244 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1251 thread_status=MagickTrue;
1252 if (is_gray == MagickFalse)
1253 thread_status=InverseFourierTransformChannel(magnitude_image,
1254 phase_image,GreenPixelChannel,modulus,fourier_image,exception);
1255 if (thread_status == MagickFalse)
1256 status=thread_status;
1258 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1265 thread_status=MagickTrue;
1266 if (is_gray == MagickFalse)
1267 thread_status=InverseFourierTransformChannel(magnitude_image,
1268 phase_image,BluePixelChannel,modulus,fourier_image,exception);
1269 if (thread_status == MagickFalse)
1270 status=thread_status;
1272 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1279 thread_status=MagickTrue;
1280 if (magnitude_image->colorspace == CMYKColorspace)
1281 thread_status=InverseFourierTransformChannel(magnitude_image,
1282 phase_image,BlackPixelChannel,modulus,fourier_image,exception);
1283 if (thread_status == MagickFalse)
1284 status=thread_status;
1286 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1293 thread_status=MagickTrue;
1294 if (magnitude_image->alpha_trait == BlendPixelTrait)
1295 thread_status=InverseFourierTransformChannel(magnitude_image,
1296 phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
1297 if (thread_status == MagickFalse)
1298 status=thread_status;
1301 if (status == MagickFalse)
1302 fourier_image=DestroyImage(fourier_image);
1307 return(fourier_image);