Changeset 1983 for trunk/psLib/src/image/psImageConvolve.c

Timestamp:

Oct 6, 2004, 11:31:30 AM (22 years ago)

Author:

desonia

Message:

work on FFT and convolution. Fourier-domain convolution is still
broken.

File:

• trunk/psLib/src/image/psImageConvolve.c (modified) (7 diffs)

trunk/psLib/src/image/psImageConvolve.c

@@ -5 +5 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-01 20:58:32 $
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-06 21:31:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 #include "psLogMsg.h"
 #include "psError.h"
+#include "psImageIO.h"
 
 #include "psImageErrors.h"
@@ -315 +316 @@
         int kRows = kernel->image->numRows;
         int kCols = kernel->image->numCols;
-        int x0 = (paddedCols - kCols) / 2;
-        int y0 = (paddedRows - kRows) / 2;
-        if (x0 < 0 || y0 < 0) {
+        if (kRows >= numRows || kCols >= numCols) {
             psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
                        PS_ERR_BAD_PARAMETER_SIZE, true,
                        PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE,
-                       kernel->image->numCols,kernel->image->numRows,
+                       kCols,kRows,
                        numCols, numRows);
             psFree(out);
@@ -350 +349 @@
         }
 
+        psImageWriteSection(paddedImage,0,0,0,NULL,0,"paddedImage.fits");
+
         // pad the kernel to the same size of paddedImage
         psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL);
         memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image
-        for (int row=0;row<kRows;row++) {
-            psKernelType* padRow = paddedKernel->data.PS_TYPE_KERNEL_DATA[row+y0];
-            psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row];
-            for (int col=0;col<kCols;col++) {
-                padRow[col+x0] = kernelRow[col];
+        int yMax = kernel->yMax;
+        int xMax = kernel->xMax;
+        for (int row = kernel->yMin; row <= yMax;row++) {
+            int padRow = row;
+            if (padRow < 0) {
+                padRow += paddedRows;
             }
-        }
+            psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow];
+            psKernelType* kernelRow = kernel->kernel[row];
+            for (int col = kernel->xMin; col <= xMax; col++) {
+                if (col < 0) {
+                    padData[col+paddedCols] = kernelRow[col];
+                } else {
+                    padData[col] = kernelRow[col];
+                }
+            }
+        }
+
+        psImageWriteSection(paddedKernel,0,0,0,NULL,0,"paddedKernel.fits");
 
         psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD);
@@ -380 +393 @@
 
         // convolution in fourier domain is just a pixel-wise multiplication
-        psImage* outFourier = psImageAlloc(paddedCols,paddedRows,inFourier->type.type);
-        psBinaryOp(outFourier,inFourier,"*",kernelFourier);
-
-        psImage* complexOut = psImageFFT(NULL, outFourier, PS_FFT_REVERSE);
+        psBinaryOp(inFourier,inFourier,"*",kernelFourier);
+
+        psImage* complexOut = psImageFFT(NULL, inFourier,
+                                         PS_FFT_REVERSE);
         if (complexOut == NULL) {
             psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve",
@@ -392 +405 @@
         }
 
+        {
+            psImage* tmp = psImageReal(NULL,complexOut);
+            psImageWriteSection(tmp,0,0,0,NULL,0,"complexOut.fits");
+            psFree(tmp);
+        }
+
         // subset out the padded area now.
         psImage* complexOutSansPad = psImageSubset(complexOut,
@@ -397 +416 @@
                                      FOURIER_PADDING+numCols,FOURIER_PADDING+numRows);
 
-        // since our image was real, the result is to be real as well
-        out = psImageReal(out,complexOutSansPad);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
+        float factor = 1.0f/numCols/numRows;
+        for (int row = 0; row < numRows; row++) {
+            psF32* outRow = out->data.F32[row];
+            psC32* resultRow = complexOutSansPad->data.C32[row];
+            for (int col = 0; col < numCols; col++) {
+                outRow[col] = crealf(resultRow[col])*factor;
+            }
+        }
 
         psFree(complexOut); // frees complexOutSansPad, as it is a child.