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