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

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