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