Changeset 26898 for trunk/ppStack/src/ppStackMatch.c

Timestamp:

Feb 10, 2010, 7:41:23 PM (16 years ago)

Author:

eugene

Message:

updates from eam_branches/20091201

File:

• trunk/ppStack/src/ppStackMatch.c (modified) (11 diffs)

trunk/ppStack/src/ppStackMatch.c

@@ -144 +144 @@
     psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
 
-    psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
+    psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model
     psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis,
                                                   "PSPHOT.BACKGROUND.BINNING"); // Binning for model
@@ -274 +274 @@
                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);
 
-            pmSubtractionAnalysis(readout->analysis, NULL, kernels, region,
+            pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
                                   readout->image->numCols, readout->image->numRows);
 
             psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
+            bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
             psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
+            psImageCovarianceSetThreads(oldThreads);
             psFree(readout->covariance);
             readout->covariance = covar;
@@ -298 +300 @@
             float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
             float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel
+            float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw
             float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images
+            float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky
+            float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation
+
             const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
             pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
@@ -315 +321 @@
             float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
 
+            bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE");        // Scale kernel parameters?
+            float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling
+            float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling
+            float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling
+            if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                        "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
+                        scaleRef, scaleMin, scaleMax);
+                return false;
+            }
+
+
             // These values are specified specifically for stacking
             const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
@@ -358 +376 @@
             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
 
+#if 1
             // Add the background into the target image
             psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
+#endif
 
 #ifdef TESTING
@@ -395 +415 @@
             if (kernel) {
                 if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
-                                               stride, kernelError, maskVal, maskBad, maskPoor,
+                                               stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
                                                poorFrac, badFrac)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
@@ -408 +428 @@
                 }
             } else {
+                // Scale the input parameters
+                psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
+                if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy,
+                                                       options->inputSeeing->data.F32[index],
+                                                       options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to scale kernel parameters");
+                    psFree(fake);
+                    psFree(optWidths);
+                    psFree(stampSources);
+                    psFree(conv);
+                    psFree(widthsCopy);
+                    if (threads > 0) {
+                        pmSubtractionThreadsFinalize(readout, fake);
+                    }
+                    return false;
+                }
+
                 if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
-                                        threshold, stampSources, stampsName, type, size, order, widths,
+                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
                                         orders, inner, ringsOrder, binning, penalty,
-                                        optimum, optWidths, optOrder, optThresh, iter, rej, sysError,
-                                        kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac,
-                                        PM_SUBTRACTION_MODE_2)) {
+                                        optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
+                                        sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
+                                        poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
                     psFree(fake);
@@ -419 +456 @@
                     psFree(stampSources);
                     psFree(conv);
+                    psFree(widthsCopy);
                     if (threads > 0) {
                         pmSubtractionThreadsFinalize(readout, fake);
@@ -424 +462 @@
                     return false;
                 }
+                psFree(widthsCopy);
             }
 
@@ -495 +534 @@
             while ((item = psMetadataGetAndIncrement(iter))) {
                 assert(item->type == PS_DATA_UNKNOWN);
-                // Set the normalisation dimensions, since these will be otherwise unavailable when reading
-                // the images by scans.
                 pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
-                kernel->numCols = readout->image->numCols;
-                kernel->numRows = readout->image->numRows;
-
                 kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
             }
@@ -526 +560 @@
         }
 
+        // Kernel normalisation
+        {
+            double sum = 0.0;           // Sum of chi^2
+            int num = 0;                // Number of measurements of chi^2
+            psString regex = NULL;      // Regular expression
+            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
+            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
+            psFree(regex);
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_TYPE_F32);
+                float norm = item->data.F32; // Normalisation
+                sum += norm;
+                num++;
+            }
+            psFree(iter);
+            float conv = sum/num;       // Mean normalisation from convolution
+            float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
+            float renorm =  stars / conv; // Renormalisation to apply
+            psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
+                     index, renorm, conv, stars);
+            psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
+            psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
+        }
+
         // Reject image completely if the maximum deconvolution fraction exceeds the limit
         float deconv = psMetadataLookupF32(NULL, conv->analysis,
                                            PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
         if (deconv > deconvLimit) {
-            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
-                      deconv, deconvLimit);
+            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n",
+                      deconv, deconvLimit, index);
             psFree(conv);
             return NULL;