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