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-2014 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/artifact.h"
47 #include "MagickCore/attribute.h"
48 #include "MagickCore/blob.h"
49 #include "MagickCore/cache.h"
50 #include "MagickCore/image.h"
51 #include "MagickCore/image-private.h"
52 #include "MagickCore/list.h"
53 #include "MagickCore/fourier.h"
54 #include "MagickCore/log.h"
55 #include "MagickCore/memory_.h"
56 #include "MagickCore/monitor.h"
57 #include "MagickCore/monitor-private.h"
58 #include "MagickCore/pixel-accessor.h"
59 #include "MagickCore/pixel-private.h"
60 #include "MagickCore/property.h"
61 #include "MagickCore/quantum-private.h"
62 #include "MagickCore/resource_.h"
63 #include "MagickCore/string-private.h"
64 #include "MagickCore/thread-private.h"
65 #if defined(MAGICKCORE_FFTW_DELEGATE)
66 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
70 #if !defined(MAGICKCORE_HAVE_CABS)
71 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
73 #if !defined(MAGICKCORE_HAVE_CARG)
74 #define carg(z) (atan2(cimag(z),creal(z)))
76 #if !defined(MAGICKCORE_HAVE_CIMAG)
77 #define cimag(z) (z[1])
79 #if !defined(MAGICKCORE_HAVE_CREAL)
80 #define creal(z) (z[0])
87 typedef struct _FourierInfo
104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 % C o m p l e x I m a g e s %
112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
114 % ComplexImages() performs complex mathematics on an image sequence.
116 % The format of the ComplexImages method is:
118 % MagickBooleanType ComplexImages(Image *images,
119 % const ComplexOperator op,ExceptionInfo *exception)
121 % A description of each parameter follows:
123 % o image: the image.
125 % o op: A complex op.
127 % o exception: return any errors or warnings in this structure.
130 MagickExport Image *ComplexImages(const Image *images,const ComplexOperator op,
131 ExceptionInfo *exception)
133 #define ComplexImageTag "Complex/Image"
170 assert(images != (Image *) NULL);
171 assert(images->signature == MagickSignature);
172 if (images->debug != MagickFalse)
173 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
174 assert(exception != (ExceptionInfo *) NULL);
175 assert(exception->signature == MagickSignature);
176 if (images->next == (Image *) NULL)
178 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
179 "ImageSequenceRequired","`%s'",images->filename);
180 return((Image *) NULL);
182 image=CloneImage(images,images->columns,images->rows,MagickTrue,exception);
183 if (image == (Image *) NULL)
184 return((Image *) NULL);
185 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
187 image=DestroyImageList(image);
190 complex_images=NewImageList();
191 AppendImageToList(&complex_images,image);
192 image=CloneImage(images,images->columns,images->rows,MagickTrue,exception);
193 if (image == (Image *) NULL)
195 complex_images=DestroyImageList(complex_images);
196 return(complex_images);
198 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
200 complex_images=DestroyImageList(complex_images);
201 return(complex_images);
203 AppendImageToList(&complex_images,image);
205 Apply complex mathematics to image pixels.
207 artifact=GetImageArtifact(image,"complex:snr");
209 if (artifact != (const char *) NULL)
210 snr=StringToDouble(artifact,(char **) NULL);
212 Ai_image=images->next;
214 Bi_image=images->next;
215 if ((images->next->next != (Image *) NULL) &&
216 (images->next->next->next != (Image *) NULL))
218 Br_image=images->next->next;
219 Bi_image=images->next->next->next;
221 Cr_image=complex_images;
222 Ci_image=complex_images->next;
223 Ar_view=AcquireVirtualCacheView(Ar_image,exception);
224 Ai_view=AcquireVirtualCacheView(Ai_image,exception);
225 Br_view=AcquireVirtualCacheView(Br_image,exception);
226 Bi_view=AcquireVirtualCacheView(Bi_image,exception);
227 Cr_view=AcquireAuthenticCacheView(Cr_image,exception);
228 Ci_view=AcquireAuthenticCacheView(Ci_image,exception);
231 #if defined(MAGICKCORE_OPENMP_SUPPORT)
232 #pragma omp parallel for schedule(static,4) shared(progress,status) \
233 magick_threads(images,complex_images,images->rows,1L)
235 for (y=0; y < (ssize_t) images->rows; y++)
237 register const Quantum
250 if (status == MagickFalse)
252 Ar=GetCacheViewVirtualPixels(Ar_view,0,y,images->columns,1,exception);
253 Ai=GetCacheViewVirtualPixels(Ai_view,0,y,images->columns,1,exception);
254 Br=GetCacheViewVirtualPixels(Br_view,0,y,images->columns,1,exception);
255 Bi=GetCacheViewVirtualPixels(Bi_view,0,y,images->columns,1,exception);
256 Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,images->columns,1,exception);
257 Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,images->columns,1,exception);
258 if ((Ar == (const Quantum *) NULL) || (Ai == (const Quantum *) NULL) ||
259 (Br == (const Quantum *) NULL) || (Bi == (const Quantum *) NULL) ||
260 (Cr == (Quantum *) NULL) || (Ci == (Quantum *) NULL))
265 for (x=0; x < (ssize_t) images->columns; x++)
270 for (i=0; i < (ssize_t) GetPixelChannels(images); i++)
274 case AddComplexOperator:
280 case ConjugateComplexOperator:
287 case DivideComplexOperator:
292 gamma=PerceptibleReciprocal(Br[i]*Br[i]+Bi[i]*Bi[i]+snr);
293 Cr[i]=gamma*(Ar[i]*Br[i]+Ai[i]*Bi[i]);
294 Ci[i]=gamma*(Ai[i]*Br[i]-Ar[i]*Bi[i]);
297 case MagnitudePhaseComplexOperator:
299 Cr[i]=sqrt(Ar[i]*Ar[i]+Ai[i]*Ai[i]);
300 Ci[i]=atan2(Ai[i],Ar[i])/(2.0*MagickPI)+0.5;
303 case MultiplyComplexOperator:
305 Cr[i]=QuantumScale*(Ar[i]*Br[i]-Ai[i]*Bi[i]);
306 Ci[i]=QuantumScale*(Ai[i]*Br[i]+Ar[i]*Bi[i]);
309 case RealImaginaryComplexOperator:
311 Cr[i]=Ar[i]*cos(2.0*MagickPI*(Ai[i]-0.5));
312 Ci[i]=Ar[i]*sin(2.0*MagickPI*(Ai[i]-0.5));
315 case SubtractComplexOperator:
322 Ar+=GetPixelChannels(Ar_image);
323 Ai+=GetPixelChannels(Ai_image);
324 Br+=GetPixelChannels(Br_image);
325 Bi+=GetPixelChannels(Bi_image);
326 Cr+=GetPixelChannels(Cr_image);
327 Ci+=GetPixelChannels(Ci_image);
330 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
332 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
334 if (images->progress_monitor != (MagickProgressMonitor) NULL)
339 #if defined(MAGICKCORE_OPENMP_SUPPORT)
340 #pragma omp critical (MagickCore_ComplexImages)
342 proceed=SetImageProgress(images,ComplexImageTag,progress++,
344 if (proceed == MagickFalse)
348 Cr_view=DestroyCacheView(Cr_view);
349 Ci_view=DestroyCacheView(Ci_view);
350 Br_view=DestroyCacheView(Br_view);
351 Bi_view=DestroyCacheView(Bi_view);
352 Ar_view=DestroyCacheView(Ar_view);
353 Ai_view=DestroyCacheView(Ai_view);
354 if (status == MagickFalse)
355 complex_images=DestroyImageList(complex_images);
356 return(complex_images);
360 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
364 % 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 %
368 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
370 % ForwardFourierTransformImage() implements the discrete Fourier transform
371 % (DFT) of the image either as a magnitude / phase or real / imaginary image
374 % The format of the ForwadFourierTransformImage method is:
376 % Image *ForwardFourierTransformImage(const Image *image,
377 % const MagickBooleanType modulus,ExceptionInfo *exception)
379 % A description of each parameter follows:
381 % o image: the image.
383 % o modulus: if true, return as transform as a magnitude / phase pair
384 % otherwise a real / imaginary image pair.
386 % o exception: return any errors or warnings in this structure.
390 #if defined(MAGICKCORE_FFTW_DELEGATE)
392 static MagickBooleanType RollFourier(const size_t width,const size_t height,
393 const ssize_t x_offset,const ssize_t y_offset,double *roll_pixels)
411 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
413 source_info=AcquireVirtualMemory(height,width*sizeof(*source_pixels));
414 if (source_info == (MemoryInfo *) NULL)
416 source_pixels=(double *) GetVirtualMemoryBlob(source_info);
418 for (y=0L; y < (ssize_t) height; y++)
421 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
423 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
425 for (x=0L; x < (ssize_t) width; x++)
428 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
430 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
432 source_pixels[v*width+u]=roll_pixels[i++];
435 (void) CopyMagickMemory(roll_pixels,source_pixels,height*width*
436 sizeof(*source_pixels));
437 source_info=RelinquishVirtualMemory(source_info);
441 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
442 const size_t height,double *source_pixels,double *forward_pixels)
457 center=(ssize_t) floor((double) width/2L)+1L;
458 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,
460 if (status == MagickFalse)
462 for (y=0L; y < (ssize_t) height; y++)
463 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
464 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
465 for (y=1; y < (ssize_t) height; y++)
466 for (x=0L; x < (ssize_t) (width/2L-1L); x++)
467 forward_pixels[(height-y)*width+width/2L-x-1L]=
468 source_pixels[y*center+x+1L];
469 for (x=0L; x < (ssize_t) (width/2L); x++)
470 forward_pixels[-x+width/2L-1L]=forward_pixels[x+width/2L+1L];
474 static void CorrectPhaseLHS(const size_t width,const size_t height,
475 double *fourier_pixels)
483 for (y=0L; y < (ssize_t) height; y++)
484 for (x=0L; x < (ssize_t) (width/2L); x++)
485 fourier_pixels[y*width+x]*=(-1.0);
488 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
489 Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
520 magnitude_image=GetFirstImageInList(image);
521 phase_image=GetNextImageInList(image);
522 if (phase_image == (Image *) NULL)
524 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
525 "TwoOrMoreImagesRequired","`%s'",image->filename);
529 Create "Fourier Transform" image from constituent arrays.
531 magnitude_info=AcquireVirtualMemory((size_t) fourier_info->height,
532 fourier_info->width*sizeof(*magnitude_pixels));
533 phase_info=AcquireVirtualMemory((size_t) fourier_info->height,
534 fourier_info->width*sizeof(*phase_pixels));
535 if ((magnitude_info == (MemoryInfo *) NULL) ||
536 (phase_info == (MemoryInfo *) NULL))
538 if (phase_info != (MemoryInfo *) NULL)
539 phase_info=RelinquishVirtualMemory(phase_info);
540 if (magnitude_info != (MemoryInfo *) NULL)
541 magnitude_info=RelinquishVirtualMemory(magnitude_info);
542 (void) ThrowMagickException(exception,GetMagickModule(),
543 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
546 magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
547 (void) ResetMagickMemory(magnitude_pixels,0,fourier_info->height*
548 fourier_info->width*sizeof(*magnitude_pixels));
549 phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
550 (void) ResetMagickMemory(phase_pixels,0,fourier_info->height*
551 fourier_info->width*sizeof(*phase_pixels));
552 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
553 magnitude,magnitude_pixels);
554 if (status != MagickFalse)
555 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
557 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
558 if (fourier_info->modulus != MagickFalse)
561 for (y=0L; y < (ssize_t) fourier_info->height; y++)
562 for (x=0L; x < (ssize_t) fourier_info->width; x++)
564 phase_pixels[i]/=(2.0*MagickPI);
565 phase_pixels[i]+=0.5;
569 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
571 for (y=0L; y < (ssize_t) fourier_info->height; y++)
573 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
575 if (q == (Quantum *) NULL)
577 for (x=0L; x < (ssize_t) fourier_info->width; x++)
579 switch (fourier_info->channel)
581 case RedPixelChannel:
584 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
585 magnitude_pixels[i]),q);
588 case GreenPixelChannel:
590 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
591 magnitude_pixels[i]),q);
594 case BluePixelChannel:
596 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
597 magnitude_pixels[i]),q);
600 case BlackPixelChannel:
602 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
603 magnitude_pixels[i]),q);
606 case AlphaPixelChannel:
608 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
609 magnitude_pixels[i]),q);
614 q+=GetPixelChannels(magnitude_image);
616 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
617 if (status == MagickFalse)
620 magnitude_view=DestroyCacheView(magnitude_view);
622 phase_view=AcquireAuthenticCacheView(phase_image,exception);
623 for (y=0L; y < (ssize_t) fourier_info->height; y++)
625 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
627 if (q == (Quantum *) NULL)
629 for (x=0L; x < (ssize_t) fourier_info->width; x++)
631 switch (fourier_info->channel)
633 case RedPixelChannel:
636 SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
640 case GreenPixelChannel:
642 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
646 case BluePixelChannel:
648 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
652 case BlackPixelChannel:
654 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
658 case AlphaPixelChannel:
660 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
666 q+=GetPixelChannels(phase_image);
668 status=SyncCacheViewAuthenticPixels(phase_view,exception);
669 if (status == MagickFalse)
672 phase_view=DestroyCacheView(phase_view);
673 phase_info=RelinquishVirtualMemory(phase_info);
674 magnitude_info=RelinquishVirtualMemory(magnitude_info);
678 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
679 const Image *image,double *magnitude_pixels,double *phase_pixels,
680 ExceptionInfo *exception)
701 register const Quantum
712 Generate the forward Fourier transform.
714 source_info=AcquireVirtualMemory((size_t) fourier_info->height,
715 fourier_info->width*sizeof(*source_pixels));
716 if (source_info == (MemoryInfo *) NULL)
718 (void) ThrowMagickException(exception,GetMagickModule(),
719 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
722 source_pixels=(double *) GetVirtualMemoryBlob(source_info);
723 ResetMagickMemory(source_pixels,0,fourier_info->height*fourier_info->width*
724 sizeof(*source_pixels));
726 image_view=AcquireVirtualCacheView(image,exception);
727 for (y=0L; y < (ssize_t) fourier_info->height; y++)
729 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
731 if (p == (const Quantum *) NULL)
733 for (x=0L; x < (ssize_t) fourier_info->width; x++)
735 switch (fourier_info->channel)
737 case RedPixelChannel:
740 source_pixels[i]=QuantumScale*GetPixelRed(image,p);
743 case GreenPixelChannel:
745 source_pixels[i]=QuantumScale*GetPixelGreen(image,p);
748 case BluePixelChannel:
750 source_pixels[i]=QuantumScale*GetPixelBlue(image,p);
753 case BlackPixelChannel:
755 source_pixels[i]=QuantumScale*GetPixelBlack(image,p);
758 case AlphaPixelChannel:
760 source_pixels[i]=QuantumScale*GetPixelAlpha(image,p);
765 p+=GetPixelChannels(image);
768 image_view=DestroyCacheView(image_view);
769 forward_info=AcquireVirtualMemory((size_t) fourier_info->height,
770 fourier_info->center*sizeof(*forward_pixels));
771 if (forward_info == (MemoryInfo *) NULL)
773 (void) ThrowMagickException(exception,GetMagickModule(),
774 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
775 source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
778 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
779 #if defined(MAGICKCORE_OPENMP_SUPPORT)
780 #pragma omp critical (MagickCore_ForwardFourierTransform)
782 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->height,
783 source_pixels,forward_pixels,FFTW_ESTIMATE);
784 fftw_execute(fftw_r2c_plan);
785 fftw_destroy_plan(fftw_r2c_plan);
786 source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
787 value=GetImageArtifact(image,"fourier:normalize");
788 if ((value == (const char *) NULL) || (LocaleCompare(value,"forward") == 0))
794 Normalize Fourier transform.
797 gamma=PerceptibleReciprocal((double) fourier_info->width*
798 fourier_info->height);
799 for (y=0L; y < (ssize_t) fourier_info->height; y++)
800 for (x=0L; x < (ssize_t) fourier_info->center; x++)
802 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
803 forward_pixels[i]*=gamma;
805 forward_pixels[i][0]*=gamma;
806 forward_pixels[i][1]*=gamma;
812 Generate magnitude and phase (or real and imaginary).
815 if (fourier_info->modulus != MagickFalse)
816 for (y=0L; y < (ssize_t) fourier_info->height; y++)
817 for (x=0L; x < (ssize_t) fourier_info->center; x++)
819 magnitude_pixels[i]=cabs(forward_pixels[i]);
820 phase_pixels[i]=carg(forward_pixels[i]);
824 for (y=0L; y < (ssize_t) fourier_info->height; y++)
825 for (x=0L; x < (ssize_t) fourier_info->center; x++)
827 magnitude_pixels[i]=creal(forward_pixels[i]);
828 phase_pixels[i]=cimag(forward_pixels[i]);
831 forward_info=(MemoryInfo *) RelinquishVirtualMemory(forward_info);
835 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
836 const PixelChannel channel,const MagickBooleanType modulus,
837 Image *fourier_image,ExceptionInfo *exception)
856 fourier_info.width=image->columns;
857 fourier_info.height=image->rows;
858 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
859 ((image->rows % 2) != 0))
861 extent=image->columns < image->rows ? image->rows : image->columns;
862 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
864 fourier_info.height=fourier_info.width;
865 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
866 fourier_info.channel=channel;
867 fourier_info.modulus=modulus;
868 magnitude_info=AcquireVirtualMemory((size_t) fourier_info.height,
869 fourier_info.center*sizeof(*magnitude_pixels));
870 phase_info=AcquireVirtualMemory((size_t) fourier_info.height,
871 fourier_info.center*sizeof(*phase_pixels));
872 if ((magnitude_info == (MemoryInfo *) NULL) ||
873 (phase_info == (MemoryInfo *) NULL))
875 if (phase_info != (MemoryInfo *) NULL)
876 phase_info=RelinquishVirtualMemory(phase_info);
877 if (magnitude_info == (MemoryInfo *) NULL)
878 magnitude_info=RelinquishVirtualMemory(magnitude_info);
879 (void) ThrowMagickException(exception,GetMagickModule(),
880 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
883 magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
884 phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
885 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
886 phase_pixels,exception);
887 if (status != MagickFalse)
888 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
889 phase_pixels,exception);
890 phase_info=RelinquishVirtualMemory(phase_info);
891 magnitude_info=RelinquishVirtualMemory(magnitude_info);
896 MagickExport Image *ForwardFourierTransformImage(const Image *image,
897 const MagickBooleanType modulus,ExceptionInfo *exception)
902 fourier_image=NewImageList();
903 #if !defined(MAGICKCORE_FFTW_DELEGATE)
905 (void) ThrowMagickException(exception,GetMagickModule(),
906 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
918 width=image->columns;
920 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
921 ((image->rows % 2) != 0))
923 extent=image->columns < image->rows ? image->rows : image->columns;
924 width=(extent & 0x01) == 1 ? extent+1UL : extent;
927 magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
928 if (magnitude_image != (Image *) NULL)
933 magnitude_image->storage_class=DirectClass;
934 magnitude_image->depth=32UL;
935 phase_image=CloneImage(image,width,height,MagickTrue,exception);
936 if (phase_image == (Image *) NULL)
937 magnitude_image=DestroyImage(magnitude_image);
944 phase_image->storage_class=DirectClass;
945 phase_image->depth=32UL;
946 AppendImageToList(&fourier_image,magnitude_image);
947 AppendImageToList(&fourier_image,phase_image);
949 is_gray=IsImageGray(image,exception);
950 #if defined(MAGICKCORE_OPENMP_SUPPORT)
951 #pragma omp parallel sections
954 #if defined(MAGICKCORE_OPENMP_SUPPORT)
961 if (is_gray != MagickFalse)
962 thread_status=ForwardFourierTransformChannel(image,
963 GrayPixelChannel,modulus,fourier_image,exception);
965 thread_status=ForwardFourierTransformChannel(image,
966 RedPixelChannel,modulus,fourier_image,exception);
967 if (thread_status == MagickFalse)
968 status=thread_status;
970 #if defined(MAGICKCORE_OPENMP_SUPPORT)
977 thread_status=MagickTrue;
978 if (is_gray == MagickFalse)
979 thread_status=ForwardFourierTransformChannel(image,
980 GreenPixelChannel,modulus,fourier_image,exception);
981 if (thread_status == MagickFalse)
982 status=thread_status;
984 #if defined(MAGICKCORE_OPENMP_SUPPORT)
991 thread_status=MagickTrue;
992 if (is_gray == MagickFalse)
993 thread_status=ForwardFourierTransformChannel(image,
994 BluePixelChannel,modulus,fourier_image,exception);
995 if (thread_status == MagickFalse)
996 status=thread_status;
998 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1005 thread_status=MagickTrue;
1006 if (image->colorspace == CMYKColorspace)
1007 thread_status=ForwardFourierTransformChannel(image,
1008 BlackPixelChannel,modulus,fourier_image,exception);
1009 if (thread_status == MagickFalse)
1010 status=thread_status;
1012 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1019 thread_status=MagickTrue;
1020 if (image->alpha_trait == BlendPixelTrait)
1021 thread_status=ForwardFourierTransformChannel(image,
1022 AlphaPixelChannel,modulus,fourier_image,exception);
1023 if (thread_status == MagickFalse)
1024 status=thread_status;
1027 if (status == MagickFalse)
1028 fourier_image=DestroyImageList(fourier_image);
1034 return(fourier_image);
1038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1042 % 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 %
1046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1048 % InverseFourierTransformImage() implements the inverse discrete Fourier
1049 % transform (DFT) of the image either as a magnitude / phase or real /
1050 % imaginary image pair.
1052 % The format of the InverseFourierTransformImage method is:
1054 % Image *InverseFourierTransformImage(const Image *magnitude_image,
1055 % const Image *phase_image,const MagickBooleanType modulus,
1056 % ExceptionInfo *exception)
1058 % A description of each parameter follows:
1060 % o magnitude_image: the magnitude or real image.
1062 % o phase_image: the phase or imaginary image.
1064 % o modulus: if true, return transform as a magnitude / phase pair
1065 % otherwise a real / imaginary image pair.
1067 % o exception: return any errors or warnings in this structure.
1071 #if defined(MAGICKCORE_FFTW_DELEGATE)
1072 static MagickBooleanType InverseQuadrantSwap(const size_t width,
1073 const size_t height,const double *source,double *destination)
1085 center=(ssize_t) floor((double) width/2L)+1L;
1086 for (y=1L; y < (ssize_t) height; y++)
1087 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1088 destination[(height-y)*center-x+width/2L]=source[y*width+x];
1089 for (y=0L; y < (ssize_t) height; y++)
1090 destination[y*center]=source[y*width+width/2L];
1091 for (x=0L; x < center; x++)
1092 destination[x]=source[center-x-1L];
1093 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1096 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
1097 const Image *magnitude_image,const Image *phase_image,
1098 fftw_complex *fourier_pixels,ExceptionInfo *exception)
1117 register const Quantum
1128 Inverse fourier - read image and break down into a double array.
1130 magnitude_info=AcquireVirtualMemory((size_t)fourier_info->height,
1131 fourier_info->width*sizeof(*magnitude_pixels));
1132 phase_info=AcquireVirtualMemory((size_t) fourier_info->height,
1133 fourier_info->width*sizeof(*phase_pixels));
1134 inverse_info=AcquireVirtualMemory((size_t) fourier_info->height,
1135 fourier_info->center*sizeof(*inverse_pixels));
1136 if ((magnitude_info == (MemoryInfo *) NULL) ||
1137 (phase_info == (MemoryInfo *) NULL) ||
1138 (inverse_info == (MemoryInfo *) NULL))
1140 if (magnitude_info != (MemoryInfo *) NULL)
1141 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1142 if (phase_info != (MemoryInfo *) NULL)
1143 phase_info=RelinquishVirtualMemory(phase_info);
1144 if (inverse_info != (MemoryInfo *) NULL)
1145 inverse_info=RelinquishVirtualMemory(inverse_info);
1146 (void) ThrowMagickException(exception,GetMagickModule(),
1147 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1148 magnitude_image->filename);
1149 return(MagickFalse);
1151 magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
1152 phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
1153 inverse_pixels=(double *) GetVirtualMemoryBlob(inverse_info);
1155 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1156 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1158 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1160 if (p == (const Quantum *) NULL)
1162 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1164 switch (fourier_info->channel)
1166 case RedPixelChannel:
1169 magnitude_pixels[i]=QuantumScale*GetPixelRed(magnitude_image,p);
1172 case GreenPixelChannel:
1174 magnitude_pixels[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
1177 case BluePixelChannel:
1179 magnitude_pixels[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
1182 case BlackPixelChannel:
1184 magnitude_pixels[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
1187 case AlphaPixelChannel:
1189 magnitude_pixels[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
1194 p+=GetPixelChannels(magnitude_image);
1197 magnitude_view=DestroyCacheView(magnitude_view);
1198 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1199 magnitude_pixels,inverse_pixels);
1200 (void) CopyMagickMemory(magnitude_pixels,inverse_pixels,fourier_info->height*
1201 fourier_info->center*sizeof(*magnitude_pixels));
1203 phase_view=AcquireVirtualCacheView(phase_image,exception);
1204 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1206 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1208 if (p == (const Quantum *) NULL)
1210 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1212 switch (fourier_info->channel)
1214 case RedPixelChannel:
1217 phase_pixels[i]=QuantumScale*GetPixelRed(phase_image,p);
1220 case GreenPixelChannel:
1222 phase_pixels[i]=QuantumScale*GetPixelGreen(phase_image,p);
1225 case BluePixelChannel:
1227 phase_pixels[i]=QuantumScale*GetPixelBlue(phase_image,p);
1230 case BlackPixelChannel:
1232 phase_pixels[i]=QuantumScale*GetPixelBlack(phase_image,p);
1235 case AlphaPixelChannel:
1237 phase_pixels[i]=QuantumScale*GetPixelAlpha(phase_image,p);
1242 p+=GetPixelChannels(phase_image);
1245 if (fourier_info->modulus != MagickFalse)
1248 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1249 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1251 phase_pixels[i]-=0.5;
1252 phase_pixels[i]*=(2.0*MagickPI);
1256 phase_view=DestroyCacheView(phase_view);
1257 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1258 if (status != MagickFalse)
1259 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1260 phase_pixels,inverse_pixels);
1261 (void) CopyMagickMemory(phase_pixels,inverse_pixels,fourier_info->height*
1262 fourier_info->center*sizeof(*phase_pixels));
1263 inverse_info=RelinquishVirtualMemory(inverse_info);
1268 if (fourier_info->modulus != MagickFalse)
1269 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1270 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1272 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1273 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1274 magnitude_pixels[i]*sin(phase_pixels[i]);
1276 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1277 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1282 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1283 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1285 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1286 fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1288 fourier_pixels[i][0]=magnitude_pixels[i];
1289 fourier_pixels[i][1]=phase_pixels[i];
1293 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1294 phase_info=RelinquishVirtualMemory(phase_info);
1298 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1299 fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception)
1326 source_info=AcquireVirtualMemory((size_t) fourier_info->height,
1327 fourier_info->width*sizeof(*source_pixels));
1328 if (source_info == (MemoryInfo *) NULL)
1330 (void) ThrowMagickException(exception,GetMagickModule(),
1331 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1332 return(MagickFalse);
1334 source_pixels=(double *) GetVirtualMemoryBlob(source_info);
1335 value=GetImageArtifact(image,"fourier:normalize");
1336 if (LocaleCompare(value,"inverse") == 0)
1342 Normalize inverse transform.
1345 gamma=PerceptibleReciprocal((double) fourier_info->width*
1346 fourier_info->height);
1347 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1348 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1350 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1351 fourier_pixels[i]*=gamma;
1353 fourier_pixels[i][0]*=gamma;
1354 fourier_pixels[i][1]*=gamma;
1359 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1360 #pragma omp critical (MagickCore_InverseFourierTransform)
1363 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1364 fourier_pixels,source_pixels,FFTW_ESTIMATE);
1365 fftw_execute(fftw_c2r_plan);
1366 fftw_destroy_plan(fftw_c2r_plan);
1369 image_view=AcquireAuthenticCacheView(image,exception);
1370 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1372 if (y >= (ssize_t) image->rows)
1374 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1375 image->columns ? image->columns : fourier_info->width,1UL,exception);
1376 if (q == (Quantum *) NULL)
1378 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1380 if (x < (ssize_t) image->columns)
1381 switch (fourier_info->channel)
1383 case RedPixelChannel:
1386 SetPixelRed(image,ClampToQuantum(QuantumRange*source_pixels[i]),q);
1389 case GreenPixelChannel:
1391 SetPixelGreen(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1395 case BluePixelChannel:
1397 SetPixelBlue(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1401 case BlackPixelChannel:
1403 SetPixelBlack(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1407 case AlphaPixelChannel:
1409 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1415 q+=GetPixelChannels(image);
1417 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1420 image_view=DestroyCacheView(image_view);
1421 source_info=RelinquishVirtualMemory(source_info);
1425 static MagickBooleanType InverseFourierTransformChannel(
1426 const Image *magnitude_image,const Image *phase_image,
1427 const PixelChannel channel,const MagickBooleanType modulus,
1428 Image *fourier_image,ExceptionInfo *exception)
1445 fourier_info.width=magnitude_image->columns;
1446 fourier_info.height=magnitude_image->rows;
1447 if ((magnitude_image->columns != magnitude_image->rows) ||
1448 ((magnitude_image->columns % 2) != 0) ||
1449 ((magnitude_image->rows % 2) != 0))
1451 extent=magnitude_image->columns < magnitude_image->rows ?
1452 magnitude_image->rows : magnitude_image->columns;
1453 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1455 fourier_info.height=fourier_info.width;
1456 fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
1457 fourier_info.channel=channel;
1458 fourier_info.modulus=modulus;
1459 inverse_info=AcquireVirtualMemory((size_t) fourier_info.height,
1460 fourier_info.center*sizeof(*inverse_pixels));
1461 if (inverse_info == (MemoryInfo *) NULL)
1463 (void) ThrowMagickException(exception,GetMagickModule(),
1464 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1465 magnitude_image->filename);
1466 return(MagickFalse);
1468 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1469 status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1470 inverse_pixels,exception);
1471 if (status != MagickFalse)
1472 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1474 inverse_info=RelinquishVirtualMemory(inverse_info);
1479 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1480 const Image *phase_image,const MagickBooleanType modulus,
1481 ExceptionInfo *exception)
1486 assert(magnitude_image != (Image *) NULL);
1487 assert(magnitude_image->signature == MagickSignature);
1488 if (magnitude_image->debug != MagickFalse)
1489 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1490 magnitude_image->filename);
1491 if (phase_image == (Image *) NULL)
1493 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1494 "TwoOrMoreImagesRequired","`%s'",magnitude_image->filename);
1495 return((Image *) NULL);
1497 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1498 fourier_image=(Image *) NULL;
1500 (void) ThrowMagickException(exception,GetMagickModule(),
1501 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
1502 magnitude_image->filename);
1505 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1506 magnitude_image->rows,MagickTrue,exception);
1507 if (fourier_image != (Image *) NULL)
1514 is_gray=IsImageGray(magnitude_image,exception);
1515 if (is_gray != MagickFalse)
1516 is_gray=IsImageGray(phase_image,exception);
1517 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1518 #pragma omp parallel sections
1521 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1528 if (is_gray != MagickFalse)
1529 thread_status=InverseFourierTransformChannel(magnitude_image,
1530 phase_image,GrayPixelChannel,modulus,fourier_image,exception);
1532 thread_status=InverseFourierTransformChannel(magnitude_image,
1533 phase_image,RedPixelChannel,modulus,fourier_image,exception);
1534 if (thread_status == MagickFalse)
1535 status=thread_status;
1537 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1544 thread_status=MagickTrue;
1545 if (is_gray == MagickFalse)
1546 thread_status=InverseFourierTransformChannel(magnitude_image,
1547 phase_image,GreenPixelChannel,modulus,fourier_image,exception);
1548 if (thread_status == MagickFalse)
1549 status=thread_status;
1551 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1558 thread_status=MagickTrue;
1559 if (is_gray == MagickFalse)
1560 thread_status=InverseFourierTransformChannel(magnitude_image,
1561 phase_image,BluePixelChannel,modulus,fourier_image,exception);
1562 if (thread_status == MagickFalse)
1563 status=thread_status;
1565 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1572 thread_status=MagickTrue;
1573 if (magnitude_image->colorspace == CMYKColorspace)
1574 thread_status=InverseFourierTransformChannel(magnitude_image,
1575 phase_image,BlackPixelChannel,modulus,fourier_image,exception);
1576 if (thread_status == MagickFalse)
1577 status=thread_status;
1579 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1586 thread_status=MagickTrue;
1587 if (magnitude_image->alpha_trait == BlendPixelTrait)
1588 thread_status=InverseFourierTransformChannel(magnitude_image,
1589 phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
1590 if (thread_status == MagickFalse)
1591 status=thread_status;
1594 if (status == MagickFalse)
1595 fourier_image=DestroyImage(fourier_image);
1600 return(fourier_image);