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