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