granicus.if.org Git - imagemagick/blob - MagickCore/fourier.c

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