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