MagickCore  7.1.0
fourier.c
Go to the documentation of this file.
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % FFFFF OOO U U RRRR IIIII EEEEE RRRR %
7 % F O O U U R R I E R R %
8 % FFF O O U U RRRR I EEE RRRR %
9 % F O O U U R R I E R R %
10 % F OOO UUU R R IIIII EEEEE R R %
11 % %
12 % %
13 % MagickCore Discrete Fourier Transform Methods %
14 % %
15 % Software Design %
16 % Sean Burke %
17 % Fred Weinhaus %
18 % Cristy %
19 % July 2009 %
20 % %
21 % %
22 % Copyright 1999-2021 ImageMagick Studio LLC, a non-profit organization %
23 % dedicated to making software imaging solutions freely available. %
24 % %
25 % You may not use this file except in compliance with the License. You may %
26 % obtain a copy of the License at %
27 % %
28 % https://imagemagick.org/script/license.php %
29 % %
30 % Unless required by applicable law or agreed to in writing, software %
31 % distributed under the License is distributed on an "AS IS" BASIS, %
32 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
33 % See the License for the specific language governing permissions and %
34 % limitations under the License. %
35 % %
36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 %
38 %
39 %
40 */
41 
42 /*
43  Include declarations.
44 */
45 #include "MagickCore/studio.h"
46 #include "MagickCore/artifact.h"
47 #include "MagickCore/attribute.h"
48 #include "MagickCore/blob.h"
49 #include "MagickCore/cache.h"
50 #include "MagickCore/image.h"
52 #include "MagickCore/list.h"
53 #include "MagickCore/fourier.h"
54 #include "MagickCore/log.h"
55 #include "MagickCore/memory_.h"
56 #include "MagickCore/monitor.h"
59 #include "MagickCore/property.h"
61 #include "MagickCore/resource_.h"
64 #if defined(MAGICKCORE_FFTW_DELEGATE)
65 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
66 #include <complex.h>
67 #endif
68 #include <fftw3.h>
69 #if !defined(MAGICKCORE_HAVE_CABS)
70 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
71 #endif
72 #if !defined(MAGICKCORE_HAVE_CARG)
73 #define carg(z) (atan2(cimag(z),creal(z)))
74 #endif
75 #if !defined(MAGICKCORE_HAVE_CIMAG)
76 #define cimag(z) (z[1])
77 #endif
78 #if !defined(MAGICKCORE_HAVE_CREAL)
79 #define creal(z) (z[0])
80 #endif
81 #endif
82 
83 /*
84  Typedef declarations.
85 */
86 typedef struct _FourierInfo
87 {
90 
93 
94  size_t
96  height;
97 
98  ssize_t
100 } FourierInfo;
101 
102 /*
103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 % %
105 % %
106 % %
107 % C o m p l e x I m a g e s %
108 % %
109 % %
110 % %
111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112 %
113 % ComplexImages() performs complex mathematics on an image sequence.
114 %
115 % The format of the ComplexImages method is:
116 %
117 % MagickBooleanType ComplexImages(Image *images,const ComplexOperator op,
118 % ExceptionInfo *exception)
119 %
120 % A description of each parameter follows:
121 %
122 % o image: the image.
123 %
124 % o op: A complex operator.
125 %
126 % o exception: return any errors or warnings in this structure.
127 %
128 */
130  ExceptionInfo *exception)
131 {
132 #define ComplexImageTag "Complex/Image"
133 
134  CacheView
135  *Ai_view,
136  *Ar_view,
137  *Bi_view,
138  *Br_view,
139  *Ci_view,
140  *Cr_view;
141 
142  const char
143  *artifact;
144 
145  const Image
146  *Ai_image,
147  *Ar_image,
148  *Bi_image,
149  *Br_image;
150 
151  double
152  snr;
153 
154  Image
155  *Ci_image,
156  *complex_images,
157  *Cr_image,
158  *image;
159 
161  status;
162 
164  progress;
165 
166  size_t
167  columns,
168  number_channels,
169  rows;
170 
171  ssize_t
172  y;
173 
174  assert(images != (Image *) NULL);
175  assert(images->signature == MagickCoreSignature);
176  if (images->debug != MagickFalse)
177  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
178  assert(exception != (ExceptionInfo *) NULL);
179  assert(exception->signature == MagickCoreSignature);
180  if (images->next == (Image *) NULL)
181  {
183  "ImageSequenceRequired","`%s'",images->filename);
184  return((Image *) NULL);
185  }
186  image=CloneImage(images,0,0,MagickTrue,exception);
187  if (image == (Image *) NULL)
188  return((Image *) NULL);
189  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
190  {
191  image=DestroyImageList(image);
192  return(image);
193  }
194  image->depth=32UL;
195  complex_images=NewImageList();
196  AppendImageToList(&complex_images,image);
197  image=CloneImage(images->next,0,0,MagickTrue,exception);
198  if (image == (Image *) NULL)
199  {
200  complex_images=DestroyImageList(complex_images);
201  return(complex_images);
202  }
203  image->depth=32UL;
204  AppendImageToList(&complex_images,image);
205  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
206  {
207  complex_images=DestroyImageList(complex_images);
208  return(complex_images);
209  }
210  /*
211  Apply complex mathematics to image pixels.
212  */
213  artifact=GetImageArtifact(images,"complex:snr");
214  snr=0.0;
215  if (artifact != (const char *) NULL)
216  snr=StringToDouble(artifact,(char **) NULL);
217  Ar_image=images;
218  Ai_image=images->next;
219  Br_image=images;
220  Bi_image=images->next;
221  if ((images->next->next != (Image *) NULL) &&
222  (images->next->next->next != (Image *) NULL))
223  {
224  Br_image=images->next->next;
225  Bi_image=images->next->next->next;
226  }
227  Cr_image=complex_images;
228  Ci_image=complex_images->next;
229  number_channels=MagickMin(MagickMin(MagickMin(
230  Ar_image->number_channels,Ai_image->number_channels),MagickMin(
231  Br_image->number_channels,Bi_image->number_channels)),MagickMin(
232  Cr_image->number_channels,Ci_image->number_channels));
233  Ar_view=AcquireVirtualCacheView(Ar_image,exception);
234  Ai_view=AcquireVirtualCacheView(Ai_image,exception);
235  Br_view=AcquireVirtualCacheView(Br_image,exception);
236  Bi_view=AcquireVirtualCacheView(Bi_image,exception);
237  Cr_view=AcquireAuthenticCacheView(Cr_image,exception);
238  Ci_view=AcquireAuthenticCacheView(Ci_image,exception);
239  status=MagickTrue;
240  progress=0;
241  columns=MagickMin(Cr_image->columns,Ci_image->columns);
242  rows=MagickMin(Cr_image->rows,Ci_image->rows);
243 #if defined(MAGICKCORE_OPENMP_SUPPORT)
244  #pragma omp parallel for schedule(static) shared(progress,status) \
245  magick_number_threads(Cr_image,complex_images,rows,1L)
246 #endif
247  for (y=0; y < (ssize_t) rows; y++)
248  {
249  const Quantum
250  *magick_restrict Ai,
251  *magick_restrict Ar,
252  *magick_restrict Bi,
253  *magick_restrict Br;
254 
255  Quantum
256  *magick_restrict Ci,
257  *magick_restrict Cr;
258 
259  ssize_t
260  x;
261 
262  if (status == MagickFalse)
263  continue;
264  Ar=GetCacheViewVirtualPixels(Ar_view,0,y,columns,1,exception);
265  Ai=GetCacheViewVirtualPixels(Ai_view,0,y,columns,1,exception);
266  Br=GetCacheViewVirtualPixels(Br_view,0,y,columns,1,exception);
267  Bi=GetCacheViewVirtualPixels(Bi_view,0,y,columns,1,exception);
268  Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,columns,1,exception);
269  Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,columns,1,exception);
270  if ((Ar == (const Quantum *) NULL) || (Ai == (const Quantum *) NULL) ||
271  (Br == (const Quantum *) NULL) || (Bi == (const Quantum *) NULL) ||
272  (Cr == (Quantum *) NULL) || (Ci == (Quantum *) NULL))
273  {
274  status=MagickFalse;
275  continue;
276  }
277  for (x=0; x < (ssize_t) columns; x++)
278  {
279  ssize_t
280  i;
281 
282  for (i=0; i < (ssize_t) number_channels; i++)
283  {
284  switch (op)
285  {
286  case AddComplexOperator:
287  {
288  Cr[i]=Ar[i]+Br[i];
289  Ci[i]=Ai[i]+Bi[i];
290  break;
291  }
293  default:
294  {
295  Cr[i]=Ar[i];
296  Ci[i]=(-Ai[i]);
297  break;
298  }
300  {
301  double
302  gamma;
303 
305  QuantumScale*Bi[i]*Bi[i]+snr);
306  Cr[i]=gamma*(QuantumScale*Ar[i]*Br[i]+QuantumScale*Ai[i]*Bi[i]);
307  Ci[i]=gamma*(QuantumScale*Ai[i]*Br[i]-QuantumScale*Ar[i]*Bi[i]);
308  break;
309  }
311  {
312  Cr[i]=sqrt(QuantumScale*Ar[i]*Ar[i]+QuantumScale*Ai[i]*Ai[i]);
313  Ci[i]=atan2((double) Ai[i],(double) Ar[i])/(2.0*MagickPI)+0.5;
314  break;
315  }
317  {
318  Cr[i]=(QuantumScale*Ar[i]*Br[i]-QuantumScale*Ai[i]*Bi[i]);
319  Ci[i]=(QuantumScale*Ai[i]*Br[i]+QuantumScale*Ar[i]*Bi[i]);
320  break;
321  }
323  {
324  Cr[i]=Ar[i]*cos(2.0*MagickPI*(Ai[i]-0.5));
325  Ci[i]=Ar[i]*sin(2.0*MagickPI*(Ai[i]-0.5));
326  break;
327  }
329  {
330  Cr[i]=Ar[i]-Br[i];
331  Ci[i]=Ai[i]-Bi[i];
332  break;
333  }
334  }
335  }
336  Ar+=GetPixelChannels(Ar_image);
337  Ai+=GetPixelChannels(Ai_image);
338  Br+=GetPixelChannels(Br_image);
339  Bi+=GetPixelChannels(Bi_image);
340  Cr+=GetPixelChannels(Cr_image);
341  Ci+=GetPixelChannels(Ci_image);
342  }
343  if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
344  status=MagickFalse;
345  if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
346  status=MagickFalse;
347  if (images->progress_monitor != (MagickProgressMonitor) NULL)
348  {
350  proceed;
351 
352 #if defined(MAGICKCORE_OPENMP_SUPPORT)
353  #pragma omp atomic
354 #endif
355  progress++;
356  proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
357  if (proceed == MagickFalse)
358  status=MagickFalse;
359  }
360  }
361  Cr_view=DestroyCacheView(Cr_view);
362  Ci_view=DestroyCacheView(Ci_view);
363  Br_view=DestroyCacheView(Br_view);
364  Bi_view=DestroyCacheView(Bi_view);
365  Ar_view=DestroyCacheView(Ar_view);
366  Ai_view=DestroyCacheView(Ai_view);
367  if (status == MagickFalse)
368  complex_images=DestroyImageList(complex_images);
369  return(complex_images);
370 }
371 
372 /*
373 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
374 % %
375 % %
376 % %
377 % F o r w a r d F o u r i e r T r a n s f o r m I m a g e %
378 % %
379 % %
380 % %
381 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
382 %
383 % ForwardFourierTransformImage() implements the discrete Fourier transform
384 % (DFT) of the image either as a magnitude / phase or real / imaginary image
385 % pair.
386 %
387 % The format of the ForwadFourierTransformImage method is:
388 %
389 % Image *ForwardFourierTransformImage(const Image *image,
390 % const MagickBooleanType modulus,ExceptionInfo *exception)
391 %
392 % A description of each parameter follows:
393 %
394 % o image: the image.
395 %
396 % o modulus: if true, return as transform as a magnitude / phase pair
397 % otherwise a real / imaginary image pair.
398 %
399 % o exception: return any errors or warnings in this structure.
400 %
401 */
402 
403 #if defined(MAGICKCORE_FFTW_DELEGATE)
404 
405 static MagickBooleanType RollFourier(const size_t width,const size_t height,
406  const ssize_t x_offset,const ssize_t y_offset,double *roll_pixels)
407 {
408  double
409  *source_pixels;
410 
411  MemoryInfo
412  *source_info;
413 
414  ssize_t
415  i,
416  x;
417 
418  ssize_t
419  u,
420  v,
421  y;
422 
423  /*
424  Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
425  */
426  source_info=AcquireVirtualMemory(width,height*sizeof(*source_pixels));
427  if (source_info == (MemoryInfo *) NULL)
428  return(MagickFalse);
429  source_pixels=(double *) GetVirtualMemoryBlob(source_info);
430  i=0L;
431  for (y=0L; y < (ssize_t) height; y++)
432  {
433  if (y_offset < 0L)
434  v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
435  else
436  v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
437  y+y_offset;
438  for (x=0L; x < (ssize_t) width; x++)
439  {
440  if (x_offset < 0L)
441  u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
442  else
443  u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
444  x+x_offset;
445  source_pixels[v*width+u]=roll_pixels[i++];
446  }
447  }
448  (void) memcpy(roll_pixels,source_pixels,height*width*sizeof(*source_pixels));
449  source_info=RelinquishVirtualMemory(source_info);
450  return(MagickTrue);
451 }
452 
453 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
454  const size_t height,double *source_pixels,double *forward_pixels)
455 {
457  status;
458 
459  ssize_t
460  x;
461 
462  ssize_t
463  center,
464  y;
465 
466  /*
467  Swap quadrants.
468  */
469  center=(ssize_t) (width/2L)+1L;
470  status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,
471  source_pixels);
472  if (status == MagickFalse)
473  return(MagickFalse);
474  for (y=0L; y < (ssize_t) height; y++)
475  for (x=0L; x < (ssize_t) (width/2L); x++)
476  forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
477  for (y=1; y < (ssize_t) height; y++)
478  for (x=0L; x < (ssize_t) (width/2L); x++)
479  forward_pixels[(height-y)*width+width/2L-x-1L]=
480  source_pixels[y*center+x+1L];
481  for (x=0L; x < (ssize_t) (width/2L); x++)
482  forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
483  return(MagickTrue);
484 }
485 
486 static void CorrectPhaseLHS(const size_t width,const size_t height,
487  double *fourier_pixels)
488 {
489  ssize_t
490  x;
491 
492  ssize_t
493  y;
494 
495  for (y=0L; y < (ssize_t) height; y++)
496  for (x=0L; x < (ssize_t) (width/2L); x++)
497  fourier_pixels[y*width+x]*=(-1.0);
498 }
499 
500 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
501  Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
502 {
503  CacheView
504  *magnitude_view,
505  *phase_view;
506 
507  double
508  *magnitude_pixels,
509  *phase_pixels;
510 
511  Image
512  *magnitude_image,
513  *phase_image;
514 
516  status;
517 
518  MemoryInfo
519  *magnitude_info,
520  *phase_info;
521 
522  Quantum
523  *q;
524 
525  ssize_t
526  x;
527 
528  ssize_t
529  i,
530  y;
531 
532  magnitude_image=GetFirstImageInList(image);
533  phase_image=GetNextImageInList(image);
534  if (phase_image == (Image *) NULL)
535  {
537  "ImageSequenceRequired","`%s'",image->filename);
538  return(MagickFalse);
539  }
540  /*
541  Create "Fourier Transform" image from constituent arrays.
542  */
543  magnitude_info=AcquireVirtualMemory((size_t) fourier_info->width,
544  fourier_info->height*sizeof(*magnitude_pixels));
545  phase_info=AcquireVirtualMemory((size_t) fourier_info->width,
546  fourier_info->height*sizeof(*phase_pixels));
547  if ((magnitude_info == (MemoryInfo *) NULL) ||
548  (phase_info == (MemoryInfo *) NULL))
549  {
550  if (phase_info != (MemoryInfo *) NULL)
551  phase_info=RelinquishVirtualMemory(phase_info);
552  if (magnitude_info != (MemoryInfo *) NULL)
553  magnitude_info=RelinquishVirtualMemory(magnitude_info);
554  (void) ThrowMagickException(exception,GetMagickModule(),
555  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
556  return(MagickFalse);
557  }
558  magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
559  (void) memset(magnitude_pixels,0,fourier_info->width*
560  fourier_info->height*sizeof(*magnitude_pixels));
561  phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
562  (void) memset(phase_pixels,0,fourier_info->width*
563  fourier_info->height*sizeof(*phase_pixels));
564  status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
565  magnitude,magnitude_pixels);
566  if (status != MagickFalse)
567  status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
568  phase_pixels);
569  CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
570  if (fourier_info->modulus != MagickFalse)
571  {
572  i=0L;
573  for (y=0L; y < (ssize_t) fourier_info->height; y++)
574  for (x=0L; x < (ssize_t) fourier_info->width; x++)
575  {
576  phase_pixels[i]/=(2.0*MagickPI);
577  phase_pixels[i]+=0.5;
578  i++;
579  }
580  }
581  magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
582  i=0L;
583  for (y=0L; y < (ssize_t) fourier_info->height; y++)
584  {
585  q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
586  exception);
587  if (q == (Quantum *) NULL)
588  break;
589  for (x=0L; x < (ssize_t) fourier_info->width; x++)
590  {
591  switch (fourier_info->channel)
592  {
593  case RedPixelChannel:
594  default:
595  {
596  SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
597  magnitude_pixels[i]),q);
598  break;
599  }
600  case GreenPixelChannel:
601  {
602  SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
603  magnitude_pixels[i]),q);
604  break;
605  }
606  case BluePixelChannel:
607  {
608  SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
609  magnitude_pixels[i]),q);
610  break;
611  }
612  case BlackPixelChannel:
613  {
614  SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
615  magnitude_pixels[i]),q);
616  break;
617  }
618  case AlphaPixelChannel:
619  {
620  SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
621  magnitude_pixels[i]),q);
622  break;
623  }
624  }
625  i++;
626  q+=GetPixelChannels(magnitude_image);
627  }
628  status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
629  if (status == MagickFalse)
630  break;
631  }
632  magnitude_view=DestroyCacheView(magnitude_view);
633  i=0L;
634  phase_view=AcquireAuthenticCacheView(phase_image,exception);
635  for (y=0L; y < (ssize_t) fourier_info->height; y++)
636  {
637  q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
638  exception);
639  if (q == (Quantum *) NULL)
640  break;
641  for (x=0L; x < (ssize_t) fourier_info->width; x++)
642  {
643  switch (fourier_info->channel)
644  {
645  case RedPixelChannel:
646  default:
647  {
649  phase_pixels[i]),q);
650  break;
651  }
652  case GreenPixelChannel:
653  {
655  phase_pixels[i]),q);
656  break;
657  }
658  case BluePixelChannel:
659  {
661  phase_pixels[i]),q);
662  break;
663  }
664  case BlackPixelChannel:
665  {
667  phase_pixels[i]),q);
668  break;
669  }
670  case AlphaPixelChannel:
671  {
673  phase_pixels[i]),q);
674  break;
675  }
676  }
677  i++;
678  q+=GetPixelChannels(phase_image);
679  }
680  status=SyncCacheViewAuthenticPixels(phase_view,exception);
681  if (status == MagickFalse)
682  break;
683  }
684  phase_view=DestroyCacheView(phase_view);
685  phase_info=RelinquishVirtualMemory(phase_info);
686  magnitude_info=RelinquishVirtualMemory(magnitude_info);
687  return(status);
688 }
689 
690 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
691  const Image *image,double *magnitude_pixels,double *phase_pixels,
692  ExceptionInfo *exception)
693 {
694  CacheView
695  *image_view;
696 
697  const char
698  *value;
699 
700  double
701  *source_pixels;
702 
703  fftw_complex
704  *forward_pixels;
705 
706  fftw_plan
707  fftw_r2c_plan;
708 
709  MemoryInfo
710  *forward_info,
711  *source_info;
712 
713  const Quantum
714  *p;
715 
716  ssize_t
717  i,
718  x;
719 
720  ssize_t
721  y;
722 
723  /*
724  Generate the forward Fourier transform.
725  */
726  source_info=AcquireVirtualMemory((size_t) fourier_info->width,
727  fourier_info->height*sizeof(*source_pixels));
728  if (source_info == (MemoryInfo *) NULL)
729  {
730  (void) ThrowMagickException(exception,GetMagickModule(),
731  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
732  return(MagickFalse);
733  }
734  source_pixels=(double *) GetVirtualMemoryBlob(source_info);
735  memset(source_pixels,0,fourier_info->width*fourier_info->height*
736  sizeof(*source_pixels));
737  i=0L;
738  image_view=AcquireVirtualCacheView(image,exception);
739  for (y=0L; y < (ssize_t) fourier_info->height; y++)
740  {
741  p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
742  exception);
743  if (p == (const Quantum *) NULL)
744  break;
745  for (x=0L; x < (ssize_t) fourier_info->width; x++)
746  {
747  switch (fourier_info->channel)
748  {
749  case RedPixelChannel:
750  default:
751  {
752  source_pixels[i]=QuantumScale*GetPixelRed(image,p);
753  break;
754  }
755  case GreenPixelChannel:
756  {
757  source_pixels[i]=QuantumScale*GetPixelGreen(image,p);
758  break;
759  }
760  case BluePixelChannel:
761  {
762  source_pixels[i]=QuantumScale*GetPixelBlue(image,p);
763  break;
764  }
765  case BlackPixelChannel:
766  {
767  source_pixels[i]=QuantumScale*GetPixelBlack(image,p);
768  break;
769  }
770  case AlphaPixelChannel:
771  {
772  source_pixels[i]=QuantumScale*GetPixelAlpha(image,p);
773  break;
774  }
775  }
776  i++;
777  p+=GetPixelChannels(image);
778  }
779  }
780  image_view=DestroyCacheView(image_view);
781  forward_info=AcquireVirtualMemory((size_t) fourier_info->width,
782  (fourier_info->height/2+1)*sizeof(*forward_pixels));
783  if (forward_info == (MemoryInfo *) NULL)
784  {
785  (void) ThrowMagickException(exception,GetMagickModule(),
786  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
787  source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
788  return(MagickFalse);
789  }
790  forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
791 #if defined(MAGICKCORE_OPENMP_SUPPORT)
792  #pragma omp critical (MagickCore_ForwardFourierTransform)
793 #endif
794  fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->height,
795  source_pixels,forward_pixels,FFTW_ESTIMATE);
796  fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
797  fftw_destroy_plan(fftw_r2c_plan);
798  source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
799  value=GetImageArtifact(image,"fourier:normalize");
800  if ((value == (const char *) NULL) || (LocaleCompare(value,"forward") == 0))
801  {
802  double
803  gamma;
804 
805  /*
806  Normalize forward transform.
807  */
808  i=0L;
809  gamma=PerceptibleReciprocal((double) fourier_info->width*
810  fourier_info->height);
811  for (y=0L; y < (ssize_t) fourier_info->height; y++)
812  for (x=0L; x < (ssize_t) fourier_info->center; x++)
813  {
814 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
815  forward_pixels[i]*=gamma;
816 #else
817  forward_pixels[i][0]*=gamma;
818  forward_pixels[i][1]*=gamma;
819 #endif
820  i++;
821  }
822  }
823  /*
824  Generate magnitude and phase (or real and imaginary).
825  */
826  i=0L;
827  if (fourier_info->modulus != MagickFalse)
828  for (y=0L; y < (ssize_t) fourier_info->height; y++)
829  for (x=0L; x < (ssize_t) fourier_info->center; x++)
830  {
831  magnitude_pixels[i]=cabs(forward_pixels[i]);
832  phase_pixels[i]=carg(forward_pixels[i]);
833  i++;
834  }
835  else
836  for (y=0L; y < (ssize_t) fourier_info->height; y++)
837  for (x=0L; x < (ssize_t) fourier_info->center; x++)
838  {
839  magnitude_pixels[i]=creal(forward_pixels[i]);
840  phase_pixels[i]=cimag(forward_pixels[i]);
841  i++;
842  }
843  forward_info=(MemoryInfo *) RelinquishVirtualMemory(forward_info);
844  return(MagickTrue);
845 }
846 
847 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
848  const PixelChannel channel,const MagickBooleanType modulus,
849  Image *fourier_image,ExceptionInfo *exception)
850 {
851  double
852  *magnitude_pixels,
853  *phase_pixels;
854 
856  fourier_info;
857 
859  status;
860 
861  MemoryInfo
862  *magnitude_info,
863  *phase_info;
864 
865  fourier_info.width=image->columns;
866  fourier_info.height=image->rows;
867  if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
868  ((image->rows % 2) != 0))
869  {
870  size_t extent=image->columns < image->rows ? image->rows : image->columns;
871  fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
872  }
873  fourier_info.height=fourier_info.width;
874  fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
875  fourier_info.channel=channel;
876  fourier_info.modulus=modulus;
877  magnitude_info=AcquireVirtualMemory((size_t) fourier_info.width,
878  (fourier_info.height/2+1)*sizeof(*magnitude_pixels));
879  phase_info=AcquireVirtualMemory((size_t) fourier_info.width,
880  (fourier_info.height/2+1)*sizeof(*phase_pixels));
881  if ((magnitude_info == (MemoryInfo *) NULL) ||
882  (phase_info == (MemoryInfo *) NULL))
883  {
884  if (phase_info != (MemoryInfo *) NULL)
885  phase_info=RelinquishVirtualMemory(phase_info);
886  if (magnitude_info == (MemoryInfo *) NULL)
887  magnitude_info=RelinquishVirtualMemory(magnitude_info);
888  (void) ThrowMagickException(exception,GetMagickModule(),
889  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
890  return(MagickFalse);
891  }
892  magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
893  phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
894  status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
895  phase_pixels,exception);
896  if (status != MagickFalse)
897  status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
898  phase_pixels,exception);
899  phase_info=RelinquishVirtualMemory(phase_info);
900  magnitude_info=RelinquishVirtualMemory(magnitude_info);
901  return(status);
902 }
903 #endif
904 
906  const MagickBooleanType modulus,ExceptionInfo *exception)
907 {
908  Image
909  *fourier_image;
910 
911  fourier_image=NewImageList();
912 #if !defined(MAGICKCORE_FFTW_DELEGATE)
913  (void) modulus;
914  (void) ThrowMagickException(exception,GetMagickModule(),
915  MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
916  image->filename);
917 #else
918  {
919  Image
920  *magnitude_image;
921 
922  size_t
923  height,
924  width;
925 
926  width=image->columns;
927  height=image->rows;
928  if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
929  ((image->rows % 2) != 0))
930  {
931  size_t extent=image->columns < image->rows ? image->rows :
932  image->columns;
933  width=(extent & 0x01) == 1 ? extent+1UL : extent;
934  }
935  height=width;
936  magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
937  if (magnitude_image != (Image *) NULL)
938  {
939  Image
940  *phase_image;
941 
942  magnitude_image->storage_class=DirectClass;
943  magnitude_image->depth=32UL;
944  phase_image=CloneImage(image,width,height,MagickTrue,exception);
945  if (phase_image == (Image *) NULL)
946  magnitude_image=DestroyImage(magnitude_image);
947  else
948  {
950  is_gray,
951  status;
952 
953  phase_image->storage_class=DirectClass;
954  phase_image->depth=32UL;
955  AppendImageToList(&fourier_image,magnitude_image);
956  AppendImageToList(&fourier_image,phase_image);
957  status=MagickTrue;
958  is_gray=IsImageGray(image);
959 #if defined(MAGICKCORE_OPENMP_SUPPORT)
960  #pragma omp parallel sections
961 #endif
962  {
963 #if defined(MAGICKCORE_OPENMP_SUPPORT)
964  #pragma omp section
965 #endif
966  {
968  thread_status;
969 
970  if (is_gray != MagickFalse)
971  thread_status=ForwardFourierTransformChannel(image,
972  GrayPixelChannel,modulus,fourier_image,exception);
973  else
974  thread_status=ForwardFourierTransformChannel(image,
975  RedPixelChannel,modulus,fourier_image,exception);
976  if (thread_status == MagickFalse)
977  status=thread_status;
978  }
979 #if defined(MAGICKCORE_OPENMP_SUPPORT)
980  #pragma omp section
981 #endif
982  {
984  thread_status;
985 
986  thread_status=MagickTrue;
987  if (is_gray == MagickFalse)
988  thread_status=ForwardFourierTransformChannel(image,
989  GreenPixelChannel,modulus,fourier_image,exception);
990  if (thread_status == MagickFalse)
991  status=thread_status;
992  }
993 #if defined(MAGICKCORE_OPENMP_SUPPORT)
994  #pragma omp section
995 #endif
996  {
998  thread_status;
999 
1000  thread_status=MagickTrue;
1001  if (is_gray == MagickFalse)
1002  thread_status=ForwardFourierTransformChannel(image,
1003  BluePixelChannel,modulus,fourier_image,exception);
1004  if (thread_status == MagickFalse)
1005  status=thread_status;
1006  }
1007 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1008  #pragma omp section
1009 #endif
1010  {
1012  thread_status;
1013 
1014  thread_status=MagickTrue;
1015  if (image->colorspace == CMYKColorspace)
1016  thread_status=ForwardFourierTransformChannel(image,
1017  BlackPixelChannel,modulus,fourier_image,exception);
1018  if (thread_status == MagickFalse)
1019  status=thread_status;
1020  }
1021 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1022  #pragma omp section
1023 #endif
1024  {
1026  thread_status;
1027 
1028  thread_status=MagickTrue;
1029  if (image->alpha_trait != UndefinedPixelTrait)
1030  thread_status=ForwardFourierTransformChannel(image,
1031  AlphaPixelChannel,modulus,fourier_image,exception);
1032  if (thread_status == MagickFalse)
1033  status=thread_status;
1034  }
1035  }
1036  if (status == MagickFalse)
1037  fourier_image=DestroyImageList(fourier_image);
1038  fftw_cleanup();
1039  }
1040  }
1041  }
1042 #endif
1043  return(fourier_image);
1044 }
1045 
1046 /*
1047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1048 % %
1049 % %
1050 % %
1051 % I n v e r s e F o u r i e r T r a n s f o r m I m a g e %
1052 % %
1053 % %
1054 % %
1055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1056 %
1057 % InverseFourierTransformImage() implements the inverse discrete Fourier
1058 % transform (DFT) of the image either as a magnitude / phase or real /
1059 % imaginary image pair.
1060 %
1061 % The format of the InverseFourierTransformImage method is:
1062 %
1063 % Image *InverseFourierTransformImage(const Image *magnitude_image,
1064 % const Image *phase_image,const MagickBooleanType modulus,
1065 % ExceptionInfo *exception)
1066 %
1067 % A description of each parameter follows:
1068 %
1069 % o magnitude_image: the magnitude or real image.
1070 %
1071 % o phase_image: the phase or imaginary image.
1072 %
1073 % o modulus: if true, return transform as a magnitude / phase pair
1074 % otherwise a real / imaginary image pair.
1075 %
1076 % o exception: return any errors or warnings in this structure.
1077 %
1078 */
1079 
1080 #if defined(MAGICKCORE_FFTW_DELEGATE)
1081 static MagickBooleanType InverseQuadrantSwap(const size_t width,
1082  const size_t height,const double *source,double *destination)
1083 {
1084  ssize_t
1085  x;
1086 
1087  ssize_t
1088  center,
1089  y;
1090 
1091  /*
1092  Swap quadrants.
1093  */
1094  center=(ssize_t) (width/2L)+1L;
1095  for (y=1L; y < (ssize_t) height; y++)
1096  for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1097  destination[(height-y)*center-x+width/2L]=source[y*width+x];
1098  for (y=0L; y < (ssize_t) height; y++)
1099  destination[y*center]=source[y*width+width/2L];
1100  for (x=0L; x < center; x++)
1101  destination[x]=source[center-x-1L];
1102  return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1103 }
1104 
1105 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
1106  const Image *magnitude_image,const Image *phase_image,
1107  fftw_complex *fourier_pixels,ExceptionInfo *exception)
1108 {
1109  CacheView
1110  *magnitude_view,
1111  *phase_view;
1112 
1113  double
1114  *inverse_pixels,
1115  *magnitude_pixels,
1116  *phase_pixels;
1117 
1119  status;
1120 
1121  MemoryInfo
1122  *inverse_info,
1123  *magnitude_info,
1124  *phase_info;
1125 
1126  const Quantum
1127  *p;
1128 
1129  ssize_t
1130  i,
1131  x;
1132 
1133  ssize_t
1134  y;
1135 
1136  /*
1137  Inverse fourier - read image and break down into a double array.
1138  */
1139  magnitude_info=AcquireVirtualMemory((size_t) fourier_info->width,
1140  fourier_info->height*sizeof(*magnitude_pixels));
1141  phase_info=AcquireVirtualMemory((size_t) fourier_info->width,
1142  fourier_info->height*sizeof(*phase_pixels));
1143  inverse_info=AcquireVirtualMemory((size_t) fourier_info->width,
1144  (fourier_info->height/2+1)*sizeof(*inverse_pixels));
1145  if ((magnitude_info == (MemoryInfo *) NULL) ||
1146  (phase_info == (MemoryInfo *) NULL) ||
1147  (inverse_info == (MemoryInfo *) NULL))
1148  {
1149  if (magnitude_info != (MemoryInfo *) NULL)
1150  magnitude_info=RelinquishVirtualMemory(magnitude_info);
1151  if (phase_info != (MemoryInfo *) NULL)
1152  phase_info=RelinquishVirtualMemory(phase_info);
1153  if (inverse_info != (MemoryInfo *) NULL)
1154  inverse_info=RelinquishVirtualMemory(inverse_info);
1155  (void) ThrowMagickException(exception,GetMagickModule(),
1156  ResourceLimitError,"MemoryAllocationFailed","`%s'",
1157  magnitude_image->filename);
1158  return(MagickFalse);
1159  }
1160  magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
1161  phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
1162  inverse_pixels=(double *) GetVirtualMemoryBlob(inverse_info);
1163  i=0L;
1164  magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1165  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1166  {
1167  p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1168  exception);
1169  if (p == (const Quantum *) NULL)
1170  break;
1171  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1172  {
1173  switch (fourier_info->channel)
1174  {
1175  case RedPixelChannel:
1176  default:
1177  {
1178  magnitude_pixels[i]=QuantumScale*GetPixelRed(magnitude_image,p);
1179  break;
1180  }
1181  case GreenPixelChannel:
1182  {
1183  magnitude_pixels[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
1184  break;
1185  }
1186  case BluePixelChannel:
1187  {
1188  magnitude_pixels[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
1189  break;
1190  }
1191  case BlackPixelChannel:
1192  {
1193  magnitude_pixels[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
1194  break;
1195  }
1196  case AlphaPixelChannel:
1197  {
1198  magnitude_pixels[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
1199  break;
1200  }
1201  }
1202  i++;
1203  p+=GetPixelChannels(magnitude_image);
1204  }
1205  }
1206  magnitude_view=DestroyCacheView(magnitude_view);
1207  status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1208  magnitude_pixels,inverse_pixels);
1209  (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
1210  fourier_info->center*sizeof(*magnitude_pixels));
1211  i=0L;
1212  phase_view=AcquireVirtualCacheView(phase_image,exception);
1213  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1214  {
1215  p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1216  exception);
1217  if (p == (const Quantum *) NULL)
1218  break;
1219  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1220  {
1221  switch (fourier_info->channel)
1222  {
1223  case RedPixelChannel:
1224  default:
1225  {
1226  phase_pixels[i]=QuantumScale*GetPixelRed(phase_image,p);
1227  break;
1228  }
1229  case GreenPixelChannel:
1230  {
1231  phase_pixels[i]=QuantumScale*GetPixelGreen(phase_image,p);
1232  break;
1233  }
1234  case BluePixelChannel:
1235  {
1236  phase_pixels[i]=QuantumScale*GetPixelBlue(phase_image,p);
1237  break;
1238  }
1239  case BlackPixelChannel:
1240  {
1241  phase_pixels[i]=QuantumScale*GetPixelBlack(phase_image,p);
1242  break;
1243  }
1244  case AlphaPixelChannel:
1245  {
1246  phase_pixels[i]=QuantumScale*GetPixelAlpha(phase_image,p);
1247  break;
1248  }
1249  }
1250  i++;
1251  p+=GetPixelChannels(phase_image);
1252  }
1253  }
1254  if (fourier_info->modulus != MagickFalse)
1255  {
1256  i=0L;
1257  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1258  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1259  {
1260  phase_pixels[i]-=0.5;
1261  phase_pixels[i]*=(2.0*MagickPI);
1262  i++;
1263  }
1264  }
1265  phase_view=DestroyCacheView(phase_view);
1266  CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1267  if (status != MagickFalse)
1268  status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1269  phase_pixels,inverse_pixels);
1270  (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
1271  fourier_info->center*sizeof(*phase_pixels));
1272  inverse_info=RelinquishVirtualMemory(inverse_info);
1273  /*
1274  Merge two sets.
1275  */
1276  i=0L;
1277  if (fourier_info->modulus != MagickFalse)
1278  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1279  for (x=0L; x < (ssize_t) fourier_info->center; x++)
1280  {
1281 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1282  fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1283  magnitude_pixels[i]*sin(phase_pixels[i]);
1284 #else
1285  fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1286  fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1287 #endif
1288  i++;
1289  }
1290  else
1291  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1292  for (x=0L; x < (ssize_t) fourier_info->center; x++)
1293  {
1294 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1295  fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1296 #else
1297  fourier_pixels[i][0]=magnitude_pixels[i];
1298  fourier_pixels[i][1]=phase_pixels[i];
1299 #endif
1300  i++;
1301  }
1302  magnitude_info=RelinquishVirtualMemory(magnitude_info);
1303  phase_info=RelinquishVirtualMemory(phase_info);
1304  return(status);
1305 }
1306 
1307 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1308  fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception)
1309 {
1310  CacheView
1311  *image_view;
1312 
1313  const char
1314  *value;
1315 
1316  double
1317  *source_pixels;
1318 
1319  fftw_plan
1320  fftw_c2r_plan;
1321 
1322  MemoryInfo
1323  *source_info;
1324 
1325  Quantum
1326  *q;
1327 
1328  ssize_t
1329  i,
1330  x;
1331 
1332  ssize_t
1333  y;
1334 
1335  source_info=AcquireVirtualMemory((size_t) fourier_info->width,
1336  fourier_info->height*sizeof(*source_pixels));
1337  if (source_info == (MemoryInfo *) NULL)
1338  {
1339  (void) ThrowMagickException(exception,GetMagickModule(),
1340  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1341  return(MagickFalse);
1342  }
1343  source_pixels=(double *) GetVirtualMemoryBlob(source_info);
1344  value=GetImageArtifact(image,"fourier:normalize");
1345  if (LocaleCompare(value,"inverse") == 0)
1346  {
1347  double
1348  gamma;
1349 
1350  /*
1351  Normalize inverse transform.
1352  */
1353  i=0L;
1354  gamma=PerceptibleReciprocal((double) fourier_info->width*
1355  fourier_info->height);
1356  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1357  for (x=0L; x < (ssize_t) fourier_info->center; x++)
1358  {
1359 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1360  fourier_pixels[i]*=gamma;
1361 #else
1362  fourier_pixels[i][0]*=gamma;
1363  fourier_pixels[i][1]*=gamma;
1364 #endif
1365  i++;
1366  }
1367  }
1368 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1369  #pragma omp critical (MagickCore_InverseFourierTransform)
1370 #endif
1371  fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1372  fourier_pixels,source_pixels,FFTW_ESTIMATE);
1373  fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
1374  fftw_destroy_plan(fftw_c2r_plan);
1375  i=0L;
1376  image_view=AcquireAuthenticCacheView(image,exception);
1377  for (y=0L; y < (ssize_t) fourier_info->height; y++)
1378  {
1379  if (y >= (ssize_t) image->rows)
1380  break;
1381  q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1382  image->columns ? image->columns : fourier_info->width,1UL,exception);
1383  if (q == (Quantum *) NULL)
1384  break;
1385  for (x=0L; x < (ssize_t) fourier_info->width; x++)
1386  {
1387  if (x < (ssize_t) image->columns)
1388  switch (fourier_info->channel)
1389  {
1390  case RedPixelChannel:
1391  default:
1392  {
1393  SetPixelRed(image,ClampToQuantum(QuantumRange*source_pixels[i]),q);
1394  break;
1395  }
1396  case GreenPixelChannel:
1397  {
1398  SetPixelGreen(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1399  q);
1400  break;
1401  }
1402  case BluePixelChannel:
1403  {
1404  SetPixelBlue(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1405  q);
1406  break;
1407  }
1408  case BlackPixelChannel:
1409  {
1410  SetPixelBlack(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1411  q);
1412  break;
1413  }
1414  case AlphaPixelChannel:
1415  {
1416  SetPixelAlpha(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1417  q);
1418  break;
1419  }
1420  }
1421  i++;
1422  q+=GetPixelChannels(image);
1423  }
1424  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1425  break;
1426  }
1427  image_view=DestroyCacheView(image_view);
1428  source_info=RelinquishVirtualMemory(source_info);
1429  return(MagickTrue);
1430 }
1431 
1432 static MagickBooleanType InverseFourierTransformChannel(
1433  const Image *magnitude_image,const Image *phase_image,
1434  const PixelChannel channel,const MagickBooleanType modulus,
1435  Image *fourier_image,ExceptionInfo *exception)
1436 {
1437  fftw_complex
1438  *inverse_pixels;
1439 
1440  FourierInfo
1441  fourier_info;
1442 
1444  status;
1445 
1446  MemoryInfo
1447  *inverse_info;
1448 
1449  fourier_info.width=magnitude_image->columns;
1450  fourier_info.height=magnitude_image->rows;
1451  if ((magnitude_image->columns != magnitude_image->rows) ||
1452  ((magnitude_image->columns % 2) != 0) ||
1453  ((magnitude_image->rows % 2) != 0))
1454  {
1455  size_t extent=magnitude_image->columns < magnitude_image->rows ?
1456  magnitude_image->rows : magnitude_image->columns;
1457  fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1458  }
1459  fourier_info.height=fourier_info.width;
1460  fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
1461  fourier_info.channel=channel;
1462  fourier_info.modulus=modulus;
1463  inverse_info=AcquireVirtualMemory((size_t) fourier_info.width,
1464  (fourier_info.height/2+1)*sizeof(*inverse_pixels));
1465  if (inverse_info == (MemoryInfo *) NULL)
1466  {
1467  (void) ThrowMagickException(exception,GetMagickModule(),
1468  ResourceLimitError,"MemoryAllocationFailed","`%s'",
1469  magnitude_image->filename);
1470  return(MagickFalse);
1471  }
1472  inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1473  status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1474  inverse_pixels,exception);
1475  if (status != MagickFalse)
1476  status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1477  exception);
1478  inverse_info=RelinquishVirtualMemory(inverse_info);
1479  return(status);
1480 }
1481 #endif
1482 
1484  const Image *phase_image,const MagickBooleanType modulus,
1485  ExceptionInfo *exception)
1486 {
1487  Image
1488  *fourier_image;
1489 
1490  assert(magnitude_image != (Image *) NULL);
1491  assert(magnitude_image->signature == MagickCoreSignature);
1492  if (magnitude_image->debug != MagickFalse)
1494  magnitude_image->filename);
1495  if (phase_image == (Image *) NULL)
1496  {
1497  (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1498  "ImageSequenceRequired","`%s'",magnitude_image->filename);
1499  return((Image *) NULL);
1500  }
1501 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1502  fourier_image=(Image *) NULL;
1503  (void) modulus;
1504  (void) ThrowMagickException(exception,GetMagickModule(),
1505  MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1506  magnitude_image->filename);
1507 #else
1508  {
1509  fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1510  magnitude_image->rows,MagickTrue,exception);
1511  if (fourier_image != (Image *) NULL)
1512  {
1514  is_gray,
1515  status;
1516 
1517  status=MagickTrue;
1518  is_gray=IsImageGray(magnitude_image);
1519  if (is_gray != MagickFalse)
1520  is_gray=IsImageGray(phase_image);
1521 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1522  #pragma omp parallel sections
1523 #endif
1524  {
1525 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1526  #pragma omp section
1527 #endif
1528  {
1530  thread_status;
1531 
1532  if (is_gray != MagickFalse)
1533  thread_status=InverseFourierTransformChannel(magnitude_image,
1534  phase_image,GrayPixelChannel,modulus,fourier_image,exception);
1535  else
1536  thread_status=InverseFourierTransformChannel(magnitude_image,
1537  phase_image,RedPixelChannel,modulus,fourier_image,exception);
1538  if (thread_status == MagickFalse)
1539  status=thread_status;
1540  }
1541 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1542  #pragma omp section
1543 #endif
1544  {
1546  thread_status;
1547 
1548  thread_status=MagickTrue;
1549  if (is_gray == MagickFalse)
1550  thread_status=InverseFourierTransformChannel(magnitude_image,
1551  phase_image,GreenPixelChannel,modulus,fourier_image,exception);
1552  if (thread_status == MagickFalse)
1553  status=thread_status;
1554  }
1555 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1556  #pragma omp section
1557 #endif
1558  {
1560  thread_status;
1561 
1562  thread_status=MagickTrue;
1563  if (is_gray == MagickFalse)
1564  thread_status=InverseFourierTransformChannel(magnitude_image,
1565  phase_image,BluePixelChannel,modulus,fourier_image,exception);
1566  if (thread_status == MagickFalse)
1567  status=thread_status;
1568  }
1569 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1570  #pragma omp section
1571 #endif
1572  {
1574  thread_status;
1575 
1576  thread_status=MagickTrue;
1577  if (magnitude_image->colorspace == CMYKColorspace)
1578  thread_status=InverseFourierTransformChannel(magnitude_image,
1579  phase_image,BlackPixelChannel,modulus,fourier_image,exception);
1580  if (thread_status == MagickFalse)
1581  status=thread_status;
1582  }
1583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1584  #pragma omp section
1585 #endif
1586  {
1588  thread_status;
1589 
1590  thread_status=MagickTrue;
1591  if (magnitude_image->alpha_trait != UndefinedPixelTrait)
1592  thread_status=InverseFourierTransformChannel(magnitude_image,
1593  phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
1594  if (thread_status == MagickFalse)
1595  status=thread_status;
1596  }
1597  }
1598  if (status == MagickFalse)
1599  fourier_image=DestroyImage(fourier_image);
1600  }
1601  fftw_cleanup();
1602  }
1603 #endif
1604  return(fourier_image);
1605 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
size_t height
Definition: fourier.c:95
MagickExport CacheView * DestroyCacheView(CacheView *cache_view)
Definition: cache-view.c:252
MagickExport MemoryInfo * RelinquishVirtualMemory(MemoryInfo *memory_info)
Definition: memory.c:1229
MagickProgressMonitor progress_monitor
Definition: image.h:303
static Quantum GetPixelAlpha(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
static Quantum GetPixelRed(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport MemoryInfo * AcquireVirtualMemory(const size_t count, const size_t quantum)
Definition: memory.c:705
size_t signature
Definition: exception.h:123
PixelChannel channel
Definition: fourier.c:89
MagickExport const char * GetImageArtifact(const Image *image, const char *artifact)
Definition: artifact.c:273
#define ComplexImageTag
#define MagickPI
Definition: image-private.h:42
MagickExport const Quantum * GetCacheViewVirtualPixels(const CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:651
ClassType storage_class
Definition: image.h:154
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:133
static Quantum ClampToQuantum(const MagickRealType quantum)
Definition: quantum.h:85
Definition: image.h:151
MagickExport Image * ForwardFourierTransformImage(const Image *image, const MagickBooleanType modulus, ExceptionInfo *exception)
Definition: fourier.c:905
#define MagickCoreSignature
MagickExport Quantum * GetCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:299
MagickExport Image * GetFirstImageInList(const Image *images)
Definition: list.c:576
MagickBooleanType
Definition: magick-type.h:161
MagickExport Image * NewImageList(void)
Definition: list.c:953
static double PerceptibleReciprocal(const double x)
size_t width
Definition: fourier.c:95
static Quantum GetPixelGreen(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
PixelTrait alpha_trait
Definition: image.h:280
static Quantum GetPixelBlack(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport Quantum * QueueCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:977
MagickExport MagickBooleanType ThrowMagickException(ExceptionInfo *exception, const char *module, const char *function, const size_t line, const ExceptionType severity, const char *tag, const char *format,...)
Definition: exception.c:1145
MagickExport MagickBooleanType LogMagickEvent(const LogEventType type, const char *module, const char *function, const size_t line, const char *format,...)
Definition: log.c:1662
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:119
size_t columns
Definition: image.h:172
ssize_t center
Definition: fourier.c:99
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
struct _Image * next
Definition: image.h:348
struct _FourierInfo FourierInfo
static void SetPixelBlue(const Image *magick_restrict image, const Quantum blue, Quantum *magick_restrict pixel)
MagickExport MagickBooleanType SetImageStorageClass(Image *image, const ClassType storage_class, ExceptionInfo *exception)
Definition: image.c:2616
MagickExport Image * DestroyImageList(Image *images)
Definition: list.c:477
PixelChannel
Definition: pixel.h:70
static size_t GetPixelChannels(const Image *magick_restrict image)
MagickBooleanType modulus
Definition: fourier.c:92
MagickExport int LocaleCompare(const char *p, const char *q)
Definition: locale.c:1401
char filename[MagickPathExtent]
Definition: image.h:319
#define GetMagickModule()
Definition: log.h:28
MagickExport CacheView * AcquireVirtualCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:149
MagickExport MagickBooleanType IsImageGray(const Image *image)
Definition: attribute.c:1776
unsigned short Quantum
Definition: magick-type.h:86
MagickExport Image * GetNextImageInList(const Image *images)
Definition: list.c:786
size_t number_channels
Definition: image.h:283
MagickExport void AppendImageToList(Image **images, const Image *append)
Definition: list.c:80
#define MagickMin(x, y)
Definition: image-private.h:39
static void SetPixelAlpha(const Image *magick_restrict image, const Quantum alpha, Quantum *magick_restrict pixel)
static double StringToDouble(const char *magick_restrict string, char *magick_restrict *sentinal)
static void SetPixelRed(const Image *magick_restrict image, const Quantum red, Quantum *magick_restrict pixel)
ComplexOperator
Definition: fourier.h:25
#define MagickExport
MagickExport MagickBooleanType SyncCacheViewAuthenticPixels(CacheView *magick_restrict cache_view, ExceptionInfo *exception)
Definition: cache-view.c:1100
MagickExport CacheView * AcquireAuthenticCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:112
static void SetPixelBlack(const Image *magick_restrict image, const Quantum black, Quantum *magick_restrict pixel)
static Quantum GetPixelBlue(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport void * GetVirtualMemoryBlob(const MemoryInfo *memory_info)
Definition: memory.c:1090
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1179
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:787
ColorspaceType colorspace
Definition: image.h:157
#define QuantumRange
Definition: magick-type.h:87
MagickExport Image * ComplexImages(const Image *images, const ComplexOperator op, ExceptionInfo *exception)
Definition: fourier.c:129
MagickExport MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
Definition: monitor.c:136
MagickBooleanType debug
Definition: image.h:334
static void SetPixelGreen(const Image *magick_restrict image, const Quantum green, Quantum *magick_restrict pixel)
size_t depth
Definition: image.h:172
MagickExport Image * InverseFourierTransformImage(const Image *magnitude_image, const Image *phase_image, const MagickBooleanType modulus, ExceptionInfo *exception)
Definition: fourier.c:1483