45 #include "magick/studio.h"
46 #include "magick/artifact.h"
47 #include "magick/attribute.h"
48 #include "magick/blob.h"
49 #include "magick/cache.h"
50 #include "magick/image.h"
51 #include "magick/image-private.h"
52 #include "magick/list.h"
53 #include "magick/fourier.h"
54 #include "magick/log.h"
55 #include "magick/memory_.h"
56 #include "magick/monitor.h"
57 #include "magick/monitor-private.h"
58 #include "magick/pixel-accessor.h"
59 #include "magick/pixel-private.h"
60 #include "magick/property.h"
61 #include "magick/quantum-private.h"
62 #include "magick/resource_.h"
63 #include "magick/string-private.h"
64 #include "magick/thread-private.h"
65 #if defined(MAGICKCORE_FFTW_DELEGATE)
67 #define ENABLE_FFTW_DELEGATE
68 #elif !defined(__cplusplus) && !defined(c_plusplus)
69 #define ENABLE_FFTW_DELEGATE
72 #if defined(ENABLE_FFTW_DELEGATE)
73 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
77 #if !defined(MAGICKCORE_HAVE_CABS)
78 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
80 #if !defined(MAGICKCORE_HAVE_CARG)
81 #define carg(z) (atan2(cimag(z),creal(z)))
83 #if !defined(MAGICKCORE_HAVE_CIMAG)
84 #define cimag(z) (z[1])
86 #if !defined(MAGICKCORE_HAVE_CREAL)
87 #define creal(z) (z[0])
137 MagickExport
Image *ComplexImages(
const Image *images,
const ComplexOperator op,
140 #define ComplexImageTag "Complex/Image"
181 assert(images != (
Image *) NULL);
182 assert(images->signature == MagickCoreSignature);
184 assert(exception->signature == MagickCoreSignature);
185 if (IsEventLogging() != MagickFalse)
186 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",images->filename);
187 if (images->next == (
Image *) NULL)
189 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
190 "ImageSequenceRequired",
"`%s'",images->filename);
191 return((
Image *) NULL);
193 image=CloneImage(images,0,0,MagickTrue,exception);
194 if (image == (
Image *) NULL)
195 return((
Image *) NULL);
196 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
198 image=DestroyImageList(image);
202 complex_images=NewImageList();
203 AppendImageToList(&complex_images,image);
204 image=CloneImage(images->next,0,0,MagickTrue,exception);
205 if (image == (
Image *) NULL)
207 complex_images=DestroyImageList(complex_images);
208 return(complex_images);
210 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
212 image=DestroyImageList(image);
216 AppendImageToList(&complex_images,image);
220 artifact=GetImageArtifact(image,
"complex:snr");
222 if (artifact != (
const char *) NULL)
223 snr=StringToDouble(artifact,(
char **) NULL);
225 Ai_image=images->next;
227 Bi_image=images->next;
228 if ((images->next->next != (
Image *) NULL) &&
229 (images->next->next->next != (
Image *) NULL))
231 Br_image=images->next->next;
232 Bi_image=images->next->next->next;
234 Cr_image=complex_images;
235 Ci_image=complex_images->next;
236 Ar_view=AcquireVirtualCacheView(Ar_image,exception);
237 Ai_view=AcquireVirtualCacheView(Ai_image,exception);
238 Br_view=AcquireVirtualCacheView(Br_image,exception);
239 Bi_view=AcquireVirtualCacheView(Bi_image,exception);
240 Cr_view=AcquireAuthenticCacheView(Cr_image,exception);
241 Ci_view=AcquireAuthenticCacheView(Ci_image,exception);
244 columns=MagickMin(Cr_image->columns,Ci_image->columns);
245 rows=MagickMin(Cr_image->rows,Ci_image->rows);
246 #if defined(MAGICKCORE_OPENMP_SUPPORT)
247 #pragma omp parallel for schedule(static) shared(progress,status) \
248 magick_number_threads(Cr_image,complex_images,rows,1L)
250 for (y=0; y < (ssize_t) rows; y++)
265 if (status == MagickFalse)
267 Ar=GetCacheViewVirtualPixels(Ar_view,0,y,columns,1,exception);
268 Ai=GetCacheViewVirtualPixels(Ai_view,0,y,columns,1,exception);
269 Br=GetCacheViewVirtualPixels(Br_view,0,y,columns,1,exception);
270 Bi=GetCacheViewVirtualPixels(Bi_view,0,y,columns,1,exception);
271 Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,columns,1,exception);
272 Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,columns,1,exception);
282 for (x=0; x < (ssize_t) columns; x++)
285 ai = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Ai)),
286 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Ai)),
287 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Ai)),
288 (Quantum) (image->matte != MagickFalse ? QuantumScale*
289 (MagickRealType) GetPixelOpacity(Ai) : (MagickRealType)
291 ar = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Ar)),
292 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Ar)),
293 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Ar)),
294 (Quantum) (image->matte != MagickFalse ? QuantumScale*
295 (MagickRealType) GetPixelOpacity(Ar) : (MagickRealType)
297 bi = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Bi)),
298 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Bi)),
299 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Bi)),
300 (Quantum) (image->matte != MagickFalse ? QuantumScale*
301 (MagickRealType) GetPixelOpacity(Bi) : (MagickRealType)
303 br = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Br)),
304 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Br)),
305 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Br)),
306 (Quantum) (image->matte != MagickFalse ? QuantumScale*
307 (MagickRealType) GetPixelOpacity(Br) : (MagickRealType)
314 case AddComplexOperator:
316 cr.red=ar.red+br.red;
317 ci.red=ai.red+bi.red;
318 cr.green=ar.green+br.green;
319 ci.green=ai.green+bi.green;
320 cr.blue=ar.blue+br.blue;
321 ci.blue=ai.blue+bi.blue;
322 cr.opacity=ar.opacity+br.opacity;
323 ci.opacity=ai.opacity+bi.opacity;
326 case ConjugateComplexOperator:
332 ci.green=(-ai.green);
335 cr.opacity=ar.opacity;
336 ci.opacity=(-ai.opacity);
339 case DivideComplexOperator:
344 gamma=PerceptibleReciprocal((
double) br.red*(
double) br.red+
345 (
double) bi.red*(
double) bi.red+snr);
346 cr.red=gamma*((double) ar.red*(
double) br.red+(double) ai.red*
348 ci.red=gamma*((double) ai.red*(
double) br.red-(double) ar.red*
350 gamma=PerceptibleReciprocal((
double) br.green*(
double) br.green+
351 (
double) bi.green*(
double) bi.green+snr);
352 cr.green=gamma*((double) ar.green*(
double) br.green+
353 (double) ai.green*(
double) bi.green);
354 ci.green=gamma*((double) ai.green*(
double) br.green-
355 (double) ar.green*(
double) bi.green);
356 gamma=PerceptibleReciprocal((
double) br.blue*(
double) br.blue+
357 (
double) bi.blue*(
double) bi.blue+snr);
358 cr.blue=gamma*((double) ar.blue*(
double) br.blue+
359 (double) ai.blue*(
double) bi.blue);
360 ci.blue=gamma*((double) ai.blue*(
double) br.blue-
361 (double) ar.blue*(
double) bi.blue);
362 gamma=PerceptibleReciprocal((
double) br.opacity*(
double) br.opacity+
363 (
double) bi.opacity*(
double) bi.opacity+snr);
364 cr.opacity=gamma*((double) ar.opacity*(
double) br.opacity+
365 (double) ai.opacity*(
double) bi.opacity);
366 ci.opacity=gamma*((double) ai.opacity*(
double) br.opacity-
367 (double) ar.opacity*(
double) bi.opacity);
370 case MagnitudePhaseComplexOperator:
372 cr.red=sqrt((
double) ar.red*(
double) ar.red+(
double) ai.red*
374 ci.red=atan2((
double) ai.red,(
double) ar.red)/(2.0*MagickPI)+0.5;
375 cr.green=sqrt((
double) ar.green*(
double) ar.green+(
double) ai.green*
377 ci.green=atan2((
double) ai.green,(
double) ar.green)/(2.0*MagickPI)+
379 cr.blue=sqrt((
double) ar.blue*(
double) ar.blue+(
double) ai.blue*
381 ci.blue=atan2((
double) ai.blue,(
double) ar.blue)/(2.0*MagickPI)+0.5;
382 cr.opacity=sqrt((
double) ar.opacity*(
double) ar.opacity+
383 (
double) ai.opacity*(
double) ai.opacity);
384 ci.opacity=atan2((
double) ai.opacity,(
double) ar.opacity)/
388 case MultiplyComplexOperator:
390 cr.red=((double) ar.red*(
double) br.red-(double) ai.red*
392 ci.red=((double) ai.red*(
double) br.red+(double) ar.red*
394 cr.green=((double) ar.green*(
double) br.green-(double) ai.green*
396 ci.green=((double) ai.green*(
double) br.green+(double) ar.green*
398 cr.blue=((double) ar.blue*(
double) br.blue-(double) ai.blue*
400 ci.blue=((double) ai.blue*(
double) br.blue+(double) ar.blue*
402 cr.opacity=((double) ar.opacity*(
double) br.opacity-
403 (double) ai.opacity*(
double) bi.opacity);
404 ci.opacity=((double) ai.opacity*(
double) br.opacity+
405 (double) ar.opacity*(
double) bi.opacity);
408 case RealImaginaryComplexOperator:
410 cr.red=(double) ar.red*cos(2.0*MagickPI*((
double) ai.red-0.5));
411 ci.red=(double) ar.red*sin(2.0*MagickPI*((
double) ai.red-0.5));
412 cr.green=(double) ar.green*cos(2.0*MagickPI*((
double) ai.green-0.5));
413 ci.green=(double) ar.green*sin(2.0*MagickPI*((
double) ai.green-0.5));
414 cr.blue=(double) ar.blue*cos(2.0*MagickPI*((
double) ai.blue-0.5));
415 ci.blue=(double) ar.blue*sin(2.0*MagickPI*((
double) ai.blue-0.5));
416 cr.opacity=(double) ar.opacity*cos(2.0*MagickPI*((
double) ai.opacity-
418 ci.opacity=(double) ar.opacity*sin(2.0*MagickPI*((
double) ai.opacity-
422 case SubtractComplexOperator:
424 cr.red=ar.red-br.red;
425 ci.red=ai.red-bi.red;
426 cr.green=ar.green-br.green;
427 ci.green=ai.green-bi.green;
428 cr.blue=ar.blue-br.blue;
429 ci.blue=ai.blue-bi.blue;
430 cr.opacity=ar.opacity-br.opacity;
431 ci.opacity=ai.opacity-bi.opacity;
435 Cr->red=(MagickRealType) QuantumRange*(MagickRealType) cr.red;
436 Ci->red=(MagickRealType) QuantumRange*(MagickRealType) ci.red;
437 Cr->green=(MagickRealType) QuantumRange*(MagickRealType) cr.green;
438 Ci->green=(MagickRealType) QuantumRange*(MagickRealType) ci.green;
439 Cr->blue=(MagickRealType) QuantumRange*(MagickRealType) cr.blue;
440 Ci->blue=(MagickRealType) QuantumRange*(MagickRealType) ci.blue;
441 if (images->matte != MagickFalse)
443 Cr->opacity=(MagickRealType) QuantumRange*(MagickRealType) cr.opacity;
444 Ci->opacity=(MagickRealType) QuantumRange*(MagickRealType) ci.opacity;
453 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
455 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
457 if (images->progress_monitor != (MagickProgressMonitor) NULL)
462 #if defined(MAGICKCORE_OPENMP_SUPPORT)
466 proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
467 if (proceed == MagickFalse)
471 Cr_view=DestroyCacheView(Cr_view);
472 Ci_view=DestroyCacheView(Ci_view);
473 Br_view=DestroyCacheView(Br_view);
474 Bi_view=DestroyCacheView(Bi_view);
475 Ar_view=DestroyCacheView(Ar_view);
476 Ai_view=DestroyCacheView(Ai_view);
477 if (status == MagickFalse)
478 complex_images=DestroyImageList(complex_images);
479 return(complex_images);
513 #if defined(ENABLE_FFTW_DELEGATE)
515 static MagickBooleanType RollFourier(
const size_t width,
const size_t height,
516 const ssize_t x_offset,
const ssize_t y_offset,
double *roll_pixels)
534 source_info=AcquireVirtualMemory(width,height*
sizeof(*source_pixels));
537 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
539 for (y=0L; y < (ssize_t) height; y++)
542 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
544 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
546 for (x=0L; x < (ssize_t) width; x++)
549 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
551 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
553 source_pixels[v*width+u]=roll_pixels[i++];
556 (void) memcpy(roll_pixels,source_pixels,height*width*
sizeof(*source_pixels));
557 source_info=RelinquishVirtualMemory(source_info);
561 static MagickBooleanType ForwardQuadrantSwap(
const size_t width,
562 const size_t height,
double *source_pixels,
double *forward_pixels)
575 center=(ssize_t) (width/2L)+1L;
576 status=RollFourier((
size_t) center,height,0L,(ssize_t) height/2L,
578 if (status == MagickFalse)
580 for (y=0L; y < (ssize_t) height; y++)
581 for (x=0L; x < (ssize_t) (width/2L); x++)
582 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
583 for (y=1; y < (ssize_t) height; y++)
584 for (x=0L; x < (ssize_t) (width/2L); x++)
585 forward_pixels[(height-y)*width+width/2L-x-1L]=
586 source_pixels[y*center+x+1L];
587 for (x=0L; x < (ssize_t) (width/2L); x++)
588 forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
592 static void CorrectPhaseLHS(
const size_t width,
const size_t height,
593 double *fourier_pixels)
599 for (y=0L; y < (ssize_t) height; y++)
600 for (x=0L; x < (ssize_t) (width/2L); x++)
601 fourier_pixels[y*width+x]*=(-1.0);
604 static MagickBooleanType ForwardFourier(
const FourierInfo *fourier_info,
637 magnitude_image=GetFirstImageInList(image);
638 phase_image=GetNextImageInList(image);
639 if (phase_image == (
Image *) NULL)
641 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
642 "ImageSequenceRequired",
"`%s'",image->filename);
648 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
649 sizeof(*magnitude_pixels));
650 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
651 sizeof(*phase_pixels));
652 if ((magnitude_info == (
MemoryInfo *) NULL) ||
656 phase_info=RelinquishVirtualMemory(phase_info);
658 magnitude_info=RelinquishVirtualMemory(magnitude_info);
659 (void) ThrowMagickException(exception,GetMagickModule(),
660 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
663 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
664 (void) memset(magnitude_pixels,0,fourier_info->width*
665 fourier_info->height*
sizeof(*magnitude_pixels));
666 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
667 (void) memset(phase_pixels,0,fourier_info->width*
668 fourier_info->height*
sizeof(*phase_pixels));
669 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
670 magnitude,magnitude_pixels);
671 if (status != MagickFalse)
672 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
674 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
675 if (fourier_info->modulus != MagickFalse)
678 for (y=0L; y < (ssize_t) fourier_info->height; y++)
679 for (x=0L; x < (ssize_t) fourier_info->width; x++)
681 phase_pixels[i]/=(2.0*MagickPI);
682 phase_pixels[i]+=0.5;
686 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
688 for (y=0L; y < (ssize_t) fourier_info->height; y++)
690 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
694 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
695 for (x=0L; x < (ssize_t) fourier_info->width; x++)
697 switch (fourier_info->channel)
702 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
707 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
712 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
717 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
722 SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType) QuantumRange*
723 magnitude_pixels[i]));
728 SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
735 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
736 if (status == MagickFalse)
739 magnitude_view=DestroyCacheView(magnitude_view);
741 phase_view=AcquireAuthenticCacheView(phase_image,exception);
742 for (y=0L; y < (ssize_t) fourier_info->height; y++)
744 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
748 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
749 for (x=0L; x < (ssize_t) fourier_info->width; x++)
751 switch (fourier_info->channel)
756 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
761 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
766 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
771 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
776 SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
781 SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
788 status=SyncCacheViewAuthenticPixels(phase_view,exception);
789 if (status == MagickFalse)
792 phase_view=DestroyCacheView(phase_view);
793 phase_info=RelinquishVirtualMemory(phase_info);
794 magnitude_info=RelinquishVirtualMemory(magnitude_info);
798 static MagickBooleanType ForwardFourierTransform(
FourierInfo *fourier_info,
799 const Image *image,
double *magnitude_pixels,
double *phase_pixels,
836 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
838 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
839 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
842 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
843 sizeof(*source_pixels));
846 (void) ThrowMagickException(exception,GetMagickModule(),
847 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
850 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
851 memset(source_pixels,0,fourier_info->width*fourier_info->height*
852 sizeof(*source_pixels));
854 image_view=AcquireVirtualCacheView(image,exception);
855 for (y=0L; y < (ssize_t) fourier_info->height; y++)
857 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
861 indexes=GetCacheViewVirtualIndexQueue(image_view);
862 for (x=0L; x < (ssize_t) fourier_info->width; x++)
864 switch (fourier_info->channel)
869 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
874 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
879 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
884 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
889 source_pixels[i]=QuantumScale*(MagickRealType)
890 GetPixelIndex(indexes+x);
895 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
903 image_view=DestroyCacheView(image_view);
904 forward_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/2+
905 1)*
sizeof(*forward_pixels));
908 (void) ThrowMagickException(exception,GetMagickModule(),
909 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
910 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
913 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
914 #if defined(MAGICKCORE_OPENMP_SUPPORT)
915 #pragma omp critical (MagickCore_ForwardFourierTransform)
917 fftw_r2c_plan=fftw_plan_dft_r2c_2d((
int) fourier_info->width,
918 (
int) fourier_info->height,source_pixels,forward_pixels,FFTW_ESTIMATE);
919 fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
920 fftw_destroy_plan(fftw_r2c_plan);
921 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
922 value=GetImageArtifact(image,
"fourier:normalize");
923 if ((value == (
const char *) NULL) || (LocaleCompare(value,
"forward") == 0))
932 gamma=PerceptibleReciprocal((
double) fourier_info->width*
933 fourier_info->height);
934 for (y=0L; y < (ssize_t) fourier_info->height; y++)
935 for (x=0L; x < (ssize_t) fourier_info->center; x++)
937 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
938 forward_pixels[i]*=gamma;
940 forward_pixels[i][0]*=gamma;
941 forward_pixels[i][1]*=gamma;
950 if (fourier_info->modulus != MagickFalse)
951 for (y=0L; y < (ssize_t) fourier_info->height; y++)
952 for (x=0L; x < (ssize_t) fourier_info->center; x++)
954 magnitude_pixels[i]=cabs(forward_pixels[i]);
955 phase_pixels[i]=carg(forward_pixels[i]);
959 for (y=0L; y < (ssize_t) fourier_info->height; y++)
960 for (x=0L; x < (ssize_t) fourier_info->center; x++)
962 magnitude_pixels[i]=creal(forward_pixels[i]);
963 phase_pixels[i]=cimag(forward_pixels[i]);
966 forward_info=(
MemoryInfo *) RelinquishVirtualMemory(forward_info);
970 static MagickBooleanType ForwardFourierTransformChannel(
const Image *image,
971 const ChannelType channel,
const MagickBooleanType modulus,
988 fourier_info.width=image->columns;
989 fourier_info.height=image->rows;
990 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
991 ((image->rows % 2) != 0))
993 size_t extent=image->columns < image->rows ? image->rows : image->columns;
994 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
996 fourier_info.height=fourier_info.width;
997 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
998 fourier_info.channel=channel;
999 fourier_info.modulus=modulus;
1000 magnitude_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1001 1)*
sizeof(*magnitude_pixels));
1002 phase_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+1)*
1003 sizeof(*phase_pixels));
1004 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1008 phase_info=RelinquishVirtualMemory(phase_info);
1010 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1011 (void) ThrowMagickException(exception,GetMagickModule(),
1012 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1013 return(MagickFalse);
1015 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1016 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1017 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
1018 phase_pixels,exception);
1019 if (status != MagickFalse)
1020 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
1021 phase_pixels,exception);
1022 phase_info=RelinquishVirtualMemory(phase_info);
1023 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1028 MagickExport
Image *ForwardFourierTransformImage(
const Image *image,
1034 fourier_image=NewImageList();
1035 #if !defined(ENABLE_FFTW_DELEGATE)
1037 (void) ThrowMagickException(exception,GetMagickModule(),
1038 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1049 width=image->columns;
1051 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
1052 ((image->rows % 2) != 0))
1054 size_t extent=image->columns < image->rows ? image->rows :
1056 width=(extent & 0x01) == 1 ? extent+1UL : extent;
1059 magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
1060 if (magnitude_image != (
Image *) NULL)
1065 magnitude_image->storage_class=DirectClass;
1066 magnitude_image->depth=32UL;
1067 phase_image=CloneImage(image,width,height,MagickTrue,exception);
1068 if (phase_image == (
Image *) NULL)
1069 magnitude_image=DestroyImage(magnitude_image);
1076 phase_image->storage_class=DirectClass;
1077 phase_image->depth=32UL;
1078 AppendImageToList(&fourier_image,magnitude_image);
1079 AppendImageToList(&fourier_image,phase_image);
1081 is_gray=IsGrayImage(image,exception);
1082 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1083 #pragma omp parallel sections
1086 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1093 if (is_gray != MagickFalse)
1094 thread_status=ForwardFourierTransformChannel(image,
1095 GrayChannels,modulus,fourier_image,exception);
1097 thread_status=ForwardFourierTransformChannel(image,RedChannel,
1098 modulus,fourier_image,exception);
1099 if (thread_status == MagickFalse)
1100 status=thread_status;
1102 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1109 thread_status=MagickTrue;
1110 if (is_gray == MagickFalse)
1111 thread_status=ForwardFourierTransformChannel(image,
1112 GreenChannel,modulus,fourier_image,exception);
1113 if (thread_status == MagickFalse)
1114 status=thread_status;
1116 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1123 thread_status=MagickTrue;
1124 if (is_gray == MagickFalse)
1125 thread_status=ForwardFourierTransformChannel(image,
1126 BlueChannel,modulus,fourier_image,exception);
1127 if (thread_status == MagickFalse)
1128 status=thread_status;
1130 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1137 thread_status=MagickTrue;
1138 if (image->matte != MagickFalse)
1139 thread_status=ForwardFourierTransformChannel(image,
1140 OpacityChannel,modulus,fourier_image,exception);
1141 if (thread_status == MagickFalse)
1142 status=thread_status;
1144 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1151 thread_status=MagickTrue;
1152 if (image->colorspace == CMYKColorspace)
1153 thread_status=ForwardFourierTransformChannel(image,
1154 IndexChannel,modulus,fourier_image,exception);
1155 if (thread_status == MagickFalse)
1156 status=thread_status;
1159 if (status == MagickFalse)
1160 fourier_image=DestroyImageList(fourier_image);
1166 return(fourier_image);
1203 #if defined(ENABLE_FFTW_DELEGATE)
1204 static MagickBooleanType InverseQuadrantSwap(
const size_t width,
1205 const size_t height,
const double *source,
double *destination)
1215 center=(ssize_t) (width/2L)+1L;
1216 for (y=1L; y < (ssize_t) height; y++)
1217 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1218 destination[(height-y)*center-x+width/2L]=source[y*width+x];
1219 for (y=0L; y < (ssize_t) height; y++)
1220 destination[y*center]=source[y*width+width/2L];
1221 for (x=0L; x < center; x++)
1222 destination[x]=source[center-x-1L];
1223 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1226 static MagickBooleanType InverseFourier(
FourierInfo *fourier_info,
1227 const Image *magnitude_image,
const Image *phase_image,
1261 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1262 sizeof(*magnitude_pixels));
1263 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1264 sizeof(*phase_pixels));
1265 inverse_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/
1266 2+1)*
sizeof(*inverse_pixels));
1267 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1272 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1274 phase_info=RelinquishVirtualMemory(phase_info);
1276 inverse_info=RelinquishVirtualMemory(inverse_info);
1277 (void) ThrowMagickException(exception,GetMagickModule(),
1278 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1279 magnitude_image->filename);
1280 return(MagickFalse);
1282 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1283 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1284 inverse_pixels=(
double *) GetVirtualMemoryBlob(inverse_info);
1286 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1287 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1289 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1293 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
1294 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1296 switch (fourier_info->channel)
1301 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
1306 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
1311 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
1314 case OpacityChannel:
1316 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
1321 magnitude_pixels[i]=QuantumScale*(MagickRealType)
1322 GetPixelIndex(indexes+x);
1327 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
1335 magnitude_view=DestroyCacheView(magnitude_view);
1336 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1337 magnitude_pixels,inverse_pixels);
1338 (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1339 fourier_info->center*
sizeof(*magnitude_pixels));
1341 phase_view=AcquireVirtualCacheView(phase_image,exception);
1342 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1344 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1348 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
1349 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1351 switch (fourier_info->channel)
1356 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
1361 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
1366 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
1369 case OpacityChannel:
1371 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
1376 phase_pixels[i]=QuantumScale*(MagickRealType)
1377 GetPixelIndex(indexes+x);
1382 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
1390 if (fourier_info->modulus != MagickFalse)
1393 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1394 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1396 phase_pixels[i]-=0.5;
1397 phase_pixels[i]*=(2.0*MagickPI);
1401 phase_view=DestroyCacheView(phase_view);
1402 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1403 if (status != MagickFalse)
1404 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1405 phase_pixels,inverse_pixels);
1406 (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1407 fourier_info->center*
sizeof(*phase_pixels));
1408 inverse_info=RelinquishVirtualMemory(inverse_info);
1413 if (fourier_info->modulus != MagickFalse)
1414 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1415 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1417 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1418 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1419 magnitude_pixels[i]*sin(phase_pixels[i]);
1421 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1422 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1427 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1428 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1430 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1431 fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1433 fourier_pixels[i][0]=magnitude_pixels[i];
1434 fourier_pixels[i][1]=phase_pixels[i];
1438 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1439 phase_info=RelinquishVirtualMemory(phase_info);
1443 static MagickBooleanType InverseFourierTransform(
FourierInfo *fourier_info,
1475 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
1477 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1478 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
1479 return(MagickFalse);
1481 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1482 sizeof(*source_pixels));
1485 (void) ThrowMagickException(exception,GetMagickModule(),
1486 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1487 return(MagickFalse);
1489 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
1490 value=GetImageArtifact(image,
"fourier:normalize");
1491 if (LocaleCompare(value,
"inverse") == 0)
1500 gamma=PerceptibleReciprocal((
double) fourier_info->width*
1501 fourier_info->height);
1502 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1503 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1505 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1506 fourier_pixels[i]*=gamma;
1508 fourier_pixels[i][0]*=gamma;
1509 fourier_pixels[i][1]*=gamma;
1514 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1515 #pragma omp critical (MagickCore_InverseFourierTransform)
1517 fftw_c2r_plan=fftw_plan_dft_c2r_2d((
int) fourier_info->width,
1518 (
int) fourier_info->height,fourier_pixels,source_pixels,FFTW_ESTIMATE);
1519 fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1520 fftw_destroy_plan(fftw_c2r_plan);
1522 image_view=AcquireAuthenticCacheView(image,exception);
1523 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1525 if (y >= (ssize_t) image->rows)
1527 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1528 image->columns ? image->columns : fourier_info->width,1UL,exception);
1531 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1532 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1534 if (x < (ssize_t) image->columns)
1535 switch (fourier_info->channel)
1540 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
1546 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
1552 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
1556 case OpacityChannel:
1558 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*
1564 SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType)
1565 QuantumRange*source_pixels[i]));
1570 SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*
1578 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1581 image_view=DestroyCacheView(image_view);
1582 source_info=RelinquishVirtualMemory(source_info);
1586 static MagickBooleanType InverseFourierTransformChannel(
1587 const Image *magnitude_image,
const Image *phase_image,
1588 const ChannelType channel,
const MagickBooleanType modulus,
1603 fourier_info.width=magnitude_image->columns;
1604 fourier_info.height=magnitude_image->rows;
1605 if ((magnitude_image->columns != magnitude_image->rows) ||
1606 ((magnitude_image->columns % 2) != 0) ||
1607 ((magnitude_image->rows % 2) != 0))
1609 size_t extent=magnitude_image->columns < magnitude_image->rows ?
1610 magnitude_image->rows : magnitude_image->columns;
1611 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1613 fourier_info.height=fourier_info.width;
1614 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1615 fourier_info.channel=channel;
1616 fourier_info.modulus=modulus;
1617 inverse_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1618 1)*
sizeof(*inverse_pixels));
1621 (void) ThrowMagickException(exception,GetMagickModule(),
1622 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1623 magnitude_image->filename);
1624 return(MagickFalse);
1626 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1627 status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1628 inverse_pixels,exception);
1629 if (status != MagickFalse)
1630 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1632 inverse_info=RelinquishVirtualMemory(inverse_info);
1637 MagickExport
Image *InverseFourierTransformImage(
const Image *magnitude_image,
1638 const Image *phase_image,
const MagickBooleanType modulus,
1644 assert(magnitude_image != (
Image *) NULL);
1645 assert(magnitude_image->signature == MagickCoreSignature);
1646 if (IsEventLogging() != MagickFalse)
1647 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1648 magnitude_image->filename);
1649 if (phase_image == (
Image *) NULL)
1651 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1652 "ImageSequenceRequired",
"`%s'",magnitude_image->filename);
1653 return((
Image *) NULL);
1655 #if !defined(ENABLE_FFTW_DELEGATE)
1656 fourier_image=(
Image *) NULL;
1658 (void) ThrowMagickException(exception,GetMagickModule(),
1659 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1660 magnitude_image->filename);
1663 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1664 magnitude_image->rows,MagickTrue,exception);
1665 if (fourier_image != (
Image *) NULL)
1672 is_gray=IsGrayImage(magnitude_image,exception);
1673 if (is_gray != MagickFalse)
1674 is_gray=IsGrayImage(phase_image,exception);
1675 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1676 #pragma omp parallel sections
1679 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1686 if (is_gray != MagickFalse)
1687 thread_status=InverseFourierTransformChannel(magnitude_image,
1688 phase_image,GrayChannels,modulus,fourier_image,exception);
1690 thread_status=InverseFourierTransformChannel(magnitude_image,
1691 phase_image,RedChannel,modulus,fourier_image,exception);
1692 if (thread_status == MagickFalse)
1693 status=thread_status;
1695 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1702 thread_status=MagickTrue;
1703 if (is_gray == MagickFalse)
1704 thread_status=InverseFourierTransformChannel(magnitude_image,
1705 phase_image,GreenChannel,modulus,fourier_image,exception);
1706 if (thread_status == MagickFalse)
1707 status=thread_status;
1709 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1716 thread_status=MagickTrue;
1717 if (is_gray == MagickFalse)
1718 thread_status=InverseFourierTransformChannel(magnitude_image,
1719 phase_image,BlueChannel,modulus,fourier_image,exception);
1720 if (thread_status == MagickFalse)
1721 status=thread_status;
1723 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1730 thread_status=MagickTrue;
1731 if (magnitude_image->matte != MagickFalse)
1732 thread_status=InverseFourierTransformChannel(magnitude_image,
1733 phase_image,OpacityChannel,modulus,fourier_image,exception);
1734 if (thread_status == MagickFalse)
1735 status=thread_status;
1737 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1744 thread_status=MagickTrue;
1745 if (magnitude_image->colorspace == CMYKColorspace)
1746 thread_status=InverseFourierTransformChannel(magnitude_image,
1747 phase_image,IndexChannel,modulus,fourier_image,exception);
1748 if (thread_status == MagickFalse)
1749 status=thread_status;
1752 if (status == MagickFalse)
1753 fourier_image=DestroyImage(fourier_image);
1758 return(fourier_image);