Changeset 14480 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Aug 13, 2007, 4:14:30 PM (19 years ago)

Author:

Paul Price

Message:

Adding enhancement to extract the pixels for each stamp from the image, so that pmSubtractionCalculateEquation doesn't need to look at the input image (so I can later put in a fake target PSF). The stamps are extracted as 'kernels' so that offsets need not be performed.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (14 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-09 21:20:35 $
+ *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-14 02:14:30 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -217 +217 @@
                                    int index, // Kernel basis function index
                                    int x, int y, // Pixel around which to convolve
-                                   const psImage *image // Image to convolve
+                                   const psKernel *image // Image to convolve (a kernel for convenience)
                                    )
 {
@@ -225 +225 @@
           int u = kernels->u->data.S32[index]; // Offset in x
           int v = kernels->v->data.S32[index]; // Offset in y
-          float value = image->data.F32[y + v][x + u]; // Value of convolution
+          float value = image->kernel[y + v][x + u]; // Value of convolution
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
               // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-              value -= image->data.F32[y][x];
+              value -= image->kernel[y][x];
           }
           return value;
@@ -242 +242 @@
           for (int v = vStart; v <= vStop; v++) {
               for (int u = uStart; u <= uStop; u++) {
-                  sum += image->data.F32[y + v][x + u];
+                  sum += image->kernel[y + v][x + u];
               }
           }
@@ -248 +248 @@
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
               // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-              sum -= image->data.F32[y][x];
+              sum -= image->kernel[y][x];
           }
           return sum;
@@ -260 +260 @@
               for (int v = -size; v <= size; v++) {
                   for (int u = -size; u <= size; u++) {
-                      sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
+                      sum += kernel->kernel[v][u] * image->kernel[y + v][x + u];
                   }
               }
@@ -268 +268 @@
           int u = kernels->u->data.S32[index]; // Offset in x
           int v = kernels->v->data.S32[index]; // Offset in y
-          float value = image->data.F32[y + v][x + u]; // Value of convolution
+          float value = image->kernel[y + v][x + u]; // Value of convolution
           // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              value -= image->data.F32[y][x];
+              value -= image->kernel[y][x];
           }
           return value;
@@ -281 +281 @@
           for (int v = -size; v <= size; v++) {
               for (int u = -size; u <= size; u++) {
-                  sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u];
+                  sum += kernel->kernel[v][u] * image->kernel[y + v][x + u];
                   // Photometric scaling is already built in to the precalculated kernel
               }
@@ -289 +289 @@
       case PM_SUBTRACTION_KERNEL_RINGS: {
           if (index == kernels->subIndex) {
-              return image->data.F32[y][x];
+              return image->kernel[y][x];
           }
           psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
@@ -299 +299 @@
           for (int j = 0; j < num; j++) {
               int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
-              sum += image->data.F32[y + v][x + u] * poly->data.F32[j];
+              sum += image->kernel[y + v][x + u] * poly->data.F32[j];
           }
           // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
-          return sum - num * image->data.F32[y][x];
+          return sum - num * image->kernel[y][x];
       }
       default:
@@ -436 +436 @@
 
 
-bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input,
-                                    const psImage *weight, const pmSubtractionKernels *kernels, int footprint)
+bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint)
 {
     PS_ASSERT_ARRAY_NON_NULL(stamps, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
-    PS_ASSERT_IMAGE_NON_NULL(input, false);
-    PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(reference, input, false);
     PS_ASSERT_PTR_NON_NULL(kernels, false);
     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
-    if (weight) {
-        PS_ASSERT_IMAGE_NON_NULL(weight, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, input, false);
-        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
-    }
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
-
-    // Image size
-    int numCols = reference->numCols;
-    int numRows = reference->numRows;
 
     int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation
@@ -491 +476 @@
             psVectorInit(stampVector, 0.0);
 
-            float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord
-            float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord
-            psImage *polyValues = spatialPolyValues(spatialOrder, xNorm, yNorm); // Spatial polynomial terms
-
-            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
-                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
-                    float invNoise2 = 1.0 / (weight ? weight->data.F32[y][x] :
-                                             input->data.F32[y][x]); // Inverse square noise
+            // Spatial polynomial terms
+            psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm);
+
+            psKernel *reference = stamp->reference; // Reference postage stamp
+            psKernel *input = stamp->input; // Input postage stamp
+            psKernel *weight = stamp->weight; // Weight map postage stamp
+
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    float invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse square noise
 
                     // Generate the convolutions
@@ -519 +506 @@
                                 convolutions->data.F64[j] * invNoise2;
                         }
-                        stampVector->data.F64[i] += input->data.F32[y][x] * convolutions->data.F64[i] *
+                        stampVector->data.F64[i] += input->kernel[y][x] * convolutions->data.F64[i] *
                             invNoise2;
 
@@ -527 +514 @@
                     // Background only terms
                     stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2;
-                    stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2;
+                    stampVector->data.F64[bgIndex] += input->kernel[y][x] * invNoise2;
                 }
             }