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