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