% Software Design %
% Sean Burke %
% Fred Weinhaus %
-% John Cristy %
+% Cristy %
% July 2009 %
% %
% %
-% Copyright 1999-2013 ImageMagick Studio LLC, a non-profit organization %
+% Copyright 1999-2014 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 %
Include declarations.
*/
#include "MagickCore/studio.h"
+#include "MagickCore/artifact.h"
#include "MagickCore/attribute.h"
#include "MagickCore/blob.h"
#include "MagickCore/cache.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)
% %
% %
% %
+% 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 op.
+%
+% 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 == MagickSignature);
+ if (images->debug != MagickFalse)
+ (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
+ assert(exception != (ExceptionInfo *) NULL);
+ assert(exception->signature == MagickSignature);
+ if (images->next == (Image *) NULL)
+ {
+ (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
+ "ImageSequenceRequired","`%s'",images->filename);
+ return((Image *) NULL);
+ }
+ image=CloneImage(images,images->columns,images->rows,MagickTrue,exception);
+ if (image == (Image *) NULL)
+ return((Image *) NULL);
+ if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
+ {
+ image=DestroyImageList(image);
+ return(image);
+ }
+ complex_images=NewImageList();
+ AppendImageToList(&complex_images,image);
+ image=CloneImage(images,images->columns,images->rows,MagickTrue,exception);
+ if (image == (Image *) NULL)
+ {
+ complex_images=DestroyImageList(complex_images);
+ return(complex_images);
+ }
+ if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
+ {
+ 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,4) shared(progress,status) \
+ magick_threads(images,complex_images,images->rows,1L)
+#endif
+ for (y=0; y < (ssize_t) images->rows; y++)
+ {
+ register const Quantum
+ *restrict Ai,
+ *restrict Ar,
+ *restrict Bi,
+ *restrict Br;
+
+ register Quantum
+ *restrict Ci,
+ *restrict Cr;
+
+ register ssize_t
+ x;
+
+ if (status == MagickFalse)
+ continue;
+ Ar=GetCacheViewVirtualPixels(Ar_view,0,y,images->columns,1,exception);
+ Ai=GetCacheViewVirtualPixels(Ai_view,0,y,images->columns,1,exception);
+ Br=GetCacheViewVirtualPixels(Br_view,0,y,images->columns,1,exception);
+ Bi=GetCacheViewVirtualPixels(Bi_view,0,y,images->columns,1,exception);
+ Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,images->columns,1,exception);
+ Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,images->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 critical (MagickCore_ComplexImages)
+#endif
+ 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
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% %
+% %
+% %
% 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 %
% %
% %
#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 *fourier_pixels)
+ const ssize_t x_offset,const ssize_t y_offset,double *roll_pixels)
{
double
- *roll_pixels;
+ *source_pixels;
MemoryInfo
- *roll_info;
+ *source_info;
register ssize_t
i,
/*
Move zero frequency (DC, average color) from (0,0) to (width/2,height/2).
*/
- roll_info=AcquireVirtualMemory(height,width*sizeof(*roll_pixels));
- if (roll_info == (MemoryInfo *) NULL)
+ source_info=AcquireVirtualMemory(height,width*sizeof(*source_pixels));
+ if (source_info == (MemoryInfo *) NULL)
return(MagickFalse);
- roll_pixels=(double *) GetVirtualMemoryBlob(roll_info);
+ source_pixels=(double *) GetVirtualMemoryBlob(source_info);
i=0L;
for (y=0L; y < (ssize_t) height; y++)
{
else
u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width :
x+x_offset;
- roll_pixels[v*width+u]=fourier_pixels[i++];
+ source_pixels[v*width+u]=roll_pixels[i++];
}
}
- (void) CopyMagickMemory(fourier_pixels,roll_pixels,height*width*
- sizeof(*roll_pixels));
- roll_info=RelinquishVirtualMemory(roll_info);
+ (void) CopyMagickMemory(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,double *destination)
+ const size_t height,double *source_pixels,double *forward_pixels)
{
MagickBooleanType
status;
Swap quadrants.
*/
center=(ssize_t) floor((double) width/2L)+1L;
- status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L,source);
+ 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-1L); x++)
- destination[width*y+x+width/2L]=source[center*y+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-1L); x++)
- destination[width*(height-y)+width/2L-x-1L]=source[center*y+x+1L];
+ forward_pixels[(height-y)*width+width/2L-x-1L]=
+ source_pixels[y*center+x+1L];
for (x=0L; x < (ssize_t) (width/2L); x++)
- destination[-x+width/2L-1L]=destination[x+width/2L+1L];
+ forward_pixels[-x+width/2L-1L]=forward_pixels[x+width/2L+1L];
return(MagickTrue);
}
static void CorrectPhaseLHS(const size_t width,const size_t height,
- double *fourier)
+ double *fourier_pixels)
{
register ssize_t
x;
for (y=0L; y < (ssize_t) height; y++)
for (x=0L; x < (ssize_t) (width/2L); x++)
- fourier[y*width+x]*=(-1.0);
+ fourier_pixels[y*width+x]*=(-1.0);
}
static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info,
phase_pixels=(double *) GetVirtualMemoryBlob(phase_info);
(void) ResetMagickMemory(phase_pixels,0,fourier_info->height*
fourier_info->width*sizeof(*phase_pixels));
- status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,
+ status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,
magnitude,magnitude_pixels);
if (status != MagickFalse)
- status=ForwardQuadrantSwap(fourier_info->height,fourier_info->height,phase,
+ status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase,
phase_pixels);
- CorrectPhaseLHS(fourier_info->height,fourier_info->height,phase_pixels);
+ CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
if (fourier_info->modulus != MagickFalse)
{
i=0L;
}
static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info,
- const Image *image,double *magnitude,double *phase,ExceptionInfo *exception)
+ const Image *image,double *magnitude_pixels,double *phase_pixels,
+ ExceptionInfo *exception)
{
CacheView
*image_view;
+ const char
+ *value;
+
double
- n,
*source_pixels;
fftw_complex
- *destination_pixels;
+ *forward_pixels;
fftw_plan
fftw_r2c_plan;
MemoryInfo
- *destination_info,
+ *forward_info,
*source_info;
register const Quantum
}
}
image_view=DestroyCacheView(image_view);
- destination_info=AcquireVirtualMemory((size_t) fourier_info->height,
- fourier_info->center*sizeof(*destination_pixels));
- if (destination_info == (MemoryInfo *) NULL)
+ forward_info=AcquireVirtualMemory((size_t) fourier_info->height,
+ fourier_info->center*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);
}
- destination_pixels=(fftw_complex *) GetVirtualMemoryBlob(destination_info);
+ 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,destination_pixels,FFTW_ESTIMATE);
+ source_pixels,forward_pixels,FFTW_ESTIMATE);
fftw_execute(fftw_r2c_plan);
fftw_destroy_plan(fftw_r2c_plan);
source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
- /*
- Normalize Fourier transform.
- */
- n=(double) fourier_info->width*(double) fourier_info->width;
- i=0L;
- for (y=0L; y < (ssize_t) fourier_info->height; y++)
- for (x=0L; x < (ssize_t) fourier_info->center; x++)
+ 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)
- destination_pixels[i]/=n;
+ forward_pixels[i]*=gamma;
#else
- destination_pixels[i][0]/=n;
- destination_pixels[i][1]/=n;
+ forward_pixels[i][0]*=gamma;
+ forward_pixels[i][1]*=gamma;
#endif
- i++;
+ i++;
+ }
}
/*
Generate magnitude and phase (or real and imaginary).
for (y=0L; y < (ssize_t) fourier_info->height; y++)
for (x=0L; x < (ssize_t) fourier_info->center; x++)
{
- magnitude[i]=cabs(destination_pixels[i]);
- phase[i]=carg(destination_pixels[i]);
+ 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[i]=creal(destination_pixels[i]);
- phase[i]=cimag(destination_pixels[i]);
+ magnitude_pixels[i]=creal(forward_pixels[i]);
+ phase_pixels[i]=cimag(forward_pixels[i]);
i++;
}
- destination_info=(MemoryInfo *) RelinquishVirtualMemory(destination_info);
+ forward_info=(MemoryInfo *) RelinquishVirtualMemory(forward_info);
return(MagickTrue);
}
extent;
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,
+ height,
width;
width=image->columns;
+ height=image->rows;
if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
((image->rows % 2) != 0))
{
extent=image->columns < image->rows ? image->rows : image->columns;
width=(extent & 0x01) == 1 ? extent+1UL : extent;
}
- magnitude_image=CloneImage(image,width,width,MagickFalse,exception);
+ height=width;
+ magnitude_image=CloneImage(image,width,height,MagickTrue,exception);
if (magnitude_image != (Image *) NULL)
{
Image
magnitude_image->storage_class=DirectClass;
magnitude_image->depth=32UL;
- phase_image=CloneImage(image,width,width,MagickFalse,exception);
+ phase_image=CloneImage(image,width,height,MagickTrue,exception);
if (phase_image == (Image *) NULL)
magnitude_image=DestroyImage(magnitude_image);
else
center=(ssize_t) floor((double) width/2L)+1L;
for (y=1L; y < (ssize_t) height; y++)
for (x=0L; x < (ssize_t) (width/2L+1L); x++)
- destination[center*(height-y)-x+width/2L]=source[y*width+x];
+ destination[(height-y)*center-x+width/2L]=source[y*width+x];
for (y=0L; y < (ssize_t) height; y++)
- destination[center*y]=source[y*width+width/2L];
+ 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,
- ExceptionInfo *exception)
+ const Image *magnitude_image,const Image *phase_image,
+ fftw_complex *fourier_pixels,ExceptionInfo *exception)
{
CacheView
*magnitude_view,
*phase_view;
double
- *magnitude,
- *phase,
+ *inverse_pixels,
*magnitude_pixels,
*phase_pixels;
MagickBooleanType
status;
+ MemoryInfo
+ *inverse_info,
+ *magnitude_info,
+ *phase_info;
+
register const Quantum
*p;
/*
Inverse fourier - read image and break down into a double array.
*/
- magnitude_pixels=(double *) AcquireQuantumMemory((size_t)
- fourier_info->height,fourier_info->width*sizeof(*magnitude_pixels));
- if (magnitude_pixels == (double *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- return(MagickFalse);
- }
- phase_pixels=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
+ magnitude_info=AcquireVirtualMemory((size_t)fourier_info->height,
+ fourier_info->width*sizeof(*magnitude_pixels));
+ phase_info=AcquireVirtualMemory((size_t) fourier_info->height,
fourier_info->width*sizeof(*phase_pixels));
- if (phase_pixels == (double *) NULL)
+ inverse_info=AcquireVirtualMemory((size_t) fourier_info->height,
+ fourier_info->center*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);
- magnitude_pixels=(double *) RelinquishMagickMemory(magnitude_pixels);
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+=GetPixelChannels(magnitude_image);
}
}
+ magnitude_view=DestroyCacheView(magnitude_view);
+ status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
+ magnitude_pixels,inverse_pixels);
+ (void) CopyMagickMemory(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++)
i++;
}
}
- magnitude_view=DestroyCacheView(magnitude_view);
phase_view=DestroyCacheView(phase_view);
- magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
- fourier_info->center*sizeof(*magnitude));
- if (magnitude == (double *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- magnitude_pixels=(double *) RelinquishMagickMemory(magnitude_pixels);
- phase_pixels=(double *) RelinquishMagickMemory(phase_pixels);
- return(MagickFalse);
- }
- status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
- magnitude_pixels,magnitude);
- magnitude_pixels=(double *) RelinquishMagickMemory(magnitude_pixels);
- phase=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
- fourier_info->width*sizeof(*phase));
- if (phase == (double *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- phase_pixels=(double *) RelinquishMagickMemory(phase_pixels);
- return(MagickFalse);
- }
- CorrectPhaseLHS(fourier_info->width,fourier_info->width,phase_pixels);
+ CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels);
if (status != MagickFalse)
status=InverseQuadrantSwap(fourier_info->width,fourier_info->height,
- phase_pixels,phase);
- phase_pixels=(double *) RelinquishMagickMemory(phase_pixels);
+ phase_pixels,inverse_pixels);
+ (void) CopyMagickMemory(phase_pixels,inverse_pixels,fourier_info->height*
+ fourier_info->center*sizeof(*phase_pixels));
+ inverse_info=RelinquishVirtualMemory(inverse_info);
/*
Merge two sets.
*/
for (x=0L; x < (ssize_t) fourier_info->center; x++)
{
#if defined(MAGICKCORE_HAVE_COMPLEX_H)
- fourier[i]=magnitude[i]*cos(phase[i])+I*magnitude[i]*sin(phase[i]);
+ fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I*
+ magnitude_pixels[i]*sin(phase_pixels[i]);
#else
- fourier[i][0]=magnitude[i]*cos(phase[i]);
- fourier[i][1]=magnitude[i]*sin(phase[i]);
+ fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]);
+ fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]);
#endif
i++;
}
for (x=0L; x < (ssize_t) fourier_info->center; x++)
{
#if defined(MAGICKCORE_HAVE_COMPLEX_H)
- fourier[i]=magnitude[i]+I*phase[i];
+ fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i];
#else
- fourier[i][0]=magnitude[i];
- fourier[i][1]=phase[i];
+ fourier_pixels[i][0]=magnitude_pixels[i];
+ fourier_pixels[i][1]=phase_pixels[i];
#endif
i++;
}
- phase=(double *) RelinquishMagickMemory(phase);
- magnitude=(double *) RelinquishMagickMemory(magnitude);
+ magnitude_info=RelinquishVirtualMemory(magnitude_info);
+ phase_info=RelinquishVirtualMemory(phase_info);
return(status);
}
static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info,
- fftw_complex *fourier,Image *image,ExceptionInfo *exception)
+ fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception)
{
CacheView
*image_view;
+ const char
+ *value;
+
double
- *source;
+ *source_pixels;
fftw_plan
fftw_c2r_plan;
+ MemoryInfo
+ *source_info;
+
register Quantum
*q;
ssize_t
y;
- source=(double *) AcquireQuantumMemory((size_t) fourier_info->height,
- fourier_info->width*sizeof(*source));
- if (source == (double *) NULL)
+ source_info=AcquireVirtualMemory((size_t) fourier_info->height,
+ fourier_info->width*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,source,FFTW_ESTIMATE);
+ fourier_pixels,source_pixels,FFTW_ESTIMATE);
fftw_execute(fftw_c2r_plan);
fftw_destroy_plan(fftw_c2r_plan);
}
case RedPixelChannel:
default:
{
- SetPixelRed(image,ClampToQuantum(QuantumRange*source[i]),q);
+ SetPixelRed(image,ClampToQuantum(QuantumRange*source_pixels[i]),q);
break;
}
case GreenPixelChannel:
{
- SetPixelGreen(image,ClampToQuantum(QuantumRange*source[i]),q);
+ SetPixelGreen(image,ClampToQuantum(QuantumRange*source_pixels[i]),
+ q);
break;
}
case BluePixelChannel:
{
- SetPixelBlue(image,ClampToQuantum(QuantumRange*source[i]),q);
+ SetPixelBlue(image,ClampToQuantum(QuantumRange*source_pixels[i]),
+ q);
break;
}
case BlackPixelChannel:
{
- SetPixelBlack(image,ClampToQuantum(QuantumRange*source[i]),q);
+ SetPixelBlack(image,ClampToQuantum(QuantumRange*source_pixels[i]),
+ q);
break;
}
case AlphaPixelChannel:
{
- SetPixelAlpha(image,ClampToQuantum(QuantumRange*source[i]),q);
+ SetPixelAlpha(image,ClampToQuantum(QuantumRange*source_pixels[i]),
+ q);
break;
}
}
break;
}
image_view=DestroyCacheView(image_view);
- source=(double *) RelinquishMagickMemory(source);
+ source_info=RelinquishVirtualMemory(source_info);
return(MagickTrue);
}
const PixelChannel channel,const MagickBooleanType modulus,
Image *fourier_image,ExceptionInfo *exception)
{
- double
- *magnitude,
- *phase;
-
fftw_complex
- *fourier;
+ *inverse_pixels;
FourierInfo
fourier_info;
MagickBooleanType
status;
+ MemoryInfo
+ *inverse_info;
+
size_t
extent;
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))
fourier_info.center=(ssize_t) floor((double) fourier_info.width/2L)+1L;
fourier_info.channel=channel;
fourier_info.modulus=modulus;
- magnitude=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
- fourier_info.center*sizeof(*magnitude));
- if (magnitude == (double *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- return(MagickFalse);
- }
- phase=(double *) AcquireQuantumMemory((size_t) fourier_info.height,
- fourier_info.center*sizeof(*phase));
- if (phase == (double *) NULL)
- {
- (void) ThrowMagickException(exception,GetMagickModule(),
- ResourceLimitError,"MemoryAllocationFailed","`%s'",
- magnitude_image->filename);
- magnitude=(double *) RelinquishMagickMemory(magnitude);
- return(MagickFalse);
- }
- fourier=(fftw_complex *) AcquireQuantumMemory((size_t) fourier_info.height,
- fourier_info.center*sizeof(*fourier));
- if (fourier == (fftw_complex *) NULL)
+ inverse_info=AcquireVirtualMemory((size_t) fourier_info.height,
+ fourier_info.center*sizeof(*inverse_pixels));
+ if (inverse_info == (MemoryInfo *) NULL)
{
(void) ThrowMagickException(exception,GetMagickModule(),
ResourceLimitError,"MemoryAllocationFailed","`%s'",
magnitude_image->filename);
- phase=(double *) RelinquishMagickMemory(phase);
- magnitude=(double *) RelinquishMagickMemory(magnitude);
return(MagickFalse);
}
- status=InverseFourier(&fourier_info,magnitude_image,phase_image,fourier,
- exception);
+ 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,fourier,fourier_image,
+ status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image,
exception);
- fourier=(fftw_complex *) RelinquishMagickMemory(fourier);
- phase=(double *) RelinquishMagickMemory(phase);
- magnitude=(double *) RelinquishMagickMemory(magnitude);
+ inverse_info=RelinquishVirtualMemory(inverse_info);
return(status);
}
#endif
#else
{
fourier_image=CloneImage(magnitude_image,magnitude_image->columns,
- magnitude_image->rows,MagickFalse,exception);
+ magnitude_image->rows,MagickTrue,exception);
if (fourier_image != (Image *) NULL)
{
MagickBooleanType