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 = { QuantumScale*GetPixelRed(Ai), QuantumScale*GetPixelGreen(Ai),
279 QuantumScale*GetPixelBlue(Ai), image->matte != MagickFalse ?
280 QuantumScale*GetPixelOpacity(Ai) : OpaqueOpacity, 0 },
281 ar = { QuantumScale*GetPixelRed(Ar), QuantumScale*GetPixelGreen(Ar),
282 QuantumScale*GetPixelBlue(Ar), image->matte != MagickFalse ?
283 QuantumScale*GetPixelOpacity(Ar) : OpaqueOpacity, 0 },
284 bi = { QuantumScale*GetPixelRed(Bi), QuantumScale*GetPixelGreen(Bi),
285 QuantumScale*GetPixelBlue(Bi), image->matte != MagickFalse ?
286 QuantumScale*GetPixelOpacity(Bi) : OpaqueOpacity, 0 },
287 br = { QuantumScale*GetPixelRed(Br), QuantumScale*GetPixelGreen(Br),
288 QuantumScale*GetPixelBlue(Br), image->matte != MagickFalse ?
289 QuantumScale*GetPixelOpacity(Br) : OpaqueOpacity, 0 },
295 case AddComplexOperator:
297 cr.red=ar.red+br.red;
298 ci.red=ai.red+bi.red;
299 cr.green=ar.green+br.green;
300 ci.green=ai.green+bi.green;
301 cr.blue=ar.blue+br.blue;
302 ci.blue=ai.blue+bi.blue;
303 cr.opacity=ar.opacity+br.opacity;
304 ci.opacity=ai.opacity+bi.opacity;
307 case ConjugateComplexOperator:
313 ci.green=(-ai.green);
316 cr.opacity=ar.opacity;
317 ci.opacity=(-ai.opacity);
320 case DivideComplexOperator:
325 gamma=PerceptibleReciprocal(br.red*br.red+bi.red*bi.red+snr);
326 cr.red=gamma*(ar.red*br.red+ai.red*bi.red);
327 ci.red=gamma*(ai.red*br.red-ar.red*bi.red);
328 gamma=PerceptibleReciprocal(br.green*br.green+bi.green*bi.green+snr);
329 cr.green=gamma*(ar.green*br.green+ai.green*bi.green);
330 ci.green=gamma*(ai.green*br.green-ar.green*bi.green);
331 gamma=PerceptibleReciprocal(br.blue*br.blue+bi.blue*bi.blue+snr);
332 cr.blue=gamma*(ar.blue*br.blue+ai.blue*bi.blue);
333 ci.blue=gamma*(ai.blue*br.blue-ar.blue*bi.blue);
334 gamma=PerceptibleReciprocal(br.opacity*br.opacity+bi.opacity*
336 cr.opacity=gamma*(ar.opacity*br.opacity+ai.opacity*bi.opacity);
337 ci.opacity=gamma*(ai.opacity*br.opacity-ar.opacity*bi.opacity);
340 case MagnitudePhaseComplexOperator:
342 cr.red=sqrt((
double) ar.red*ar.red+ai.red*ai.red);
343 ci.red=atan2((
double) ai.red,(
double) ar.red)/(2.0*MagickPI)+0.5;
344 cr.green=sqrt((
double) ar.green*ar.green+ai.green*ai.green);
345 ci.green=atan2((
double) ai.green,(
double) ar.green)/
347 cr.blue=sqrt((
double) ar.blue*ar.blue+ai.blue*ai.blue);
348 ci.blue=atan2((
double) ai.blue,(
double) ar.blue)/(2.0*MagickPI)+0.5;
349 cr.opacity=sqrt((
double) ar.opacity*ar.opacity+ai.opacity*ai.opacity);
350 ci.opacity=atan2((
double) ai.opacity,(
double) ar.opacity)/
354 case MultiplyComplexOperator:
356 cr.red=(ar.red*br.red-(double) ai.red*bi.red);
357 ci.red=(ai.red*br.red+(double) ar.red*bi.red);
358 cr.green=(ar.green*br.green-(double) ai.green*bi.green);
359 ci.green=(ai.green*br.green+(double) ar.green*bi.green);
360 cr.blue=(ar.blue*br.blue-(double) ai.blue*bi.blue);
361 ci.blue=(ai.blue*br.blue+(double) ar.blue*bi.blue);
362 cr.opacity=(ar.opacity*br.opacity-(double) ai.opacity*bi.opacity);
363 ci.opacity=(ai.opacity*br.opacity+(double) ar.opacity*bi.opacity);
366 case RealImaginaryComplexOperator:
368 cr.red=ar.red*cos(2.0*MagickPI*(ai.red-0.5));
369 ci.red=ar.red*sin(2.0*MagickPI*(ai.red-0.5));
370 cr.green=ar.green*cos(2.0*MagickPI*(ai.green-0.5));
371 ci.green=ar.green*sin(2.0*MagickPI*(ai.green-0.5));
372 cr.blue=ar.blue*cos(2.0*MagickPI*(ai.blue-0.5));
373 ci.blue=ar.blue*sin(2.0*MagickPI*(ai.blue-0.5));
374 cr.opacity=ar.opacity*cos(2.0*MagickPI*(ai.opacity-0.5));
375 ci.opacity=ar.opacity*sin(2.0*MagickPI*(ai.opacity-0.5));
378 case SubtractComplexOperator:
380 cr.red=ar.red-br.red;
381 ci.red=ai.red-bi.red;
382 cr.green=ar.green-br.green;
383 ci.green=ai.green-bi.green;
384 cr.blue=ar.blue-br.blue;
385 ci.blue=ai.blue-bi.blue;
386 cr.opacity=ar.opacity-br.opacity;
387 ci.opacity=ai.opacity-bi.opacity;
391 Cr->red=QuantumRange*cr.red;
392 Ci->red=QuantumRange*ci.red;
393 Cr->green=QuantumRange*cr.green;
394 Ci->green=QuantumRange*ci.green;
395 Cr->blue=QuantumRange*cr.blue;
396 Ci->blue=QuantumRange*ci.blue;
397 if (images->matte != MagickFalse)
399 Cr->opacity=QuantumRange*cr.opacity;
400 Ci->opacity=QuantumRange*ci.opacity;
409 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
411 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
413 if (images->progress_monitor != (MagickProgressMonitor) NULL)
418 #if defined(MAGICKCORE_OPENMP_SUPPORT)
422 proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
423 if (proceed == MagickFalse)
427 Cr_view=DestroyCacheView(Cr_view);
428 Ci_view=DestroyCacheView(Ci_view);
429 Br_view=DestroyCacheView(Br_view);
430 Bi_view=DestroyCacheView(Bi_view);
431 Ar_view=DestroyCacheView(Ar_view);
432 Ai_view=DestroyCacheView(Ai_view);
433 if (status == MagickFalse)
434 complex_images=DestroyImageList(complex_images);
435 return(complex_images);
469 #if defined(MAGICKCORE_FFTW_DELEGATE)
471 static MagickBooleanType RollFourier(
const size_t width,
const size_t height,
472 const ssize_t x_offset,
const ssize_t y_offset,
double *roll_pixels)
492 source_info=AcquireVirtualMemory(width,height*
sizeof(*source_pixels));
495 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
497 for (y=0L; y < (ssize_t) height; y++)
500 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
502 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
504 for (x=0L; x < (ssize_t) width; x++)
507 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
509 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
511 source_pixels[v*width+u]=roll_pixels[i++];
514 (void) memcpy(roll_pixels,source_pixels,height*width*
sizeof(*source_pixels));
515 source_info=RelinquishVirtualMemory(source_info);
519 static MagickBooleanType ForwardQuadrantSwap(
const size_t width,
520 const size_t height,
double *source_pixels,
double *forward_pixels)
535 center=(ssize_t) (width/2L)+1L;
536 status=RollFourier((
size_t) center,height,0L,(ssize_t) height/2L,
538 if (status == MagickFalse)
540 for (y=0L; y < (ssize_t) height; y++)
541 for (x=0L; x < (ssize_t) (width/2L); x++)
542 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
543 for (y=1; y < (ssize_t) height; y++)
544 for (x=0L; x < (ssize_t) (width/2L); x++)
545 forward_pixels[(height-y)*width+width/2L-x-1L]=
546 source_pixels[y*center+x+1L];
547 for (x=0L; x < (ssize_t) (width/2L); x++)
548 forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
552 static void CorrectPhaseLHS(
const size_t width,
const size_t height,
553 double *fourier_pixels)
561 for (y=0L; y < (ssize_t) height; y++)
562 for (x=0L; x < (ssize_t) (width/2L); x++)
563 fourier_pixels[y*width+x]*=(-1.0);
566 static MagickBooleanType ForwardFourier(
const FourierInfo *fourier_info,
601 magnitude_image=GetFirstImageInList(image);
602 phase_image=GetNextImageInList(image);
603 if (phase_image == (
Image *) NULL)
605 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
606 "ImageSequenceRequired",
"`%s'",image->filename);
612 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
613 sizeof(*magnitude_pixels));
614 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
615 sizeof(*phase_pixels));
616 if ((magnitude_info == (
MemoryInfo *) NULL) ||
620 phase_info=RelinquishVirtualMemory(phase_info);
622 magnitude_info=RelinquishVirtualMemory(magnitude_info);
623 (void) ThrowMagickException(exception,GetMagickModule(),
624 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
627 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
628 (void) memset(magnitude_pixels,0,fourier_info->width*
629 fourier_info->height*
sizeof(*magnitude_pixels));
630 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
631 (void) memset(phase_pixels,0,fourier_info->width*
632 fourier_info->height*
sizeof(*phase_pixels));
633 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
634 magnitude,magnitude_pixels);
635 if (status != MagickFalse)
636 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
638 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
639 if (fourier_info->modulus != MagickFalse)
642 for (y=0L; y < (ssize_t) fourier_info->height; y++)
643 for (x=0L; x < (ssize_t) fourier_info->width; x++)
645 phase_pixels[i]/=(2.0*MagickPI);
646 phase_pixels[i]+=0.5;
650 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
652 for (y=0L; y < (ssize_t) fourier_info->height; y++)
654 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
658 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
659 for (x=0L; x < (ssize_t) fourier_info->width; x++)
661 switch (fourier_info->channel)
666 SetPixelRed(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
671 SetPixelGreen(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
676 SetPixelBlue(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
681 SetPixelOpacity(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
686 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*
687 magnitude_pixels[i]));
692 SetPixelGray(q,ClampToQuantum(QuantumRange*magnitude_pixels[i]));
699 status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
700 if (status == MagickFalse)
703 magnitude_view=DestroyCacheView(magnitude_view);
705 phase_view=AcquireAuthenticCacheView(phase_image,exception);
706 for (y=0L; y < (ssize_t) fourier_info->height; y++)
708 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
712 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
713 for (x=0L; x < (ssize_t) fourier_info->width; x++)
715 switch (fourier_info->channel)
720 SetPixelRed(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
725 SetPixelGreen(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
730 SetPixelBlue(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
735 SetPixelOpacity(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
740 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*phase_pixels[i]));
745 SetPixelGray(q,ClampToQuantum(QuantumRange*phase_pixels[i]));
752 status=SyncCacheViewAuthenticPixels(phase_view,exception);
753 if (status == MagickFalse)
756 phase_view=DestroyCacheView(phase_view);
757 phase_info=RelinquishVirtualMemory(phase_info);
758 magnitude_info=RelinquishVirtualMemory(magnitude_info);
762 static MagickBooleanType ForwardFourierTransform(
FourierInfo *fourier_info,
763 const Image *image,
double *magnitude_pixels,
double *phase_pixels,
802 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
804 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
805 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
808 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
809 sizeof(*source_pixels));
812 (void) ThrowMagickException(exception,GetMagickModule(),
813 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
816 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
817 memset(source_pixels,0,fourier_info->width*fourier_info->height*
818 sizeof(*source_pixels));
820 image_view=AcquireVirtualCacheView(image,exception);
821 for (y=0L; y < (ssize_t) fourier_info->height; y++)
823 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
827 indexes=GetCacheViewVirtualIndexQueue(image_view);
828 for (x=0L; x < (ssize_t) fourier_info->width; x++)
830 switch (fourier_info->channel)
835 source_pixels[i]=QuantumScale*GetPixelRed(p);
840 source_pixels[i]=QuantumScale*GetPixelGreen(p);
845 source_pixels[i]=QuantumScale*GetPixelBlue(p);
850 source_pixels[i]=QuantumScale*GetPixelOpacity(p);
855 source_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
860 source_pixels[i]=QuantumScale*GetPixelGray(p);
868 image_view=DestroyCacheView(image_view);
869 forward_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/2+
870 1)*
sizeof(*forward_pixels));
873 (void) ThrowMagickException(exception,GetMagickModule(),
874 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
875 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
878 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
879 #if defined(MAGICKCORE_OPENMP_SUPPORT)
880 #pragma omp critical (MagickCore_ForwardFourierTransform)
882 fftw_r2c_plan=fftw_plan_dft_r2c_2d((
int) fourier_info->width,
883 (
int) fourier_info->height,source_pixels,forward_pixels,FFTW_ESTIMATE);
884 fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
885 fftw_destroy_plan(fftw_r2c_plan);
886 source_info=(
MemoryInfo *) RelinquishVirtualMemory(source_info);
887 value=GetImageArtifact(image,
"fourier:normalize");
888 if ((value == (
const char *) NULL) || (LocaleCompare(value,
"forward") == 0))
897 gamma=PerceptibleReciprocal((
double) fourier_info->width*
898 fourier_info->height);
899 for (y=0L; y < (ssize_t) fourier_info->height; y++)
900 for (x=0L; x < (ssize_t) fourier_info->center; x++)
902 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
903 forward_pixels[i]*=gamma;
905 forward_pixels[i][0]*=gamma;
906 forward_pixels[i][1]*=gamma;
915 if (fourier_info->modulus != MagickFalse)
916 for (y=0L; y < (ssize_t) fourier_info->height; y++)
917 for (x=0L; x < (ssize_t) fourier_info->center; x++)
919 magnitude_pixels[i]=cabs(forward_pixels[i]);
920 phase_pixels[i]=carg(forward_pixels[i]);
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]=creal(forward_pixels[i]);
928 phase_pixels[i]=cimag(forward_pixels[i]);
931 forward_info=(
MemoryInfo *) RelinquishVirtualMemory(forward_info);
935 static MagickBooleanType ForwardFourierTransformChannel(
const Image *image,
936 const ChannelType channel,
const MagickBooleanType modulus,
953 fourier_info.width=image->columns;
954 fourier_info.height=image->rows;
955 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
956 ((image->rows % 2) != 0))
958 size_t extent=image->columns < image->rows ? image->rows : image->columns;
959 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
961 fourier_info.height=fourier_info.width;
962 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
963 fourier_info.channel=channel;
964 fourier_info.modulus=modulus;
965 magnitude_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
966 1)*
sizeof(*magnitude_pixels));
967 phase_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+1)*
968 sizeof(*phase_pixels));
969 if ((magnitude_info == (
MemoryInfo *) NULL) ||
973 phase_info=RelinquishVirtualMemory(phase_info);
975 magnitude_info=RelinquishVirtualMemory(magnitude_info);
976 (void) ThrowMagickException(exception,GetMagickModule(),
977 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
980 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
981 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
982 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
983 phase_pixels,exception);
984 if (status != MagickFalse)
985 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
986 phase_pixels,exception);
987 phase_info=RelinquishVirtualMemory(phase_info);
988 magnitude_info=RelinquishVirtualMemory(magnitude_info);
993 MagickExport
Image *ForwardFourierTransformImage(
const Image *image,
999 fourier_image=NewImageList();
1000 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1002 (void) ThrowMagickException(exception,GetMagickModule(),
1003 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1014 width=image->columns;
1016 if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
1017 ((image->rows % 2) != 0))
1019 size_t extent=image->columns < image->rows ? image->rows :
1021 width=(extent & 0x01) == 1 ? extent+1UL : extent;
1024 magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
1025 if (magnitude_image != (
Image *) NULL)
1030 magnitude_image->storage_class=DirectClass;
1031 magnitude_image->depth=32UL;
1032 phase_image=CloneImage(image,width,height,MagickTrue,exception);
1033 if (phase_image == (
Image *) NULL)
1034 magnitude_image=DestroyImage(magnitude_image);
1041 phase_image->storage_class=DirectClass;
1042 phase_image->depth=32UL;
1043 AppendImageToList(&fourier_image,magnitude_image);
1044 AppendImageToList(&fourier_image,phase_image);
1046 is_gray=IsGrayImage(image,exception);
1047 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1048 #pragma omp parallel sections
1051 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1058 if (is_gray != MagickFalse)
1059 thread_status=ForwardFourierTransformChannel(image,
1060 GrayChannels,modulus,fourier_image,exception);
1062 thread_status=ForwardFourierTransformChannel(image,RedChannel,
1063 modulus,fourier_image,exception);
1064 if (thread_status == MagickFalse)
1065 status=thread_status;
1067 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1074 thread_status=MagickTrue;
1075 if (is_gray == MagickFalse)
1076 thread_status=ForwardFourierTransformChannel(image,
1077 GreenChannel,modulus,fourier_image,exception);
1078 if (thread_status == MagickFalse)
1079 status=thread_status;
1081 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1088 thread_status=MagickTrue;
1089 if (is_gray == MagickFalse)
1090 thread_status=ForwardFourierTransformChannel(image,
1091 BlueChannel,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 (image->matte != MagickFalse)
1104 thread_status=ForwardFourierTransformChannel(image,
1105 OpacityChannel,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 (image->colorspace == CMYKColorspace)
1118 thread_status=ForwardFourierTransformChannel(image,
1119 IndexChannel,modulus,fourier_image,exception);
1120 if (thread_status == MagickFalse)
1121 status=thread_status;
1124 if (status == MagickFalse)
1125 fourier_image=DestroyImageList(fourier_image);
1131 return(fourier_image);
1168 #if defined(MAGICKCORE_FFTW_DELEGATE)
1169 static MagickBooleanType InverseQuadrantSwap(
const size_t width,
1170 const size_t height,
const double *source,
double *destination)
1182 center=(ssize_t) (width/2L)+1L;
1183 for (y=1L; y < (ssize_t) height; y++)
1184 for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1185 destination[(height-y)*center-x+width/2L]=source[y*width+x];
1186 for (y=0L; y < (ssize_t) height; y++)
1187 destination[y*center]=source[y*width+width/2L];
1188 for (x=0L; x < center; x++)
1189 destination[x]=source[center-x-1L];
1190 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1193 static MagickBooleanType InverseFourier(
FourierInfo *fourier_info,
1194 const Image *magnitude_image,
const Image *phase_image,
1230 magnitude_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1231 sizeof(*magnitude_pixels));
1232 phase_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1233 sizeof(*phase_pixels));
1234 inverse_info=AcquireVirtualMemory(fourier_info->width,(fourier_info->height/
1235 2+1)*
sizeof(*inverse_pixels));
1236 if ((magnitude_info == (
MemoryInfo *) NULL) ||
1241 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1243 phase_info=RelinquishVirtualMemory(phase_info);
1245 inverse_info=RelinquishVirtualMemory(inverse_info);
1246 (void) ThrowMagickException(exception,GetMagickModule(),
1247 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1248 magnitude_image->filename);
1249 return(MagickFalse);
1251 magnitude_pixels=(
double *) GetVirtualMemoryBlob(magnitude_info);
1252 phase_pixels=(
double *) GetVirtualMemoryBlob(phase_info);
1253 inverse_pixels=(
double *) GetVirtualMemoryBlob(inverse_info);
1255 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1256 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1258 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1262 indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
1263 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1265 switch (fourier_info->channel)
1270 magnitude_pixels[i]=QuantumScale*GetPixelRed(p);
1275 magnitude_pixels[i]=QuantumScale*GetPixelGreen(p);
1280 magnitude_pixels[i]=QuantumScale*GetPixelBlue(p);
1283 case OpacityChannel:
1285 magnitude_pixels[i]=QuantumScale*GetPixelOpacity(p);
1290 magnitude_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
1295 magnitude_pixels[i]=QuantumScale*GetPixelGray(p);
1303 magnitude_view=DestroyCacheView(magnitude_view);
1304 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1305 magnitude_pixels,inverse_pixels);
1306 (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1307 fourier_info->center*
sizeof(*magnitude_pixels));
1309 phase_view=AcquireVirtualCacheView(phase_image,exception);
1310 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1312 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1316 indexes=GetCacheViewAuthenticIndexQueue(phase_view);
1317 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1319 switch (fourier_info->channel)
1324 phase_pixels[i]=QuantumScale*GetPixelRed(p);
1329 phase_pixels[i]=QuantumScale*GetPixelGreen(p);
1334 phase_pixels[i]=QuantumScale*GetPixelBlue(p);
1337 case OpacityChannel:
1339 phase_pixels[i]=QuantumScale*GetPixelOpacity(p);
1344 phase_pixels[i]=QuantumScale*GetPixelIndex(indexes+x);
1349 phase_pixels[i]=QuantumScale*GetPixelGray(p);
1357 if (fourier_info->modulus != MagickFalse)
1360 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1361 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1363 phase_pixels[i]-=0.5;
1364 phase_pixels[i]*=(2.0*MagickPI);
1368 phase_view=DestroyCacheView(phase_view);
1369 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1370 if (status != MagickFalse)
1371 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1372 phase_pixels,inverse_pixels);
1373 (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1374 fourier_info->center*
sizeof(*phase_pixels));
1375 inverse_info=RelinquishVirtualMemory(inverse_info);
1380 if (fourier_info->modulus != MagickFalse)
1381 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1382 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1384 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1385 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1386 magnitude_pixels[i]*sin(phase_pixels[i]);
1388 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1389 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1394 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1395 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1397 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1398 fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1400 fourier_pixels[i][0]=magnitude_pixels[i];
1401 fourier_pixels[i][1]=phase_pixels[i];
1405 magnitude_info=RelinquishVirtualMemory(magnitude_info);
1406 phase_info=RelinquishVirtualMemory(phase_info);
1410 static MagickBooleanType InverseFourierTransform(
FourierInfo *fourier_info,
1444 if ((fourier_info->width >= INT_MAX) || (fourier_info->height >= INT_MAX))
1446 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1447 "WidthOrHeightExceedsLimit",
"`%s'",image->filename);
1448 return(MagickFalse);
1450 source_info=AcquireVirtualMemory(fourier_info->width,fourier_info->height*
1451 sizeof(*source_pixels));
1454 (void) ThrowMagickException(exception,GetMagickModule(),
1455 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1456 return(MagickFalse);
1458 source_pixels=(
double *) GetVirtualMemoryBlob(source_info);
1459 value=GetImageArtifact(image,
"fourier:normalize");
1460 if (LocaleCompare(value,
"inverse") == 0)
1469 gamma=PerceptibleReciprocal((
double) fourier_info->width*
1470 fourier_info->height);
1471 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1472 for (x=0L; x < (ssize_t) fourier_info->center; x++)
1474 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1475 fourier_pixels[i]*=gamma;
1477 fourier_pixels[i][0]*=gamma;
1478 fourier_pixels[i][1]*=gamma;
1483 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1484 #pragma omp critical (MagickCore_InverseFourierTransform)
1486 fftw_c2r_plan=fftw_plan_dft_c2r_2d((
int) fourier_info->width,
1487 (
int) fourier_info->height,fourier_pixels,source_pixels,FFTW_ESTIMATE);
1488 fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1489 fftw_destroy_plan(fftw_c2r_plan);
1491 image_view=AcquireAuthenticCacheView(image,exception);
1492 for (y=0L; y < (ssize_t) fourier_info->height; y++)
1494 if (y >= (ssize_t) image->rows)
1496 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1497 image->columns ? image->columns : fourier_info->width,1UL,exception);
1500 indexes=GetCacheViewAuthenticIndexQueue(image_view);
1501 for (x=0L; x < (ssize_t) fourier_info->width; x++)
1503 if (x < (ssize_t) image->columns)
1504 switch (fourier_info->channel)
1509 SetPixelRed(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1514 SetPixelGreen(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1519 SetPixelBlue(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1522 case OpacityChannel:
1524 SetPixelOpacity(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1529 SetPixelIndex(indexes+x,ClampToQuantum(QuantumRange*
1535 SetPixelGray(q,ClampToQuantum(QuantumRange*source_pixels[i]));
1542 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1545 image_view=DestroyCacheView(image_view);
1546 source_info=RelinquishVirtualMemory(source_info);
1550 static MagickBooleanType InverseFourierTransformChannel(
1551 const Image *magnitude_image,
const Image *phase_image,
1552 const ChannelType channel,
const MagickBooleanType modulus,
1567 fourier_info.width=magnitude_image->columns;
1568 fourier_info.height=magnitude_image->rows;
1569 if ((magnitude_image->columns != magnitude_image->rows) ||
1570 ((magnitude_image->columns % 2) != 0) ||
1571 ((magnitude_image->rows % 2) != 0))
1573 size_t extent=magnitude_image->columns < magnitude_image->rows ?
1574 magnitude_image->rows : magnitude_image->columns;
1575 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1577 fourier_info.height=fourier_info.width;
1578 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1579 fourier_info.channel=channel;
1580 fourier_info.modulus=modulus;
1581 inverse_info=AcquireVirtualMemory(fourier_info.width,(fourier_info.height/2+
1582 1)*
sizeof(*inverse_pixels));
1585 (void) ThrowMagickException(exception,GetMagickModule(),
1586 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",
1587 magnitude_image->filename);
1588 return(MagickFalse);
1590 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1591 status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1592 inverse_pixels,exception);
1593 if (status != MagickFalse)
1594 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1596 inverse_info=RelinquishVirtualMemory(inverse_info);
1601 MagickExport
Image *InverseFourierTransformImage(
const Image *magnitude_image,
1602 const Image *phase_image,
const MagickBooleanType modulus,
1608 assert(magnitude_image != (
Image *) NULL);
1609 assert(magnitude_image->signature == MagickCoreSignature);
1610 if (IsEventLogging() != MagickFalse)
1611 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1612 magnitude_image->filename);
1613 if (phase_image == (
Image *) NULL)
1615 (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1616 "ImageSequenceRequired",
"`%s'",magnitude_image->filename);
1617 return((
Image *) NULL);
1619 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1620 fourier_image=(
Image *) NULL;
1622 (void) ThrowMagickException(exception,GetMagickModule(),
1623 MissingDelegateWarning,
"DelegateLibrarySupportNotBuiltIn",
"`%s' (FFTW)",
1624 magnitude_image->filename);
1627 fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1628 magnitude_image->rows,MagickTrue,exception);
1629 if (fourier_image != (
Image *) NULL)
1636 is_gray=IsGrayImage(magnitude_image,exception);
1637 if (is_gray != MagickFalse)
1638 is_gray=IsGrayImage(phase_image,exception);
1639 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1640 #pragma omp parallel sections
1643 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1650 if (is_gray != MagickFalse)
1651 thread_status=InverseFourierTransformChannel(magnitude_image,
1652 phase_image,GrayChannels,modulus,fourier_image,exception);
1654 thread_status=InverseFourierTransformChannel(magnitude_image,
1655 phase_image,RedChannel,modulus,fourier_image,exception);
1656 if (thread_status == MagickFalse)
1657 status=thread_status;
1659 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1666 thread_status=MagickTrue;
1667 if (is_gray == MagickFalse)
1668 thread_status=InverseFourierTransformChannel(magnitude_image,
1669 phase_image,GreenChannel,modulus,fourier_image,exception);
1670 if (thread_status == MagickFalse)
1671 status=thread_status;
1673 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1680 thread_status=MagickTrue;
1681 if (is_gray == MagickFalse)
1682 thread_status=InverseFourierTransformChannel(magnitude_image,
1683 phase_image,BlueChannel,modulus,fourier_image,exception);
1684 if (thread_status == MagickFalse)
1685 status=thread_status;
1687 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1694 thread_status=MagickTrue;
1695 if (magnitude_image->matte != MagickFalse)
1696 thread_status=InverseFourierTransformChannel(magnitude_image,
1697 phase_image,OpacityChannel,modulus,fourier_image,exception);
1698 if (thread_status == MagickFalse)
1699 status=thread_status;
1701 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1708 thread_status=MagickTrue;
1709 if (magnitude_image->colorspace == CMYKColorspace)
1710 thread_status=InverseFourierTransformChannel(magnitude_image,
1711 phase_image,IndexChannel,modulus,fourier_image,exception);
1712 if (thread_status == MagickFalse)
1713 status=thread_status;
1716 if (status == MagickFalse)
1717 fourier_image=DestroyImage(fourier_image);
1722 return(fourier_image);