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