Changeset 28667 for trunk/psLib/src/imageops

Timestamp:

Jul 13, 2010, 2:34:04 PM (16 years ago)

Author:

Paul Price

Message:

It's bothered me for a while that the noise model that we have been using is only kind of spatially variable. The variance map is intended to track the change in the noise properties as a function of position, and while it does keep track of different noise levels due to, e.g., different detectors, vignetting, etc., it has not tracked the different noise that comes from convolution kernels changing over an image. This had been an intentional oversight --- the prescription we were using for calculating the variance and covariance removed the spatial variation (by prescribing that the variance level under a convolution remain unchanged) and we were sweeping it under the carpet so we could calculate the covariance pseudo-matrix and use that for setting the absolute level of the noise.

I've now figured out how to have our cake and eat it too. By changing where the scalings are applied, the variance map can track the changing noise introduced by a varying convolution kernel, while the covariance pseudo-matrix tracks an absolute change. This fix is purely algorithmic --- it does not affect any of the variances/covariances stored on disk, nor is it affected by them. I hope that this means that any Magic recertification required is either minimal or completely unnecessary.

Again, all that has changed is where the scalings are applied.

Location:

trunk/psLib/src/imageops

Files:

• psImageCovariance.c (modified) (15 diffs)
• psImageInterpolate.c (modified) (4 diffs)

trunk/psLib/src/imageops/psImageCovariance.c

@@ -34 +34 @@
 static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
                                       const psKernel *kernel, // Convolution kernel
-                                      int x, int y            // Coordinates in output covariance matrix
+                                      int x, int y,           // Coordinates in output covariance matrix
+                                      float scale             // Scale to apply
                                       )
 {
@@ -83 +84 @@
     }
 
-    return sum;
+    return scale * sum;
 }
 
@@ -91 +92 @@
     PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     psAssert(job->args, "No job arguments");
-    psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
+    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
 
     psKernel *out = job->args->data[0]; // Output covariance matrix
@@ -98 +99 @@
     int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
     int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
-
-    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+    float scale = PS_SCALAR_VALUE(job->args->data[5], F32); // Scaling to apply
+
+    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
 
     return true;
@@ -127 +129 @@
 
     // Check for non-finite elements
+    double sumKernel = 0.0, sumKernel2 = 0.0; // Sum of the kernel
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
@@ -135 +138 @@
                 return NULL;
             }
+            sumKernel += kernel->kernel[y][x];
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -155 +160 @@
     int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
     psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+    float scale = 1.0 / sumKernel2;          // Scaling to apply
 
     for (int y = yMin; y <= yMax; y++) {
@@ -165 +171 @@
                 PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
                 PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, scale, PS_TYPE_F32);
                 if (!psThreadJobAddPending(job)) {
                     psFree(covar);
@@ -170 +177 @@
                 }
             } else {
-                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
-            }
-        }
-    }
-    psFree(covar);
+                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y, scale);
+            }
+        }
+    }
 
     if (threaded && !psThreadPoolWait(true)) {
@@ -180 +186 @@
         return false;
     }
+
+    psFree(covar);
 
     return out;
@@ -194 +202 @@
 
     // Check for non-finite elements
+    double sumKernel2 = 0.0; // Sum of the squared kernel
     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
@@ -202 +211 @@
                 return NAN;
             }
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
         }
     }
@@ -215 +225 @@
     }
 
-    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
+    float scale = 1.0 / sumKernel2;     // Scale to apply
+    float factor = imageCovarianceCalculate(covar, kernel, 0, 0, scale); // Covariance factor
     psFree(covar);
     return factor;
@@ -338 +349 @@
         }
     }
-    psFree(covar);
-
     if (threaded && !psThreadPoolWait(true)) {
         psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
         return false;
     }
+    psFree(covar);
 
     return out;
@@ -616 +626 @@
     if (set && !threaded) {
         {
-            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 6);
             task->function = &imageCovarianceCalculateThread;
             psThreadTaskAdd(task);

trunk/psLib/src/imageops/psImageInterpolate.c

@@ -427 +427 @@
 
 // Determine the result of the interpolation after all the math has been done
-#define INTERPOLATE_RESULT() \
-    psImageInterpolateStatus status = PS_INTERPOLATE_STATUS_ERROR; /* Status of interpolation */ \
-    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage; \
-    if (wantVariance) { \
-        *varianceValue = sumVariance / (sumKernel2 - sumBad); \
-    } \
-    if (sumKernel == 0.0) { \
-        /* No kernel contributions */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->badMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_BAD; \
-    } else if (sumBad == 0) { \
-        /* Completely good pixel */ \
-        status = PS_INTERPOLATE_STATUS_GOOD; \
-    } else if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) { \
-        /* Some pixels masked: poor pixel */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->poorMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_POOR; \
-    } else { \
-        /* Many pixels (or a few important pixels) masked: bad pixel */ \
-        if (haveMask && maskValue) { \
-            *maskValue |= interp->badMask; \
-        } \
-        status = PS_INTERPOLATE_STATUS_BAD; \
-    }
+static psImageInterpolateStatus interpolateResult(const psImageInterpolation *interp,
+                                                  double *imageValue, double *varianceValue,
+                                                  psImageMaskType *maskValue,
+                                                  double sumImage, double sumVariance, double sumBad,
+                                                  double sumKernel, double sumKernel2,
+                                                  bool wantVariance, bool haveMask)
+{
+    *imageValue = sumKernel > 0 ? sumImage / sumKernel : interp->badImage;
+    if (wantVariance) {
+        if (sumBad > 0) {
+            sumVariance *= sumKernel2 / (sumKernel2 - sumBad);
+        }
+        *varianceValue = sumVariance / PS_SQR(sumKernel);
+    }
+    if (sumKernel == 0.0) {
+        // No kernel contributions at all
+        if (haveMask && maskValue) {
+            *maskValue |= interp->badMask;
+        }
+        return PS_INTERPOLATE_STATUS_BAD;
+    }
+    if (sumBad == 0) {
+        // Completely good pixel
+        return PS_INTERPOLATE_STATUS_GOOD;
+    }
+    if (sumBad < PS_SQR(interp->poorFrac) * sumKernel2) {
+        // Some pixels masked: poor pixel
+        if (haveMask && maskValue) {
+            *maskValue |= interp->poorMask;
+        }
+        return PS_INTERPOLATE_STATUS_POOR;
+    }
+    // Many pixels (or a few important pixels) masked: bad pixel
+    if (haveMask && maskValue) {
+        *maskValue |= interp->badMask;
+    }
+    return PS_INTERPOLATE_STATUS_BAD;
+}
 
 // Interpolation engine for separable interpolation kernels
@@ -703 +713 @@
     }
 
-    INTERPOLATE_RESULT();
-
     psFree(xKernelNew);
     psFree(yKernelNew);
@@ -710 +718 @@
     psFree(yKernel2New);
 
-    return status;
+    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
+                             sumKernel, sumKernel2, wantVariance, haveMask);
 }
 
@@ -861 +870 @@
     }
 
-    INTERPOLATE_RESULT();
-
-    return status;
+    return interpolateResult(interp, imageValue, varianceValue, maskValue, sumImage, sumVariance, sumBad,
+                             sumKernel, sumKernel2, wantVariance, haveMask);
 }