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