Changeset 34226 for trunk/psphot/src

Timestamp:

Jul 31, 2012, 9:31:06 AM (14 years ago)

Author:

bills

Message:

If number of objects is above a recipe limit reduce the number of radial
apertures processed. Also add recipe value to cause final linear fit
for detected sources be performed first ommitting negative flux sources
followed by the fit for matched sources. Set to false for now pending
further testing

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotRadialBins.c (modified) (1 diff)
• psphotStackReadout.c (modified) (3 diffs)

trunk/psphot/src/psphot.h

@@ -264 +264 @@
 bool            psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS);
 bool            psphotEllipticalContour (pmSource *source);
+bool            psphotLimitRadialApertures(psMetadata *recipe, long nSources);
 
 // psphotVisual functions

trunk/psphot/src/psphotRadialBins.c

@@ -202 +202 @@
 }
 
+// If the number of sources is large, edit the recipe to remove radial bins beyond
+// a value specified in the recipe
+bool psphotLimitRadialApertures(psMetadata *recipe, long nSources) {
+
+    bool status = false;
+    long sourceLimit = psMetadataLookupS32 (&status, recipe, "RADIAL.NUM.SOURCES.LIMIT");
+    if (!status) {
+        sourceLimit = 50000;
+    }
+    if (nSources < sourceLimit) {
+        return true;
+    }
+    psF32 maxRadius = psMetadataLookupF32 (&status, recipe, "RADIAL.SOURCES.OVER.LIMIT.RADIUS");
+    if (!status) {
+        maxRadius = 20.0;
+    }
+    psLogMsg ("psphot", PS_LOG_INFO, "Number of objects: %ld is greater than limit %ld. Limiting radial annular bins to %.3f pixels\n", 
+            nSources, sourceLimit, maxRadius);
+
+    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    if (!radMin || !radMin->n) {
+        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)");
+        return false;
+    }
+    if (!radMax || !radMax->n) {
+        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)");
+        return false;
+    }
+    if (radMax->n != radMin->n) {
+        psError (PSPHOT_ERR_CONFIG, true, "length of radMin %ld and radMax %ld not equal)", radMin->n, radMax->n);
+        return false;
+    }
+    int i = 0;
+    psF32 lastRadius = 0;
+    for (; i < radMax->n; i++) {
+        if (radMax->data.F32[i] > maxRadius) {
+            break;
+        }
+        lastRadius = radMax->data.F32[i];
+    }
+    if (i == radMax->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "Radius of all bins is within the limit. lastRadius: %.3f\n", lastRadius);
+        return true;
+    }
+    psLogMsg ("psphot", PS_LOG_INFO, "  radius limit exceeded at bin %d of %ld. New max radius is %.3f\n",
+            i, radMax->n, lastRadius);
+
+    psVector *radMinNew = psVectorRealloc(radMin, i);
+    if (radMinNew != radMin) {
+        psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.LOWER", PS_META_REPLACE, "", radMinNew);
+        // XXX: I don't need this so I?
+        // psFree(radMinNew);
+    }
+
+    psVector *radMaxNew = psVectorRealloc(radMax, i);
+    if (radMaxNew != radMax) {
+        psMetadataAddVector(recipe, PS_LIST_TAIL, "RADIAL.ANNULAR.BINS.UPPER", PS_META_REPLACE, "", radMaxNew);
+        // XXX: I don't need to free this do I?
+        // psFree(radMaxNew);
+    }
+
+    return true;
+}
+
 // the area-weighted mean radius is given by:
 

trunk/psphot/src/psphotStackReadout.c

@@ -286 +286 @@
     }
 
+    // gcc doesn't like the label to refer to a declaration so we declare this here
+    bool splitLinearFit;
+
 pass1finish:
 
-#ifdef notyet
-    // Split the fit of detected sources from matched sources. But not yet.
-    // NOTE: apply to ALL sources. Only include sources with postitive flux in the fit
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
-#endif
+    splitLinearFit = psMetadataLookupBool(NULL, recipe, "PSPHOT.STACK.SPLIT.LINEAR.FIT");
+    if (splitLinearFit) {
+        psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
+        // Fit the detected sources separately from matched that wea are about to create.
+        // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
+        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
+    }
 
     // generate the objects (objects unify the sources from the different images) NOTE: could
@@ -300 +305 @@
     psMemDump("matchsources");
 
+    // check the source density. If it too high change the number of radial bins
+    // in the recipe.
+    psphotLimitRadialApertures(recipe, objects->n);
+
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
@@ -308 +317 @@
     psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
 
-#ifdef notyet
-    // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
-    // not be included in this fit.
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
-#else
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
-#endif
+    if (splitLinearFit) {
+        // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will
+        // not be included in this fit.
+        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
+    } else {
+        // Fit all sources together
+        psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources)
+    }
 
     // measure the radial profiles to the sky (only measures new objects)