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