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

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