Changeset 27840 for branches/simtest_nebulous_branches/ppStack/src/ppStackMatch.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• ppStack/src (modified) (1 prop)
• ppStack/src/ppStackMatch.c (modified) (22 diffs)

branches/simtest_nebulous_branches/ppStack/src

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h

branches/simtest_nebulous_branches/ppStack/src/ppStackMatch.c

@@ -14 +14 @@
 #define FAKE_SIZE 1                     // Size of fake convolution kernel
 #define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
-                     PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
-#define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux
+                     PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
+#define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
 
-//#define TESTING                         // Enable debugging output
+// #define TESTING                         // Enable debugging output
 
 #ifdef TESTING
@@ -31 +31 @@
     psFree(resolved);
     if (!fits) {
-        psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name);
+        psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name);
         return false;
     }
     psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
     if (!image) {
-        psError(PS_ERR_IO, false, "Unable to read previously produced image: %s", name);
+        psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name);
         psFitsClose(fits);
         return false;
@@ -87 +87 @@
     x->n = y->n = numGood;
 
-    psTree *tree = psTreePlant(2, 2, x, y); // kd tree
+    psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree
 
     psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
@@ -115 +115 @@
     psFree(coords);
     psFree(tree);
+    psFree(x);
+    psFree(y);
 
     psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
@@ -144 +146 @@
     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
@@ -150 +152 @@
     psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
     if (!psImageUnbin(unbinned, binned, binning)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model");
+        psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model");
+        psFree(binned);
         psFree(unbinned);
         return NULL;
     }
-
-    // XXX should these really be here?? (probably not...)
-    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
-    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
-
-     return unbinned;
+    psFree(binned);
+
+    return unbinned;
 }
+
+// Renormalise a readout's variance map
+bool stackRenormaliseReadout(const pmConfig *config, // Configuration
+                             pmReadout *readout      // Readout to renormalise
+    )
+{
+#if 1
+    bool mdok; // Status of metadata lookups
+
+    psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack
+    psAssert(recipe, "Need PPSTACK recipe");
+
+    if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
+
+    int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
+    if (!mdok) {
+        psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
+        return false;
+    }
+    float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
+    if (!mdok) {
+        psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
+        return false;
+    }
+    float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
+    if (!mdok) {
+        psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
+        return false;
+    }
+
+    psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
+
+    psImageCovarianceTransfer(readout->variance, readout->covariance);
+    return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
+#else
+    return true;
+#endif
+}
+
 
 
@@ -178 +217 @@
     int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
 
-    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in
     psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
     psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
@@ -190 +229 @@
 
     if (!pmReadoutMaskNonfinite(readout, maskVal)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in readout.");
+        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout.");
         return false;
     }
@@ -217 +256 @@
             psFree(resolved);
             if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
+                psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel");
                 psFitsClose(fits);
                 return false;
@@ -223 +262 @@
             psFitsClose(fits);
 
-            if (!readImage(&readout->image, options->imageNames->data[index], config) ||
-                !readImage(&readout->mask, options->maskNames->data[index], config) ||
-                !readImage(&readout->variance, options->varianceNames->data[index], config)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced image.");
+            if (!readImage(&readout->image, options->convImages->data[index], config) ||
+                !readImage(&readout->mask, options->convMasks->data[index], config) ||
+                !readImage(&readout->variance, options->convVariances->data[index], config)) {
+                psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image.");
                 return false;
             }
@@ -232 +271 @@
             psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
                                                    PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
-
-            pmSubtractionAnalysis(readout->analysis, kernels, region,
+            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
+                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);
+
+            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;
@@ -256 +299 @@
             int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
             float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
-            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
+            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
@@ -273 +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(PPSTACK_ERR_CONFIG, 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
@@ -283 +343 @@
             pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
 
+            psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
+            psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+            if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
+                psError(PPSTACK_ERR_DATA, false, "Can't measure background for image.");
+                psFree(fake);
+                psFree(optWidths);
+                psFree(conv);
+                psFree(bg);
+                psFree(rng);
+                return false;
+            }
+            float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
+            psFree(rng);
+            psFree(bg);
+
             // For the sake of stamps, remove nearby sources
             psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
                                                        footprint); // Filtered list of sources
 
+            bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, NULL, NULL, options->psf, NAN, footprint + size,
-                                          false, true)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+                                          stampSources, SOURCE_MASK, NULL, NULL, options->psf,
+                                          minFlux, footprint + size, false, true)) {
+                psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
                 psFree(optWidths);
@@ -296 +372 @@
                 return false;
             }
+            pmReadoutFakeThreads(oldThreads);
 
             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
@@ -328 +407 @@
 
             if (threads > 0) {
-                pmSubtractionThreadsInit(readout, fake);
+                pmSubtractionThreadsInit();
             }
 
@@ -335 +414 @@
                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
             if (kernel) {
-                if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis,
-                                               stride, sysError, maskVal, maskBad, maskPoor,
+                if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
+                                               stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
                                                poorFrac, badFrac)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to convolve images.");
+                    psError(psErrorCodeLast(), false, "Unable to convolve images.");
                     psFree(fake);
                     psFree(optWidths);
@@ -344 +423 @@
                     psFree(conv);
                     if (threads > 0) {
-                        pmSubtractionThreadsFinalize(readout, fake);
+                        pmSubtractionThreadsFinalize();
                     }
                     return false;
                 }
             } else {
-                if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing,
-                                        threshold, stampSources, stampsName, type, size, order, widths,
-                                        orders, inner, ringsOrder, binning, penalty,
-                                        optimum, optWidths, optOrder, optThresh, iter, rej, sysError,
-                                        maskVal, maskBad, maskPoor, poorFrac, badFrac,
-                                        PM_SUBTRACTION_MODE_1)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+                // 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(psErrorCodeLast(), false, "Unable to scale kernel parameters");
                     psFree(fake);
                     psFree(optWidths);
                     psFree(stampSources);
                     psFree(conv);
+                    psFree(widthsCopy);
                     if (threads > 0) {
-                        pmSubtractionThreadsFinalize(readout, fake);
+                        pmSubtractionThreadsFinalize();
                     }
                     return false;
                 }
-            }
+
+                if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
+                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
+                                        orders, inner, ringsOrder, binning, penalty,
+                                        optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
+                                        sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
+                                        poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
+                    psError(psErrorCodeLast(), false, "Unable to match images.");
+                    psFree(fake);
+                    psFree(optWidths);
+                    psFree(stampSources);
+                    psFree(conv);
+                    psFree(widthsCopy);
+                    if (threads > 0) {
+                        pmSubtractionThreadsFinalize();
+                    }
+                    return false;
+                }
+                psFree(widthsCopy);
+            }
+
 
 #ifdef TESTING
@@ -396 +495 @@
 
             if (threads > 0) {
-                pmSubtractionThreadsFinalize(readout, fake);
+                pmSubtractionThreadsFinalize();
             }
 
@@ -436 +535 @@
             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);
             }
@@ -467 +561 @@
         }
 
+        // 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;
@@ -491 +610 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
-      psWarning("Can't measure background for image.");
-      psErrorClear();
+        psWarning("Can't measure background for image.");
+        psErrorClear();
     } else {
-      if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
-        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
-                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
-        (void)psBinaryOp(readout->image, readout->image, "-",
-                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
-      }
-    }
-
-
-    // Measure the variance level for the weighting
-    if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
+        if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) {
+            psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
+                     psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
+            (void)psBinaryOp(readout->image, readout->image, "-",
+                             psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
+        }
+    }
+
+    if (!stackRenormaliseReadout(config, readout)) {
         psFree(rng);
         psFree(bg);
         return false;
     }
-    options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
-                                                  psImageCovarianceFactor(readout->covariance));
-    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
-                     "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]);
+
+    // Measure the variance level for the weighting
+    if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) {
+        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
+            psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
+            psFree(rng);
+            psFree(bg);
+            return false;
+        }
+        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
+                                                      psImageCovarianceFactor(readout->covariance));
+    } else {
+        options->weightings->data.F32[index] = 1.0;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n",
+             index, options->weightings->data.F32[index]);
 
     psFree(rng);

@@ -10 +10 @@
 stamp-h1
 ppStackVersionDefinitions.h
+ppStackErrorCodes.c
+ppStackErrorCodes.h