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]-Ai[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]*exp(Ai[i]);
 301             Ci[i]=Ar[i]*(cos(Ai[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   if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
 847       ((image->rows % 2) != 0))
 848     {
 849       extent=image->columns < image->rows ? image->rows : image->columns;
 850       fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
 851     }
 852   fourier_info.height=fourier_info.width;
 853   fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
 854   fourier_info.channel=channel;
 855   fourier_info.modulus=modulus;
 856   magnitude_info=AcquireVirtualMemory((size_t) fourier_info.height,
 857     fourier_info.center*sizeof(*magnitude_pixels));
 858   phase_info=AcquireVirtualMemory((size_t) fourier_info.height,
 859     fourier_info.center*sizeof(*phase_pixels));
 860   if ((magnitude_info == (MemoryInfo *) NULL) ||
 861       (phase_info == (MemoryInfo *) NULL))
 862     {
 863       if (phase_info != (MemoryInfo *) NULL)
 864         phase_info=RelinquishVirtualMemory(phase_info);
 865       if (magnitude_info == (MemoryInfo *) NULL)
 866         magnitude_info=RelinquishVirtualMemory(magnitude_info);
 867       (void) ThrowMagickException(exception,GetMagickModule(),
 868         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 869       return(MagickFalse);
 870     }
 871   magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
 872   phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
 873   status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
 874     phase_pixels,exception);
 875   if (status != MagickFalse)
 876     status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
 877       phase_pixels,exception);
 878   phase_info=RelinquishVirtualMemory(phase_info);
 879   magnitude_info=RelinquishVirtualMemory(magnitude_info);
 880   return(status);
 881 }
 882 #endif
 883
 884 MagickExport Image *ForwardFourierTransformImage(const Image *image,
 885   const MagickBooleanType modulus,ExceptionInfo *exception)
 886 {
 887   Image
 888     *fourier_image;
 889
 890   fourier_image=NewImageList();
 891 #if !defined(MAGICKCORE_FFTW_DELEGATE)
 892   (void) modulus;
 893   (void) ThrowMagickException(exception,GetMagickModule(),
 894     MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
 895     image->filename);
 896 #else
 897   {
 898     Image
 899       *magnitude_image;
 900
 901     size_t
 902       extent,
 903       width;
 904
 905     width=image->columns;
 906     if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
 907         ((image->rows % 2) != 0))
 908       {
 909         extent=image->columns < image->rows ? image->rows : image->columns;
 910         width=(extent & 0x01) == 1 ? extent+1UL : extent;
 911       }
 912     magnitude_image=CloneImage(image,width,width,MagickTrue,exception);
 913     if (magnitude_image != (Image *) NULL)
 914       {
 915         Image
 916           *phase_image;
 917
 918         magnitude_image->storage_class=DirectClass;
 919         magnitude_image->depth=32UL;
 920         phase_image=CloneImage(image,width,width,MagickTrue,exception);
 921         if (phase_image == (Image *) NULL)
 922           magnitude_image=DestroyImage(magnitude_image);
 923         else
 924           {
 925             MagickBooleanType
 926               is_gray,
 927               status;
 928
 929             phase_image->storage_class=DirectClass;
 930             phase_image->depth=32UL;
 931             AppendImageToList(&fourier_image,magnitude_image);
 932             AppendImageToList(&fourier_image,phase_image);
 933             status=MagickTrue;
 934             is_gray=IsImageGray(image,exception);
 935 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 936             #pragma omp parallel sections
 937 #endif
 938             {
 939 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 940               #pragma omp section
 941 #endif
 942               {
 943                 MagickBooleanType
 944                   thread_status;
 945
 946                 if (is_gray != MagickFalse)
 947                   thread_status=ForwardFourierTransformChannel(image,
 948                     GrayPixelChannel,modulus,fourier_image,exception);
 949                 else
 950                   thread_status=ForwardFourierTransformChannel(image,
 951                     RedPixelChannel,modulus,fourier_image,exception);
 952                 if (thread_status == MagickFalse)
 953                   status=thread_status;
 954               }
 955 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 956               #pragma omp section
 957 #endif
 958               {
 959                 MagickBooleanType
 960                   thread_status;
 961
 962                 thread_status=MagickTrue;
 963                 if (is_gray == MagickFalse)
 964                   thread_status=ForwardFourierTransformChannel(image,
 965                     GreenPixelChannel,modulus,fourier_image,exception);
 966                 if (thread_status == MagickFalse)
 967                   status=thread_status;
 968               }
 969 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 970               #pragma omp section
 971 #endif
 972               {
 973                 MagickBooleanType
 974                   thread_status;
 975
 976                 thread_status=MagickTrue;
 977                 if (is_gray == MagickFalse)
 978                   thread_status=ForwardFourierTransformChannel(image,
 979                     BluePixelChannel,modulus,fourier_image,exception);
 980                 if (thread_status == MagickFalse)
 981                   status=thread_status;
 982               }
 983 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 984               #pragma omp section
 985 #endif
 986               {
 987                 MagickBooleanType
 988                   thread_status;
 989
 990                 thread_status=MagickTrue;
 991                 if (image->colorspace == CMYKColorspace)
 992                   thread_status=ForwardFourierTransformChannel(image,
 993                     BlackPixelChannel,modulus,fourier_image,exception);
 994                 if (thread_status == MagickFalse)
 995                   status=thread_status;
 996               }
 997 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 998               #pragma omp section
 999 #endif
1000               {
1001                 MagickBooleanType
1002                   thread_status;
1003
1004                 thread_status=MagickTrue;
1005                 if (image->alpha_trait == BlendPixelTrait)
1006                   thread_status=ForwardFourierTransformChannel(image,
1007                     AlphaPixelChannel,modulus,fourier_image,exception);
1008                 if (thread_status == MagickFalse)
1009                   status=thread_status;
1010               }
1011             }
1012             if (status == MagickFalse)
1013               fourier_image=DestroyImageList(fourier_image);
1014             fftw_cleanup();
1015           }
1016       }
1017   }
1018 #endif
1019   return(fourier_image);
1020 }
1021 \f
1022 /*
1023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1024 %                                                                             %
1025 %                                                                             %
1026 %                                                                             %
1027 %     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                 %
1028 %                                                                             %
1029 %                                                                             %
1030 %                                                                             %
1031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1032 %
1033 %  InverseFourierTransformImage() implements the inverse discrete Fourier
1034 %  transform (DFT) of the image either as a magnitude / phase or real /
1035 %  imaginary image pair.
1036 %
1037 %  The format of the InverseFourierTransformImage method is:
1038 %
1039 %      Image *InverseFourierTransformImage(const Image *magnitude_image,
1040 %        const Image *phase_image,const MagickBooleanType modulus,
1041 %        ExceptionInfo *exception)
1042 %
1043 %  A description of each parameter follows:
1044 %
1045 %    o magnitude_image: the magnitude or real image.
1046 %
1047 %    o phase_image: the phase or imaginary image.
1048 %
1049 %    o modulus: if true, return transform as a magnitude / phase pair
1050 %      otherwise a real / imaginary image pair.
1051 %
1052 %    o exception: return any errors or warnings in this structure.
1053 %
1054 */
1055
1056 #if defined(MAGICKCORE_FFTW_DELEGATE)
1057 static MagickBooleanType InverseQuadrantSwap(const size_t width,
1058   const size_t height,const double *source,double *destination)
1059 {
1060   register ssize_t
1061     x;
1062
1063   ssize_t
1064     center,
1065     y;
1066
1067   /*
1068     Swap quadrants.
1069   */
1070   center=(ssize_t) floor((double) width/2L)+1L;
1071   for (y=1L; y < (ssize_t) height; y++)
1072     for (x=0L; x < (ssize_t) (width/2L+1L); x++)
1073       destination[(height-y)*center-x+width/2L]=source[y*width+x];
1074   for (y=0L; y < (ssize_t) height; y++)
1075     destination[y*center]=source[y*width+width/2L];
1076   for (x=0L; x < center; x++)
1077     destination[x]=source[center-x-1L];
1078   return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
1079 }
1080
1081 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
1082   const Image *magnitude_image,const Image *phase_image,
1083   fftw_complex *fourier_pixels,ExceptionInfo *exception)
1084 {
1085   CacheView
1086     *magnitude_view,
1087     *phase_view;
1088
1089   double
1090     *inverse_pixels,
1091     *magnitude_pixels,
1092     *phase_pixels;
1093
1094   MagickBooleanType
1095     status;
1096
1097   MemoryInfo
1098     *inverse_info,
1099     *magnitude_info,
1100     *phase_info;
1101
1102   register const Quantum
1103     *p;
1104
1105   register ssize_t
1106     i,
1107     x;
1108
1109   ssize_t
1110     y;
1111
1112   /*
1113     Inverse fourier - read image and break down into a double array.
1114   */
1115   magnitude_info=AcquireVirtualMemory((size_t)fourier_info->height,
1116     fourier_info->width*sizeof(*magnitude_pixels));
1117   phase_info=AcquireVirtualMemory((size_t) fourier_info->height,
1118     fourier_info->width*sizeof(*phase_pixels));
1119   inverse_info=AcquireVirtualMemory((size_t) fourier_info->height,
1120     fourier_info->center*sizeof(*inverse_pixels));
1121   if ((magnitude_info == (MemoryInfo *) NULL) ||
1122       (phase_info == (MemoryInfo *) NULL) ||
1123       (inverse_info == (MemoryInfo *) NULL))
1124     {
1125       if (magnitude_info != (MemoryInfo *) NULL)
1126         magnitude_info=RelinquishVirtualMemory(magnitude_info);
1127       if (phase_info != (MemoryInfo *) NULL)
1128         phase_info=RelinquishVirtualMemory(phase_info);
1129       if (inverse_info != (MemoryInfo *) NULL)
1130         inverse_info=RelinquishVirtualMemory(inverse_info);
1131       (void) ThrowMagickException(exception,GetMagickModule(),
1132         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1133         magnitude_image->filename);
1134       return(MagickFalse);
1135     }
1136   magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
1137   phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
1138   inverse_pixels=(double *) GetVirtualMemoryBlob(inverse_info);
1139   i=0L;
1140   magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
1141   for (y=0L; y < (ssize_t) fourier_info->height; y++)
1142   {
1143     p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
1144       exception);
1145     if (p == (const Quantum *) NULL)
1146       break;
1147     for (x=0L; x < (ssize_t) fourier_info->width; x++)
1148     {
1149       switch (fourier_info->channel)
1150       {
1151         case RedPixelChannel:
1152         default:
1153         {
1154           magnitude_pixels[i]=QuantumScale*GetPixelRed(magnitude_image,p);
1155           break;
1156         }
1157         case GreenPixelChannel:
1158         {
1159           magnitude_pixels[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
1160           break;
1161         }
1162         case BluePixelChannel:
1163         {
1164           magnitude_pixels[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
1165           break;
1166         }
1167         case BlackPixelChannel:
1168         {
1169           magnitude_pixels[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
1170           break;
1171         }
1172         case AlphaPixelChannel:
1173         {
1174           magnitude_pixels[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
1175           break;
1176         }
1177       }
1178       i++;
1179       p+=GetPixelChannels(magnitude_image);
1180     }
1181   }
1182   magnitude_view=DestroyCacheView(magnitude_view);
1183   status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1184     magnitude_pixels,inverse_pixels);
1185   (void) CopyMagickMemory(magnitude_pixels,inverse_pixels,fourier_info->height*
1186     fourier_info->center*sizeof(*magnitude_pixels));
1187   i=0L;
1188   phase_view=AcquireVirtualCacheView(phase_image,exception);
1189   for (y=0L; y < (ssize_t) fourier_info->height; y++)
1190   {
1191     p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
1192       exception);
1193     if (p == (const Quantum *) NULL)
1194       break;
1195     for (x=0L; x < (ssize_t) fourier_info->width; x++)
1196     {
1197       switch (fourier_info->channel)
1198       {
1199         case RedPixelChannel:
1200         default:
1201         {
1202           phase_pixels[i]=QuantumScale*GetPixelRed(phase_image,p);
1203           break;
1204         }
1205         case GreenPixelChannel:
1206         {
1207           phase_pixels[i]=QuantumScale*GetPixelGreen(phase_image,p);
1208           break;
1209         }
1210         case BluePixelChannel:
1211         {
1212           phase_pixels[i]=QuantumScale*GetPixelBlue(phase_image,p);
1213           break;
1214         }
1215         case BlackPixelChannel:
1216         {
1217           phase_pixels[i]=QuantumScale*GetPixelBlack(phase_image,p);
1218           break;
1219         }
1220         case AlphaPixelChannel:
1221         {
1222           phase_pixels[i]=QuantumScale*GetPixelAlpha(phase_image,p);
1223           break;
1224         }
1225       }
1226       i++;
1227       p+=GetPixelChannels(phase_image);
1228     }
1229   }
1230   if (fourier_info->modulus != MagickFalse)
1231     {
1232       i=0L;
1233       for (y=0L; y < (ssize_t) fourier_info->height; y++)
1234         for (x=0L; x < (ssize_t) fourier_info->width; x++)
1235         {
1236           phase_pixels[i]-=0.5;
1237           phase_pixels[i]*=(2.0*MagickPI);
1238           i++;
1239         }
1240     }
1241   phase_view=DestroyCacheView(phase_view);
1242   CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
1243   if (status != MagickFalse)
1244     status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1245       phase_pixels,inverse_pixels);
1246   (void) CopyMagickMemory(phase_pixels,inverse_pixels,fourier_info->height*
1247     fourier_info->center*sizeof(*phase_pixels));
1248   inverse_info=RelinquishVirtualMemory(inverse_info);
1249   /*
1250     Merge two sets.
1251   */
1252   i=0L;
1253   if (fourier_info->modulus != MagickFalse)
1254     for (y=0L; y < (ssize_t) fourier_info->height; y++)
1255        for (x=0L; x < (ssize_t) fourier_info->center; x++)
1256        {
1257 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1258          fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
1259            magnitude_pixels[i]*sin(phase_pixels[i]);
1260 #else
1261          fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
1262          fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
1263 #endif
1264          i++;
1265       }
1266   else
1267     for (y=0L; y < (ssize_t) fourier_info->height; y++)
1268       for (x=0L; x < (ssize_t) fourier_info->center; x++)
1269       {
1270 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1271         fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
1272 #else
1273         fourier_pixels[i][0]=magnitude_pixels[i];
1274         fourier_pixels[i][1]=phase_pixels[i];
1275 #endif
1276         i++;
1277       }
1278   magnitude_info=RelinquishVirtualMemory(magnitude_info);
1279   phase_info=RelinquishVirtualMemory(phase_info);
1280   return(status);
1281 }
1282
1283 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1284   fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception)
1285 {
1286   CacheView
1287     *image_view;
1288
1289   const char
1290     *value;
1291
1292   double
1293     *source_pixels;
1294
1295   fftw_plan
1296     fftw_c2r_plan;
1297
1298   MemoryInfo
1299     *source_info;
1300
1301   register Quantum
1302     *q;
1303
1304   register ssize_t
1305     i,
1306     x;
1307
1308   ssize_t
1309     y;
1310
1311   source_info=AcquireVirtualMemory((size_t) fourier_info->height,
1312     fourier_info->width*sizeof(*source_pixels));
1313   if (source_info == (MemoryInfo *) NULL)
1314     {
1315       (void) ThrowMagickException(exception,GetMagickModule(),
1316         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1317       return(MagickFalse);
1318     }
1319   source_pixels=(double *) GetVirtualMemoryBlob(source_info);
1320   value=GetImageArtifact(image,"fourier:normalize");
1321   if (LocaleCompare(value,"inverse") == 0)
1322     {
1323       double
1324         gamma;
1325
1326       /*
1327         Normalize inverse transform.
1328       */
1329       i=0L;
1330       gamma=PerceptibleReciprocal((double) fourier_info->width*
1331         fourier_info->height);
1332       for (y=0L; y < (ssize_t) fourier_info->height; y++)
1333         for (x=0L; x < (ssize_t) fourier_info->center; x++)
1334         {
1335 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1336           fourier_pixels[i]*=gamma;
1337 #else
1338           fourier_pixels[i][0]*=gamma;
1339           fourier_pixels[i][1]*=gamma;
1340 #endif
1341           i++;
1342         }
1343     }
1344 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1345   #pragma omp critical (MagickCore_InverseFourierTransform)
1346 #endif
1347   {
1348     fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1349       fourier_pixels,source_pixels,FFTW_ESTIMATE);
1350     fftw_execute(fftw_c2r_plan);
1351     fftw_destroy_plan(fftw_c2r_plan);
1352   }
1353   i=0L;
1354   image_view=AcquireAuthenticCacheView(image,exception);
1355   for (y=0L; y < (ssize_t) fourier_info->height; y++)
1356   {
1357     if (y >= (ssize_t) image->rows)
1358       break;
1359     q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1360       image->columns ? image->columns : fourier_info->width,1UL,exception);
1361     if (q == (Quantum *) NULL)
1362       break;
1363     for (x=0L; x < (ssize_t) fourier_info->width; x++)
1364     {
1365       if (x < (ssize_t) image->columns)
1366         switch (fourier_info->channel)
1367         {
1368           case RedPixelChannel:
1369           default:
1370           {
1371             SetPixelRed(image,ClampToQuantum(QuantumRange*source_pixels[i]),q);
1372             break;
1373           }
1374           case GreenPixelChannel:
1375           {
1376             SetPixelGreen(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1377               q);
1378             break;
1379           }
1380           case BluePixelChannel:
1381           {
1382             SetPixelBlue(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1383               q);
1384             break;
1385           }
1386           case BlackPixelChannel:
1387           {
1388             SetPixelBlack(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1389               q);
1390             break;
1391           }
1392           case AlphaPixelChannel:
1393           {
1394             SetPixelAlpha(image,ClampToQuantum(QuantumRange*source_pixels[i]),
1395               q);
1396             break;
1397           }
1398         }
1399       i++;
1400       q+=GetPixelChannels(image);
1401     }
1402     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1403       break;
1404   }
1405   image_view=DestroyCacheView(image_view);
1406   source_info=RelinquishVirtualMemory(source_info);
1407   return(MagickTrue);
1408 }
1409
1410 static MagickBooleanType InverseFourierTransformChannel(
1411   const Image *magnitude_image,const Image *phase_image,
1412   const PixelChannel channel,const MagickBooleanType modulus,
1413   Image *fourier_image,ExceptionInfo *exception)
1414 {
1415   fftw_complex
1416     *inverse_pixels;
1417
1418   FourierInfo
1419     fourier_info;
1420
1421   MagickBooleanType
1422     status;
1423
1424   MemoryInfo
1425     *inverse_info;
1426
1427   size_t
1428     extent;
1429
1430   fourier_info.width=magnitude_image->columns;
1431   if ((magnitude_image->columns != magnitude_image->rows) ||
1432       ((magnitude_image->columns % 2) != 0) ||
1433       ((magnitude_image->rows % 2) != 0))
1434     {
1435       extent=magnitude_image->columns < magnitude_image->rows ?
1436         magnitude_image->rows : magnitude_image->columns;
1437       fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1438     }
1439   fourier_info.height=fourier_info.width;
1440   fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
1441   fourier_info.channel=channel;
1442   fourier_info.modulus=modulus;
1443   inverse_info=AcquireVirtualMemory((size_t) fourier_info.height,
1444     fourier_info.center*sizeof(*inverse_pixels));
1445   if (inverse_info == (MemoryInfo *) NULL)
1446     {
1447       (void) ThrowMagickException(exception,GetMagickModule(),
1448         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1449         magnitude_image->filename);
1450       return(MagickFalse);
1451     }
1452   inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
1453   status=InverseFourier(&fourier_info,magnitude_image,phase_image,
1454     inverse_pixels,exception);
1455   if (status != MagickFalse)
1456     status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
1457       exception);
1458   inverse_info=RelinquishVirtualMemory(inverse_info);
1459   return(status);
1460 }
1461 #endif
1462
1463 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1464   const Image *phase_image,const MagickBooleanType modulus,
1465   ExceptionInfo *exception)
1466 {
1467   Image
1468     *fourier_image;
1469
1470   assert(magnitude_image != (Image *) NULL);
1471   assert(magnitude_image->signature == MagickSignature);
1472   if (magnitude_image->debug != MagickFalse)
1473     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1474       magnitude_image->filename);
1475   if (phase_image == (Image *) NULL)
1476     {
1477       (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1478         "TwoOrMoreImagesRequired","`%s'",magnitude_image->filename);
1479       return((Image *) NULL);
1480     }
1481 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1482   fourier_image=(Image *) NULL;
1483   (void) modulus;
1484   (void) ThrowMagickException(exception,GetMagickModule(),
1485     MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","'%s' (FFTW)",
1486     magnitude_image->filename);
1487 #else
1488   {
1489     fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1490       magnitude_image->rows,MagickTrue,exception);
1491     if (fourier_image != (Image *) NULL)
1492       {
1493         MagickBooleanType
1494           is_gray,
1495           status;
1496
1497         status=MagickTrue;
1498         is_gray=IsImageGray(magnitude_image,exception);
1499         if (is_gray != MagickFalse)
1500           is_gray=IsImageGray(phase_image,exception);
1501 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1502         #pragma omp parallel sections
1503 #endif
1504         {
1505 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1506           #pragma omp section
1507 #endif
1508           {
1509             MagickBooleanType
1510               thread_status;
1511
1512             if (is_gray != MagickFalse)
1513               thread_status=InverseFourierTransformChannel(magnitude_image,
1514                 phase_image,GrayPixelChannel,modulus,fourier_image,exception);
1515             else
1516               thread_status=InverseFourierTransformChannel(magnitude_image,
1517                 phase_image,RedPixelChannel,modulus,fourier_image,exception);
1518             if (thread_status == MagickFalse)
1519               status=thread_status;
1520           }
1521 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1522           #pragma omp section
1523 #endif
1524           {
1525             MagickBooleanType
1526               thread_status;
1527
1528             thread_status=MagickTrue;
1529             if (is_gray == MagickFalse)
1530               thread_status=InverseFourierTransformChannel(magnitude_image,
1531                 phase_image,GreenPixelChannel,modulus,fourier_image,exception);
1532             if (thread_status == MagickFalse)
1533               status=thread_status;
1534           }
1535 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1536           #pragma omp section
1537 #endif
1538           {
1539             MagickBooleanType
1540               thread_status;
1541
1542             thread_status=MagickTrue;
1543             if (is_gray == MagickFalse)
1544               thread_status=InverseFourierTransformChannel(magnitude_image,
1545                 phase_image,BluePixelChannel,modulus,fourier_image,exception);
1546             if (thread_status == MagickFalse)
1547               status=thread_status;
1548           }
1549 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1550           #pragma omp section
1551 #endif
1552           {
1553             MagickBooleanType
1554               thread_status;
1555
1556             thread_status=MagickTrue;
1557             if (magnitude_image->colorspace == CMYKColorspace)
1558               thread_status=InverseFourierTransformChannel(magnitude_image,
1559                 phase_image,BlackPixelChannel,modulus,fourier_image,exception);
1560             if (thread_status == MagickFalse)
1561               status=thread_status;
1562           }
1563 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1564           #pragma omp section
1565 #endif
1566           {
1567             MagickBooleanType
1568               thread_status;
1569
1570             thread_status=MagickTrue;
1571             if (magnitude_image->alpha_trait == BlendPixelTrait)
1572               thread_status=InverseFourierTransformChannel(magnitude_image,
1573                 phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
1574             if (thread_status == MagickFalse)
1575               status=thread_status;
1576           }
1577         }
1578         if (status == MagickFalse)
1579           fourier_image=DestroyImage(fourier_image);
1580       }
1581     fftw_cleanup();
1582   }
1583 #endif
1584   return(fourier_image);
1585 }