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