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