]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/fourier.c
(no commit message)
[imagemagick] / MagickCore / fourier.c
index 6e18dcb5a848d7faa172758e7743a9a8f04fa127..a649eca2f53c54d77987d1ba80e19eca51fcdeea 100644 (file)
 %                              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  %
@@ -60,6 +60,7 @@
 #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)
@@ -115,19 +116,19 @@ typedef struct _FourierInfo
 %  The format of the ComplexImages method is:
 %
 %      MagickBooleanType ComplexImages(Image *images,
-%        const ComplexOperator operator,ExceptionInfo *exception)
+%        const ComplexOperator op,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o operator: A complex operator.
+%    o op: A complex op.
 %
 %    o exception: return any errors or warnings in this structure.
 %
 */
-MagickExport Image *ComplexImages(const Image *images,
-  const ComplexOperator operator,ExceptionInfo *exception)
+MagickExport Image *ComplexImages(const Image *images,const ComplexOperator op,
+  ExceptionInfo *exception)
 {
 #define ComplexImageTag  "Complex/Image"
 
@@ -139,16 +140,23 @@ MagickExport Image *ComplexImages(const Image *images,
     *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;
+    *Cr_image,
+    *image;
 
   MagickBooleanType
     status;
@@ -171,29 +179,41 @@ MagickExport Image *ComplexImages(const Image *images,
         "ImageSequenceRequired","`%s'",images->filename);
       return((Image *) NULL);
     }
-  complex_images=CloneImage(images,images->columns,images->rows,MagickTrue,
-    exception);
-  if (complex_images == (Image *) NULL)
+  image=CloneImage(images,images->columns,images->rows,MagickTrue,exception);
+  if (image == (Image *) NULL)
     return((Image *) NULL);
-  complex_images->next=CloneImage(images,images->columns,images->rows,
-    MagickTrue,exception);
-  if (complex_images->next == (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;
-  if ((images->next->next == (Image *) NULL) ||
-      (images->next->next->next == (Image *) NULL))
-    {
-      Br_image=images;
-      Bi_image=images->next;
-    }
-  else
+  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;
@@ -227,6 +247,8 @@ MagickExport Image *ComplexImages(const Image *images,
     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);
@@ -242,6 +264,68 @@ MagickExport Image *ComplexImages(const Image *images,
       }
     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;
@@ -770,6 +854,7 @@ static MagickBooleanType ForwardFourierTransformChannel(const Image *image,
     extent;
 
   fourier_info.width=image->columns;
+  fourier_info.height=image->rows;
   if ((image->columns != image->rows) || ((image->columns % 2) != 0) ||
       ((image->rows % 2) != 0))
     {
@@ -827,16 +912,19 @@ MagickExport Image *ForwardFourierTransformImage(const Image *image,
 
     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
@@ -844,7 +932,7 @@ MagickExport Image *ForwardFourierTransformImage(const Image *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
@@ -1355,6 +1443,7 @@ static MagickBooleanType InverseFourierTransformChannel(
     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))
@@ -1414,7 +1503,7 @@ MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image,
 #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