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