/MagickCore/fourier.c
C | 1589 lines | 1319 code | 101 blank | 169 comment | 251 complexity | d681ab508d88c7a9ea2f173b17e95104 MD5 | raw file
Possible License(s): MPL-2.0-no-copyleft-exception
- /*
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % %
- % %
- % FFFFF OOO U U RRRR IIIII EEEEE RRRR %
- % F O O U U R R I E R R %
- % FFF O O U U RRRR I EEE RRRR %
- % F O O U U R R I E R R %
- % F OOO UUU R R IIIII EEEEE R R %
- % %
- % %
- % MagickCore Discrete Fourier Transform Methods %
- % %
- % Software Design %
- % Sean Burke %
- % Fred Weinhaus %
- % Cristy %
- % July 2009 %
- % %
- % %
- % Copyright 1999-2019 ImageMagick Studio LLC, a non-profit organization %
- % dedicated to making software imaging solutions freely available. %
- % %
- % You may not use this file except in compliance with the License. You may %
- % obtain a copy of the License at %
- % %
- % https://imagemagick.org/script/license.php %
- % %
- % Unless required by applicable law or agreed to in writing, software %
- % distributed under the License is distributed on an "AS IS" BASIS, %
- % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
- % See the License for the specific language governing permissions and %
- % limitations under the License. %
- % %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- %
- %
- */
- /*
- Include declarations.
- */
- #include "MagickCore/studio.h"
- #include "MagickCore/artifact.h"
- #include "MagickCore/attribute.h"
- #include "MagickCore/blob.h"
- #include "MagickCore/cache.h"
- #include "MagickCore/image.h"
- #include "MagickCore/image-private.h"
- #include "MagickCore/list.h"
- #include "MagickCore/fourier.h"
- #include "MagickCore/log.h"
- #include "MagickCore/memory_.h"
- #include "MagickCore/monitor.h"
- #include "MagickCore/monitor-private.h"
- #include "MagickCore/pixel-accessor.h"
- #include "MagickCore/pixel-private.h"
- #include "MagickCore/property.h"
- #include "MagickCore/quantum-private.h"
- #include "MagickCore/resource_.h"
- #include "MagickCore/string-private.h"
- #include "MagickCore/thread-private.h"
- #if defined(MAGICKCORE_FFTW_DELEGATE)
- #if defined(MAGICKCORE_HAVE_COMPLEX_H)
- #include <complex.h>
- #endif
- #include <fftw3.h>
- #if !defined(MAGICKCORE_HAVE_CABS)
- #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1]))
- #endif
- #if !defined(MAGICKCORE_HAVE_CARG)
- #define carg(z) (atan2(cimag(z),creal(z)))
- #endif
- #if !defined(MAGICKCORE_HAVE_CIMAG)
- #define cimag(z) (z[1])
- #endif
- #if !defined(MAGICKCORE_HAVE_CREAL)
- #define creal(z) (z[0])
- #endif
- #endif
- /*
- Typedef declarations.
- */
- typedef struct _FourierInfo
- {
- PixelChannel
- channel;
- MagickBooleanType
- modulus;
- size_t
- width,
- height;
- ssize_t
- center;
- } FourierInfo;
- /*
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % %
- % %
- % C o m p l e x I m a g e s %
- % %
- % %
- % %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % ComplexImages() performs complex mathematics on an image sequence.
- %
- % The format of the ComplexImages method is:
- %
- % MagickBooleanType ComplexImages(Image *images,const ComplexOperator op,
- % ExceptionInfo *exception)
- %
- % A description of each parameter follows:
- %
- % o image: the image.
- %
- % o op: A complex operator.
- %
- % o exception: return any errors or warnings in this structure.
- %
- */
- MagickExport Image *ComplexImages(const Image *images,const ComplexOperator op,
- ExceptionInfo *exception)
- {
- #define ComplexImageTag "Complex/Image"
- CacheView
- *Ai_view,
- *Ar_view,
- *Bi_view,
- *Br_view,
- *Ci_view,
- *Cr_view;
- const char
- *artifact;
- const Image
- *Ai_image,
- *Ar_image,
- *Bi_image,
- *Br_image;
- double
- snr;
- Image
- *Ci_image,
- *complex_images,
- *Cr_image,
- *image;
- MagickBooleanType
- status;
- MagickOffsetType
- progress;
- ssize_t
- y;
- assert(images != (Image *) NULL);
- assert(images->signature == MagickCoreSignature);
- if (images->debug != MagickFalse)
- (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
- assert(exception != (ExceptionInfo *) NULL);
- assert(exception->signature == MagickCoreSignature);
- if (images->next == (Image *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
- "ImageSequenceRequired","`%s'",images->filename);
- return((Image *) NULL);
- }
- image=CloneImage(images,0,0,MagickTrue,exception);
- if (image == (Image *) NULL)
- return((Image *) NULL);
- if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
- {
- image=DestroyImageList(image);
- return(image);
- }
- image->depth=32UL;
- complex_images=NewImageList();
- AppendImageToList(&complex_images,image);
- image=CloneImage(images,0,0,MagickTrue,exception);
- if (image == (Image *) NULL)
- {
- complex_images=DestroyImageList(complex_images);
- return(complex_images);
- }
- AppendImageToList(&complex_images,image);
- /*
- Apply complex mathematics to image pixels.
- */
- artifact=GetImageArtifact(image,"complex:snr");
- snr=0.0;
- if (artifact != (const char *) NULL)
- snr=StringToDouble(artifact,(char **) NULL);
- Ar_image=images;
- Ai_image=images->next;
- Br_image=images;
- Bi_image=images->next;
- if ((images->next->next != (Image *) NULL) &&
- (images->next->next->next != (Image *) NULL))
- {
- Br_image=images->next->next;
- Bi_image=images->next->next->next;
- }
- Cr_image=complex_images;
- Ci_image=complex_images->next;
- Ar_view=AcquireVirtualCacheView(Ar_image,exception);
- Ai_view=AcquireVirtualCacheView(Ai_image,exception);
- Br_view=AcquireVirtualCacheView(Br_image,exception);
- Bi_view=AcquireVirtualCacheView(Bi_image,exception);
- Cr_view=AcquireAuthenticCacheView(Cr_image,exception);
- Ci_view=AcquireAuthenticCacheView(Ci_image,exception);
- status=MagickTrue;
- progress=0;
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp parallel for schedule(static) shared(progress,status) \
- magick_number_threads(images,complex_images,images->rows,1L)
- #endif
- for (y=0; y < (ssize_t) images->rows; y++)
- {
- register const Quantum
- *magick_restrict Ai,
- *magick_restrict Ar,
- *magick_restrict Bi,
- *magick_restrict Br;
- register Quantum
- *magick_restrict Ci,
- *magick_restrict Cr;
- register ssize_t
- x;
- if (status == MagickFalse)
- continue;
- Ar=GetCacheViewVirtualPixels(Ar_view,0,y,Ar_image->columns,1,exception);
- Ai=GetCacheViewVirtualPixels(Ai_view,0,y,Ai_image->columns,1,exception);
- Br=GetCacheViewVirtualPixels(Br_view,0,y,Br_image->columns,1,exception);
- Bi=GetCacheViewVirtualPixels(Bi_view,0,y,Bi_image->columns,1,exception);
- Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,Cr_image->columns,1,exception);
- Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,Ci_image->columns,1,exception);
- if ((Ar == (const Quantum *) NULL) || (Ai == (const Quantum *) NULL) ||
- (Br == (const Quantum *) NULL) || (Bi == (const Quantum *) NULL) ||
- (Cr == (Quantum *) NULL) || (Ci == (Quantum *) NULL))
- {
- status=MagickFalse;
- continue;
- }
- for (x=0; x < (ssize_t) images->columns; x++)
- {
- register ssize_t
- i;
- for (i=0; i < (ssize_t) GetPixelChannels(images); i++)
- {
- switch (op)
- {
- case AddComplexOperator:
- {
- Cr[i]=Ar[i]+Br[i];
- Ci[i]=Ai[i]+Bi[i];
- break;
- }
- case ConjugateComplexOperator:
- default:
- {
- Cr[i]=Ar[i];
- Ci[i]=(-Bi[i]);
- break;
- }
- case DivideComplexOperator:
- {
- double
- gamma;
- gamma=PerceptibleReciprocal(Br[i]*Br[i]+Bi[i]*Bi[i]+snr);
- Cr[i]=gamma*(Ar[i]*Br[i]+Ai[i]*Bi[i]);
- Ci[i]=gamma*(Ai[i]*Br[i]-Ar[i]*Bi[i]);
- break;
- }
- case MagnitudePhaseComplexOperator:
- {
- Cr[i]=sqrt(Ar[i]*Ar[i]+Ai[i]*Ai[i]);
- Ci[i]=atan2(Ai[i],Ar[i])/(2.0*MagickPI)+0.5;
- break;
- }
- case MultiplyComplexOperator:
- {
- Cr[i]=QuantumScale*(Ar[i]*Br[i]-Ai[i]*Bi[i]);
- Ci[i]=QuantumScale*(Ai[i]*Br[i]+Ar[i]*Bi[i]);
- break;
- }
- case RealImaginaryComplexOperator:
- {
- Cr[i]=Ar[i]*cos(2.0*MagickPI*(Ai[i]-0.5));
- Ci[i]=Ar[i]*sin(2.0*MagickPI*(Ai[i]-0.5));
- break;
- }
- case SubtractComplexOperator:
- {
- Cr[i]=Ar[i]-Br[i];
- Ci[i]=Ai[i]-Bi[i];
- break;
- }
- }
- }
- Ar+=GetPixelChannels(Ar_image);
- Ai+=GetPixelChannels(Ai_image);
- Br+=GetPixelChannels(Br_image);
- Bi+=GetPixelChannels(Bi_image);
- Cr+=GetPixelChannels(Cr_image);
- Ci+=GetPixelChannels(Ci_image);
- }
- if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse)
- status=MagickFalse;
- if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse)
- status=MagickFalse;
- if (images->progress_monitor != (MagickProgressMonitor) NULL)
- {
- MagickBooleanType
- proceed;
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp atomic
- #endif
- progress++;
- proceed=SetImageProgress(images,ComplexImageTag,progress,images->rows);
- if (proceed == MagickFalse)
- status=MagickFalse;
- }
- }
- Cr_view=DestroyCacheView(Cr_view);
- Ci_view=DestroyCacheView(Ci_view);
- Br_view=DestroyCacheView(Br_view);
- Bi_view=DestroyCacheView(Bi_view);
- Ar_view=DestroyCacheView(Ar_view);
- Ai_view=DestroyCacheView(Ai_view);
- if (status == MagickFalse)
- complex_images=DestroyImageList(complex_images);
- return(complex_images);
- }
- /*
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % %
- % %
- % 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 %
- % %
- % %
- % %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % ForwardFourierTransformImage() implements the discrete Fourier transform
- % (DFT) of the image either as a magnitude / phase or real / imaginary image
- % pair.
- %
- % The format of the ForwadFourierTransformImage method is:
- %
- % Image *ForwardFourierTransformImage(const Image *image,
- % const MagickBooleanType modulus,ExceptionInfo *exception)
- %
- % A description of each parameter follows:
- %
- % o image: the image.
- %
- % o modulus: if true, return as transform as a magnitude / phase pair
- % otherwise a real / imaginary image pair.
- %
- % o exception: return any errors or warnings in this structure.
- %
- */
- #if defined(MAGICKCORE_FFTW_DELEGATE)
- static MagickBooleanType RollFourier(const size_t width,const size_t height,
- const ssize_t x_offset,const ssize_t y_offset,double *roll_pixels)
- {
- double
- *source_pixels;
- MemoryInfo
- *source_info;
- register ssize_t
- i,
- x;
- ssize_t
- u,
- v,
- y;
- /*
- Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
- */
- source_info=AcquireVirtualMemory(width,height*sizeof(*source_pixels));
- if (source_info == (MemoryInfo *) NULL)
- return(MagickFalse);
- source_pixels=(double *) GetVirtualMemoryBlob(source_info);
- i=0L;
- for (y=0L; y < (ssize_t) height; y++)
- {
- if (y_offset < 0L)
- v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset;
- else
- v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height :
- y+y_offset;
- for (x=0L; x < (ssize_t) width; x++)
- {
- if (x_offset < 0L)
- u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset;
- else
- u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
- x+x_offset;
- source_pixels[v*width+u]=roll_pixels[i++];
- }
- }
- (void) memcpy(roll_pixels,source_pixels,height*width*
- sizeof(*source_pixels));
- source_info=RelinquishVirtualMemory(source_info);
- return(MagickTrue);
- }
- static MagickBooleanType ForwardQuadrantSwap(const size_t width,
- const size_t height,double *source_pixels,double *forward_pixels)
- {
- MagickBooleanType
- status;
- register ssize_t
- x;
- ssize_t
- center,
- y;
- /*
- Swap quadrants.
- */
- center=(ssize_t) (width/2L)+1L;
- status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,
- source_pixels);
- if (status == MagickFalse)
- return(MagickFalse);
- for (y=0L; y < (ssize_t) height; y++)
- for (x=0L; x < (ssize_t) (width/2L); x++)
- forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x];
- for (y=1; y < (ssize_t) height; y++)
- for (x=0L; x < (ssize_t) (width/2L); x++)
- forward_pixels[(height-y)*width+width/2L-x-1L]=
- source_pixels[y*center+x+1L];
- for (x=0L; x < (ssize_t) (width/2L); x++)
- forward_pixels[width/2L-x-1L]=source_pixels[x+1L];
- return(MagickTrue);
- }
- static void CorrectPhaseLHS(const size_t width,const size_t height,
- double *fourier_pixels)
- {
- register ssize_t
- x;
- ssize_t
- y;
- for (y=0L; y < (ssize_t) height; y++)
- for (x=0L; x < (ssize_t) (width/2L); x++)
- fourier_pixels[y*width+x]*=(-1.0);
- }
- static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
- Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
- {
- CacheView
- *magnitude_view,
- *phase_view;
- double
- *magnitude_pixels,
- *phase_pixels;
- Image
- *magnitude_image,
- *phase_image;
- MagickBooleanType
- status;
- MemoryInfo
- *magnitude_info,
- *phase_info;
- register Quantum
- *q;
- register ssize_t
- x;
- ssize_t
- i,
- y;
- magnitude_image=GetFirstImageInList(image);
- phase_image=GetNextImageInList(image);
- if (phase_image == (Image *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
- "ImageSequenceRequired","`%s'",image->filename);
- return(MagickFalse);
- }
- /*
- Create "Fourier Transform" image from constituent arrays.
- */
- magnitude_info=AcquireVirtualMemory((size_t) fourier_info->width,
- fourier_info->height*sizeof(*magnitude_pixels));
- phase_info=AcquireVirtualMemory((size_t) fourier_info->width,
- fourier_info->height*sizeof(*phase_pixels));
- if ((magnitude_info == (MemoryInfo *) NULL) ||
- (phase_info == (MemoryInfo *) NULL))
- {
- if (phase_info != (MemoryInfo *) NULL)
- phase_info=RelinquishVirtualMemory(phase_info);
- if (magnitude_info != (MemoryInfo *) NULL)
- magnitude_info=RelinquishVirtualMemory(magnitude_info);
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
- return(MagickFalse);
- }
- magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
- (void) memset(magnitude_pixels,0,fourier_info->width*
- fourier_info->height*sizeof(*magnitude_pixels));
- phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
- (void) memset(phase_pixels,0,fourier_info->width*
- fourier_info->height*sizeof(*phase_pixels));
- status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
- magnitude,magnitude_pixels);
- if (status != MagickFalse)
- status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
- phase_pixels);
- CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
- if (fourier_info->modulus != MagickFalse)
- {
- i=0L;
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- phase_pixels[i]/=(2.0*MagickPI);
- phase_pixels[i]+=0.5;
- i++;
- }
- }
- magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception);
- i=0L;
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- {
- q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL,
- exception);
- if (q == (Quantum *) NULL)
- break;
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- switch (fourier_info->channel)
- {
- case RedPixelChannel:
- default:
- {
- SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange*
- magnitude_pixels[i]),q);
- break;
- }
- case GreenPixelChannel:
- {
- SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange*
- magnitude_pixels[i]),q);
- break;
- }
- case BluePixelChannel:
- {
- SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange*
- magnitude_pixels[i]),q);
- break;
- }
- case BlackPixelChannel:
- {
- SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange*
- magnitude_pixels[i]),q);
- break;
- }
- case AlphaPixelChannel:
- {
- SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange*
- magnitude_pixels[i]),q);
- break;
- }
- }
- i++;
- q+=GetPixelChannels(magnitude_image);
- }
- status=SyncCacheViewAuthenticPixels(magnitude_view,exception);
- if (status == MagickFalse)
- break;
- }
- magnitude_view=DestroyCacheView(magnitude_view);
- i=0L;
- phase_view=AcquireAuthenticCacheView(phase_image,exception);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- {
- q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL,
- exception);
- if (q == (Quantum *) NULL)
- break;
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- switch (fourier_info->channel)
- {
- case RedPixelChannel:
- default:
- {
- SetPixelRed(phase_image,ClampToQuantum(QuantumRange*
- phase_pixels[i]),q);
- break;
- }
- case GreenPixelChannel:
- {
- SetPixelGreen(phase_image,ClampToQuantum(QuantumRange*
- phase_pixels[i]),q);
- break;
- }
- case BluePixelChannel:
- {
- SetPixelBlue(phase_image,ClampToQuantum(QuantumRange*
- phase_pixels[i]),q);
- break;
- }
- case BlackPixelChannel:
- {
- SetPixelBlack(phase_image,ClampToQuantum(QuantumRange*
- phase_pixels[i]),q);
- break;
- }
- case AlphaPixelChannel:
- {
- SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange*
- phase_pixels[i]),q);
- break;
- }
- }
- i++;
- q+=GetPixelChannels(phase_image);
- }
- status=SyncCacheViewAuthenticPixels(phase_view,exception);
- if (status == MagickFalse)
- break;
- }
- phase_view=DestroyCacheView(phase_view);
- phase_info=RelinquishVirtualMemory(phase_info);
- magnitude_info=RelinquishVirtualMemory(magnitude_info);
- return(status);
- }
- static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
- const Image *image,double *magnitude_pixels,double *phase_pixels,
- ExceptionInfo *exception)
- {
- CacheView
- *image_view;
- const char
- *value;
- double
- *source_pixels;
- fftw_complex
- *forward_pixels;
- fftw_plan
- fftw_r2c_plan;
- MemoryInfo
- *forward_info,
- *source_info;
- register const Quantum
- *p;
- register ssize_t
- i,
- x;
- ssize_t
- y;
- /*
- Generate the forward Fourier transform.
- */
- source_info=AcquireVirtualMemory((size_t) fourier_info->width,
- fourier_info->height*sizeof(*source_pixels));
- if (source_info == (MemoryInfo *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
- return(MagickFalse);
- }
- source_pixels=(double *) GetVirtualMemoryBlob(source_info);
- memset(source_pixels,0,fourier_info->width*fourier_info->height*
- sizeof(*source_pixels));
- i=0L;
- image_view=AcquireVirtualCacheView(image,exception);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- {
- p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL,
- exception);
- if (p == (const Quantum *) NULL)
- break;
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- switch (fourier_info->channel)
- {
- case RedPixelChannel:
- default:
- {
- source_pixels[i]=QuantumScale*GetPixelRed(image,p);
- break;
- }
- case GreenPixelChannel:
- {
- source_pixels[i]=QuantumScale*GetPixelGreen(image,p);
- break;
- }
- case BluePixelChannel:
- {
- source_pixels[i]=QuantumScale*GetPixelBlue(image,p);
- break;
- }
- case BlackPixelChannel:
- {
- source_pixels[i]=QuantumScale*GetPixelBlack(image,p);
- break;
- }
- case AlphaPixelChannel:
- {
- source_pixels[i]=QuantumScale*GetPixelAlpha(image,p);
- break;
- }
- }
- i++;
- p+=GetPixelChannels(image);
- }
- }
- image_view=DestroyCacheView(image_view);
- forward_info=AcquireVirtualMemory((size_t) fourier_info->width,
- (fourier_info->height/2+1)*sizeof(*forward_pixels));
- if (forward_info == (MemoryInfo *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
- source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
- return(MagickFalse);
- }
- forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info);
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp critical (MagickCore_ForwardFourierTransform)
- #endif
- fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->height,
- source_pixels,forward_pixels,FFTW_ESTIMATE);
- fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels);
- fftw_destroy_plan(fftw_r2c_plan);
- source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
- value=GetImageArtifact(image,"fourier:normalize");
- if ((value == (const char *) NULL) || (LocaleCompare(value,"forward") == 0))
- {
- double
- gamma;
- /*
- Normalize fourier transform.
- */
- i=0L;
- gamma=PerceptibleReciprocal((double) fourier_info->width*
- fourier_info->height);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
- {
- #if defined(MAGICKCORE_HAVE_COMPLEX_H)
- forward_pixels[i]*=gamma;
- #else
- forward_pixels[i][0]*=gamma;
- forward_pixels[i][1]*=gamma;
- #endif
- i++;
- }
- }
- /*
- Generate magnitude and phase (or real and imaginary).
- */
- i=0L;
- if (fourier_info->modulus != MagickFalse)
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
- {
- magnitude_pixels[i]=cabs(forward_pixels[i]);
- phase_pixels[i]=carg(forward_pixels[i]);
- i++;
- }
- else
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
- {
- magnitude_pixels[i]=creal(forward_pixels[i]);
- phase_pixels[i]=cimag(forward_pixels[i]);
- i++;
- }
- forward_info=(MemoryInfo *) RelinquishVirtualMemory(forward_info);
- return(MagickTrue);
- }
- static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
- const PixelChannel channel,const MagickBooleanType modulus,
- Image *fourier_image,ExceptionInfo *exception)
- {
- double
- *magnitude_pixels,
- *phase_pixels;
- FourierInfo
- fourier_info;
- MagickBooleanType
- status;
- MemoryInfo
- *magnitude_info,
- *phase_info;
- fourier_info.width=image->columns;
- fourier_info.height=image->rows;
- if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
- ((image->rows % 2) != 0))
- {
- size_t extent=image->columns < image->rows ? image->rows : image->columns;
- fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
- }
- fourier_info.height=fourier_info.width;
- fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
- fourier_info.channel=channel;
- fourier_info.modulus=modulus;
- magnitude_info=AcquireVirtualMemory((size_t) fourier_info.width,
- (fourier_info.height/2+1)*sizeof(*magnitude_pixels));
- phase_info=AcquireVirtualMemory((size_t) fourier_info.width,
- (fourier_info.height/2+1)*sizeof(*phase_pixels));
- if ((magnitude_info == (MemoryInfo *) NULL) ||
- (phase_info == (MemoryInfo *) NULL))
- {
- if (phase_info != (MemoryInfo *) NULL)
- phase_info=RelinquishVirtualMemory(phase_info);
- if (magnitude_info == (MemoryInfo *) NULL)
- magnitude_info=RelinquishVirtualMemory(magnitude_info);
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
- return(MagickFalse);
- }
- magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
- phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
- status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels,
- phase_pixels,exception);
- if (status != MagickFalse)
- status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels,
- phase_pixels,exception);
- phase_info=RelinquishVirtualMemory(phase_info);
- magnitude_info=RelinquishVirtualMemory(magnitude_info);
- return(status);
- }
- #endif
- MagickExport Image *ForwardFourierTransformImage(const Image *image,
- const MagickBooleanType modulus,ExceptionInfo *exception)
- {
- Image
- *fourier_image;
- fourier_image=NewImageList();
- #if !defined(MAGICKCORE_FFTW_DELEGATE)
- (void) modulus;
- (void) ThrowMagickException(exception,GetMagickModule(),
- MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
- image->filename);
- #else
- {
- Image
- *magnitude_image;
- size_t
- height,
- width;
- width=image->columns;
- height=image->rows;
- if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
- ((image->rows % 2) != 0))
- {
- size_t extent=image->columns < image->rows ? image->rows :
- image->columns;
- width=(extent & 0x01) == 1 ? extent+1UL : extent;
- }
- height=width;
- magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
- if (magnitude_image != (Image *) NULL)
- {
- Image
- *phase_image;
- magnitude_image->storage_class=DirectClass;
- magnitude_image->depth=32UL;
- phase_image=CloneImage(image,width,height,MagickTrue,exception);
- if (phase_image == (Image *) NULL)
- magnitude_image=DestroyImage(magnitude_image);
- else
- {
- MagickBooleanType
- is_gray,
- status;
- phase_image->storage_class=DirectClass;
- phase_image->depth=32UL;
- AppendImageToList(&fourier_image,magnitude_image);
- AppendImageToList(&fourier_image,phase_image);
- status=MagickTrue;
- is_gray=IsImageGray(image);
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp parallel sections
- #endif
- {
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- if (is_gray != MagickFalse)
- thread_status=ForwardFourierTransformChannel(image,
- GrayPixelChannel,modulus,fourier_image,exception);
- else
- thread_status=ForwardFourierTransformChannel(image,
- RedPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (is_gray == MagickFalse)
- thread_status=ForwardFourierTransformChannel(image,
- GreenPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (is_gray == MagickFalse)
- thread_status=ForwardFourierTransformChannel(image,
- BluePixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (image->colorspace == CMYKColorspace)
- thread_status=ForwardFourierTransformChannel(image,
- BlackPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (image->alpha_trait != UndefinedPixelTrait)
- thread_status=ForwardFourierTransformChannel(image,
- AlphaPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- }
- if (status == MagickFalse)
- fourier_image=DestroyImageList(fourier_image);
- fftw_cleanup();
- }
- }
- }
- #endif
- return(fourier_image);
- }
- /*
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % %
- % %
- % 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 %
- % %
- % %
- % %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % InverseFourierTransformImage() implements the inverse discrete Fourier
- % transform (DFT) of the image either as a magnitude / phase or real /
- % imaginary image pair.
- %
- % The format of the InverseFourierTransformImage method is:
- %
- % Image *InverseFourierTransformImage(const Image *magnitude_image,
- % const Image *phase_image,const MagickBooleanType modulus,
- % ExceptionInfo *exception)
- %
- % A description of each parameter follows:
- %
- % o magnitude_image: the magnitude or real image.
- %
- % o phase_image: the phase or imaginary image.
- %
- % o modulus: if true, return transform as a magnitude / phase pair
- % otherwise a real / imaginary image pair.
- %
- % o exception: return any errors or warnings in this structure.
- %
- */
- #if defined(MAGICKCORE_FFTW_DELEGATE)
- static MagickBooleanType InverseQuadrantSwap(const size_t width,
- const size_t height,const double *source,double *destination)
- {
- register ssize_t
- x;
- ssize_t
- center,
- y;
- /*
- Swap quadrants.
- */
- center=(ssize_t) (width/2L)+1L;
- for (y=1L; y < (ssize_t) height; y++)
- for (x=0L; x < (ssize_t) (width/2L+1L); x++)
- destination[(height-y)*center-x+width/2L]=source[y*width+x];
- for (y=0L; y < (ssize_t) height; y++)
- destination[y*center]=source[y*width+width/2L];
- for (x=0L; x < center; x++)
- destination[x]=source[center-x-1L];
- return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination));
- }
- static MagickBooleanType InverseFourier(FourierInfo *fourier_info,
- const Image *magnitude_image,const Image *phase_image,
- fftw_complex *fourier_pixels,ExceptionInfo *exception)
- {
- CacheView
- *magnitude_view,
- *phase_view;
- double
- *inverse_pixels,
- *magnitude_pixels,
- *phase_pixels;
- MagickBooleanType
- status;
- MemoryInfo
- *inverse_info,
- *magnitude_info,
- *phase_info;
- register const Quantum
- *p;
- register ssize_t
- i,
- x;
- ssize_t
- y;
- /*
- Inverse fourier - read image and break down into a double array.
- */
- magnitude_info=AcquireVirtualMemory((size_t) fourier_info->width,
- fourier_info->height*sizeof(*magnitude_pixels));
- phase_info=AcquireVirtualMemory((size_t) fourier_info->width,
- fourier_info->height*sizeof(*phase_pixels));
- inverse_info=AcquireVirtualMemory((size_t) fourier_info->width,
- (fourier_info->height/2+1)*sizeof(*inverse_pixels));
- if ((magnitude_info == (MemoryInfo *) NULL) ||
- (phase_info == (MemoryInfo *) NULL) ||
- (inverse_info == (MemoryInfo *) NULL))
- {
- if (magnitude_info != (MemoryInfo *) NULL)
- magnitude_info=RelinquishVirtualMemory(magnitude_info);
- if (phase_info != (MemoryInfo *) NULL)
- phase_info=RelinquishVirtualMemory(phase_info);
- if (inverse_info != (MemoryInfo *) NULL)
- inverse_info=RelinquishVirtualMemory(inverse_info);
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- return(MagickFalse);
- }
- magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info);
- phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
- inverse_pixels=(double *) GetVirtualMemoryBlob(inverse_info);
- i=0L;
- magnitude_view=AcquireVirtualCacheView(magnitude_image,exception);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- {
- p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL,
- exception);
- if (p == (const Quantum *) NULL)
- break;
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- switch (fourier_info->channel)
- {
- case RedPixelChannel:
- default:
- {
- magnitude_pixels[i]=QuantumScale*GetPixelRed(magnitude_image,p);
- break;
- }
- case GreenPixelChannel:
- {
- magnitude_pixels[i]=QuantumScale*GetPixelGreen(magnitude_image,p);
- break;
- }
- case BluePixelChannel:
- {
- magnitude_pixels[i]=QuantumScale*GetPixelBlue(magnitude_image,p);
- break;
- }
- case BlackPixelChannel:
- {
- magnitude_pixels[i]=QuantumScale*GetPixelBlack(magnitude_image,p);
- break;
- }
- case AlphaPixelChannel:
- {
- magnitude_pixels[i]=QuantumScale*GetPixelAlpha(magnitude_image,p);
- break;
- }
- }
- i++;
- p+=GetPixelChannels(magnitude_image);
- }
- }
- magnitude_view=DestroyCacheView(magnitude_view);
- status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
- magnitude_pixels,inverse_pixels);
- (void) memcpy(magnitude_pixels,inverse_pixels,fourier_info->height*
- fourier_info->center*sizeof(*magnitude_pixels));
- i=0L;
- phase_view=AcquireVirtualCacheView(phase_image,exception);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- {
- p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1,
- exception);
- if (p == (const Quantum *) NULL)
- break;
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- switch (fourier_info->channel)
- {
- case RedPixelChannel:
- default:
- {
- phase_pixels[i]=QuantumScale*GetPixelRed(phase_image,p);
- break;
- }
- case GreenPixelChannel:
- {
- phase_pixels[i]=QuantumScale*GetPixelGreen(phase_image,p);
- break;
- }
- case BluePixelChannel:
- {
- phase_pixels[i]=QuantumScale*GetPixelBlue(phase_image,p);
- break;
- }
- case BlackPixelChannel:
- {
- phase_pixels[i]=QuantumScale*GetPixelBlack(phase_image,p);
- break;
- }
- case AlphaPixelChannel:
- {
- phase_pixels[i]=QuantumScale*GetPixelAlpha(phase_image,p);
- break;
- }
- }
- i++;
- p+=GetPixelChannels(phase_image);
- }
- }
- if (fourier_info->modulus != MagickFalse)
- {
- i=0L;
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- phase_pixels[i]-=0.5;
- phase_pixels[i]*=(2.0*MagickPI);
- i++;
- }
- }
- phase_view=DestroyCacheView(phase_view);
- CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
- if (status != MagickFalse)
- status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
- phase_pixels,inverse_pixels);
- (void) memcpy(phase_pixels,inverse_pixels,fourier_info->height*
- fourier_info->center*sizeof(*phase_pixels));
- inverse_info=RelinquishVirtualMemory(inverse_info);
- /*
- Merge two sets.
- */
- i=0L;
- if (fourier_info->modulus != MagickFalse)
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
- {
- #if defined(MAGICKCORE_HAVE_COMPLEX_H)
- fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
- magnitude_pixels[i]*sin(phase_pixels[i]);
- #else
- fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
- fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
- #endif
- i++;
- }
- else
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
- {
- #if defined(MAGICKCORE_HAVE_COMPLEX_H)
- fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
- #else
- fourier_pixels[i][0]=magnitude_pixels[i];
- fourier_pixels[i][1]=phase_pixels[i];
- #endif
- i++;
- }
- magnitude_info=RelinquishVirtualMemory(magnitude_info);
- phase_info=RelinquishVirtualMemory(phase_info);
- return(status);
- }
- static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
- fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception)
- {
- CacheView
- *image_view;
- const char
- *value;
- double
- *source_pixels;
- fftw_plan
- fftw_c2r_plan;
- MemoryInfo
- *source_info;
- register Quantum
- *q;
- register ssize_t
- i,
- x;
- ssize_t
- y;
- source_info=AcquireVirtualMemory((size_t) fourier_info->width,
- fourier_info->height*sizeof(*source_pixels));
- if (source_info == (MemoryInfo *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
- return(MagickFalse);
- }
- source_pixels=(double *) GetVirtualMemoryBlob(source_info);
- value=GetImageArtifact(image,"fourier:normalize");
- if (LocaleCompare(value,"inverse") == 0)
- {
- double
- gamma;
- /*
- Normalize inverse transform.
- */
- i=0L;
- gamma=PerceptibleReciprocal((double) fourier_info->width*
- fourier_info->height);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
- {
- #if defined(MAGICKCORE_HAVE_COMPLEX_H)
- fourier_pixels[i]*=gamma;
- #else
- fourier_pixels[i][0]*=gamma;
- fourier_pixels[i][1]*=gamma;
- #endif
- i++;
- }
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp critical (MagickCore_InverseFourierTransform)
- #endif
- fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height,
- fourier_pixels,source_pixels,FFTW_ESTIMATE);
- fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels);
- fftw_destroy_plan(fftw_c2r_plan);
- i=0L;
- image_view=AcquireAuthenticCacheView(image,exception);
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- {
- if (y >= (ssize_t) image->rows)
- break;
- q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width >
- image->columns ? image->columns : fourier_info->width,1UL,exception);
- if (q == (Quantum *) NULL)
- break;
- for (x=0L; x < (ssize_t) fourier_info->width; x++)
- {
- if (x < (ssize_t) image->columns)
- switch (fourier_info->channel)
- {
- case RedPixelChannel:
- default:
- {
- SetPixelRed(image,ClampToQuantum(QuantumRange*source_pixels[i]),q);
- break;
- }
- case GreenPixelChannel:
- {
- SetPixelGreen(image,ClampToQuantum(QuantumRange*source_pixels[i]),
- q);
- break;
- }
- case BluePixelChannel:
- {
- SetPixelBlue(image,ClampToQuantum(QuantumRange*source_pixels[i]),
- q);
- break;
- }
- case BlackPixelChannel:
- {
- SetPixelBlack(image,ClampToQuantum(QuantumRange*source_pixels[i]),
- q);
- break;
- }
- case AlphaPixelChannel:
- {
- SetPixelAlpha(image,ClampToQuantum(QuantumRange*source_pixels[i]),
- q);
- break;
- }
- }
- i++;
- q+=GetPixelChannels(image);
- }
- if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
- break;
- }
- image_view=DestroyCacheView(image_view);
- source_info=RelinquishVirtualMemory(source_info);
- return(MagickTrue);
- }
- static MagickBooleanType InverseFourierTransformChannel(
- const Image *magnitude_image,const Image *phase_image,
- const PixelChannel channel,const MagickBooleanType modulus,
- Image *fourier_image,ExceptionInfo *exception)
- {
- fftw_complex
- *inverse_pixels;
- FourierInfo
- fourier_info;
- MagickBooleanType
- status;
- MemoryInfo
- *inverse_info;
- fourier_info.width=magnitude_image->columns;
- fourier_info.height=magnitude_image->rows;
- if ((magnitude_image->columns != magnitude_image->rows) ||
- ((magnitude_image->columns % 2) != 0) ||
- ((magnitude_image->rows % 2) != 0))
- {
- size_t extent=magnitude_image->columns < magnitude_image->rows ?
- magnitude_image->rows : magnitude_image->columns;
- fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent;
- }
- fourier_info.height=fourier_info.width;
- fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L;
- fourier_info.channel=channel;
- fourier_info.modulus=modulus;
- inverse_info=AcquireVirtualMemory((size_t) fourier_info.width,
- (fourier_info.height/2+1)*sizeof(*inverse_pixels));
- if (inverse_info == (MemoryInfo *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- return(MagickFalse);
- }
- inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info);
- status=InverseFourier(&fourier_info,magnitude_image,phase_image,
- inverse_pixels,exception);
- if (status != MagickFalse)
- status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
- exception);
- inverse_info=RelinquishVirtualMemory(inverse_info);
- return(status);
- }
- #endif
- MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
- const Image *phase_image,const MagickBooleanType modulus,
- ExceptionInfo *exception)
- {
- Image
- *fourier_image;
- assert(magnitude_image != (Image *) NULL);
- assert(magnitude_image->signature == MagickCoreSignature);
- if (magnitude_image->debug != MagickFalse)
- (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
- magnitude_image->filename);
- if (phase_image == (Image *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
- "ImageSequenceRequired","`%s'",magnitude_image->filename);
- return((Image *) NULL);
- }
- #if !defined(MAGICKCORE_FFTW_DELEGATE)
- fourier_image=(Image *) NULL;
- (void) modulus;
- (void) ThrowMagickException(exception,GetMagickModule(),
- MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)",
- magnitude_image->filename);
- #else
- {
- fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
- magnitude_image->rows,MagickTrue,exception);
- if (fourier_image != (Image *) NULL)
- {
- MagickBooleanType
- is_gray,
- status;
- status=MagickTrue;
- is_gray=IsImageGray(magnitude_image);
- if (is_gray != MagickFalse)
- is_gray=IsImageGray(phase_image);
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp parallel sections
- #endif
- {
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- if (is_gray != MagickFalse)
- thread_status=InverseFourierTransformChannel(magnitude_image,
- phase_image,GrayPixelChannel,modulus,fourier_image,exception);
- else
- thread_status=InverseFourierTransformChannel(magnitude_image,
- phase_image,RedPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (is_gray == MagickFalse)
- thread_status=InverseFourierTransformChannel(magnitude_image,
- phase_image,GreenPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (is_gray == MagickFalse)
- thread_status=InverseFourierTransformChannel(magnitude_image,
- phase_image,BluePixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (magnitude_image->colorspace == CMYKColorspace)
- thread_status=InverseFourierTransformChannel(magnitude_image,
- phase_image,BlackPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- #if defined(MAGICKCORE_OPENMP_SUPPORT)
- #pragma omp section
- #endif
- {
- MagickBooleanType
- thread_status;
- thread_status=MagickTrue;
- if (magnitude_image->alpha_trait != UndefinedPixelTrait)
- thread_status=InverseFourierTransformChannel(magnitude_image,
- phase_image,AlphaPixelChannel,modulus,fourier_image,exception);
- if (thread_status == MagickFalse)
- status=thread_status;
- }
- }
- if (status == MagickFalse)
- fourier_image=DestroyImage(fourier_image);
- }
- fftw_cleanup();
- }
- #endif
- return(fourier_image);
- }