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) && !defined(__cplusplus) && !defined(c_plusplus)
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])
130 MagickExport
Image *ComplexImages(
const Image *images,
const ComplexOperator op,
133 #define ComplexImageTag "Complex/Image"
174 assert(images != (
Image *) NULL);
175 assert(images->signature == MagickCoreSignature);
177 assert(exception->signature == MagickCoreSignature);
178 if (IsEventLogging() != MagickFalse)
179 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",images->filename);
180 if (images->next == (
Image *) NULL)
182 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
183 "ImageSequenceRequired",
"`%s'",images->filename);
184 return((
Image *) NULL);
186 image=CloneImage(images,0,0,MagickTrue,exception);
187 if (image == (
Image *) NULL)
188 return((
Image *) NULL);
189 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
191 image=DestroyImageList(image);
195 complex_images=NewImageList();
196 AppendImageToList(&complex_images,image);
197 image=CloneImage(images->next,0,0,MagickTrue,exception);
198 if (image == (
Image *) NULL)
200 complex_images=DestroyImageList(complex_images);
201 return(complex_images);
203 if (SetImageStorageClass(image,DirectClass) == MagickFalse)
205 image=DestroyImageList(image);
209 AppendImageToList(&complex_images,image);
213 artifact=GetImageArtifact(image,
"complex:snr");
215 if (artifact != (
const char *) NULL)
216 snr=StringToDouble(artifact,(
char **) NULL);
218 Ai_image=images->next;
220 Bi_image=images->next;
221 if ((images->next->next != (
Image *) NULL) &&
222 (images->next->next->next != (
Image *) NULL))
224 Br_image=images->next->next;
225 Bi_image=images->next->next->next;
227 Cr_image=complex_images;
228 Ci_image=complex_images->next;
229 Ar_view=AcquireVirtualCacheView(Ar_image,exception);
230 Ai_view=AcquireVirtualCacheView(Ai_image,exception);
231 Br_view=AcquireVirtualCacheView(Br_image,exception);
232 Bi_view=AcquireVirtualCacheView(Bi_image,exception);
233 Cr_view=AcquireAuthenticCacheView(Cr_image,exception);
234 Ci_view=AcquireAuthenticCacheView(Ci_image,exception);
237 columns=MagickMin(Cr_image->columns,Ci_image->columns);
238 rows=MagickMin(Cr_image->rows,Ci_image->rows);
239 #if defined(MAGICKCORE_OPENMP_SUPPORT)
240 #pragma omp parallel for schedule(static) shared(progress,status) \
241 magick_number_threads(Cr_image,complex_images,rows,1L)
243 for (y=0; y < (ssize_t) rows; y++)
258 if (status == MagickFalse)
260 Ar=GetCacheViewVirtualPixels(Ar_view,0,y,columns,1,exception);
261 Ai=GetCacheViewVirtualPixels(Ai_view,0,y,columns,1,exception);
262 Br=GetCacheViewVirtualPixels(Br_view,0,y,columns,1,exception);
263 Bi=GetCacheViewVirtualPixels(Bi_view,0,y,columns,1,exception);
264 Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,columns,1,exception);
265 Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,columns,1,exception);
275 for (x=0; x < (ssize_t) columns; x++)
278 ai = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Ai)),
279 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Ai)),
280 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Ai)),
281 (Quantum) (image->matte != MagickFalse ? QuantumScale*
282 (MagickRealType) GetPixelOpacity(Ai) : (MagickRealType)
284 ar = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Ar)),
285 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Ar)),
286 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Ar)),
287 (Quantum) (image->matte != MagickFalse ? QuantumScale*
288 (MagickRealType) GetPixelOpacity(Ar) : (MagickRealType)
290 bi = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Bi)),
291 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Bi)),
292 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Bi)),
293 (Quantum) (image->matte != MagickFalse ? QuantumScale*
294 (MagickRealType) GetPixelOpacity(Bi) : (MagickRealType)
296 br = { (Quantum) (QuantumScale*(MagickRealType) GetPixelRed(Br)),
297 (Quantum) (QuantumScale*(MagickRealType) GetPixelGreen(Br)),
298 (Quantum) (QuantumScale*(MagickRealType) GetPixelBlue(Br)),
299 (Quantum) (image->matte != MagickFalse ? QuantumScale*
300 (MagickRealType) GetPixelOpacity(Br) : (MagickRealType)
307 case AddComplexOperator:
309 cr.red=ar.red+br.red;
310 ci.red=ai.red+bi.red;
311 cr.green=ar.green+br.green;
312 ci.green=ai.green+bi.green;
313 cr.blue=ar.blue+br.blue;
314 ci.blue=ai.blue+bi.blue;
315 cr.opacity=ar.opacity+br.opacity;
316 ci.opacity=ai.opacity+bi.opacity;
319 case ConjugateComplexOperator:
325 ci.green=(-ai.green);
328 cr.opacity=ar.opacity;
329 ci.opacity=(-ai.opacity);
332 case DivideComplexOperator:
337 gamma=PerceptibleReciprocal((
double) br.red*(
double) br.red+
338 (
double) bi.red*(
double) bi.red+snr);
339 cr.red=gamma*((double) ar.red*(
double) br.red+(double) ai.red*
341 ci.red=gamma*((double) ai.red*(
double) br.red-(double) ar.red*
343 gamma=PerceptibleReciprocal((
double) br.green*(
double) br.green+
344 (
double) bi.green*(
double) bi.green+snr);
345 cr.green=gamma*((double) ar.green*(
double) br.green+
346 (double) ai.green*(
double) bi.green);
347 ci.green=gamma*((double) ai.green*(
double) br.green-
348 (double) ar.green*(
double) bi.green);
349 gamma=PerceptibleReciprocal((
double) br.blue*(
double) br.blue+
350 (
double) bi.blue*(
double) bi.blue+snr);
351 cr.blue=gamma*((double) ar.blue*(
double) br.blue+
352 (double) ai.blue*(
double) bi.blue);
353 ci.blue=gamma*((double) ai.blue*(
double) br.blue-
354 (double) ar.blue*(
double) bi.blue);
355 gamma=PerceptibleReciprocal((
double) br.opacity*(
double) br.opacity+
356 (
double) bi.opacity*(
double) bi.opacity+snr);
357 cr.opacity=gamma*((double) ar.opacity*(
double) br.opacity+
358 (double) ai.opacity*(
double) bi.opacity);
359 ci.opacity=gamma*((double) ai.opacity*(
double) br.opacity-
360 (double) ar.opacity*(
double) bi.opacity);
363 case MagnitudePhaseComplexOperator:
365 cr.red=sqrt((
double) ar.red*(
double) ar.red+(
double) ai.red*
367 ci.red=atan2((
double) ai.red,(
double) ar.red)/(2.0*MagickPI)+0.5;
368 cr.green=sqrt((
double) ar.green*(
double) ar.green+(
double) ai.green*
370 ci.green=atan2((
double) ai.green,(
double) ar.green)/(2.0*MagickPI)+
372 cr.blue=sqrt((
double) ar.blue*(
double) ar.blue+(
double) ai.blue*
374 ci.blue=atan2((
double) ai.blue,(
double) ar.blue)/(2.0*MagickPI)+0.5;
375 cr.opacity=sqrt((
double) ar.opacity*(
double) ar.opacity+
376 (
double) ai.opacity*(
double) ai.opacity);
377 ci.opacity=atan2((
double) ai.opacity,(
double) ar.opacity)/
381 case MultiplyComplexOperator:
383 cr.red=((double) ar.red*(
double) br.red-(double) ai.red*
385 ci.red=((double) ai.red*(
double) br.red+(double) ar.red*
387 cr.green=((double) ar.green*(
double) br.green-(double) ai.green*
389 ci.green=((double) ai.green*(
double) br.green+(double) ar.green*
391 cr.blue=((double) ar.blue*(
double) br.blue-(double) ai.blue*
393 ci.blue=((double) ai.blue*(
double) br.blue+(double) ar.blue*
395 cr.opacity=((double) ar.opacity*(
double) br.opacity-
396 (double) ai.opacity*(
double) bi.opacity);
397 ci.opacity=((double) ai.opacity*(
double) br.opacity+
398 (double) ar.opacity*(
double) bi.opacity);
401 case RealImaginaryComplexOperator:
403 cr.red=(double) ar.red*cos(2.0*MagickPI*((
double) ai.red-0.5));
404 ci.red=(double) ar.red*sin(2.0*MagickPI*((
double) ai.red-0.5));
405 cr.green=(double) ar.green*cos(2.0*MagickPI*((
double) ai.green-0.5));
406 ci.green=(double) ar.green*sin(2.0*MagickPI*((
double) ai.green-0.5));
407 cr.blue=(double) ar.blue*cos(2.0*MagickPI*((
double) ai.blue-0.5));
408 ci.blue=(double) ar.blue*sin(2.0*MagickPI*((
double) ai.blue-0.5));
409 cr.opacity=(double) ar.opacity*cos(2.0*MagickPI*((
double) ai.opacity-
411 ci.opacity=(double) ar.opacity*sin(2.0*MagickPI*((
double) ai.opacity-
415 case SubtractComplexOperator:
417 cr.red=ar.red-br.red;
418 ci.red=ai.red-bi.red;
419 cr.green=ar.green-br.green;
420 ci.green=ai.green-bi.green;
421 cr.blue=ar.blue-br.blue;
422 ci.blue=ai.blue-bi.blue;
423 cr.opacity=ar.opacity-br.opacity;
424 ci.opacity=ai.opacity-bi.opacity;
428 Cr->red=(MagickRealType) QuantumRange*(MagickRealType) cr.red;
429 Ci->red=(MagickRealType) QuantumRange*(MagickRealType) ci.red;
430 Cr->green=(MagickRealType) QuantumRange*(MagickRealType) cr.green;
431 Ci->green=(MagickRealType) QuantumRange*(MagickRealType) ci.green;
432 Cr->blue=(MagickRealType) QuantumRange*(MagickRealType) cr.blue;
433 Ci->blue=(MagickRealType) QuantumRange*(MagickRealType) ci.blue;
434 if (images->matte != MagickFalse)
436 Cr->opacity=(MagickRealType) QuantumRange*(MagickRealType) cr.opacity;
437 Ci->opacity=(MagickRealType) QuantumRange*(MagickRealType) ci.opacity;
446 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
448 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
450 if (images->progress_monitor != (MagickProgressMonitor) NULL)
455 #if defined(MAGICKCORE_OPENMP_SUPPORT)
459 proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
460 if (proceed == MagickFalse)
464 Cr_view=DestroyCacheView(Cr_view);
465 Ci_view=DestroyCacheView(Ci_view);
466 Br_view=DestroyCacheView(Br_view);
467 Bi_view=DestroyCacheView(Bi_view);
468 Ar_view=DestroyCacheView(Ar_view);
469 Ai_view=DestroyCacheView(Ai_view);
470 if (status == MagickFalse)
471 complex_images=DestroyImageList(complex_images);
472 return(complex_images);
506 #if defined(MAGICKCORE_FFTW_DELEGATE) && !defined(__cplusplus) && !defined(c_plusplus)
508 static MagickBooleanType RollFourier(
const size_t width,
const size_t height,
509 const ssize_t x_offset,
const ssize_t y_offset,
double *roll_pixels)
527 source_info=AcquireVirtualMemory(width,height*
sizeof(*source_pixels));
530 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
532 for (y=0L; y < (ssize_t) height; y++)
535 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
537 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
539 for (x=0L; x < (ssize_t) width; x++)
542 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
544 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
546 source_pixels[v*width+u]=roll_pixels[i++];
549 (void) memcpy(roll_pixels,source_pixels,height*width*
sizeof(*source_pixels));
550 source_info=RelinquishVirtualMemory(source_info);
554 static MagickBooleanType ForwardQuadrantSwap(
const size_t width,
555 const size_t height,
double *source_pixels,
double *forward_pixels)
568 center=(ssize_t) (width/2L)+1L;
569 status=RollFourier((
size_t) center,height,0L,(ssize_t) height/2L,
571 if (status == MagickFalse)
573 for (y=0L; y < (ssize_t) height; y++)
574 for (x=0L; x < (ssize_t) (width/2L); x++)
575 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
576 for (y=1; y < (ssize_t) height; y++)
577 for (x=0L; x < (ssize_t) (width/2L); x++)
578 forward_pixels[(height-y)*width+width/2L-x-1L]=
579 source_pixels[y*center+x+1L];
580 for (x=0L; x < (ssize_t) (width/2L); x++)
581 forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
585 static void CorrectPhaseLHS(
const size_t width,
const size_t height,
586 double *fourier_pixels)
592 for (y=0L; y < (ssize_t) height; y++)
593 for (x=0L; x < (ssize_t) (width/2L); x++)
594 fourier_pixels[y*width+x]*=(-1.0);
597 static MagickBooleanType ForwardFourier(
const FourierInfo *fourier_info,
630 magnitude_image=GetFirstImageInList(image);
631 phase_image=GetNextImageInList(image);
632 if (phase_image == (
Image *) NULL)
634 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
635 "ImageSequenceRequired",
"`%s'",image->filename);
641 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
642 sizeof(*magnitude_pixels));
643 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
644 sizeof(*phase_pixels));
645 if ((magnitude_info == (
MemoryInfo *) NULL) ||
649 phase_info=RelinquishVirtualMemory(phase_info);
651 magnitude_info=RelinquishVirtualMemory(magnitude_info);
652 (void) ThrowMagickException(exception,GetMagickModule(),
653 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
656 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
657 (void) memset(magnitude_pixels,0,fourier_info->width*
658 fourier_info->height*
sizeof(*magnitude_pixels));
659 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
660 (void) memset(phase_pixels,0,fourier_info->width*
661 fourier_info->height*
sizeof(*phase_pixels));
662 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
663 magnitude,magnitude_pixels);
664 if (status != MagickFalse)
665 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
667 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
668 if (fourier_info->modulus != MagickFalse)
671 for (y=0L; y < (ssize_t) fourier_info->height; y++)
672 for (x=0L; x < (ssize_t) fourier_info->width; x++)
674 phase_pixels[i]/=(2.0*MagickPI);
675 phase_pixels[i]+=0.5;
679 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
681 for (y=0L; y < (ssize_t) fourier_info->height; y++)
683 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
687 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
688 for (x=0L; x < (ssize_t) fourier_info->width; x++)
690 switch (fourier_info->channel)
695 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
700 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
705 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
710 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
715 SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType) QuantumRange*
716 magnitude_pixels[i]));
721 SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*magnitude_pixels[i]));
728 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
729 if (status == MagickFalse)
732 magnitude_view=DestroyCacheView(magnitude_view);
734 phase_view=AcquireAuthenticCacheView(phase_image,exception);
735 for (y=0L; y < (ssize_t) fourier_info->height; y++)
737 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
741 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
742 for (x=0L; x < (ssize_t) fourier_info->width; x++)
744 switch (fourier_info->channel)
749 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
754 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
759 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
764 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
769 SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
774 SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*phase_pixels[i]));
781 status=SyncCacheViewAuthenticPixels(phase_view,exception);
782 if (status == MagickFalse)
785 phase_view=DestroyCacheView(phase_view);
786 phase_info=RelinquishVirtualMemory(phase_info);
787 magnitude_info=RelinquishVirtualMemory(magnitude_info);
791 static MagickBooleanType ForwardFourierTransform(
FourierInfo *fourier_info,
792 const Image *image,
double *magnitude_pixels,
double *phase_pixels,
829 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
831 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
832 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
835 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
836 sizeof(*source_pixels));
839 (void) ThrowMagickException(exception,GetMagickModule(),
840 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
843 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
844 memset(source_pixels,0,fourier_info->width*fourier_info->height*
845 sizeof(*source_pixels));
847 image_view=AcquireVirtualCacheView(image,exception);
848 for (y=0L; y < (ssize_t) fourier_info->height; y++)
850 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
854 indexes=GetCacheViewVirtualIndexQueue(image_view);
855 for (x=0L; x < (ssize_t) fourier_info->width; x++)
857 switch (fourier_info->channel)
862 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
867 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
872 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
877 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
882 source_pixels[i]=QuantumScale*(MagickRealType)
883 GetPixelIndex(indexes+x);
888 source_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
896 image_view=DestroyCacheView(image_view);
897 forward_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/2+
898 1)*
sizeof(*forward_pixels));
901 (void) ThrowMagickException(exception,GetMagickModule(),
902 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
903 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
906 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
907 #if defined(MAGICKCORE_OPENMP_SUPPORT)
908 #pragma omp critical (MagickCore_ForwardFourierTransform)
910 fftw_r2c_plan=fftw_plan_dft_r2c_2d((
int) fourier_info->width,
911 (
int) fourier_info->height,source_pixels,forward_pixels,FFTW_ESTIMATE);
912 fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
913 fftw_destroy_plan(fftw_r2c_plan);
914 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
915 value=GetImageArtifact(image,
"fourier:normalize");
916 if ((value == (
const char *) NULL) || (LocaleCompare(value,
"forward") == 0))
925 gamma=PerceptibleReciprocal((
double) fourier_info->width*
926 fourier_info->height);
927 for (y=0L; y < (ssize_t) fourier_info->height; y++)
928 for (x=0L; x < (ssize_t) fourier_info->center; x++)
930 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
931 forward_pixels[i]*=gamma;
933 forward_pixels[i][0]*=gamma;
934 forward_pixels[i][1]*=gamma;
943 if (fourier_info->modulus != MagickFalse)
944 for (y=0L; y < (ssize_t) fourier_info->height; y++)
945 for (x=0L; x < (ssize_t) fourier_info->center; x++)
947 magnitude_pixels[i]=cabs(forward_pixels[i]);
948 phase_pixels[i]=carg(forward_pixels[i]);
952 for (y=0L; y < (ssize_t) fourier_info->height; y++)
953 for (x=0L; x < (ssize_t) fourier_info->center; x++)
955 magnitude_pixels[i]=creal(forward_pixels[i]);
956 phase_pixels[i]=cimag(forward_pixels[i]);
959 forward_info=(
MemoryInfo *) RelinquishVirtualMemory(forward_info);
963 static MagickBooleanType ForwardFourierTransformChannel(
const Image *image,
964 const ChannelType channel,
const MagickBooleanType modulus,
981 fourier_info.width=image->columns;
982 fourier_info.height=image->rows;
983 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
984 ((image->rows % 2) != 0))
986 size_t extent=image->columns < image->rows ? image->rows : image->columns;
987 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
989 fourier_info.height=fourier_info.width;
990 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
991 fourier_info.channel=channel;
992 fourier_info.modulus=modulus;
993 magnitude_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
994 1)*
sizeof(*magnitude_pixels));
995 phase_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+1)*
996 sizeof(*phase_pixels));
997 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1001 phase_info=RelinquishVirtualMemory(phase_info);
1003 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1004 (void) ThrowMagickException(exception,GetMagickModule(),
1005 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1006 return(MagickFalse);
1008 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1009 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1010 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
1011 phase_pixels,exception);
1012 if (status != MagickFalse)
1013 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
1014 phase_pixels,exception);
1015 phase_info=RelinquishVirtualMemory(phase_info);
1016 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1021 MagickExport
Image *ForwardFourierTransformImage(
const Image *image,
1027 fourier_image=NewImageList();
1028 #if !defined(MAGICKCORE_FFTW_DELEGATE) || defined(__cplusplus) || defined(c_plusplus)
1030 (void) ThrowMagickException(exception,GetMagickModule(),
1031 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1042 width=image->columns;
1044 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
1045 ((image->rows % 2) != 0))
1047 size_t extent=image->columns < image->rows ? image->rows :
1049 width=(extent & 0x01) == 1 ? extent+1UL : extent;
1052 magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
1053 if (magnitude_image != (
Image *) NULL)
1058 magnitude_image->storage_class=DirectClass;
1059 magnitude_image->depth=32UL;
1060 phase_image=CloneImage(image,width,height,MagickTrue,exception);
1061 if (phase_image == (
Image *) NULL)
1062 magnitude_image=DestroyImage(magnitude_image);
1069 phase_image->storage_class=DirectClass;
1070 phase_image->depth=32UL;
1071 AppendImageToList(&fourier_image,magnitude_image);
1072 AppendImageToList(&fourier_image,phase_image);
1074 is_gray=IsGrayImage(image,exception);
1075 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1076 #pragma omp parallel sections
1079 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1086 if (is_gray != MagickFalse)
1087 thread_status=ForwardFourierTransformChannel(image,
1088 GrayChannels,modulus,fourier_image,exception);
1090 thread_status=ForwardFourierTransformChannel(image,RedChannel,
1091 modulus,fourier_image,exception);
1092 if (thread_status == MagickFalse)
1093 status=thread_status;
1095 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1102 thread_status=MagickTrue;
1103 if (is_gray == MagickFalse)
1104 thread_status=ForwardFourierTransformChannel(image,
1105 GreenChannel,modulus,fourier_image,exception);
1106 if (thread_status == MagickFalse)
1107 status=thread_status;
1109 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1116 thread_status=MagickTrue;
1117 if (is_gray == MagickFalse)
1118 thread_status=ForwardFourierTransformChannel(image,
1119 BlueChannel,modulus,fourier_image,exception);
1120 if (thread_status == MagickFalse)
1121 status=thread_status;
1123 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1130 thread_status=MagickTrue;
1131 if (image->matte != MagickFalse)
1132 thread_status=ForwardFourierTransformChannel(image,
1133 OpacityChannel,modulus,fourier_image,exception);
1134 if (thread_status == MagickFalse)
1135 status=thread_status;
1137 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1144 thread_status=MagickTrue;
1145 if (image->colorspace == CMYKColorspace)
1146 thread_status=ForwardFourierTransformChannel(image,
1147 IndexChannel,modulus,fourier_image,exception);
1148 if (thread_status == MagickFalse)
1149 status=thread_status;
1152 if (status == MagickFalse)
1153 fourier_image=DestroyImageList(fourier_image);
1159 return(fourier_image);
1196 #if defined(MAGICKCORE_FFTW_DELEGATE) && !defined(__cplusplus) && !defined(c_plusplus)
1197 static MagickBooleanType InverseQuadrantSwap(
const size_t width,
1198 const size_t height,
const double *source,
double *destination)
1208 center=(ssize_t) (width/2L)+1L;
1209 for (y=1L; y < (ssize_t) height; y++)
1210 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1211 destination[(height-y)*center-x+width/2L]=source[y*width+x];
1212 for (y=0L; y < (ssize_t) height; y++)
1213 destination[y*center]=source[y*width+width/2L];
1214 for (x=0L; x < center; x++)
1215 destination[x]=source[center-x-1L];
1216 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1219 static MagickBooleanType InverseFourier(
FourierInfo *fourier_info,
1220 const Image *magnitude_image,
const Image *phase_image,
1254 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1255 sizeof(*magnitude_pixels));
1256 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1257 sizeof(*phase_pixels));
1258 inverse_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/
1259 2+1)*
sizeof(*inverse_pixels));
1260 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1265 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1267 phase_info=RelinquishVirtualMemory(phase_info);
1269 inverse_info=RelinquishVirtualMemory(inverse_info);
1270 (void) ThrowMagickException(exception,GetMagickModule(),
1271 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1272 magnitude_image->filename);
1273 return(MagickFalse);
1275 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1276 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1277 inverse_pixels=(
double *) GetVirtualMemoryBlob(inverse_info);
1279 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1280 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1282 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1286 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
1287 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1289 switch (fourier_info->channel)
1294 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
1299 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
1304 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
1307 case OpacityChannel:
1309 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
1314 magnitude_pixels[i]=QuantumScale*(MagickRealType)
1315 GetPixelIndex(indexes+x);
1320 magnitude_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
1328 magnitude_view=DestroyCacheView(magnitude_view);
1329 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1330 magnitude_pixels,inverse_pixels);
1331 (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1332 fourier_info->center*
sizeof(*magnitude_pixels));
1334 phase_view=AcquireVirtualCacheView(phase_image,exception);
1335 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1337 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1341 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
1342 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1344 switch (fourier_info->channel)
1349 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelRed(p);
1354 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelGreen(p);
1359 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelBlue(p);
1362 case OpacityChannel:
1364 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelOpacity(p);
1369 phase_pixels[i]=QuantumScale*(MagickRealType)
1370 GetPixelIndex(indexes+x);
1375 phase_pixels[i]=QuantumScale*(MagickRealType) GetPixelGray(p);
1383 if (fourier_info->modulus != MagickFalse)
1386 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1387 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1389 phase_pixels[i]-=0.5;
1390 phase_pixels[i]*=(2.0*MagickPI);
1394 phase_view=DestroyCacheView(phase_view);
1395 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1396 if (status != MagickFalse)
1397 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1398 phase_pixels,inverse_pixels);
1399 (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1400 fourier_info->center*
sizeof(*phase_pixels));
1401 inverse_info=RelinquishVirtualMemory(inverse_info);
1406 if (fourier_info->modulus != MagickFalse)
1407 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1408 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1410 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1411 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+
1412 (MagickRealType) I*magnitude_pixels[i]*sin(phase_pixels[i]);
1414 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1415 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1420 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1421 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1423 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1424 fourier_pixels[i]=magnitude_pixels[i]+(MagickRealType) I*
1427 fourier_pixels[i][0]=magnitude_pixels[i];
1428 fourier_pixels[i][1]=phase_pixels[i];
1432 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1433 phase_info=RelinquishVirtualMemory(phase_info);
1437 static MagickBooleanType InverseFourierTransform(
FourierInfo *fourier_info,
1469 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
1471 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1472 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
1473 return(MagickFalse);
1475 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1476 sizeof(*source_pixels));
1479 (void) ThrowMagickException(exception,GetMagickModule(),
1480 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1481 return(MagickFalse);
1483 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
1484 value=GetImageArtifact(image,
"fourier:normalize");
1485 if (LocaleCompare(value,
"inverse") == 0)
1494 gamma=PerceptibleReciprocal((
double) fourier_info->width*
1495 fourier_info->height);
1496 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1497 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1499 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1500 fourier_pixels[i]*=gamma;
1502 fourier_pixels[i][0]*=gamma;
1503 fourier_pixels[i][1]*=gamma;
1508 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1509 #pragma omp critical (MagickCore_InverseFourierTransform)
1511 fftw_c2r_plan=fftw_plan_dft_c2r_2d((
int) fourier_info->width,
1512 (
int) fourier_info->height,fourier_pixels,source_pixels,FFTW_ESTIMATE);
1513 fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1514 fftw_destroy_plan(fftw_c2r_plan);
1516 image_view=AcquireAuthenticCacheView(image,exception);
1517 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1519 if (y >= (ssize_t) image->rows)
1521 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1522 image->columns ? image->columns : fourier_info->width,1UL,exception);
1525 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1526 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1528 if (x < (ssize_t) image->columns)
1529 switch (fourier_info->channel)
1534 SetPixelRed(q,ClampToQuantum((MagickRealType) QuantumRange*
1540 SetPixelGreen(q,ClampToQuantum((MagickRealType) QuantumRange*
1546 SetPixelBlue(q,ClampToQuantum((MagickRealType) QuantumRange*
1550 case OpacityChannel:
1552 SetPixelOpacity(q,ClampToQuantum((MagickRealType) QuantumRange*
1558 SetPixelIndex(indexes+x,ClampToQuantum((MagickRealType)
1559 QuantumRange*source_pixels[i]));
1564 SetPixelGray(q,ClampToQuantum((MagickRealType) QuantumRange*
1572 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1575 image_view=DestroyCacheView(image_view);
1576 source_info=RelinquishVirtualMemory(source_info);
1580 static MagickBooleanType InverseFourierTransformChannel(
1581 const Image *magnitude_image,
const Image *phase_image,
1582 const ChannelType channel,
const MagickBooleanType modulus,
1597 fourier_info.width=magnitude_image->columns;
1598 fourier_info.height=magnitude_image->rows;
1599 if ((magnitude_image->columns != magnitude_image->rows) ||
1600 ((magnitude_image->columns % 2) != 0) ||
1601 ((magnitude_image->rows % 2) != 0))
1603 size_t extent=magnitude_image->columns < magnitude_image->rows ?
1604 magnitude_image->rows : magnitude_image->columns;
1605 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1607 fourier_info.height=fourier_info.width;
1608 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1609 fourier_info.channel=channel;
1610 fourier_info.modulus=modulus;
1611 inverse_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1612 1)*
sizeof(*inverse_pixels));
1615 (void) ThrowMagickException(exception,GetMagickModule(),
1616 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1617 magnitude_image->filename);
1618 return(MagickFalse);
1620 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1621 status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1622 inverse_pixels,exception);
1623 if (status != MagickFalse)
1624 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1626 inverse_info=RelinquishVirtualMemory(inverse_info);
1631 MagickExport
Image *InverseFourierTransformImage(
const Image *magnitude_image,
1632 const Image *phase_image,
const MagickBooleanType modulus,
1638 assert(magnitude_image != (
Image *) NULL);
1639 assert(magnitude_image->signature == MagickCoreSignature);
1640 if (IsEventLogging() != MagickFalse)
1641 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1642 magnitude_image->filename);
1643 if (phase_image == (
Image *) NULL)
1645 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1646 "ImageSequenceRequired",
"`%s'",magnitude_image->filename);
1647 return((
Image *) NULL);
1649 #if !defined(MAGICKCORE_FFTW_DELEGATE) || defined(__cplusplus) || defined(c_plusplus)
1650 fourier_image=(
Image *) NULL;
1652 (void) ThrowMagickException(exception,GetMagickModule(),
1653 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1654 magnitude_image->filename);
1657 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1658 magnitude_image->rows,MagickTrue,exception);
1659 if (fourier_image != (
Image *) NULL)
1666 is_gray=IsGrayImage(magnitude_image,exception);
1667 if (is_gray != MagickFalse)
1668 is_gray=IsGrayImage(phase_image,exception);
1669 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1670 #pragma omp parallel sections
1673 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1680 if (is_gray != MagickFalse)
1681 thread_status=InverseFourierTransformChannel(magnitude_image,
1682 phase_image,GrayChannels,modulus,fourier_image,exception);
1684 thread_status=InverseFourierTransformChannel(magnitude_image,
1685 phase_image,RedChannel,modulus,fourier_image,exception);
1686 if (thread_status == MagickFalse)
1687 status=thread_status;
1689 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1696 thread_status=MagickTrue;
1697 if (is_gray == MagickFalse)
1698 thread_status=InverseFourierTransformChannel(magnitude_image,
1699 phase_image,GreenChannel,modulus,fourier_image,exception);
1700 if (thread_status == MagickFalse)
1701 status=thread_status;
1703 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1710 thread_status=MagickTrue;
1711 if (is_gray == MagickFalse)
1712 thread_status=InverseFourierTransformChannel(magnitude_image,
1713 phase_image,BlueChannel,modulus,fourier_image,exception);
1714 if (thread_status == MagickFalse)
1715 status=thread_status;
1717 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1724 thread_status=MagickTrue;
1725 if (magnitude_image->matte != MagickFalse)
1726 thread_status=InverseFourierTransformChannel(magnitude_image,
1727 phase_image,OpacityChannel,modulus,fourier_image,exception);
1728 if (thread_status == MagickFalse)
1729 status=thread_status;
1731 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1738 thread_status=MagickTrue;
1739 if (magnitude_image->colorspace == CMYKColorspace)
1740 thread_status=InverseFourierTransformChannel(magnitude_image,
1741 phase_image,IndexChannel,modulus,fourier_image,exception);
1742 if (thread_status == MagickFalse)
1743 status=thread_status;
1746 if (status == MagickFalse)
1747 fourier_image=DestroyImage(fourier_image);
1752 return(fourier_image);