Changeset 21150 for branches/pap_branch_20090108/ppStack/src/ppStackMatch.c

Timestamp:

Jan 22, 2009, 9:06:25 AM (17 years ago)

Author:

Paul Price

Message:

Adding new input parameter: limit to the 'maximum deconvolution fraction' --- when this value is in excess of unity (e.g., 1.01 is probably OK, 1.1 is bad), the convolved image tends to be bad. We now detect these and throw them out.

File:

• branches/pap_branch_20090108/ppStack/src/ppStackMatch.c (modified) (4 diffs)

branches/pap_branch_20090108/ppStack/src/ppStackMatch.c

@@ -173 +173 @@
     psAssert(recipe, "We've thrown an error on this before.");
 
+    float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
+
     // Look up appropriate values from the ppSub recipe
     psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
@@ -204 +206 @@
         assert(outName);
         // Read convolution kernel
-        {
-            psString filename = NULL;   // Output filename
-            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
-            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-            psFree(filename);
-            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-            psFree(resolved);
-            if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
-                psFitsClose(fits);
-                return false;
-            }
+        psString filename = NULL;   // Output filename
+        psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+        psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+        psFree(filename);
+        psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+        psFree(resolved);
+        if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
             psFitsClose(fits);
-
-            // Add in variance factor
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
-            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            if (!isfinite(vf)) {
-                vf = 1.0;
-            }
-            if (isfinite(vfItem->data.F32)) {
-                vfItem->data.F32 *= vf;
-            } else {
-                vfItem->data.F32 = vf;
-            }
+            return false;
+        }
+        psFitsClose(fits);
+
+        // Add in variance factor
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
+                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+        float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
+        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+        if (!isfinite(vf)) {
+            vf = 1.0;
+        }
+        if (isfinite(vfItem->data.F32)) {
+            vfItem->data.F32 *= vf;
+        } else {
+            vfItem->data.F32 = vf;
         }
 
@@ -253 +253 @@
         psFree(maskName);
         psFree(weightName);
+
+        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
+                                               PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+
+        pmSubtractionAnalysis(readout->analysis, kernels, region,
+                              readout->image->numCols, readout->image->numRows);
     } else {
 #endif
@@ -535 +541 @@
     }
 
+    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    float deconv = psMetadataLookupF32(NULL, output->analysis,
+                                       PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
+    if (deconv > deconvLimit) {
+        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
+                  deconv, deconvLimit);
+        psFree(output);
+        return NULL;
+    }
+
     // Renormalise the variances if desired
     if (renorm) {