Changeset 14329

Timestamp:

Jul 19, 2007, 3:53:41 PM (19 years ago)

Author:

Paul Price

Message:

Not using psPolynomials to do spatial polynomials --- it's awkward,
confusing (x,y or y,x?) and I'm afraid the way it behaves will get
changed underneath me.
Adding function to generate image of each of the components of the
solution kernel.
Ensuring entire image gets convolved.

File:

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

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-19 21:42:12 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-20 01:53:41 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -41 +41 @@
 
     psImage *polyValues = psImageAlloc(spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); // Matrix to return
-    psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, spatialOrder, spatialOrder); // Polynomial
-
-    for (int j = 0; j <= spatialOrder; j++) {
-        for (int i = 0; i <= spatialOrder - j; i++) {
-            poly->coeff[i][j] = 1.0;
-            polyValues->data.F64[j][i] = psPolynomial2DEval(poly, x, y);
-            poly->coeff[i][j] = 0.0;
-        }
-    }
-    psFree(poly);
+    polyValues->data.F64[0][0] = 1.0;
+
+    double value = 1.0;
+    for (int i = 1; i <= spatialOrder; i++) {
+        value *= x;
+        polyValues->data.F64[0][i] = value;
+    }
+
+    value = 1.0;
+    for (int j = 1; j <= spatialOrder; j++) {
+        value *= y;
+        polyValues->data.F64[j][0] = value;
+    }
+
+    for (int j = 1; j <= spatialOrder; j++) {
+        for (int i = 1; i <= spatialOrder - j; i++) {
+            polyValues->data.F64[j][i] = polyValues->data.F64[j][0] * polyValues->data.F64[0][i];
+        }
+    }
 
     return polyValues;
@@ -875 +884 @@
 }
 
+psArray *pmSubtractionKernelSolutions(const psVector *solution, const pmSubtractionKernels *kernels,
+                                      float x, float y
+                                      )
+{
+    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
+    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
+    PS_ASSERT_VECTOR_SIZE(kernels->u, solution->n - 1, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->xOrder, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->yOrder, NULL);
+    if (kernels->type == PM_SUBTRACTION_KERNEL_ISIS) {
+        PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, NULL);
+        PS_ASSERT_ARRAY_SIZE(kernels->preCalc, solution->n - 1, NULL);
+    }
+
+    psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return
+    psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector
+    psVectorInit(fakeSolution, 0.0);
+
+    for (int i = 0; i < solution->n - 1; i++) {
+        fakeSolution->data.F64[i] = solution->data.F64[i];
+        images->data[i] = pmSubtractionKernelImage(fakeSolution, kernels, x, y);
+        fakeSolution->data.F64[i] = 0.0;
+    }
+
+    psFree(fakeSolution);
+
+    return images;
+}
+
+
 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
                            const psImage *inImage, const psImage *inWeight, const psImage *subMask,
@@ -953 +995 @@
     psKernel *kernelWeight = NULL;      // Kernel for the weight map
 
-    for (int j = size; j < inImage->numRows - fullSize - size; j += fullSize) {
-        for (int i = size; i < inImage->numCols - fullSize - size; i += fullSize) {
+    for (int j = size; j < inImage->numRows - size; j += fullSize) {
+        for (int i = size; i < inImage->numCols - size; i += fullSize) {
 
             // Only generate polynomial values every kernel footprint, since we have already assumed
@@ -960 +1002 @@
             psFree(polyValues);
             polyValues = spatialPolyValues(kernels->spatialOrder,
-                                           2.0 * (float)(i + size - numCols/2.0) / (float)numCols,
-                                           2.0 * (float)(j + size - numRows/2.0) / (float)numRows);
+                                           2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols,
+                                           2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);
 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
@@ -975 +1017 @@
 #endif
 
-            for (int y = j; y < j + fullSize; y++) {
-                for (int x = i; x < i + fullSize; x++) {
+            for (int y = j; y < PS_MIN(j + fullSize, numRows - size); y++) {
+                for (int x = i; x < PS_MIN(i + fullSize, numCols - size); x++) {
                     // Check and propagate the kernel footprint, if required
                     if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &