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