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