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-2011 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/attribute.h"
  47 #include "MagickCore/cache.h"
  48 #include "MagickCore/image.h"
  49 #include "MagickCore/image-private.h"
  50 #include "MagickCore/list.h"
  51 #include "MagickCore/fourier.h"
  52 #include "MagickCore/log.h"
  53 #include "MagickCore/memory_.h"
  54 #include "MagickCore/monitor.h"
  55 #include "MagickCore/pixel-accessor.h"
  56 #include "MagickCore/property.h"
  57 #include "MagickCore/quantum-private.h"
  58 #include "MagickCore/thread-private.h"
  59 #if defined(MAGICKCORE_FFTW_DELEGATE)
  60 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
  61 #include <complex.h>
  62 #endif
  63 #include <fftw3.h>
  64 #if !defined(MAGICKCORE_HAVE_CABS)
  65 #define cabs(z)  (sqrt(z[0]*z[0]+z[1]*z[1]))
  66 #endif
  67 #if !defined(MAGICKCORE_HAVE_CARG)
  68 #define carg(z)  (atan2(cimag(z),creal(z)))
  69 #endif
  70 #if !defined(MAGICKCORE_HAVE_CIMAG)
  71 #define cimag(z)  (z[1])
  72 #endif
  73 #if !defined(MAGICKCORE_HAVE_CREAL)
  74 #define creal(z)  (z[0])
  75 #endif
  76 #endif
  77 \f
  78 /*
  79   Typedef declarations.
  80 */
  81 typedef struct _FourierInfo
  82 {
  83   ChannelType
  84     channel;
  85
  86   MagickBooleanType
  87     modulus;
  88
  89   size_t
  90     width,
  91     height;
  92
  93   ssize_t
  94     center;
  95 } FourierInfo;
  96 \f
  97 /*
  98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  99 %                                                                             %
 100 %                                                                             %
 101 %                                                                             %
 102 %     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                 %
 103 %                                                                             %
 104 %                                                                             %
 105 %                                                                             %
 106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 107 %
 108 %  ForwardFourierTransformImage() implements the discrete Fourier transform
 109 %  (DFT) of the image either as a magnitude / phase or real / imaginary image
 110 %  pair.
 111 %
 112 %  The format of the ForwadFourierTransformImage method is:
 113 %
 114 %      Image *ForwardFourierTransformImage(const Image *image,
 115 %        const MagickBooleanType modulus,ExceptionInfo *exception)
 116 %
 117 %  A description of each parameter follows:
 118 %
 119 %    o image: the image.
 120 %
 121 %    o modulus: if true, return as transform as a magnitude / phase pair
 122 %      otherwise a real / imaginary image pair.
 123 %
 124 %    o exception: return any errors or warnings in this structure.
 125 %
 126 */
 127
 128 #if defined(MAGICKCORE_FFTW_DELEGATE)
 129
 130 static MagickBooleanType RollFourier(const size_t width,const size_t height,
 131   const ssize_t x_offset,const ssize_t y_offset,double *fourier)
 132 {
 133   double
 134     *roll;
 135
 136   register ssize_t
 137     i,
 138     x;
 139
 140   ssize_t
 141     u,
 142     v,
 143     y;
 144
 145   /*
 146     Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
 147   */
 148   roll=(double *) AcquireQuantumMemory((size_t) height,width*sizeof(*roll));
 149   if (roll == (double *) NULL)
 150     return(MagickFalse);
 151   i=0L;
 152   for (y=0L; y < (ssize_t) height; y++)
 153   {
 154     if (y_offset < 0L)
 155       v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
 156     else
 157       v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
 158         y+y_offset;
 159     for (x=0L; x < (ssize_t) width; x++)
 160     {
 161       if (x_offset < 0L)
 162         u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
 163       else
 164         u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
 165           x+x_offset;
 166       roll[v*width+u]=fourier[i++];
 167     }
 168   }
 169   (void) CopyMagickMemory(fourier,roll,height*width*sizeof(*roll));
 170   roll=(double *) RelinquishMagickMemory(roll);
 171   return(MagickTrue);
 172 }
 173
 174 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
 175   const size_t height,double *source,double *destination)
 176 {
 177   MagickBooleanType
 178     status;
 179
 180   register ssize_t
 181     x;
 182
 183   ssize_t
 184     center,
 185     y;
 186
 187   /*
 188     Swap quadrants.
 189   */
 190   center=(ssize_t) floor((double) width/2L)+1L;
 191   status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
 192   if (status == MagickFalse)
 193     return(MagickFalse);
 194   for (y=0L; y < (ssize_t) height; y++)
 195     for (x=0L; x < (ssize_t) (width/2L-1L); x++)
 196       destination[width*y+x+width/2L]=source[center*y+x];
 197   for (y=1; y < (ssize_t) height; y++)
 198     for (x=0L; x < (ssize_t) (width/2L-1L); x++)
 199       destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
 200   for (x=0L; x < (ssize_t) (width/2L); x++)
 201     destination[-x+width/2L-1L]=destination[x+width/2L+1L];
 202   return(MagickTrue);
 203 }
 204
 205 static void CorrectPhaseLHS(const size_t width,const size_t height,
 206   double *fourier)
 207 {
 208   register ssize_t
 209     x;
 210
 211   ssize_t
 212     y;
 213
 214   for (y=0L; y < (ssize_t) height; y++)
 215     for (x=0L; x < (ssize_t) (width/2L); x++)
 216       fourier[y*width+x]*=(-1.0);
 217 }
 218
 219 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
 220   Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
 221 {
 222   CacheView
 223     *magnitude_view,
 224     *phase_view;
 225
 226   double
 227     *magnitude_source,
 228     *phase_source;
 229
 230   Image
 231     *magnitude_image,
 232     *phase_image;
 233
 234   MagickBooleanType
 235     status;
 236
 237   register ssize_t
 238     x;
 239
 240   register Quantum
 241     *q;
 242
 243   ssize_t
 244     i,
 245     y;
 246
 247   magnitude_image=GetFirstImageInList(image);
 248   phase_image=GetNextImageInList(image);
 249   if (phase_image == (Image *) NULL)
 250     {
 251       (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
 252         "ImageSequenceRequired","`%s'",image->filename);
 253       return(MagickFalse);
 254     }
 255   /*
 256     Create "Fourier Transform" image from constituent arrays.
 257   */
 258   magnitude_source=(double *) AcquireQuantumMemory((size_t)
 259     fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
 260   if (magnitude_source == (double *) NULL)
 261     return(MagickFalse);
 262   (void) ResetMagickMemory(magnitude_source,0,fourier_info->height*
 263     fourier_info->width*sizeof(*magnitude_source));
 264   phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 265     fourier_info->width*sizeof(*phase_source));
 266   if (phase_source == (double *) NULL)
 267     {
 268       (void) ThrowMagickException(exception,GetMagickModule(),
 269         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 270       magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 271       return(MagickFalse);
 272     }
 273   status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
 274     magnitude,magnitude_source);
 275   if (status != MagickFalse)
 276     status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
 277       phase_source);
 278   CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
 279   if (fourier_info->modulus != MagickFalse)
 280     {
 281       i=0L;
 282       for (y=0L; y < (ssize_t) fourier_info->height; y++)
 283         for (x=0L; x < (ssize_t) fourier_info->width; x++)
 284         {
 285           phase_source[i]/=(2.0*MagickPI);
 286           phase_source[i]+=0.5;
 287           i++;
 288         }
 289     }
 290   magnitude_view=AcquireCacheView(magnitude_image);
 291   phase_view=AcquireCacheView(phase_image);
 292   i=0L;
 293   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 294   {
 295     q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
 296       exception);
 297     if (q == (const Quantum *) NULL)
 298       break;
 299     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 300     {
 301       switch (fourier_info->channel)
 302       {
 303         case RedChannel:
 304         default:
 305         {
 306           SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
 307             magnitude_source[i]),q);
 308           break;
 309         }
 310         case GreenChannel:
 311         {
 312           SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
 313             magnitude_source[i]),q);
 314           break;
 315         }
 316         case BlueChannel:
 317         {
 318           SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
 319             magnitude_source[i]),q);
 320           break;
 321         }
 322         case BlackChannel:
 323         {
 324           SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
 325             magnitude_source[i]),q);
 326           break;
 327         }
 328         case AlphaChannel:
 329         {
 330           SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
 331             magnitude_source[i]),q);
 332           break;
 333         }
 334         case GrayChannels:
 335         {
 336           SetPixelGray(magnitude_image,ClampToQuantum(QuantumRange*
 337             magnitude_source[i]),q);
 338           break;
 339         }
 340       }
 341       i++;
 342       q+=GetPixelChannels(magnitude_image);
 343     }
 344     status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
 345     if (status == MagickFalse)
 346       break;
 347   }
 348   i=0L;
 349   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 350   {
 351     q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
 352       exception);
 353     if (q == (const Quantum *) NULL)
 354       break;
 355     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 356     {
 357       switch (fourier_info->channel)
 358       {
 359         case RedChannel:
 360         default:
 361         {
 362           SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
 363             phase_source[i]),q);
 364           break;
 365         }
 366         case GreenChannel:
 367         {
 368           SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
 369             phase_source[i]),q);
 370           break;
 371         }
 372         case BlueChannel:
 373         {
 374           SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
 375             phase_source[i]),q);
 376           break;
 377         }
 378         case BlackChannel:
 379         {
 380           SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
 381             phase_source[i]),q);
 382           break;
 383         }
 384         case AlphaChannel:
 385         {
 386           SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
 387             phase_source[i]),q);
 388           break;
 389         }
 390         case GrayChannels:
 391         {
 392           SetPixelGray(phase_image,ClampToQuantum(QuantumRange*
 393             phase_source[i]),q);
 394           break;
 395         }
 396       }
 397       i++;
 398       q+=GetPixelChannels(phase_image);
 399     }
 400     status=SyncCacheViewAuthenticPixels(phase_view,exception);
 401     if (status == MagickFalse)
 402       break;
 403    }
 404   phase_view=DestroyCacheView(phase_view);
 405   magnitude_view=DestroyCacheView(magnitude_view);
 406   phase_source=(double *) RelinquishMagickMemory(phase_source);
 407   magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 408   return(status);
 409 }
 410
 411 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
 412   const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
 413 {
 414   CacheView
 415     *image_view;
 416
 417   double
 418     n,
 419     *source;
 420
 421   fftw_complex
 422     *fourier;
 423
 424   fftw_plan
 425     fftw_r2c_plan;
 426
 427   register const Quantum
 428     *p;
 429
 430   register ssize_t
 431     i,
 432     x;
 433
 434   ssize_t
 435     y;
 436
 437   /*
 438     Generate the forward Fourier transform.
 439   */
 440   source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 441     fourier_info->width*sizeof(*source));
 442   if (source == (double *) NULL)
 443     {
 444       (void) ThrowMagickException(exception,GetMagickModule(),
 445         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 446       return(MagickFalse);
 447     }
 448   ResetMagickMemory(source,0,fourier_info->height*fourier_info->width*
 449     sizeof(*source));
 450   i=0L;
 451   image_view=AcquireCacheView(image);
 452   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 453   {
 454     p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
 455       exception);
 456     if (p == (const Quantum *) NULL)
 457       break;
 458     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 459     {
 460       switch (fourier_info->channel)
 461       {
 462         case RedChannel:
 463         default:
 464         {
 465           source[i]=QuantumScale*GetPixelRed(image,p);
 466           break;
 467         }
 468         case GreenChannel:
 469         {
 470           source[i]=QuantumScale*GetPixelGreen(image,p);
 471           break;
 472         }
 473         case BlueChannel:
 474         {
 475           source[i]=QuantumScale*GetPixelBlue(image,p);
 476           break;
 477         }
 478         case BlackChannel:
 479         {
 480           source[i]=QuantumScale*GetPixelBlack(image,p);
 481           break;
 482         }
 483         case AlphaChannel:
 484         {
 485           source[i]=QuantumScale*GetPixelAlpha(image,p);
 486           break;
 487         }
 488         case GrayChannels:
 489         {
 490           source[i]=QuantumScale*GetPixelGray(image,p);
 491           break;
 492         }
 493       }
 494       i++;
 495       p+=GetPixelChannels(image);
 496     }
 497   }
 498   image_view=DestroyCacheView(image_view);
 499   fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info->height,
 500     fourier_info->center*sizeof(*fourier));
 501   if (fourier == (fftw_complex *) NULL)
 502     {
 503       (void) ThrowMagickException(exception,GetMagickModule(),
 504         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 505       source=(double *) RelinquishMagickMemory(source);
 506       return(MagickFalse);
 507     }
 508 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 509   #pragma omp critical (MagickCore_ForwardFourierTransform)
 510 #endif
 511   fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->width,
 512     source,fourier,FFTW_ESTIMATE);
 513   fftw_execute(fftw_r2c_plan);
 514   fftw_destroy_plan(fftw_r2c_plan);
 515   source=(double *) RelinquishMagickMemory(source);
 516   /*
 517     Normalize Fourier transform.
 518   */
 519   n=(double) fourier_info->width*(double) fourier_info->width;
 520   i=0L;
 521   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 522     for (x=0L; x < (ssize_t) fourier_info->center; x++)
 523     {
 524 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
 525       fourier[i]/=n;
 526 #else
 527       fourier[i][0]/=n;
 528       fourier[i][1]/=n;
 529 #endif
 530       i++;
 531     }
 532   /*
 533     Generate magnitude and phase (or real and imaginary).
 534   */
 535   i=0L;
 536   if (fourier_info->modulus != MagickFalse)
 537     for (y=0L; y < (ssize_t) fourier_info->height; y++)
 538       for (x=0L; x < (ssize_t) fourier_info->center; x++)
 539       {
 540         magnitude[i]=cabs(fourier[i]);
 541         phase[i]=carg(fourier[i]);
 542         i++;
 543       }
 544   else
 545     for (y=0L; y < (ssize_t) fourier_info->height; y++)
 546       for (x=0L; x < (ssize_t) fourier_info->center; x++)
 547       {
 548         magnitude[i]=creal(fourier[i]);
 549         phase[i]=cimag(fourier[i]);
 550         i++;
 551       }
 552   fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
 553   return(MagickTrue);
 554 }
 555
 556 static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
 557   const ChannelType channel,const MagickBooleanType modulus,
 558   Image *fourier_image,ExceptionInfo *exception)
 559 {
 560   double
 561     *magnitude,
 562     *phase;
 563
 564   fftw_complex
 565     *fourier;
 566
 567   FourierInfo
 568     fourier_info;
 569
 570   MagickBooleanType
 571     status;
 572
 573   size_t
 574     extent;
 575
 576   fourier_info.width=image->columns;
 577   if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
 578       ((image->rows % 2) != 0))
 579     {
 580       extent=image->columns < image->rows ? image->rows : image->columns;
 581       fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
 582     }
 583   fourier_info.height=fourier_info.width;
 584   fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
 585   fourier_info.channel=channel;
 586   fourier_info.modulus=modulus;
 587   magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
 588     fourier_info.center*sizeof(*magnitude));
 589   if (magnitude == (double *) NULL)
 590     {
 591       (void) ThrowMagickException(exception,GetMagickModule(),
 592         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 593       return(MagickFalse);
 594     }
 595   phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
 596     fourier_info.center*sizeof(*phase));
 597   if (phase == (double *) NULL)
 598     {
 599       (void) ThrowMagickException(exception,GetMagickModule(),
 600         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 601       magnitude=(double *) RelinquishMagickMemory(magnitude);
 602       return(MagickFalse);
 603     }
 604   fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
 605     fourier_info.center*sizeof(*fourier));
 606   if (fourier == (fftw_complex *) NULL)
 607     {
 608       (void) ThrowMagickException(exception,GetMagickModule(),
 609         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 610       phase=(double *) RelinquishMagickMemory(phase);
 611       magnitude=(double *) RelinquishMagickMemory(magnitude);
 612       return(MagickFalse);
 613     }
 614   status=ForwardFourierTransform(&fourier_info,image,magnitude,phase,exception);
 615   if (status != MagickFalse)
 616     status=ForwardFourier(&fourier_info,fourier_image,magnitude,phase,
 617       exception);
 618   fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
 619   phase=(double *) RelinquishMagickMemory(phase);
 620   magnitude=(double *) RelinquishMagickMemory(magnitude);
 621   return(status);
 622 }
 623 #endif
 624
 625 MagickExport Image *ForwardFourierTransformImage(const Image *image,
 626   const MagickBooleanType modulus,ExceptionInfo *exception)
 627 {
 628   Image
 629     *fourier_image;
 630
 631   fourier_image=NewImageList();
 632 #if !defined(MAGICKCORE_FFTW_DELEGATE)
 633   (void) modulus;
 634   (void) ThrowMagickException(exception,GetMagickModule(),
 635     MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
 636     image->filename);
 637 #else
 638   {
 639     Image
 640       *magnitude_image;
 641
 642     size_t
 643       extent,
 644       width;
 645
 646     width=image->columns;
 647     if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
 648         ((image->rows % 2) != 0))
 649       {
 650         extent=image->columns < image->rows ? image->rows : image->columns;
 651         width=(extent & 0x01) == 1 ? extent+1UL : extent;
 652       }
 653     magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
 654     if (magnitude_image != (Image *) NULL)
 655       {
 656         Image
 657           *phase_image;
 658
 659         magnitude_image->storage_class=DirectClass;
 660         magnitude_image->depth=32UL;
 661         phase_image=CloneImage(image,width,width,MagickFalse,exception);
 662         if (phase_image == (Image *) NULL)
 663           magnitude_image=DestroyImage(magnitude_image);
 664         else
 665           {
 666             MagickBooleanType
 667               is_gray,
 668               status;
 669
 670             phase_image->storage_class=DirectClass;
 671             phase_image->depth=32UL;
 672             AppendImageToList(&fourier_image,magnitude_image);
 673             AppendImageToList(&fourier_image,phase_image);
 674             status=MagickTrue;
 675             is_gray=IsImageGray(image,exception);
 676 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 677   #pragma omp parallel sections
 678 #endif
 679             {
 680 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 681               #pragma omp section
 682 #endif
 683               {
 684                 MagickBooleanType
 685                   thread_status;
 686
 687                 if (is_gray != MagickFalse)
 688                   thread_status=ForwardFourierTransformChannel(image,
 689                     GrayChannels,modulus,fourier_image,exception);
 690                 else
 691                   thread_status=ForwardFourierTransformChannel(image,
 692                     RedChannel,modulus,fourier_image,exception);
 693                 if (thread_status == MagickFalse)
 694                   status=thread_status;
 695               }
 696 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 697               #pragma omp section
 698 #endif
 699               {
 700                 MagickBooleanType
 701                   thread_status;
 702
 703                 thread_status=MagickTrue;
 704                 if (is_gray == MagickFalse)
 705                   thread_status=ForwardFourierTransformChannel(image,
 706                     GreenChannel,modulus,fourier_image,exception);
 707                 if (thread_status == MagickFalse)
 708                   status=thread_status;
 709               }
 710 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 711               #pragma omp section
 712 #endif
 713               {
 714                 MagickBooleanType
 715                   thread_status;
 716
 717                 thread_status=MagickTrue;
 718                 if (is_gray == MagickFalse)
 719                   thread_status=ForwardFourierTransformChannel(image,
 720                     BlueChannel,modulus,fourier_image,exception);
 721                 if (thread_status == MagickFalse)
 722                   status=thread_status;
 723               }
 724 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 725               #pragma omp section
 726 #endif
 727               {
 728                 MagickBooleanType
 729                   thread_status;
 730
 731                 thread_status=MagickTrue;
 732                 if (image->colorspace == CMYKColorspace)
 733                   thread_status=ForwardFourierTransformChannel(image,
 734                     BlackChannel,modulus,fourier_image,exception);
 735                 if (thread_status == MagickFalse)
 736                   status=thread_status;
 737               }
 738 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 739               #pragma omp section
 740 #endif
 741               {
 742                 MagickBooleanType
 743                   thread_status;
 744
 745                 thread_status=MagickTrue;
 746                 if (image->matte != MagickFalse)
 747                   thread_status=ForwardFourierTransformChannel(image,
 748                     AlphaChannel,modulus,fourier_image,exception);
 749                 if (thread_status == MagickFalse)
 750                   status=thread_status;
 751               }
 752             }
 753             if (status == MagickFalse)
 754               fourier_image=DestroyImageList(fourier_image);
 755             fftw_cleanup();
 756           }
 757       }
 758   }
 759 #endif
 760   return(fourier_image);
 761 }
 762 \f
 763 /*
 764 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 765 %                                                                             %
 766 %                                                                             %
 767 %                                                                             %
 768 %     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                 %
 769 %                                                                             %
 770 %                                                                             %
 771 %                                                                             %
 772 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 773 %
 774 %  InverseFourierTransformImage() implements the inverse discrete Fourier
 775 %  transform (DFT) of the image either as a magnitude / phase or real /
 776 %  imaginary image pair.
 777 %
 778 %  The format of the InverseFourierTransformImage method is:
 779 %
 780 %      Image *InverseFourierTransformImage(const Image *magnitude_image,
 781 %        const Image *phase_image,const MagickBooleanType modulus,
 782 %        ExceptionInfo *exception)
 783 %
 784 %  A description of each parameter follows:
 785 %
 786 %    o magnitude_image: the magnitude or real image.
 787 %
 788 %    o phase_image: the phase or imaginary image.
 789 %
 790 %    o modulus: if true, return transform as a magnitude / phase pair
 791 %      otherwise a real / imaginary image pair.
 792 %
 793 %    o exception: return any errors or warnings in this structure.
 794 %
 795 */
 796
 797 #if defined(MAGICKCORE_FFTW_DELEGATE)
 798 static MagickBooleanType InverseQuadrantSwap(const size_t width,
 799   const size_t height,const double *source,double *destination)
 800 {
 801   register ssize_t
 802     x;
 803
 804   ssize_t
 805     center,
 806     y;
 807
 808   /*
 809     Swap quadrants.
 810   */
 811   center=(ssize_t) floor((double) width/2.0)+1L;
 812   for (y=1L; y < (ssize_t) height; y++)
 813     for (x=0L; x < (ssize_t) (width/2L+1L); x++)
 814       destination[center*(height-y)-x+width/2L]=source[y*width+x];
 815   for (y=0L; y < (ssize_t) height; y++)
 816     destination[center*y]=source[y*width+width/2L];
 817   for (x=0L; x < center; x++)
 818     destination[x]=source[center-x-1L];
 819   return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
 820 }
 821
 822 static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
 823   const Image *magnitude_image,const Image *phase_image,fftw_complex *fourier,
 824   ExceptionInfo *exception)
 825 {
 826   CacheView
 827     *magnitude_view,
 828     *phase_view;
 829
 830   double
 831     *magnitude,
 832     *phase,
 833     *magnitude_source,
 834     *phase_source;
 835
 836   MagickBooleanType
 837     status;
 838
 839   register const Quantum
 840     *p;
 841
 842   register ssize_t
 843     i,
 844     x;
 845
 846   ssize_t
 847     y;
 848
 849   /*
 850     Inverse fourier - read image and break down into a double array.
 851   */
 852   magnitude_source=(double *) AcquireQuantumMemory((size_t)
 853     fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
 854   if (magnitude_source == (double *) NULL)
 855     {
 856       (void) ThrowMagickException(exception,GetMagickModule(),
 857         ResourceLimitError,"MemoryAllocationFailed","`%s'",
 858         magnitude_image->filename);
 859       return(MagickFalse);
 860     }
 861   phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 862     fourier_info->width*sizeof(*phase_source));
 863   if (phase_source == (double *) NULL)
 864     {
 865       (void) ThrowMagickException(exception,GetMagickModule(),
 866         ResourceLimitError,"MemoryAllocationFailed","`%s'",
 867         magnitude_image->filename);
 868       magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 869       return(MagickFalse);
 870     }
 871   i=0L;
 872   magnitude_view=AcquireCacheView(magnitude_image);
 873   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 874   {
 875     p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
 876       exception);
 877     if (p == (const Quantum *) NULL)
 878       break;
 879     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 880     {
 881       switch (fourier_info->channel)
 882       {
 883         case RedChannel:
 884         default:
 885         {
 886           magnitude_source[i]=QuantumScale*GetPixelRed(magnitude_image,p);
 887           break;
 888         }
 889         case GreenChannel:
 890         {
 891           magnitude_source[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
 892           break;
 893         }
 894         case BlueChannel:
 895         {
 896           magnitude_source[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
 897           break;
 898         }
 899         case BlackChannel:
 900         {
 901           magnitude_source[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
 902           break;
 903         }
 904         case AlphaChannel:
 905         {
 906           magnitude_source[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
 907           break;
 908         }
 909         case GrayChannels:
 910         {
 911           magnitude_source[i]=QuantumScale*GetPixelGray(magnitude_image,p);
 912           break;
 913         }
 914       }
 915       i++;
 916       p+=GetPixelChannels(magnitude_image);
 917     }
 918   }
 919   i=0L;
 920   phase_view=AcquireCacheView(phase_image);
 921   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 922   {
 923     p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
 924       exception);
 925     if (p == (const Quantum *) NULL)
 926       break;
 927     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 928     {
 929       switch (fourier_info->channel)
 930       {
 931         case RedChannel:
 932         default:
 933         {
 934           phase_source[i]=QuantumScale*GetPixelRed(phase_image,p);
 935           break;
 936         }
 937         case GreenChannel:
 938         {
 939           phase_source[i]=QuantumScale*GetPixelGreen(phase_image,p);
 940           break;
 941         }
 942         case BlueChannel:
 943         {
 944           phase_source[i]=QuantumScale*GetPixelBlue(phase_image,p);
 945           break;
 946         }
 947         case BlackChannel:
 948         {
 949           phase_source[i]=QuantumScale*GetPixelBlack(phase_image,p);
 950           break;
 951         }
 952         case AlphaChannel:
 953         {
 954           phase_source[i]=QuantumScale*GetPixelAlpha(phase_image,p);
 955           break;
 956         }
 957         case GrayChannels:
 958         {
 959           phase_source[i]=QuantumScale*GetPixelGray(phase_image,p);
 960           break;
 961         }
 962       }
 963       i++;
 964       p+=GetPixelChannels(phase_image);
 965     }
 966   }
 967   if (fourier_info->modulus != MagickFalse)
 968     {
 969       i=0L;
 970       for (y=0L; y < (ssize_t) fourier_info->height; y++)
 971         for (x=0L; x < (ssize_t) fourier_info->width; x++)
 972         {
 973           phase_source[i]-=0.5;
 974           phase_source[i]*=(2.0*MagickPI);
 975           i++;
 976         }
 977     }
 978   magnitude_view=DestroyCacheView(magnitude_view);
 979   phase_view=DestroyCacheView(phase_view);
 980   magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 981     fourier_info->center*sizeof(*magnitude));
 982   if (magnitude == (double *) NULL)
 983     {
 984       (void) ThrowMagickException(exception,GetMagickModule(),
 985         ResourceLimitError,"MemoryAllocationFailed","`%s'",
 986         magnitude_image->filename);
 987       magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 988       phase_source=(double *) RelinquishMagickMemory(phase_source);
 989       return(MagickFalse);
 990     }
 991   status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
 992     magnitude_source,magnitude);
 993   magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 994   phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 995     fourier_info->width*sizeof(*phase));
 996   if (phase == (double *) NULL)
 997     {
 998       (void) ThrowMagickException(exception,GetMagickModule(),
 999         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1000         magnitude_image->filename);
1001       phase_source=(double *) RelinquishMagickMemory(phase_source);
1002       return(MagickFalse);
1003     }
1004   CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1005   if (status != MagickFalse)
1006     status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1007       phase_source,phase);
1008   phase_source=(double *) RelinquishMagickMemory(phase_source);
1009   /*
1010     Merge two sets.
1011   */
1012   i=0L;
1013   if (fourier_info->modulus != MagickFalse)
1014     for (y=0L; y < (ssize_t) fourier_info->height; y++)
1015        for (x=0L; x < (ssize_t) fourier_info->center; x++)
1016        {
1017 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1018          fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
1019 #else
1020          fourier[i][0]=magnitude[i]*cos(phase[i]);
1021          fourier[i][1]=magnitude[i]*sin(phase[i]);
1022 #endif
1023          i++;
1024       }
1025   else
1026     for (y=0L; y < (ssize_t) fourier_info->height; y++)
1027       for (x=0L; x < (ssize_t) fourier_info->center; x++)
1028       {
1029 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1030         fourier[i]=magnitude[i]+I*phase[i];
1031 #else
1032         fourier[i][0]=magnitude[i];
1033         fourier[i][1]=phase[i];
1034 #endif
1035         i++;
1036       }
1037   phase=(double *) RelinquishMagickMemory(phase);
1038   magnitude=(double *) RelinquishMagickMemory(magnitude);
1039   return(status);
1040 }
1041
1042 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1043   fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1044 {
1045   CacheView
1046     *image_view;
1047
1048   double
1049     *source;
1050
1051   fftw_plan
1052     fftw_c2r_plan;
1053
1054   register Quantum
1055     *q;
1056
1057   register ssize_t
1058     i,
1059     x;
1060
1061   ssize_t
1062     y;
1063
1064   source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
1065     fourier_info->width*sizeof(*source));
1066   if (source == (double *) NULL)
1067     {
1068       (void) ThrowMagickException(exception,GetMagickModule(),
1069         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1070       return(MagickFalse);
1071     }
1072 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1073   #pragma omp critical (MagickCore_InverseFourierTransform)
1074 #endif
1075   {
1076     fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1077       fourier,source,FFTW_ESTIMATE);
1078     fftw_execute(fftw_c2r_plan);
1079     fftw_destroy_plan(fftw_c2r_plan);
1080   }
1081   i=0L;
1082   image_view=AcquireCacheView(image);
1083   for (y=0L; y < (ssize_t) fourier_info->height; y++)
1084   {
1085     if (y >= (ssize_t) image->rows)
1086       break;
1087     q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1088       image->columns ? image->columns : fourier_info->width,1UL,exception);
1089     if (q == (const Quantum *) NULL)
1090       break;
1091     for (x=0L; x < (ssize_t) fourier_info->width; x++)
1092     {
1093       switch (fourier_info->channel)
1094       {
1095         case RedChannel:
1096         default:
1097         {
1098           SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
1099           break;
1100         }
1101         case GreenChannel:
1102         {
1103           SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
1104           break;
1105         }
1106         case BlueChannel:
1107         {
1108           SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
1109           break;
1110         }
1111         case BlackChannel:
1112         {
1113           SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
1114           break;
1115         }
1116         case AlphaChannel:
1117         {
1118           SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
1119           break;
1120         }
1121         case GrayChannels:
1122         {
1123           SetPixelGray(image,ClampToQuantum(QuantumRange*source[i]),q);
1124           break;
1125         }
1126       }
1127       i++;
1128       q+=GetPixelChannels(image);
1129     }
1130     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1131       break;
1132   }
1133   image_view=DestroyCacheView(image_view);
1134   source=(double *) RelinquishMagickMemory(source);
1135   return(MagickTrue);
1136 }
1137
1138 static MagickBooleanType InverseFourierTransformChannel(
1139   const Image *magnitude_image,const Image *phase_image,
1140   const ChannelType channel,const MagickBooleanType modulus,
1141   Image *fourier_image,ExceptionInfo *exception)
1142 {
1143   double
1144     *magnitude,
1145     *phase;
1146
1147   fftw_complex
1148     *fourier;
1149
1150   FourierInfo
1151     fourier_info;
1152
1153   MagickBooleanType
1154     status;
1155
1156   size_t
1157     extent;
1158
1159   fourier_info.width=magnitude_image->columns;
1160   if ((magnitude_image->columns != magnitude_image->rows) ||
1161       ((magnitude_image->columns % 2) != 0) ||
1162       ((magnitude_image->rows % 2) != 0))
1163     {
1164       extent=magnitude_image->columns < magnitude_image->rows ?
1165         magnitude_image->rows : magnitude_image->columns;
1166       fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1167     }
1168   fourier_info.height=fourier_info.width;
1169   fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
1170   fourier_info.channel=channel;
1171   fourier_info.modulus=modulus;
1172   magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1173     fourier_info.center*sizeof(*magnitude));
1174   if (magnitude == (double *) NULL)
1175     {
1176       (void) ThrowMagickException(exception,GetMagickModule(),
1177         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1178         magnitude_image->filename);
1179       return(MagickFalse);
1180     }
1181   phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1182     fourier_info.center*sizeof(*phase));
1183   if (phase == (double *) NULL)
1184     {
1185       (void) ThrowMagickException(exception,GetMagickModule(),
1186         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1187         magnitude_image->filename);
1188       magnitude=(double *) RelinquishMagickMemory(magnitude);
1189       return(MagickFalse);
1190     }
1191   fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1192     fourier_info.center*sizeof(*fourier));
1193   if (fourier == (fftw_complex *) NULL)
1194     {
1195       (void) ThrowMagickException(exception,GetMagickModule(),
1196         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1197         magnitude_image->filename);
1198       phase=(double *) RelinquishMagickMemory(phase);
1199       magnitude=(double *) RelinquishMagickMemory(magnitude);
1200       return(MagickFalse);
1201     }
1202   status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1203    exception);
1204   if (status != MagickFalse)
1205     status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1206       exception);
1207   fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1208   phase=(double *) RelinquishMagickMemory(phase);
1209   magnitude=(double *) RelinquishMagickMemory(magnitude);
1210   return(status);
1211 }
1212 #endif
1213
1214 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1215   const Image *phase_image,const MagickBooleanType modulus,
1216   ExceptionInfo *exception)
1217 {
1218   Image
1219     *fourier_image;
1220
1221   assert(magnitude_image != (Image *) NULL);
1222   assert(magnitude_image->signature == MagickSignature);
1223   if (magnitude_image->debug != MagickFalse)
1224     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1225       magnitude_image->filename);
1226   if (phase_image == (Image *) NULL)
1227     {
1228       (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1229         "ImageSequenceRequired","`%s'",magnitude_image->filename);
1230       return((Image *) NULL);
1231     }
1232 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1233   fourier_image=(Image *) NULL;
1234   (void) modulus;
1235   (void) ThrowMagickException(exception,GetMagickModule(),
1236     MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1237     magnitude_image->filename);
1238 #else
1239   {
1240     fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1241       magnitude_image->rows,MagickFalse,exception);
1242     if (fourier_image != (Image *) NULL)
1243       {
1244         MagickBooleanType
1245           is_gray,
1246           status;
1247
1248         status=MagickTrue;
1249         is_gray=IsImageGray(magnitude_image,exception);
1250         if (is_gray != MagickFalse)
1251           is_gray=IsImageGray(phase_image,exception);
1252 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1253         #pragma omp parallel sections
1254 #endif
1255         {
1256 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1257           #pragma omp section
1258 #endif
1259           {
1260             MagickBooleanType
1261               thread_status;
1262
1263             if (is_gray != MagickFalse)
1264               thread_status=InverseFourierTransformChannel(magnitude_image,
1265                 phase_image,GrayChannels,modulus,fourier_image,exception);
1266             else
1267               thread_status=InverseFourierTransformChannel(magnitude_image,
1268                 phase_image,RedChannel,modulus,fourier_image,exception);
1269             if (thread_status == MagickFalse)
1270               status=thread_status;
1271           }
1272 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1273           #pragma omp section
1274 #endif
1275           {
1276             MagickBooleanType
1277               thread_status;
1278
1279             thread_status=MagickTrue;
1280             if (is_gray == MagickFalse)
1281               thread_status=InverseFourierTransformChannel(magnitude_image,
1282                 phase_image,GreenChannel,modulus,fourier_image,exception);
1283             if (thread_status == MagickFalse)
1284               status=thread_status;
1285           }
1286 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1287           #pragma omp section
1288 #endif
1289           {
1290             MagickBooleanType
1291               thread_status;
1292
1293             thread_status=MagickTrue;
1294             if (is_gray == MagickFalse)
1295               thread_status=InverseFourierTransformChannel(magnitude_image,
1296                 phase_image,BlueChannel,modulus,fourier_image,exception);
1297             if (thread_status == MagickFalse)
1298               status=thread_status;
1299           }
1300 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1301           #pragma omp section
1302 #endif
1303           {
1304             MagickBooleanType
1305               thread_status;
1306
1307             thread_status=MagickTrue;
1308             if (magnitude_image->colorspace == CMYKColorspace)
1309               thread_status=InverseFourierTransformChannel(magnitude_image,
1310                 phase_image,BlackChannel,modulus,fourier_image,exception);
1311             if (thread_status == MagickFalse)
1312               status=thread_status;
1313           }
1314 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1315           #pragma omp section
1316 #endif
1317           {
1318             MagickBooleanType
1319               thread_status;
1320
1321             thread_status=MagickTrue;
1322             if (magnitude_image->matte != MagickFalse)
1323               thread_status=InverseFourierTransformChannel(magnitude_image,
1324                 phase_image,AlphaChannel,modulus,fourier_image,exception);
1325             if (thread_status == MagickFalse)
1326               status=thread_status;
1327           }
1328         }
1329         if (status == MagickFalse)
1330           fourier_image=DestroyImage(fourier_image);
1331       }
1332     fftw_cleanup();
1333   }
1334 #endif
1335   return(fourier_image);
1336 }