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++)
278 ai = { (Quantum) (QuantumScale*GetPixelRed(Ai)),
279 (Quantum) (QuantumScale*GetPixelGreen(Ai)),
280 (Quantum) (QuantumScale*GetPixelBlue(Ai)),
281 (Quantum) (image->matte != MagickFalse ? QuantumScale*
282 GetPixelOpacity(Ai) : OpaqueOpacity), 0 },
283 ar = { (Quantum) (QuantumScale*GetPixelRed(Ar)),
284 (Quantum) (QuantumScale*GetPixelGreen(Ar)),
285 (Quantum) (QuantumScale*GetPixelBlue(Ar)),
286 (Quantum) (image->matte != MagickFalse ? QuantumScale*
287 GetPixelOpacity(Ar) : OpaqueOpacity), 0 },
288 bi = { (Quantum) (QuantumScale*GetPixelRed(Bi)),
289 (Quantum) (QuantumScale*GetPixelGreen(Bi)),
290 (Quantum) (QuantumScale*GetPixelBlue(Bi)),
291 (Quantum) (image->matte != MagickFalse ? QuantumScale*
292 GetPixelOpacity(Bi) : OpaqueOpacity), 0 },
293 br = { (Quantum) (QuantumScale*GetPixelRed(Br)),
294 (Quantum) (QuantumScale*GetPixelGreen(Br)),
295 (Quantum) (QuantumScale*GetPixelBlue(Br)),
296 (Quantum) (image->matte != MagickFalse ? QuantumScale*
297 GetPixelOpacity(Br) : OpaqueOpacity), 0 },
303 case AddComplexOperator:
305 cr.red=ar.red+br.red;
306 ci.red=ai.red+bi.red;
307 cr.green=ar.green+br.green;
308 ci.green=ai.green+bi.green;
309 cr.blue=ar.blue+br.blue;
310 ci.blue=ai.blue+bi.blue;
311 cr.opacity=ar.opacity+br.opacity;
312 ci.opacity=ai.opacity+bi.opacity;
315 case ConjugateComplexOperator:
321 ci.green=(-ai.green);
324 cr.opacity=ar.opacity;
325 ci.opacity=(-ai.opacity);
328 case DivideComplexOperator:
333 gamma=PerceptibleReciprocal(br.red*br.red+bi.red*bi.red+snr);
334 cr.red=gamma*(ar.red*br.red+ai.red*bi.red);
335 ci.red=gamma*(ai.red*br.red-ar.red*bi.red);
336 gamma=PerceptibleReciprocal(br.green*br.green+bi.green*bi.green+snr);
337 cr.green=gamma*(ar.green*br.green+ai.green*bi.green);
338 ci.green=gamma*(ai.green*br.green-ar.green*bi.green);
339 gamma=PerceptibleReciprocal(br.blue*br.blue+bi.blue*bi.blue+snr);
340 cr.blue=gamma*(ar.blue*br.blue+ai.blue*bi.blue);
341 ci.blue=gamma*(ai.blue*br.blue-ar.blue*bi.blue);
342 gamma=PerceptibleReciprocal(br.opacity*br.opacity+bi.opacity*
344 cr.opacity=gamma*(ar.opacity*br.opacity+ai.opacity*bi.opacity);
345 ci.opacity=gamma*(ai.opacity*br.opacity-ar.opacity*bi.opacity);
348 case MagnitudePhaseComplexOperator:
350 cr.red=sqrt((
double) ar.red*ar.red+ai.red*ai.red);
351 ci.red=atan2((
double) ai.red,(
double) ar.red)/(2.0*MagickPI)+0.5;
352 cr.green=sqrt((
double) ar.green*ar.green+ai.green*ai.green);
353 ci.green=atan2((
double) ai.green,(
double) ar.green)/
355 cr.blue=sqrt((
double) ar.blue*ar.blue+ai.blue*ai.blue);
356 ci.blue=atan2((
double) ai.blue,(
double) ar.blue)/(2.0*MagickPI)+0.5;
357 cr.opacity=sqrt((
double) ar.opacity*ar.opacity+ai.opacity*ai.opacity);
358 ci.opacity=atan2((
double) ai.opacity,(
double) ar.opacity)/
362 case MultiplyComplexOperator:
364 cr.red=(ar.red*br.red-(double) ai.red*bi.red);
365 ci.red=(ai.red*br.red+(double) ar.red*bi.red);
366 cr.green=(ar.green*br.green-(double) ai.green*bi.green);
367 ci.green=(ai.green*br.green+(double) ar.green*bi.green);
368 cr.blue=(ar.blue*br.blue-(double) ai.blue*bi.blue);
369 ci.blue=(ai.blue*br.blue+(double) ar.blue*bi.blue);
370 cr.opacity=(ar.opacity*br.opacity-(double) ai.opacity*bi.opacity);
371 ci.opacity=(ai.opacity*br.opacity+(double) ar.opacity*bi.opacity);
374 case RealImaginaryComplexOperator:
376 cr.red=ar.red*cos(2.0*MagickPI*(ai.red-0.5));
377 ci.red=ar.red*sin(2.0*MagickPI*(ai.red-0.5));
378 cr.green=ar.green*cos(2.0*MagickPI*(ai.green-0.5));
379 ci.green=ar.green*sin(2.0*MagickPI*(ai.green-0.5));
380 cr.blue=ar.blue*cos(2.0*MagickPI*(ai.blue-0.5));
381 ci.blue=ar.blue*sin(2.0*MagickPI*(ai.blue-0.5));
382 cr.opacity=ar.opacity*cos(2.0*MagickPI*(ai.opacity-0.5));
383 ci.opacity=ar.opacity*sin(2.0*MagickPI*(ai.opacity-0.5));
386 case SubtractComplexOperator:
388 cr.red=ar.red-br.red;
389 ci.red=ai.red-bi.red;
390 cr.green=ar.green-br.green;
391 ci.green=ai.green-bi.green;
392 cr.blue=ar.blue-br.blue;
393 ci.blue=ai.blue-bi.blue;
394 cr.opacity=ar.opacity-br.opacity;
395 ci.opacity=ai.opacity-bi.opacity;
399 Cr->red=QuantumRange*cr.red;
400 Ci->red=QuantumRange*ci.red;
401 Cr->green=QuantumRange*cr.green;
402 Ci->green=QuantumRange*ci.green;
403 Cr->blue=QuantumRange*cr.blue;
404 Ci->blue=QuantumRange*ci.blue;
405 if (images->matte != MagickFalse)
407 Cr->opacity=QuantumRange*cr.opacity;
408 Ci->opacity=QuantumRange*ci.opacity;
417 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
419 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
421 if (images->progress_monitor != (MagickProgressMonitor) NULL)
426 #if defined(MAGICKCORE_OPENMP_SUPPORT)
430 proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
431 if (proceed == MagickFalse)
435 Cr_view=DestroyCacheView(Cr_view);
436 Ci_view=DestroyCacheView(Ci_view);
437 Br_view=DestroyCacheView(Br_view);
438 Bi_view=DestroyCacheView(Bi_view);
439 Ar_view=DestroyCacheView(Ar_view);
440 Ai_view=DestroyCacheView(Ai_view);
441 if (status == MagickFalse)
442 complex_images=DestroyImageList(complex_images);
443 return(complex_images);
477 #if defined(MAGICKCORE_FFTW_DELEGATE)
479 static MagickBooleanType RollFourier(
const size_t width,
const size_t height,
480 const ssize_t x_offset,
const ssize_t y_offset,
double *roll_pixels)
500 source_info=AcquireVirtualMemory(width,height*
sizeof(*source_pixels));
503 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
505 for (y=0L; y < (ssize_t) height; y++)
508 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
510 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
512 for (x=0L; x < (ssize_t) width; x++)
515 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
517 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
519 source_pixels[v*width+u]=roll_pixels[i++];
522 (void) memcpy(roll_pixels,source_pixels,height*width*
sizeof(*source_pixels));
523 source_info=RelinquishVirtualMemory(source_info);
527 static MagickBooleanType ForwardQuadrantSwap(
const size_t width,
528 const size_t height,
double *source_pixels,
double *forward_pixels)
543 center=(ssize_t) (width/2L)+1L;
544 status=RollFourier((
size_t) center,height,0L,(ssize_t) height/2L,
546 if (status == MagickFalse)
548 for (y=0L; y < (ssize_t) height; y++)
549 for (x=0L; x < (ssize_t) (width/2L); x++)
550 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
551 for (y=1; y < (ssize_t) height; y++)
552 for (x=0L; x < (ssize_t) (width/2L); x++)
553 forward_pixels[(height-y)*width+width/2L-x-1L]=
554 source_pixels[y*center+x+1L];
555 for (x=0L; x < (ssize_t) (width/2L); x++)
556 forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
560 static void CorrectPhaseLHS(
const size_t width,
const size_t height,
561 double *fourier_pixels)
569 for (y=0L; y < (ssize_t) height; y++)
570 for (x=0L; x < (ssize_t) (width/2L); x++)
571 fourier_pixels[y*width+x]*=(-1.0);
574 static MagickBooleanType ForwardFourier(
const FourierInfo *fourier_info,
609 magnitude_image=GetFirstImageInList(image);
610 phase_image=GetNextImageInList(image);
611 if (phase_image == (
Image *) NULL)
613 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
614 "ImageSequenceRequired",
"`%s'",image->filename);
620 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
621 sizeof(*magnitude_pixels));
622 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
623 sizeof(*phase_pixels));
624 if ((magnitude_info == (
MemoryInfo *) NULL) ||
628 phase_info=RelinquishVirtualMemory(phase_info);
630 magnitude_info=RelinquishVirtualMemory(magnitude_info);
631 (void) ThrowMagickException(exception,GetMagickModule(),
632 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
635 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
636 (void) memset(magnitude_pixels,0,fourier_info->width*
637 fourier_info->height*
sizeof(*magnitude_pixels));
638 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
639 (void) memset(phase_pixels,0,fourier_info->width*
640 fourier_info->height*
sizeof(*phase_pixels));
641 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
642 magnitude,magnitude_pixels);
643 if (status != MagickFalse)
644 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
646 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
647 if (fourier_info->modulus != MagickFalse)
650 for (y=0L; y < (ssize_t) fourier_info->height; y++)
651 for (x=0L; x < (ssize_t) fourier_info->width; x++)
653 phase_pixels[i]/=(2.0*MagickPI);
654 phase_pixels[i]+=0.5;
658 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
660 for (y=0L; y < (ssize_t) fourier_info->height; y++)
662 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
666 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
667 for (x=0L; x < (ssize_t) fourier_info->width; x++)
669 switch (fourier_info->channel)
674 SetPixelRed(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
679 SetPixelGreen(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
684 SetPixelBlue(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
689 SetPixelOpacity(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
694 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*
695 magnitude_pixels[i]));
700 SetPixelGray(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
707 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
708 if (status == MagickFalse)
711 magnitude_view=DestroyCacheView(magnitude_view);
713 phase_view=AcquireAuthenticCacheView(phase_image,exception);
714 for (y=0L; y < (ssize_t) fourier_info->height; y++)
716 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
720 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
721 for (x=0L; x < (ssize_t) fourier_info->width; x++)
723 switch (fourier_info->channel)
728 SetPixelRed(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
733 SetPixelGreen(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
738 SetPixelBlue(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
743 SetPixelOpacity(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
748 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*phase_pixels[i]));
753 SetPixelGray(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
760 status=SyncCacheViewAuthenticPixels(phase_view,exception);
761 if (status == MagickFalse)
764 phase_view=DestroyCacheView(phase_view);
765 phase_info=RelinquishVirtualMemory(phase_info);
766 magnitude_info=RelinquishVirtualMemory(magnitude_info);
770 static MagickBooleanType ForwardFourierTransform(
FourierInfo *fourier_info,
771 const Image *image,
double *magnitude_pixels,
double *phase_pixels,
810 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
812 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
813 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
816 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
817 sizeof(*source_pixels));
820 (void) ThrowMagickException(exception,GetMagickModule(),
821 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
824 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
825 memset(source_pixels,0,fourier_info->width*fourier_info->height*
826 sizeof(*source_pixels));
828 image_view=AcquireVirtualCacheView(image,exception);
829 for (y=0L; y < (ssize_t) fourier_info->height; y++)
831 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
835 indexes=GetCacheViewVirtualIndexQueue(image_view);
836 for (x=0L; x < (ssize_t) fourier_info->width; x++)
838 switch (fourier_info->channel)
843 source_pixels[i]=QuantumScale*GetPixelRed(p);
848 source_pixels[i]=QuantumScale*GetPixelGreen(p);
853 source_pixels[i]=QuantumScale*GetPixelBlue(p);
858 source_pixels[i]=QuantumScale*GetPixelOpacity(p);
863 source_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
868 source_pixels[i]=QuantumScale*GetPixelGray(p);
876 image_view=DestroyCacheView(image_view);
877 forward_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/2+
878 1)*
sizeof(*forward_pixels));
881 (void) ThrowMagickException(exception,GetMagickModule(),
882 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
883 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
886 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
887 #if defined(MAGICKCORE_OPENMP_SUPPORT)
888 #pragma omp critical (MagickCore_ForwardFourierTransform)
890 fftw_r2c_plan=fftw_plan_dft_r2c_2d((
int) fourier_info->width,
891 (
int) fourier_info->height,source_pixels,forward_pixels,FFTW_ESTIMATE);
892 fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
893 fftw_destroy_plan(fftw_r2c_plan);
894 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
895 value=GetImageArtifact(image,
"fourier:normalize");
896 if ((value == (
const char *) NULL) || (LocaleCompare(value,
"forward") == 0))
905 gamma=PerceptibleReciprocal((
double) fourier_info->width*
906 fourier_info->height);
907 for (y=0L; y < (ssize_t) fourier_info->height; y++)
908 for (x=0L; x < (ssize_t) fourier_info->center; x++)
910 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
911 forward_pixels[i]*=gamma;
913 forward_pixels[i][0]*=gamma;
914 forward_pixels[i][1]*=gamma;
923 if (fourier_info->modulus != MagickFalse)
924 for (y=0L; y < (ssize_t) fourier_info->height; y++)
925 for (x=0L; x < (ssize_t) fourier_info->center; x++)
927 magnitude_pixels[i]=cabs(forward_pixels[i]);
928 phase_pixels[i]=carg(forward_pixels[i]);
932 for (y=0L; y < (ssize_t) fourier_info->height; y++)
933 for (x=0L; x < (ssize_t) fourier_info->center; x++)
935 magnitude_pixels[i]=creal(forward_pixels[i]);
936 phase_pixels[i]=cimag(forward_pixels[i]);
939 forward_info=(
MemoryInfo *) RelinquishVirtualMemory(forward_info);
943 static MagickBooleanType ForwardFourierTransformChannel(
const Image *image,
944 const ChannelType channel,
const MagickBooleanType modulus,
961 fourier_info.width=image->columns;
962 fourier_info.height=image->rows;
963 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
964 ((image->rows % 2) != 0))
966 size_t extent=image->columns < image->rows ? image->rows : image->columns;
967 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
969 fourier_info.height=fourier_info.width;
970 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
971 fourier_info.channel=channel;
972 fourier_info.modulus=modulus;
973 magnitude_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
974 1)*
sizeof(*magnitude_pixels));
975 phase_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+1)*
976 sizeof(*phase_pixels));
977 if ((magnitude_info == (
MemoryInfo *) NULL) ||
981 phase_info=RelinquishVirtualMemory(phase_info);
983 magnitude_info=RelinquishVirtualMemory(magnitude_info);
984 (void) ThrowMagickException(exception,GetMagickModule(),
985 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
988 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
989 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
990 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
991 phase_pixels,exception);
992 if (status != MagickFalse)
993 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
994 phase_pixels,exception);
995 phase_info=RelinquishVirtualMemory(phase_info);
996 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1001 MagickExport
Image *ForwardFourierTransformImage(
const Image *image,
1007 fourier_image=NewImageList();
1008 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1010 (void) ThrowMagickException(exception,GetMagickModule(),
1011 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1022 width=image->columns;
1024 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
1025 ((image->rows % 2) != 0))
1027 size_t extent=image->columns < image->rows ? image->rows :
1029 width=(extent & 0x01) == 1 ? extent+1UL : extent;
1032 magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
1033 if (magnitude_image != (
Image *) NULL)
1038 magnitude_image->storage_class=DirectClass;
1039 magnitude_image->depth=32UL;
1040 phase_image=CloneImage(image,width,height,MagickTrue,exception);
1041 if (phase_image == (
Image *) NULL)
1042 magnitude_image=DestroyImage(magnitude_image);
1049 phase_image->storage_class=DirectClass;
1050 phase_image->depth=32UL;
1051 AppendImageToList(&fourier_image,magnitude_image);
1052 AppendImageToList(&fourier_image,phase_image);
1054 is_gray=IsGrayImage(image,exception);
1055 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1056 #pragma omp parallel sections
1059 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1066 if (is_gray != MagickFalse)
1067 thread_status=ForwardFourierTransformChannel(image,
1068 GrayChannels,modulus,fourier_image,exception);
1070 thread_status=ForwardFourierTransformChannel(image,RedChannel,
1071 modulus,fourier_image,exception);
1072 if (thread_status == MagickFalse)
1073 status=thread_status;
1075 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1082 thread_status=MagickTrue;
1083 if (is_gray == MagickFalse)
1084 thread_status=ForwardFourierTransformChannel(image,
1085 GreenChannel,modulus,fourier_image,exception);
1086 if (thread_status == MagickFalse)
1087 status=thread_status;
1089 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1096 thread_status=MagickTrue;
1097 if (is_gray == MagickFalse)
1098 thread_status=ForwardFourierTransformChannel(image,
1099 BlueChannel,modulus,fourier_image,exception);
1100 if (thread_status == MagickFalse)
1101 status=thread_status;
1103 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1110 thread_status=MagickTrue;
1111 if (image->matte != MagickFalse)
1112 thread_status=ForwardFourierTransformChannel(image,
1113 OpacityChannel,modulus,fourier_image,exception);
1114 if (thread_status == MagickFalse)
1115 status=thread_status;
1117 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1124 thread_status=MagickTrue;
1125 if (image->colorspace == CMYKColorspace)
1126 thread_status=ForwardFourierTransformChannel(image,
1127 IndexChannel,modulus,fourier_image,exception);
1128 if (thread_status == MagickFalse)
1129 status=thread_status;
1132 if (status == MagickFalse)
1133 fourier_image=DestroyImageList(fourier_image);
1139 return(fourier_image);
1176 #if defined(MAGICKCORE_FFTW_DELEGATE)
1177 static MagickBooleanType InverseQuadrantSwap(
const size_t width,
1178 const size_t height,
const double *source,
double *destination)
1190 center=(ssize_t) (width/2L)+1L;
1191 for (y=1L; y < (ssize_t) height; y++)
1192 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1193 destination[(height-y)*center-x+width/2L]=source[y*width+x];
1194 for (y=0L; y < (ssize_t) height; y++)
1195 destination[y*center]=source[y*width+width/2L];
1196 for (x=0L; x < center; x++)
1197 destination[x]=source[center-x-1L];
1198 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1201 static MagickBooleanType InverseFourier(
FourierInfo *fourier_info,
1202 const Image *magnitude_image,
const Image *phase_image,
1238 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1239 sizeof(*magnitude_pixels));
1240 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1241 sizeof(*phase_pixels));
1242 inverse_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/
1243 2+1)*
sizeof(*inverse_pixels));
1244 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1249 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1251 phase_info=RelinquishVirtualMemory(phase_info);
1253 inverse_info=RelinquishVirtualMemory(inverse_info);
1254 (void) ThrowMagickException(exception,GetMagickModule(),
1255 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1256 magnitude_image->filename);
1257 return(MagickFalse);
1259 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1260 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1261 inverse_pixels=(
double *) GetVirtualMemoryBlob(inverse_info);
1263 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1264 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1266 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1270 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
1271 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1273 switch (fourier_info->channel)
1278 magnitude_pixels[i]=QuantumScale*GetPixelRed(p);
1283 magnitude_pixels[i]=QuantumScale*GetPixelGreen(p);
1288 magnitude_pixels[i]=QuantumScale*GetPixelBlue(p);
1291 case OpacityChannel:
1293 magnitude_pixels[i]=QuantumScale*GetPixelOpacity(p);
1298 magnitude_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
1303 magnitude_pixels[i]=QuantumScale*GetPixelGray(p);
1311 magnitude_view=DestroyCacheView(magnitude_view);
1312 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1313 magnitude_pixels,inverse_pixels);
1314 (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1315 fourier_info->center*
sizeof(*magnitude_pixels));
1317 phase_view=AcquireVirtualCacheView(phase_image,exception);
1318 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1320 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1324 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
1325 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1327 switch (fourier_info->channel)
1332 phase_pixels[i]=QuantumScale*GetPixelRed(p);
1337 phase_pixels[i]=QuantumScale*GetPixelGreen(p);
1342 phase_pixels[i]=QuantumScale*GetPixelBlue(p);
1345 case OpacityChannel:
1347 phase_pixels[i]=QuantumScale*GetPixelOpacity(p);
1352 phase_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
1357 phase_pixels[i]=QuantumScale*GetPixelGray(p);
1365 if (fourier_info->modulus != MagickFalse)
1368 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1369 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1371 phase_pixels[i]-=0.5;
1372 phase_pixels[i]*=(2.0*MagickPI);
1376 phase_view=DestroyCacheView(phase_view);
1377 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1378 if (status != MagickFalse)
1379 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1380 phase_pixels,inverse_pixels);
1381 (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1382 fourier_info->center*
sizeof(*phase_pixels));
1383 inverse_info=RelinquishVirtualMemory(inverse_info);
1388 if (fourier_info->modulus != MagickFalse)
1389 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1390 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1392 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1393 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1394 magnitude_pixels[i]*sin(phase_pixels[i]);
1396 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1397 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1402 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1403 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1405 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1406 fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1408 fourier_pixels[i][0]=magnitude_pixels[i];
1409 fourier_pixels[i][1]=phase_pixels[i];
1413 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1414 phase_info=RelinquishVirtualMemory(phase_info);
1418 static MagickBooleanType InverseFourierTransform(
FourierInfo *fourier_info,
1452 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
1454 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1455 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
1456 return(MagickFalse);
1458 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1459 sizeof(*source_pixels));
1462 (void) ThrowMagickException(exception,GetMagickModule(),
1463 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1464 return(MagickFalse);
1466 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
1467 value=GetImageArtifact(image,
"fourier:normalize");
1468 if (LocaleCompare(value,
"inverse") == 0)
1477 gamma=PerceptibleReciprocal((
double) fourier_info->width*
1478 fourier_info->height);
1479 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1480 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1482 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1483 fourier_pixels[i]*=gamma;
1485 fourier_pixels[i][0]*=gamma;
1486 fourier_pixels[i][1]*=gamma;
1491 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1492 #pragma omp critical (MagickCore_InverseFourierTransform)
1494 fftw_c2r_plan=fftw_plan_dft_c2r_2d((
int) fourier_info->width,
1495 (
int) fourier_info->height,fourier_pixels,source_pixels,FFTW_ESTIMATE);
1496 fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1497 fftw_destroy_plan(fftw_c2r_plan);
1499 image_view=AcquireAuthenticCacheView(image,exception);
1500 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1502 if (y >= (ssize_t) image->rows)
1504 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1505 image->columns ? image->columns : fourier_info->width,1UL,exception);
1508 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1509 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1511 if (x < (ssize_t) image->columns)
1512 switch (fourier_info->channel)
1517 SetPixelRed(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1522 SetPixelGreen(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1527 SetPixelBlue(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1530 case OpacityChannel:
1532 SetPixelOpacity(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1537 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*
1543 SetPixelGray(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1550 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1553 image_view=DestroyCacheView(image_view);
1554 source_info=RelinquishVirtualMemory(source_info);
1558 static MagickBooleanType InverseFourierTransformChannel(
1559 const Image *magnitude_image,
const Image *phase_image,
1560 const ChannelType channel,
const MagickBooleanType modulus,
1575 fourier_info.width=magnitude_image->columns;
1576 fourier_info.height=magnitude_image->rows;
1577 if ((magnitude_image->columns != magnitude_image->rows) ||
1578 ((magnitude_image->columns % 2) != 0) ||
1579 ((magnitude_image->rows % 2) != 0))
1581 size_t extent=magnitude_image->columns < magnitude_image->rows ?
1582 magnitude_image->rows : magnitude_image->columns;
1583 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1585 fourier_info.height=fourier_info.width;
1586 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1587 fourier_info.channel=channel;
1588 fourier_info.modulus=modulus;
1589 inverse_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1590 1)*
sizeof(*inverse_pixels));
1593 (void) ThrowMagickException(exception,GetMagickModule(),
1594 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1595 magnitude_image->filename);
1596 return(MagickFalse);
1598 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1599 status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1600 inverse_pixels,exception);
1601 if (status != MagickFalse)
1602 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1604 inverse_info=RelinquishVirtualMemory(inverse_info);
1609 MagickExport
Image *InverseFourierTransformImage(
const Image *magnitude_image,
1610 const Image *phase_image,
const MagickBooleanType modulus,
1616 assert(magnitude_image != (
Image *) NULL);
1617 assert(magnitude_image->signature == MagickCoreSignature);
1618 if (IsEventLogging() != MagickFalse)
1619 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1620 magnitude_image->filename);
1621 if (phase_image == (
Image *) NULL)
1623 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1624 "ImageSequenceRequired",
"`%s'",magnitude_image->filename);
1625 return((
Image *) NULL);
1627 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1628 fourier_image=(
Image *) NULL;
1630 (void) ThrowMagickException(exception,GetMagickModule(),
1631 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1632 magnitude_image->filename);
1635 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1636 magnitude_image->rows,MagickTrue,exception);
1637 if (fourier_image != (
Image *) NULL)
1644 is_gray=IsGrayImage(magnitude_image,exception);
1645 if (is_gray != MagickFalse)
1646 is_gray=IsGrayImage(phase_image,exception);
1647 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1648 #pragma omp parallel sections
1651 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1658 if (is_gray != MagickFalse)
1659 thread_status=InverseFourierTransformChannel(magnitude_image,
1660 phase_image,GrayChannels,modulus,fourier_image,exception);
1662 thread_status=InverseFourierTransformChannel(magnitude_image,
1663 phase_image,RedChannel,modulus,fourier_image,exception);
1664 if (thread_status == MagickFalse)
1665 status=thread_status;
1667 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1674 thread_status=MagickTrue;
1675 if (is_gray == MagickFalse)
1676 thread_status=InverseFourierTransformChannel(magnitude_image,
1677 phase_image,GreenChannel,modulus,fourier_image,exception);
1678 if (thread_status == MagickFalse)
1679 status=thread_status;
1681 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1688 thread_status=MagickTrue;
1689 if (is_gray == MagickFalse)
1690 thread_status=InverseFourierTransformChannel(magnitude_image,
1691 phase_image,BlueChannel,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 (magnitude_image->matte != MagickFalse)
1704 thread_status=InverseFourierTransformChannel(magnitude_image,
1705 phase_image,OpacityChannel,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 (magnitude_image->colorspace == CMYKColorspace)
1718 thread_status=InverseFourierTransformChannel(magnitude_image,
1719 phase_image,IndexChannel,modulus,fourier_image,exception);
1720 if (thread_status == MagickFalse)
1721 status=thread_status;
1724 if (status == MagickFalse)
1725 fourier_image=DestroyImage(fourier_image);
1730 return(fourier_image);