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

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