granicus.if.org Git - imagemagick/blob - magick/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 "magick/studio.h"
  46 #include "magick/cache.h"
  47 #include "magick/image.h"
  48 #include "magick/image-private.h"
  49 #include "magick/list.h"
  50 #include "magick/fourier.h"
  51 #include "magick/log.h"
  52 #include "magick/memory_.h"
  53 #include "magick/monitor.h"
  54 #include "magick/property.h"
  55 #include "magick/quantum-private.h"
  56 #include "magick/thread-private.h"
  57 #if defined(MAGICKCORE_FFTW_DELEGATE)
  58 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
  59 #include <complex.h>
  60 #endif
  61 #include <fftw3.h>
  62 #if !defined(MAGICKCORE_HAVE_CABS)
  63 #define cabs(z)  (sqrt(z[0]*z[0]+z[1]*z[1]))
  64 #endif
  65 #if !defined(MAGICKCORE_HAVE_CARG)
  66 #define carg(z)  (atan2(cimag(z[1]),creal(z[0])))
  67 #endif
  68 #if !defined(MAGICKCORE_HAVE_CIMAG)
  69 #define cimag(z)  (z[1])
  70 #endif
  71 #if !defined(MAGICKCORE_HAVE_CREAL)
  72 #define creal(z)  (z[0])
  73 #endif
  74 #endif
  75 \f
  76 /*
  77   Typedef declarations.
  78 */
  79 typedef struct _FourierInfo
  80 {
  81   ChannelType
  82     channel;
  83
  84   MagickBooleanType
  85     modulus;
  86
  87   size_t
  88     width,
  89     height;
  90
  91   ssize_t
  92     center;
  93 } FourierInfo;
  94 \f
  95 /*
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %                                                                             %
  98 %                                                                             %
  99 %                                                                             %
 100 %     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                 %
 101 %                                                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 105 %
 106 %  ForwardFourierTransformImage() implements the discrete Fourier transform
 107 %  (DFT) of the image either as a magnitude / phase or real / imaginary image
 108 %  pair.
 109 %
 110 %  The format of the ForwadFourierTransformImage method is:
 111 %
 112 %      Image *ForwardFourierTransformImage(const Image *image,
 113 %        const MagickBooleanType modulus,ExceptionInfo *exception)
 114 %
 115 %  A description of each parameter follows:
 116 %
 117 %    o image: the image.
 118 %
 119 %    o modulus: if true, return as transform as a magnitude / phase pair
 120 %      otherwise a real / imaginary image pair.
 121 %
 122 %    o exception: return any errors or warnings in this structure.
 123 %
 124 */
 125
 126 #if defined(MAGICKCORE_FFTW_DELEGATE)
 127
 128 static MagickBooleanType RollFourier(const size_t width,
 129   const size_t height,const ssize_t x_offset,const ssize_t y_offset,
 130   double *fourier)
 131 {
 132   double
 133     *roll;
 134
 135   ssize_t
 136     u,
 137     v,
 138     y;
 139
 140   register ssize_t
 141     i,
 142     x;
 143
 144   /*
 145     Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
 146   */
 147   roll=(double *) AcquireQuantumMemory((size_t) width,height*sizeof(*roll));
 148   if (roll == (double *) NULL)
 149     return(MagickFalse);
 150   i=0L;
 151   for (y=0L; y < (ssize_t) height; y++)
 152   {
 153     if (y_offset < 0L)
 154       v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
 155     else
 156       v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
 157         y+y_offset;
 158     for (x=0L; x < (ssize_t) width; x++)
 159     {
 160       if (x_offset < 0L)
 161         u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
 162       else
 163         u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
 164           x+x_offset;
 165       roll[v*width+u]=fourier[i++];
 166    }
 167   }
 168   (void) CopyMagickMemory(fourier,roll,width*height*sizeof(*roll));
 169   roll=(double *) RelinquishMagickMemory(roll);
 170   return(MagickTrue);
 171 }
 172
 173 static MagickBooleanType ForwardQuadrantSwap(const size_t width,
 174   const size_t height,double *source,double *destination)
 175 {
 176   ssize_t
 177     center,
 178     y;
 179
 180   MagickBooleanType
 181     status;
 182
 183   register ssize_t
 184     x;
 185
 186   /*
 187     Swap quadrants.
 188   */
 189   center=(ssize_t) floor((double) width/2L)+1L;
 190   status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
 191   if (status == MagickFalse)
 192     return(MagickFalse);
 193   for (y=0L; y < (ssize_t) height; y++)
 194     for (x=0L; x < (ssize_t) (width/2L-1L); x++)
 195       destination[width*y+x+width/2L]=source[center*y+x];
 196   for (y=1; y < (ssize_t) height; y++)
 197     for (x=0L; x < (ssize_t) (width/2L-1L); x++)
 198       destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
 199   for (x=0L; x < (ssize_t) (width/2L); x++)
 200     destination[-x+width/2L-1L]=destination[x+width/2L+1L];
 201   return(MagickTrue);
 202 }
 203
 204 static void CorrectPhaseLHS(const size_t width,
 205   const size_t height,double *fourier)
 206 {
 207   ssize_t
 208     y;
 209
 210   register ssize_t
 211     x;
 212
 213   for (y=0L; y < (ssize_t) height; y++)
 214     for (x=0L; x < (ssize_t) (width/2L); x++)
 215       fourier[y*width+x]*=(-1.0);
 216 }
 217
 218 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
 219   Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
 220 {
 221   CacheView
 222     *magnitude_view,
 223     *phase_view;
 224
 225   double
 226     *magnitude_source,
 227     *phase_source;
 228
 229   Image
 230     *magnitude_image,
 231     *phase_image;
 232
 233   ssize_t
 234     i,
 235     y;
 236
 237   MagickBooleanType
 238     status;
 239
 240   register IndexPacket
 241     *indexes;
 242
 243   register ssize_t
 244     x;
 245
 246   register PixelPacket
 247     *q;
 248
 249   magnitude_image=GetFirstImageInList(image);
 250   phase_image=GetNextImageInList(image);
 251   if (phase_image == (Image *) NULL)
 252     {
 253       (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
 254         "ImageSequenceRequired","`%s'",image->filename);
 255       return(MagickFalse);
 256     }
 257   /*
 258     Create "Fourier Transform" image from constituent arrays.
 259   */
 260   magnitude_source=(double *) AcquireQuantumMemory((size_t)
 261     fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
 262   if (magnitude_source == (double *) NULL)
 263     return(MagickFalse);
 264   (void) ResetMagickMemory(magnitude_source,0,fourier_info->width*
 265     fourier_info->height*sizeof(*magnitude_source));
 266   phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 267     fourier_info->width*sizeof(*phase_source));
 268   if (magnitude_source == (double *) NULL)
 269     {
 270       (void) ThrowMagickException(exception,GetMagickModule(),
 271         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 272       magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 273       return(MagickFalse);
 274     }
 275   status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
 276     magnitude,magnitude_source);
 277   if (status != MagickFalse)
 278     status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
 279       phase_source);
 280   CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_source);
 281   if (fourier_info->modulus != MagickFalse)
 282     {
 283       i=0L;
 284       for (y=0L; y < (ssize_t) fourier_info->height; y++)
 285         for (x=0L; x < (ssize_t) fourier_info->width; x++)
 286         {
 287           phase_source[i]/=(2.0*MagickPI);
 288           phase_source[i]+=0.5;
 289           i++;
 290         }
 291     }
 292   magnitude_view=AcquireCacheView(magnitude_image);
 293   phase_view=AcquireCacheView(phase_image);
 294   i=0L;
 295   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 296   {
 297     q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->height,1UL,
 298       exception);
 299     if (q == (PixelPacket *) NULL)
 300       break;
 301     indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
 302     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 303     {
 304       switch (fourier_info->channel)
 305       {
 306         case RedChannel:
 307         default:
 308         {
 309           q->red=ClampToQuantum(QuantumRange*magnitude_source[i]);
 310           break;
 311         }
 312         case GreenChannel:
 313         {
 314           q->green=ClampToQuantum(QuantumRange*magnitude_source[i]);
 315           break;
 316         }
 317         case BlueChannel:
 318         {
 319           q->blue=ClampToQuantum(QuantumRange*magnitude_source[i]);
 320           break;
 321         }
 322         case OpacityChannel:
 323         {
 324           q->opacity=ClampToQuantum(QuantumRange*magnitude_source[i]);
 325           break;
 326         }
 327         case IndexChannel:
 328         {
 329           indexes[x]=ClampToQuantum(QuantumRange*magnitude_source[i]);
 330           break;
 331         }
 332         case GrayChannels:
 333         {
 334           q->red=ClampToQuantum(QuantumRange*magnitude_source[i]);
 335           q->green=q->red;
 336           q->blue=q->red;
 337           break;
 338         }
 339       }
 340       i++;
 341       q++;
 342     }
 343     status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
 344     if (status == MagickFalse)
 345       break;
 346   }
 347   i=0L;
 348   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 349   {
 350     q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->height,1UL,
 351       exception);
 352     if (q == (PixelPacket *) NULL)
 353       break;
 354     indexes=GetCacheViewAuthenticIndexQueue(phase_view);
 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           q->red=ClampToQuantum(QuantumRange*phase_source[i]);
 363           break;
 364         }
 365         case GreenChannel:
 366         {
 367           q->green=ClampToQuantum(QuantumRange*phase_source[i]);
 368           break;
 369         }
 370         case BlueChannel:
 371         {
 372           q->blue=ClampToQuantum(QuantumRange*phase_source[i]);
 373           break;
 374         }
 375         case OpacityChannel:
 376         {
 377           q->opacity=ClampToQuantum(QuantumRange*phase_source[i]);
 378           break;
 379         }
 380         case IndexChannel:
 381         {
 382           indexes[x]=ClampToQuantum(QuantumRange*phase_source[i]);
 383           break;
 384         }
 385         case GrayChannels:
 386         {
 387           q->red=ClampToQuantum(QuantumRange*phase_source[i]);
 388           q->green=q->red;
 389           q->blue=q->red;
 390           break;
 391         }
 392       }
 393       i++;
 394       q++;
 395     }
 396     status=SyncCacheViewAuthenticPixels(phase_view,exception);
 397     if (status == MagickFalse)
 398       break;
 399    }
 400   phase_view=DestroyCacheView(phase_view);
 401   magnitude_view=DestroyCacheView(magnitude_view);
 402   phase_source=(double *) RelinquishMagickMemory(phase_source);
 403   magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 404   return(status);
 405 }
 406
 407 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
 408   const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
 409 {
 410   CacheView
 411     *image_view;
 412
 413   double
 414     n,
 415     *source;
 416
 417   fftw_complex
 418     *fourier;
 419
 420   fftw_plan
 421     fftw_r2c_plan;
 422
 423   ssize_t
 424     y;
 425
 426   register const IndexPacket
 427     *indexes;
 428
 429   register const PixelPacket
 430     *p;
 431
 432   register ssize_t
 433     i,
 434     x;
 435
 436   /*
 437     Generate the forward Fourier transform.
 438   */
 439   source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 440     fourier_info->width*sizeof(*source));
 441   if (source == (double *) NULL)
 442     {
 443       (void) ThrowMagickException(exception,GetMagickModule(),
 444         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 445       return(MagickFalse);
 446     }
 447   ResetMagickMemory(source,0,fourier_info->width*fourier_info->height*
 448     sizeof(*source));
 449   i=0L;
 450   image_view=AcquireCacheView(image);
 451   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 452   {
 453     p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
 454       exception);
 455     if (p == (const PixelPacket *) NULL)
 456       break;
 457     indexes=GetCacheViewVirtualIndexQueue(image_view);
 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*GetRedPixelComponent(p);
 466           break;
 467         }
 468         case GreenChannel:
 469         {
 470           source[i]=QuantumScale*GetGreenPixelComponent(p);
 471           break;
 472         }
 473         case BlueChannel:
 474         {
 475           source[i]=QuantumScale*GetBluePixelComponent(p);
 476           break;
 477         }
 478         case OpacityChannel:
 479         {
 480           source[i]=QuantumScale*GetOpacityPixelComponent(p);
 481           break;
 482         }
 483         case IndexChannel:
 484         {
 485           source[i]=QuantumScale*indexes[x];
 486           break;
 487         }
 488         case GrayChannels:
 489         {
 490           source[i]=QuantumScale*GetRedPixelComponent(p);
 491           break;
 492         }
 493       }
 494       i++;
 495       p++;
 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   MagickBooleanType
 568     status;
 569
 570   FourierInfo
 571     fourier_info;
 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=IsGrayImage(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->matte != MagickFalse)
 733                   thread_status=ForwardFourierTransformChannel(image,
 734                     OpacityChannel,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->colorspace == CMYKColorspace)
 747                   thread_status=ForwardFourierTransformChannel(image,
 748                     IndexChannel,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   ssize_t
 802     center,
 803     y;
 804
 805   register ssize_t
 806     x;
 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   ssize_t
 837     y;
 838
 839   MagickBooleanType
 840     status;
 841
 842   register const IndexPacket
 843     *indexes;
 844
 845   register const PixelPacket
 846     *p;
 847
 848   register ssize_t
 849     i,
 850     x;
 851
 852   /*
 853     Inverse fourier - read image and break down into a double array.
 854   */
 855   magnitude_source=(double *) AcquireQuantumMemory((size_t)
 856     fourier_info->height,fourier_info->width*sizeof(*magnitude_source));
 857   if (magnitude_source == (double *) NULL)
 858     {
 859       (void) ThrowMagickException(exception,GetMagickModule(),
 860         ResourceLimitError,"MemoryAllocationFailed","`%s'",
 861         magnitude_image->filename);
 862       return(MagickFalse);
 863     }
 864   phase_source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 865     fourier_info->height*sizeof(*phase_source));
 866   if (phase_source == (double *) NULL)
 867     {
 868       (void) ThrowMagickException(exception,GetMagickModule(),
 869         ResourceLimitError,"MemoryAllocationFailed","`%s'",
 870         magnitude_image->filename);
 871       magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 872       return(MagickFalse);
 873     }
 874   i=0L;
 875   magnitude_view=AcquireCacheView(magnitude_image);
 876   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 877   {
 878     p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
 879       exception);
 880     if (p == (const PixelPacket *) NULL)
 881       break;
 882     indexes=GetCacheViewAuthenticIndexQueue(magnitude_view);
 883     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 884     {
 885       switch (fourier_info->channel)
 886       {
 887         case RedChannel:
 888         default:
 889         {
 890           magnitude_source[i]=QuantumScale*GetRedPixelComponent(p);
 891           break;
 892         }
 893         case GreenChannel:
 894         {
 895           magnitude_source[i]=QuantumScale*GetGreenPixelComponent(p);
 896           break;
 897         }
 898         case BlueChannel:
 899         {
 900           magnitude_source[i]=QuantumScale*GetBluePixelComponent(p);
 901           break;
 902         }
 903         case OpacityChannel:
 904         {
 905           magnitude_source[i]=QuantumScale*GetOpacityPixelComponent(p);
 906           break;
 907         }
 908         case IndexChannel:
 909         {
 910           magnitude_source[i]=QuantumScale*indexes[x];
 911           break;
 912         }
 913         case GrayChannels:
 914         {
 915           magnitude_source[i]=QuantumScale*GetRedPixelComponent(p);
 916           break;
 917         }
 918       }
 919       i++;
 920       p++;
 921     }
 922   }
 923   i=0L;
 924   phase_view=AcquireCacheView(phase_image);
 925   for (y=0L; y < (ssize_t) fourier_info->height; y++)
 926   {
 927     p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
 928       exception);
 929     if (p == (const PixelPacket *) NULL)
 930       break;
 931     indexes=GetCacheViewAuthenticIndexQueue(phase_view);
 932     for (x=0L; x < (ssize_t) fourier_info->width; x++)
 933     {
 934       switch (fourier_info->channel)
 935       {
 936         case RedChannel:
 937         default:
 938         {
 939           phase_source[i]=QuantumScale*GetRedPixelComponent(p);
 940           break;
 941         }
 942         case GreenChannel:
 943         {
 944           phase_source[i]=QuantumScale*GetGreenPixelComponent(p);
 945           break;
 946         }
 947         case BlueChannel:
 948         {
 949           phase_source[i]=QuantumScale*GetBluePixelComponent(p);
 950           break;
 951         }
 952         case OpacityChannel:
 953         {
 954           phase_source[i]=QuantumScale*GetOpacityPixelComponent(p);
 955           break;
 956         }
 957         case IndexChannel:
 958         {
 959           phase_source[i]=QuantumScale*indexes[x];
 960           break;
 961         }
 962         case GrayChannels:
 963         {
 964           phase_source[i]=QuantumScale*GetRedPixelComponent(p);
 965           break;
 966         }
 967       }
 968       i++;
 969       p++;
 970     }
 971   }
 972   if (fourier_info->modulus != MagickFalse)
 973     {
 974       i=0L;
 975       for (y=0L; y < (ssize_t) fourier_info->height; y++)
 976         for (x=0L; x < (ssize_t) fourier_info->width; x++)
 977         {
 978           phase_source[i]-=0.5;
 979           phase_source[i]*=(2.0*MagickPI);
 980           i++;
 981         }
 982     }
 983   magnitude_view=DestroyCacheView(magnitude_view);
 984   phase_view=DestroyCacheView(phase_view);
 985   magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
 986     fourier_info->center*sizeof(*magnitude));
 987   if (magnitude == (double *) NULL)
 988     {
 989       (void) ThrowMagickException(exception,GetMagickModule(),
 990         ResourceLimitError,"MemoryAllocationFailed","`%s'",
 991         magnitude_image->filename);
 992       magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 993       phase_source=(double *) RelinquishMagickMemory(phase_source);
 994       return(MagickFalse);
 995     }
 996   status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
 997     magnitude_source,magnitude);
 998   magnitude_source=(double *) RelinquishMagickMemory(magnitude_source);
 999   phase=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
1000     fourier_info->height*sizeof(*phase));
1001   if (phase == (double *) NULL)
1002     {
1003       (void) ThrowMagickException(exception,GetMagickModule(),
1004         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1005         magnitude_image->filename);
1006       phase_source=(double *) RelinquishMagickMemory(phase_source);
1007       return(MagickFalse);
1008     }
1009   CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_source);
1010   if (status != MagickFalse)
1011     status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
1012       phase_source,phase);
1013   phase_source=(double *) RelinquishMagickMemory(phase_source);
1014   /*
1015     Merge two sets.
1016   */
1017   i=0L;
1018   if (fourier_info->modulus != MagickFalse)
1019     for (y=0L; y < (ssize_t) fourier_info->height; y++)
1020        for (x=0L; x < (ssize_t) fourier_info->center; x++)
1021        {
1022 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1023          fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
1024 #else
1025          fourier[i][0]=magnitude[i]*cos(phase[i]);
1026          fourier[i][1]=magnitude[i]*sin(phase[i]);
1027 #endif
1028          i++;
1029       }
1030   else
1031     for (y=0L; y < (ssize_t) fourier_info->height; y++)
1032       for (x=0L; x < (ssize_t) fourier_info->center; x++)
1033       {
1034 #if defined(MAGICKCORE_HAVE_COMPLEX_H)
1035         fourier[i]=magnitude[i]+I*phase[i];
1036 #else
1037         fourier[i][0]=magnitude[i];
1038         fourier[i][1]=phase[i];
1039 #endif
1040         i++;
1041       }
1042   phase=(double *) RelinquishMagickMemory(phase);
1043   magnitude=(double *) RelinquishMagickMemory(magnitude);
1044   return(status);
1045 }
1046
1047 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
1048   fftw_complex *fourier,Image *image,ExceptionInfo *exception)
1049 {
1050   CacheView
1051     *image_view;
1052
1053   double
1054     *source;
1055
1056   fftw_plan
1057     fftw_c2r_plan;
1058
1059   ssize_t
1060     y;
1061
1062   register IndexPacket
1063     *indexes;
1064
1065   register ssize_t
1066     i,
1067     x;
1068
1069   register PixelPacket
1070     *q;
1071
1072   source=(double *) AcquireQuantumMemory((size_t) fourier_info->width,
1073     fourier_info->height*sizeof(double));
1074   if (source == (double *) NULL)
1075     {
1076       (void) ThrowMagickException(exception,GetMagickModule(),
1077         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1078       return(MagickFalse);
1079     }
1080 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1081   #pragma omp critical (MagickCore_InverseFourierTransform)
1082 #endif
1083   {
1084     fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
1085       fourier,source,FFTW_ESTIMATE);
1086     fftw_execute(fftw_c2r_plan);
1087     fftw_destroy_plan(fftw_c2r_plan);
1088   }
1089   i=0L;
1090   image_view=AcquireCacheView(image);
1091   for (y=0L; y < (ssize_t) fourier_info->height; y++)
1092   {
1093     if (y >= (ssize_t) image->rows)
1094       break;
1095     q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
1096       image->columns ? image->columns : fourier_info->width,1UL,exception);
1097     if (q == (PixelPacket *) NULL)
1098       break;
1099     indexes=GetCacheViewAuthenticIndexQueue(image_view);
1100     for (x=0L; x < (ssize_t) fourier_info->width; x++)
1101     {
1102       switch (fourier_info->channel)
1103       {
1104         case RedChannel:
1105         default:
1106         {
1107           q->red=ClampToQuantum(QuantumRange*source[i]);
1108           break;
1109         }
1110         case GreenChannel:
1111         {
1112           q->green=ClampToQuantum(QuantumRange*source[i]);
1113           break;
1114         }
1115         case BlueChannel:
1116         {
1117           q->blue=ClampToQuantum(QuantumRange*source[i]);
1118           break;
1119         }
1120         case OpacityChannel:
1121         {
1122           q->opacity=ClampToQuantum(QuantumRange*source[i]);
1123           break;
1124         }
1125         case IndexChannel:
1126         {
1127           indexes[x]=ClampToQuantum(QuantumRange*source[i]);
1128           break;
1129         }
1130         case GrayChannels:
1131         {
1132           q->red=ClampToQuantum(QuantumRange*source[i]);
1133           q->green=q->red;
1134           q->blue=q->red;
1135           break;
1136         }
1137       }
1138       i++;
1139       q++;
1140     }
1141     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1142       break;
1143   }
1144   image_view=DestroyCacheView(image_view);
1145   source=(double *) RelinquishMagickMemory(source);
1146   return(MagickTrue);
1147 }
1148
1149 static MagickBooleanType InverseFourierTransformChannel(
1150   const Image *magnitude_image,const Image *phase_image,
1151   const ChannelType channel,const MagickBooleanType modulus,
1152   Image *fourier_image,ExceptionInfo *exception)
1153 {
1154   double
1155     *magnitude,
1156     *phase;
1157
1158   fftw_complex
1159     *fourier;
1160
1161   FourierInfo
1162     fourier_info;
1163
1164   MagickBooleanType
1165     status;
1166
1167   size_t
1168     extent;
1169
1170   fourier_info.width=magnitude_image->columns;
1171   if ((magnitude_image->columns != magnitude_image->rows) ||
1172       ((magnitude_image->columns % 2) != 0) ||
1173       ((magnitude_image->rows % 2) != 0))
1174     {
1175       extent=magnitude_image->columns < magnitude_image->rows ?
1176         magnitude_image->rows : magnitude_image->columns;
1177       fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
1178     }
1179   fourier_info.height=fourier_info.width;
1180   fourier_info.center=(ssize_t) floor((double) fourier_info.width/2.0)+1L;
1181   fourier_info.channel=channel;
1182   fourier_info.modulus=modulus;
1183   magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1184     fourier_info.center*sizeof(double));
1185   if (magnitude == (double *) NULL)
1186     {
1187       (void) ThrowMagickException(exception,GetMagickModule(),
1188         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1189         magnitude_image->filename);
1190       return(MagickFalse);
1191     }
1192   phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
1193     fourier_info.center*sizeof(double));
1194   if (phase == (double *) NULL)
1195     {
1196       (void) ThrowMagickException(exception,GetMagickModule(),
1197         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1198         magnitude_image->filename);
1199       magnitude=(double *) RelinquishMagickMemory(magnitude);
1200       return(MagickFalse);
1201     }
1202   fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
1203     fourier_info.center*sizeof(*fourier));
1204   if (fourier == (fftw_complex *) NULL)
1205     {
1206       (void) ThrowMagickException(exception,GetMagickModule(),
1207         ResourceLimitError,"MemoryAllocationFailed","`%s'",
1208         magnitude_image->filename);
1209       phase=(double *) RelinquishMagickMemory(phase);
1210       magnitude=(double *) RelinquishMagickMemory(magnitude);
1211       return(MagickFalse);
1212     }
1213   status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
1214    exception);
1215   if (status != MagickFalse)
1216     status=InverseFourierTransform(&fourier_info,fourier,fourier_image,
1217       exception);
1218   fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
1219   phase=(double *) RelinquishMagickMemory(phase);
1220   magnitude=(double *) RelinquishMagickMemory(magnitude);
1221   return(status);
1222 }
1223 #endif
1224
1225 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
1226   const Image *phase_image,const MagickBooleanType modulus,
1227   ExceptionInfo *exception)
1228 {
1229   Image
1230     *fourier_image;
1231
1232   assert(magnitude_image != (Image *) NULL);
1233   assert(magnitude_image->signature == MagickSignature);
1234   if (magnitude_image->debug != MagickFalse)
1235     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1236       magnitude_image->filename);
1237   if (phase_image == (Image *) NULL)
1238     {
1239       (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
1240         "ImageSequenceRequired","`%s'",magnitude_image->filename);
1241       return((Image *) NULL);
1242     }
1243 #if !defined(MAGICKCORE_FFTW_DELEGATE)
1244   fourier_image=(Image *) NULL;
1245   (void) modulus;
1246   (void) ThrowMagickException(exception,GetMagickModule(),
1247     MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
1248     magnitude_image->filename);
1249 #else
1250   {
1251     fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
1252       magnitude_image->rows,MagickFalse,exception);
1253     if (fourier_image != (Image *) NULL)
1254       {
1255         MagickBooleanType
1256           is_gray,
1257           status;
1258
1259         status=MagickTrue;
1260         is_gray=IsGrayImage(magnitude_image,exception);
1261         if (is_gray != MagickFalse)
1262           is_gray=IsGrayImage(phase_image,exception);
1263 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1264         #pragma omp parallel sections
1265 #endif
1266         {
1267 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1268           #pragma omp section
1269 #endif
1270           {
1271             MagickBooleanType
1272               thread_status;
1273
1274             if (is_gray != MagickFalse)
1275               thread_status=InverseFourierTransformChannel(magnitude_image,
1276                 phase_image,GrayChannels,modulus,fourier_image,exception);
1277             else
1278               thread_status=InverseFourierTransformChannel(magnitude_image,
1279                 phase_image,RedChannel,modulus,fourier_image,exception);
1280             if (thread_status == MagickFalse)
1281               status=thread_status;
1282           }
1283 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1284           #pragma omp section
1285 #endif
1286           {
1287             MagickBooleanType
1288               thread_status;
1289
1290             thread_status=MagickTrue;
1291             if (is_gray == MagickFalse)
1292               thread_status=InverseFourierTransformChannel(magnitude_image,
1293                 phase_image,GreenChannel,modulus,fourier_image,exception);
1294             if (thread_status == MagickFalse)
1295               status=thread_status;
1296           }
1297 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1298           #pragma omp section
1299 #endif
1300           {
1301             MagickBooleanType
1302               thread_status;
1303
1304             thread_status=MagickTrue;
1305             if (is_gray == MagickFalse)
1306               thread_status=InverseFourierTransformChannel(magnitude_image,
1307                 phase_image,BlueChannel,modulus,fourier_image,exception);
1308             if (thread_status == MagickFalse)
1309               status=thread_status;
1310           }
1311 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1312           #pragma omp section
1313 #endif
1314           {
1315             MagickBooleanType
1316               thread_status;
1317
1318             thread_status=MagickTrue;
1319             if (magnitude_image->matte != MagickFalse)
1320               thread_status=InverseFourierTransformChannel(magnitude_image,
1321                 phase_image,OpacityChannel,modulus,fourier_image,exception);
1322             if (thread_status == MagickFalse)
1323               status=thread_status;
1324           }
1325 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1326           #pragma omp section
1327 #endif
1328           {
1329             MagickBooleanType
1330               thread_status;
1331
1332             thread_status=MagickTrue;
1333             if (magnitude_image->colorspace == CMYKColorspace)
1334               thread_status=InverseFourierTransformChannel(magnitude_image,
1335                 phase_image,IndexChannel,modulus,fourier_image,exception);
1336             if (thread_status == MagickFalse)
1337               status=thread_status;
1338           }
1339         }
1340         if (status == MagickFalse)
1341           fourier_image=DestroyImage(fourier_image);
1342       }
1343     fftw_cleanup();
1344   }
1345 #endif
1346   return(fourier_image);
1347 }