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