Changeset 28253

Timestamp:

Jun 7, 2010, 4:52:17 PM (16 years ago)

Author:

Paul Price

Message:

Use a common target PSF image. This is to make everything run MUCH faster for stacks of lots of images.

Location:

trunk/ppStack/src

Files:

• Makefile.am (modified) (1 diff)
• ppStack.h (modified) (2 diffs)
• ppStackConvolve.c (modified) (10 diffs)
• ppStackFiles.c (modified) (2 diffs)
• ppStackMatch.c (modified) (10 diffs)
• ppStackTarget.c (added)

trunk/ppStack/src/Makefile.am

@@ -48 +48 @@
         ppStackPhotometry.c     \
         ppStackFinish.c         \
+        ppStackTarget.c         \
         ppStackErrorCodes.c
 

trunk/ppStack/src/ppStack.h

@@ -26 +26 @@
 typedef enum {
     PPSTACK_FILES_PREPARE,              // Files for preparation
+    PPSTACK_FILES_TARGET,               // Files for target generation
     PPSTACK_FILES_CONVOLVE,             // Files for convolution
     PPSTACK_FILES_STACK,                // Stack files
@@ -115 +116 @@
 void ppStackVersionPrint(void);
 
+/// Generate target PSF image
+psImage *ppStackTarget(ppStackOptions *options, // Options for stacking
+                       pmConfig *config         // Configuration
+    );
+
 /// Convolve image to match specified seeing
 bool ppStackMatch(pmReadout *readout,   // Readout to be convolved; replaced with output
+                  const psImage *target,   // Target PSF image
                   ppStackOptions *options, // Options for stacking
                   int index,            // Index of image to match

trunk/ppStack/src/ppStackConvolve.c

@@ -45 +45 @@
     options->convCovars = psArrayAlloc(num); // Covariance matrices
 
+    psImage *target = NULL;             // Target PSF image
+    if (options->convolve) {
+        target = ppStackTarget(options, config);
+        if (!target) {
+            psError(psErrorCodeLast(), false, "Unable to produce stack target image");
+            return false;
+        }
+    }
+
     psVector *renorms = psVectorAlloc(num, PS_TYPE_F32); // Renormalisation values for variances
     psVectorInit(renorms, NAN);
@@ -71 +80 @@
             psFree(fpaList);
             psFree(cellList);
+            psFree(target);
             return false;
         }
@@ -87 +97 @@
             psFree(fpaList);
             psFree(cellList);
+            psFree(target);
             return false;
         }
@@ -93 +104 @@
         psTimerStart("PPSTACK_MATCH");
         options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance);
-        if (!ppStackMatch(readout, options, i, config)) {
+        if (!ppStackMatch(readout, target, options, i, config)) {
             // XXX many things can cause a failure of ppStackMatch -- should some be handled differently?
             psErrorCode error = psErrorCodeLast(); // Error code
@@ -102 +113 @@
               case PPSTACK_ERR_IO:
                 psError(error, false, "Unable to match image %d due to fatal error.", i);
+                psFree(rng);
+                psFree(fpaList);
+                psFree(cellList);
+                psFree(target);
                 return false;
                 // Non-fatal errors
@@ -154 +169 @@
             psFree(cellList);
             psFree(rng);
+            psFree(target);
             return false;
         }
@@ -164 +180 @@
             psFree(rng);
             psFree(maskHeader);
+            psFree(target);
             return false;
         }
@@ -172 +189 @@
             psFree(cellList);
             psFree(rng);
+            psFree(target);
             return false;
         }
@@ -222 +240 @@
             psFree(cellList);
             psFree(rng);
+            psFree(target);
             return false;
         }
@@ -229 +248 @@
     }
     psFree(rng);
+    psFree(target);
 
     psFree(options->sourceLists); options->sourceLists = NULL;

trunk/ppStack/src/ppStackFiles.c

@@ -16 +16 @@
 static char *filesPrepare[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL };
 
+/// Files required for generating convolution target
+static char *filesTarget[] = { "PPSTACK.INPUT.VARIANCE", "PPSTACK.INPUT.MASK", NULL };
+
 /// Files required for the convolution
 static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL };
-
-//                                 "PPSTACK.CONV.KERNEL", NULL };
 
 /// Regular (convolved) stack files
@@ -41 +42 @@
     switch (list) {
       case PPSTACK_FILES_PREPARE:  return filesPrepare;
+      case PPSTACK_FILES_TARGET:   return filesTarget;
       case PPSTACK_FILES_CONVOLVE: return filesConvolve;
       case PPSTACK_FILES_STACK:    return filesStack;

trunk/ppStack/src/ppStackMatch.c

@@ -13 +13 @@
 #define MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
 #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 | 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
 
@@ -49 +46 @@
 #endif
 
-// Get coordinates from a source
-static void coordsFromSource(float *x, float *y, const pmSource *source)
-{
-    assert(x && y);
-    assert(source);
-
-    if (source->modelPSF) {
-        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-    } else {
-        *x = source->peak->xf;
-        *y = source->peak->yf;
-    }
-    return;
-}
-
-static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
-                                   int exclusion // Exclusion zone, pixels
-    )
-{
-    psAssert(sources && sources->n > 0, "Require array of sources");
-    if (exclusion <= 0) {
-        return psMemIncrRefCounter(sources);
-    }
-
-    int num = sources->n;               // Number of sources
-    psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates
-    int numGood = 0;                    // Number of good sources
-    for (int i = 0; i < num; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
-            continue;
-        }
-        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
-        numGood++;
-    }
-    x->n = y->n = numGood;
-
-    psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree
-
-    psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
-    psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source
-    int numFiltered = 0;                // Number of filtered sources
-    for (int i = 0; i < num; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
-            continue;
-        }
-        float xSource, ySource;         // Coordinates of source
-        coordsFromSource(&xSource, &ySource, source);
-
-        coords->data.F64[0] = xSource;
-        coords->data.F64[1] = ySource;
-
-        long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
-        psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
-                coords->data.F64[0], coords->data.F64[1], numWithin);
-        if (numWithin == 1) {
-            // Only itself inside the exclusion zone
-            filtered = psArrayAdd(filtered, filtered->n, source);
-        } else {
-            numFiltered++;
-        }
-    }
-    psFree(coords);
-    psFree(tree);
-    psFree(x);
-    psFree(y);
-
-    psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
-
-    return filtered;
-}
 
 // Add background into the fake image
@@ -202 +126 @@
 
 
-bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
+bool ppStackMatch(pmReadout *readout, const psImage *target, ppStackOptions *options,
+                  int index, const pmConfig *config)
 {
     assert(readout);
@@ -286 +211 @@
         } else {
 #endif
-
             // Normal operations here
-            psAssert(options->psf, "Require target PSF");
-            psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
+            psAssert(target, "Require target PSF image");
 
             int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
@@ -341 +264 @@
             }
 
-            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, 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);
-                psFree(conv);
-                return false;
-            }
-            pmReadoutFakeThreads(oldThreads);
-
+            pmReadout *fake = pmReadoutAlloc(NULL);
+            fake->image = psImageCopy(NULL, target, PS_TYPE_F32);
             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
 
@@ -420 +312 @@
                     psFree(fake);
                     psFree(optWidths);
-                    psFree(stampSources);
                     psFree(conv);
                     if (threads > 0) {
@@ -436 +327 @@
                     psFree(fake);
                     psFree(optWidths);
-                    psFree(stampSources);
                     psFree(conv);
                     psFree(widthsCopy);
@@ -446 +336 @@
 
                 if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
-                                        threshold, stampSources, stampsName, type, size, order, widthsCopy,
+                                        threshold, options->sources, stampsName, type, size, order, widthsCopy,
                                         orders, inner, ringsOrder, binning, penalty,
                                         optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
@@ -454 +344 @@
                     psFree(fake);
                     psFree(optWidths);
-                    psFree(stampSources);
                     psFree(conv);
                     psFree(widthsCopy);
@@ -492 +381 @@
             psFree(fake);
             psFree(optWidths);
-            psFree(stampSources);
 
             if (threads > 0) {