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