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