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