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

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