% 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 %
% %
% %
CacheView
*image_view;
+ const char
+ *value;
+
double
- gamma,
*source_pixels;
fftw_complex
fftw_execute(fftw_r2c_plan);
fftw_destroy_plan(fftw_r2c_plan);
source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info);
- /*
- 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++)
+ 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;
+ forward_pixels[i]*=gamma;
#else
- forward_pixels[i][0]*=gamma;
- forward_pixels[i][1]*=gamma;
+ forward_pixels[i][0]*=gamma;
+ forward_pixels[i][1]*=gamma;
#endif
- i++;
+ i++;
+ }
}
/*
Generate magnitude and phase (or real and imaginary).
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
CacheView
*image_view;
+ const char
+ *value;
+
double
*source_pixels;
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
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))
#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
projects / imagemagick / blobdiff