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