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)
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++)
279 case AddComplexOperator:
281 Cr->red=Ar->red+Br->red;
282 Ci->red=Ai->red+Bi->red;
283 Cr->green=Ar->green+Br->green;
284 Ci->green=Ai->green+Bi->green;
285 Cr->blue=Ar->blue+Br->blue;
286 Ci->blue=Ai->blue+Bi->blue;
287 if (images->matte != MagickFalse)
289 Cr->opacity=Ar->opacity+Br->opacity;
290 Ci->opacity=Ai->opacity+Bi->opacity;
294 case ConjugateComplexOperator:
300 Ci->green=(-Ai->green);
302 Ci->blue=(-Ai->blue);
303 if (images->matte != MagickFalse)
305 Cr->opacity=Ar->opacity;
306 Ci->opacity=(-Ai->opacity);
310 case DivideComplexOperator:
315 gamma=QuantumRange*PerceptibleReciprocal(QuantumScale*Br->red*Br->red+
316 QuantumScale*Bi->red*Bi->red+snr);
317 Cr->red=gamma*(QuantumScale*Ar->red*Br->red+QuantumScale*Ai->red*
319 Ci->red=gamma*(QuantumScale*Ai->red*Br->red-QuantumScale*Ar->red*
321 gamma=QuantumRange*PerceptibleReciprocal(QuantumScale*Br->green*
322 Br->green+QuantumScale*Bi->green*Bi->green+snr);
323 Cr->green=gamma*(QuantumScale*Ar->green*Br->green+QuantumScale*
324 Ai->green*Bi->green);
325 Ci->green=gamma*(QuantumScale*Ai->green*Br->green-QuantumScale*
326 Ar->green*Bi->green);
327 gamma=QuantumRange*PerceptibleReciprocal(QuantumScale*Br->blue*
328 Br->blue+QuantumScale*Bi->blue*Bi->blue+snr);
329 Cr->blue=gamma*(QuantumScale*Ar->blue*Br->blue+QuantumScale*
331 Ci->blue=gamma*(QuantumScale*Ai->blue*Br->blue-QuantumScale*
333 if (images->matte != MagickFalse)
335 gamma=QuantumRange*PerceptibleReciprocal(QuantumScale*Br->opacity*
336 Br->opacity+QuantumScale*Bi->opacity*Bi->opacity+snr);
337 Cr->opacity=gamma*(QuantumScale*Ar->opacity*Br->opacity+
338 QuantumScale*Ai->opacity*Bi->opacity);
339 Ci->opacity=gamma*(QuantumScale*Ai->opacity*Br->opacity-
340 QuantumScale*Ar->opacity*Bi->opacity);
344 case MagnitudePhaseComplexOperator:
346 Cr->red=sqrt(QuantumScale*Ar->red*Ar->red+QuantumScale*
348 Ci->red=atan2((
double) Ai->red,(
double) Ar->red)/(2.0*MagickPI)+0.5;
349 Cr->green=sqrt(QuantumScale*Ar->green*Ar->green+QuantumScale*
350 Ai->green*Ai->green);
351 Ci->green=atan2((
double) Ai->green,(
double) Ar->green)/
353 Cr->blue=sqrt(QuantumScale*Ar->blue*Ar->blue+QuantumScale*
355 Ci->blue=atan2(Ai->blue,Ar->blue)/(2.0*MagickPI)+0.5;
356 if (images->matte != MagickFalse)
358 Cr->opacity=sqrt(QuantumScale*Ar->opacity*Ar->opacity+
359 QuantumScale*Ai->opacity*Ai->opacity);
360 Ci->opacity=atan2((
double) Ai->opacity,(
double) Ar->opacity)/
365 case MultiplyComplexOperator:
367 Cr->red=(QuantumScale*Ar->red*Br->red-(double)
369 Ci->red=(QuantumScale*Ai->red*Br->red+(double)
371 Cr->green=(QuantumScale*Ar->green*Br->green-(double)
372 Ai->green*Bi->green);
373 Ci->green=(QuantumScale*Ai->green*Br->green+(double)
374 Ar->green*Bi->green);
375 Cr->blue=(QuantumScale*Ar->blue*Br->blue-(double)
377 Ci->blue=(QuantumScale*Ai->blue*Br->blue+(double)
379 if (images->matte != MagickFalse)
381 Cr->opacity=(QuantumScale*Ar->opacity*Br->opacity-
382 QuantumScale*Ai->opacity*Bi->opacity);
383 Ci->opacity=(QuantumScale*Ai->opacity*Br->opacity+
384 QuantumScale*Ar->opacity*Bi->opacity);
388 case RealImaginaryComplexOperator:
390 Cr->red=Ar->red*cos(2.0*MagickPI*(Ai->red-0.5));
391 Ci->red=Ar->red*sin(2.0*MagickPI*(Ai->red-0.5));
392 Cr->green=Ar->green*cos(2.0*MagickPI*(Ai->green-0.5));
393 Ci->green=Ar->green*sin(2.0*MagickPI*(Ai->green-0.5));
394 Cr->blue=Ar->blue*cos(2.0*MagickPI*(Ai->blue-0.5));
395 Ci->blue=Ar->blue*sin(2.0*MagickPI*(Ai->blue-0.5));
396 if (images->matte != MagickFalse)
398 Cr->opacity=Ar->opacity*cos(2.0*MagickPI*(Ai->opacity-0.5));
399 Ci->opacity=Ar->opacity*sin(2.0*MagickPI*(Ai->opacity-0.5));
403 case SubtractComplexOperator:
405 Cr->red=Ar->red-Br->red;
406 Ci->red=Ai->red-Bi->red;
407 Cr->green=Ar->green-Br->green;
408 Ci->green=Ai->green-Bi->green;
409 Cr->blue=Ar->blue-Br->blue;
410 Ci->blue=Ai->blue-Bi->blue;
411 if (Cr_image->matte != MagickFalse)
413 Cr->opacity=Ar->opacity-Br->opacity;
414 Ci->opacity=Ai->opacity-Bi->opacity;
426 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
428 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
430 if (images->progress_monitor != (MagickProgressMonitor) NULL)
435 #if defined(MAGICKCORE_OPENMP_SUPPORT)
439 proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
440 if (proceed == MagickFalse)
444 Cr_view=DestroyCacheView(Cr_view);
445 Ci_view=DestroyCacheView(Ci_view);
446 Br_view=DestroyCacheView(Br_view);
447 Bi_view=DestroyCacheView(Bi_view);
448 Ar_view=DestroyCacheView(Ar_view);
449 Ai_view=DestroyCacheView(Ai_view);
450 if (status == MagickFalse)
451 complex_images=DestroyImageList(complex_images);
452 return(complex_images);
486 #if defined(MAGICKCORE_FFTW_DELEGATE)
488 static MagickBooleanType RollFourier(
const size_t width,
const size_t height,
489 const ssize_t x_offset,
const ssize_t y_offset,
double *roll_pixels)
509 source_info=AcquireVirtualMemory(width,height*
sizeof(*source_pixels));
512 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
514 for (y=0L; y < (ssize_t) height; y++)
517 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
519 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
521 for (x=0L; x < (ssize_t) width; x++)
524 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
526 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
528 source_pixels[v*width+u]=roll_pixels[i++];
531 (void) memcpy(roll_pixels,source_pixels,height*width*
sizeof(*source_pixels));
532 source_info=RelinquishVirtualMemory(source_info);
536 static MagickBooleanType ForwardQuadrantSwap(
const size_t width,
537 const size_t height,
double *source_pixels,
double *forward_pixels)
552 center=(ssize_t) (width/2L)+1L;
553 status=RollFourier((
size_t) center,height,0L,(ssize_t) height/2L,
555 if (status == MagickFalse)
557 for (y=0L; y < (ssize_t) height; y++)
558 for (x=0L; x < (ssize_t) (width/2L); x++)
559 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
560 for (y=1; y < (ssize_t) height; y++)
561 for (x=0L; x < (ssize_t) (width/2L); x++)
562 forward_pixels[(height-y)*width+width/2L-x-1L]=
563 source_pixels[y*center+x+1L];
564 for (x=0L; x < (ssize_t) (width/2L); x++)
565 forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
569 static void CorrectPhaseLHS(
const size_t width,
const size_t height,
570 double *fourier_pixels)
578 for (y=0L; y < (ssize_t) height; y++)
579 for (x=0L; x < (ssize_t) (width/2L); x++)
580 fourier_pixels[y*width+x]*=(-1.0);
583 static MagickBooleanType ForwardFourier(
const FourierInfo *fourier_info,
618 magnitude_image=GetFirstImageInList(image);
619 phase_image=GetNextImageInList(image);
620 if (phase_image == (
Image *) NULL)
622 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
623 "ImageSequenceRequired",
"`%s'",image->filename);
629 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
630 sizeof(*magnitude_pixels));
631 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
632 sizeof(*phase_pixels));
633 if ((magnitude_info == (
MemoryInfo *) NULL) ||
637 phase_info=RelinquishVirtualMemory(phase_info);
639 magnitude_info=RelinquishVirtualMemory(magnitude_info);
640 (void) ThrowMagickException(exception,GetMagickModule(),
641 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
644 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
645 (void) memset(magnitude_pixels,0,fourier_info->width*
646 fourier_info->height*
sizeof(*magnitude_pixels));
647 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
648 (void) memset(phase_pixels,0,fourier_info->width*
649 fourier_info->height*
sizeof(*phase_pixels));
650 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
651 magnitude,magnitude_pixels);
652 if (status != MagickFalse)
653 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
655 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
656 if (fourier_info->modulus != MagickFalse)
659 for (y=0L; y < (ssize_t) fourier_info->height; y++)
660 for (x=0L; x < (ssize_t) fourier_info->width; x++)
662 phase_pixels[i]/=(2.0*MagickPI);
663 phase_pixels[i]+=0.5;
667 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
669 for (y=0L; y < (ssize_t) fourier_info->height; y++)
671 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
675 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
676 for (x=0L; x < (ssize_t) fourier_info->width; x++)
678 switch (fourier_info->channel)
683 SetPixelRed(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
688 SetPixelGreen(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
693 SetPixelBlue(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
698 SetPixelOpacity(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
703 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*
704 magnitude_pixels[i]));
709 SetPixelGray(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
716 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
717 if (status == MagickFalse)
720 magnitude_view=DestroyCacheView(magnitude_view);
722 phase_view=AcquireAuthenticCacheView(phase_image,exception);
723 for (y=0L; y < (ssize_t) fourier_info->height; y++)
725 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
729 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
730 for (x=0L; x < (ssize_t) fourier_info->width; x++)
732 switch (fourier_info->channel)
737 SetPixelRed(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
742 SetPixelGreen(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
747 SetPixelBlue(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
752 SetPixelOpacity(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
757 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*phase_pixels[i]));
762 SetPixelGray(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
769 status=SyncCacheViewAuthenticPixels(phase_view,exception);
770 if (status == MagickFalse)
773 phase_view=DestroyCacheView(phase_view);
774 phase_info=RelinquishVirtualMemory(phase_info);
775 magnitude_info=RelinquishVirtualMemory(magnitude_info);
779 static MagickBooleanType ForwardFourierTransform(
FourierInfo *fourier_info,
780 const Image *image,
double *magnitude_pixels,
double *phase_pixels,
819 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
821 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
822 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
825 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
826 sizeof(*source_pixels));
829 (void) ThrowMagickException(exception,GetMagickModule(),
830 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
833 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
834 memset(source_pixels,0,fourier_info->width*fourier_info->height*
835 sizeof(*source_pixels));
837 image_view=AcquireVirtualCacheView(image,exception);
838 for (y=0L; y < (ssize_t) fourier_info->height; y++)
840 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
844 indexes=GetCacheViewVirtualIndexQueue(image_view);
845 for (x=0L; x < (ssize_t) fourier_info->width; x++)
847 switch (fourier_info->channel)
852 source_pixels[i]=QuantumScale*GetPixelRed(p);
857 source_pixels[i]=QuantumScale*GetPixelGreen(p);
862 source_pixels[i]=QuantumScale*GetPixelBlue(p);
867 source_pixels[i]=QuantumScale*GetPixelOpacity(p);
872 source_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
877 source_pixels[i]=QuantumScale*GetPixelGray(p);
885 image_view=DestroyCacheView(image_view);
886 forward_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/2+
887 1)*
sizeof(*forward_pixels));
890 (void) ThrowMagickException(exception,GetMagickModule(),
891 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
892 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
895 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
896 #if defined(MAGICKCORE_OPENMP_SUPPORT)
897 #pragma omp critical (MagickCore_ForwardFourierTransform)
899 fftw_r2c_plan=fftw_plan_dft_r2c_2d((
int) fourier_info->width,
900 (
int) fourier_info->height,source_pixels,forward_pixels,FFTW_ESTIMATE);
901 fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
902 fftw_destroy_plan(fftw_r2c_plan);
903 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
904 value=GetImageArtifact(image,
"fourier:normalize");
905 if ((value == (
const char *) NULL) || (LocaleCompare(value,
"forward") == 0))
914 gamma=PerceptibleReciprocal((
double) fourier_info->width*
915 fourier_info->height);
916 for (y=0L; y < (ssize_t) fourier_info->height; y++)
917 for (x=0L; x < (ssize_t) fourier_info->center; x++)
919 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
920 forward_pixels[i]*=gamma;
922 forward_pixels[i][0]*=gamma;
923 forward_pixels[i][1]*=gamma;
932 if (fourier_info->modulus != MagickFalse)
933 for (y=0L; y < (ssize_t) fourier_info->height; y++)
934 for (x=0L; x < (ssize_t) fourier_info->center; x++)
936 magnitude_pixels[i]=cabs(forward_pixels[i]);
937 phase_pixels[i]=carg(forward_pixels[i]);
941 for (y=0L; y < (ssize_t) fourier_info->height; y++)
942 for (x=0L; x < (ssize_t) fourier_info->center; x++)
944 magnitude_pixels[i]=creal(forward_pixels[i]);
945 phase_pixels[i]=cimag(forward_pixels[i]);
948 forward_info=(
MemoryInfo *) RelinquishVirtualMemory(forward_info);
952 static MagickBooleanType ForwardFourierTransformChannel(
const Image *image,
953 const ChannelType channel,
const MagickBooleanType modulus,
970 fourier_info.width=image->columns;
971 fourier_info.height=image->rows;
972 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
973 ((image->rows % 2) != 0))
975 size_t extent=image->columns < image->rows ? image->rows : image->columns;
976 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
978 fourier_info.height=fourier_info.width;
979 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
980 fourier_info.channel=channel;
981 fourier_info.modulus=modulus;
982 magnitude_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
983 1)*
sizeof(*magnitude_pixels));
984 phase_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+1)*
985 sizeof(*phase_pixels));
986 if ((magnitude_info == (
MemoryInfo *) NULL) ||
990 phase_info=RelinquishVirtualMemory(phase_info);
992 magnitude_info=RelinquishVirtualMemory(magnitude_info);
993 (void) ThrowMagickException(exception,GetMagickModule(),
994 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
997 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
998 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
999 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
1000 phase_pixels,exception);
1001 if (status != MagickFalse)
1002 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
1003 phase_pixels,exception);
1004 phase_info=RelinquishVirtualMemory(phase_info);
1005 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1010 MagickExport
Image *ForwardFourierTransformImage(
const Image *image,
1016 fourier_image=NewImageList();
1017 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1019 (void) ThrowMagickException(exception,GetMagickModule(),
1020 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1031 width=image->columns;
1033 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
1034 ((image->rows % 2) != 0))
1036 size_t extent=image->columns < image->rows ? image->rows :
1038 width=(extent & 0x01) == 1 ? extent+1UL : extent;
1041 magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
1042 if (magnitude_image != (
Image *) NULL)
1047 magnitude_image->storage_class=DirectClass;
1048 magnitude_image->depth=32UL;
1049 phase_image=CloneImage(image,width,height,MagickTrue,exception);
1050 if (phase_image == (
Image *) NULL)
1051 magnitude_image=DestroyImage(magnitude_image);
1058 phase_image->storage_class=DirectClass;
1059 phase_image->depth=32UL;
1060 AppendImageToList(&fourier_image,magnitude_image);
1061 AppendImageToList(&fourier_image,phase_image);
1063 is_gray=IsGrayImage(image,exception);
1064 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1065 #pragma omp parallel sections
1068 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1075 if (is_gray != MagickFalse)
1076 thread_status=ForwardFourierTransformChannel(image,
1077 GrayChannels,modulus,fourier_image,exception);
1079 thread_status=ForwardFourierTransformChannel(image,RedChannel,
1080 modulus,fourier_image,exception);
1081 if (thread_status == MagickFalse)
1082 status=thread_status;
1084 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1091 thread_status=MagickTrue;
1092 if (is_gray == MagickFalse)
1093 thread_status=ForwardFourierTransformChannel(image,
1094 GreenChannel,modulus,fourier_image,exception);
1095 if (thread_status == MagickFalse)
1096 status=thread_status;
1098 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1105 thread_status=MagickTrue;
1106 if (is_gray == MagickFalse)
1107 thread_status=ForwardFourierTransformChannel(image,
1108 BlueChannel,modulus,fourier_image,exception);
1109 if (thread_status == MagickFalse)
1110 status=thread_status;
1112 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1119 thread_status=MagickTrue;
1120 if (image->matte != MagickFalse)
1121 thread_status=ForwardFourierTransformChannel(image,
1122 OpacityChannel,modulus,fourier_image,exception);
1123 if (thread_status == MagickFalse)
1124 status=thread_status;
1126 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1133 thread_status=MagickTrue;
1134 if (image->colorspace == CMYKColorspace)
1135 thread_status=ForwardFourierTransformChannel(image,
1136 IndexChannel,modulus,fourier_image,exception);
1137 if (thread_status == MagickFalse)
1138 status=thread_status;
1141 if (status == MagickFalse)
1142 fourier_image=DestroyImageList(fourier_image);
1148 return(fourier_image);
1185 #if defined(MAGICKCORE_FFTW_DELEGATE)
1186 static MagickBooleanType InverseQuadrantSwap(
const size_t width,
1187 const size_t height,
const double *source,
double *destination)
1199 center=(ssize_t) (width/2L)+1L;
1200 for (y=1L; y < (ssize_t) height; y++)
1201 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1202 destination[(height-y)*center-x+width/2L]=source[y*width+x];
1203 for (y=0L; y < (ssize_t) height; y++)
1204 destination[y*center]=source[y*width+width/2L];
1205 for (x=0L; x < center; x++)
1206 destination[x]=source[center-x-1L];
1207 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1210 static MagickBooleanType InverseFourier(
FourierInfo *fourier_info,
1211 const Image *magnitude_image,
const Image *phase_image,
1247 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1248 sizeof(*magnitude_pixels));
1249 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1250 sizeof(*phase_pixels));
1251 inverse_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/
1252 2+1)*
sizeof(*inverse_pixels));
1253 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1258 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1260 phase_info=RelinquishVirtualMemory(phase_info);
1262 inverse_info=RelinquishVirtualMemory(inverse_info);
1263 (void) ThrowMagickException(exception,GetMagickModule(),
1264 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1265 magnitude_image->filename);
1266 return(MagickFalse);
1268 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1269 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1270 inverse_pixels=(
double *) GetVirtualMemoryBlob(inverse_info);
1272 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1273 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1275 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1279 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
1280 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1282 switch (fourier_info->channel)
1287 magnitude_pixels[i]=QuantumScale*GetPixelRed(p);
1292 magnitude_pixels[i]=QuantumScale*GetPixelGreen(p);
1297 magnitude_pixels[i]=QuantumScale*GetPixelBlue(p);
1300 case OpacityChannel:
1302 magnitude_pixels[i]=QuantumScale*GetPixelOpacity(p);
1307 magnitude_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
1312 magnitude_pixels[i]=QuantumScale*GetPixelGray(p);
1320 magnitude_view=DestroyCacheView(magnitude_view);
1321 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1322 magnitude_pixels,inverse_pixels);
1323 (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1324 fourier_info->center*
sizeof(*magnitude_pixels));
1326 phase_view=AcquireVirtualCacheView(phase_image,exception);
1327 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1329 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1333 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
1334 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1336 switch (fourier_info->channel)
1341 phase_pixels[i]=QuantumScale*GetPixelRed(p);
1346 phase_pixels[i]=QuantumScale*GetPixelGreen(p);
1351 phase_pixels[i]=QuantumScale*GetPixelBlue(p);
1354 case OpacityChannel:
1356 phase_pixels[i]=QuantumScale*GetPixelOpacity(p);
1361 phase_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
1366 phase_pixels[i]=QuantumScale*GetPixelGray(p);
1374 if (fourier_info->modulus != MagickFalse)
1377 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1378 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1380 phase_pixels[i]-=0.5;
1381 phase_pixels[i]*=(2.0*MagickPI);
1385 phase_view=DestroyCacheView(phase_view);
1386 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1387 if (status != MagickFalse)
1388 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1389 phase_pixels,inverse_pixels);
1390 (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1391 fourier_info->center*
sizeof(*phase_pixels));
1392 inverse_info=RelinquishVirtualMemory(inverse_info);
1397 if (fourier_info->modulus != MagickFalse)
1398 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1399 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1401 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1402 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1403 magnitude_pixels[i]*sin(phase_pixels[i]);
1405 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1406 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1411 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1412 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1414 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1415 fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1417 fourier_pixels[i][0]=magnitude_pixels[i];
1418 fourier_pixels[i][1]=phase_pixels[i];
1422 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1423 phase_info=RelinquishVirtualMemory(phase_info);
1427 static MagickBooleanType InverseFourierTransform(
FourierInfo *fourier_info,
1461 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
1463 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1464 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
1465 return(MagickFalse);
1467 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1468 sizeof(*source_pixels));
1471 (void) ThrowMagickException(exception,GetMagickModule(),
1472 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1473 return(MagickFalse);
1475 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
1476 value=GetImageArtifact(image,
"fourier:normalize");
1477 if (LocaleCompare(value,
"inverse") == 0)
1486 gamma=PerceptibleReciprocal((
double) fourier_info->width*
1487 fourier_info->height);
1488 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1489 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1491 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1492 fourier_pixels[i]*=gamma;
1494 fourier_pixels[i][0]*=gamma;
1495 fourier_pixels[i][1]*=gamma;
1500 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1501 #pragma omp critical (MagickCore_InverseFourierTransform)
1503 fftw_c2r_plan=fftw_plan_dft_c2r_2d((
int) fourier_info->width,
1504 (
int) fourier_info->height,fourier_pixels,source_pixels,FFTW_ESTIMATE);
1505 fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1506 fftw_destroy_plan(fftw_c2r_plan);
1508 image_view=AcquireAuthenticCacheView(image,exception);
1509 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1511 if (y >= (ssize_t) image->rows)
1513 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1514 image->columns ? image->columns : fourier_info->width,1UL,exception);
1517 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1518 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1520 if (x < (ssize_t) image->columns)
1521 switch (fourier_info->channel)
1526 SetPixelRed(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1531 SetPixelGreen(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1536 SetPixelBlue(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1539 case OpacityChannel:
1541 SetPixelOpacity(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1546 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*
1552 SetPixelGray(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1559 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1562 image_view=DestroyCacheView(image_view);
1563 source_info=RelinquishVirtualMemory(source_info);
1567 static MagickBooleanType InverseFourierTransformChannel(
1568 const Image *magnitude_image,
const Image *phase_image,
1569 const ChannelType channel,
const MagickBooleanType modulus,
1584 fourier_info.width=magnitude_image->columns;
1585 fourier_info.height=magnitude_image->rows;
1586 if ((magnitude_image->columns != magnitude_image->rows) ||
1587 ((magnitude_image->columns % 2) != 0) ||
1588 ((magnitude_image->rows % 2) != 0))
1590 size_t extent=magnitude_image->columns < magnitude_image->rows ?
1591 magnitude_image->rows : magnitude_image->columns;
1592 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1594 fourier_info.height=fourier_info.width;
1595 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1596 fourier_info.channel=channel;
1597 fourier_info.modulus=modulus;
1598 inverse_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1599 1)*
sizeof(*inverse_pixels));
1602 (void) ThrowMagickException(exception,GetMagickModule(),
1603 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1604 magnitude_image->filename);
1605 return(MagickFalse);
1607 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1608 status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1609 inverse_pixels,exception);
1610 if (status != MagickFalse)
1611 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1613 inverse_info=RelinquishVirtualMemory(inverse_info);
1618 MagickExport
Image *InverseFourierTransformImage(
const Image *magnitude_image,
1619 const Image *phase_image,
const MagickBooleanType modulus,
1625 assert(magnitude_image != (
Image *) NULL);
1626 assert(magnitude_image->signature == MagickCoreSignature);
1627 if (IsEventLogging() != MagickFalse)
1628 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1629 magnitude_image->filename);
1630 if (phase_image == (
Image *) NULL)
1632 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1633 "ImageSequenceRequired",
"`%s'",magnitude_image->filename);
1634 return((
Image *) NULL);
1636 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1637 fourier_image=(
Image *) NULL;
1639 (void) ThrowMagickException(exception,GetMagickModule(),
1640 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1641 magnitude_image->filename);
1644 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1645 magnitude_image->rows,MagickTrue,exception);
1646 if (fourier_image != (
Image *) NULL)
1653 is_gray=IsGrayImage(magnitude_image,exception);
1654 if (is_gray != MagickFalse)
1655 is_gray=IsGrayImage(phase_image,exception);
1656 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1657 #pragma omp parallel sections
1660 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1667 if (is_gray != MagickFalse)
1668 thread_status=InverseFourierTransformChannel(magnitude_image,
1669 phase_image,GrayChannels,modulus,fourier_image,exception);
1671 thread_status=InverseFourierTransformChannel(magnitude_image,
1672 phase_image,RedChannel,modulus,fourier_image,exception);
1673 if (thread_status == MagickFalse)
1674 status=thread_status;
1676 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1683 thread_status=MagickTrue;
1684 if (is_gray == MagickFalse)
1685 thread_status=InverseFourierTransformChannel(magnitude_image,
1686 phase_image,GreenChannel,modulus,fourier_image,exception);
1687 if (thread_status == MagickFalse)
1688 status=thread_status;
1690 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1697 thread_status=MagickTrue;
1698 if (is_gray == MagickFalse)
1699 thread_status=InverseFourierTransformChannel(magnitude_image,
1700 phase_image,BlueChannel,modulus,fourier_image,exception);
1701 if (thread_status == MagickFalse)
1702 status=thread_status;
1704 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1711 thread_status=MagickTrue;
1712 if (magnitude_image->matte != MagickFalse)
1713 thread_status=InverseFourierTransformChannel(magnitude_image,
1714 phase_image,OpacityChannel,modulus,fourier_image,exception);
1715 if (thread_status == MagickFalse)
1716 status=thread_status;
1718 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1725 thread_status=MagickTrue;
1726 if (magnitude_image->colorspace == CMYKColorspace)
1727 thread_status=InverseFourierTransformChannel(magnitude_image,
1728 phase_image,IndexChannel,modulus,fourier_image,exception);
1729 if (thread_status == MagickFalse)
1730 status=thread_status;
1733 if (status == MagickFalse)
1734 fourier_image=DestroyImage(fourier_image);
1739 return(fourier_image);