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