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