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