Changeset 29060 for branches/sc_branches/trunkTest/psphot/src/psphotRadiusChecks.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotRadiusChecks.c (modified) (9 diffs)

branches/sc_branches/trunkTest/psphot/src/psphotRadiusChecks.c

@@ -8 +8 @@
                                         // and a per-object radius is calculated)
 
-bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {
+bool psphotInitRadiusPSF(psMetadata *recipe, pmReadout *readout) {
 
     bool status = true;
@@ -15 +15 @@
     PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING");
 
-    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
+    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, readout->analysis, "PSF_FIT_RADIUS");
     if (!status) {
         PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
     }
 
-    PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
+    PSF_APERTURE =  psMetadataLookupF32(&status, readout->analysis, "PSF_APERTURE");
     if (!status) {
         PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
@@ -28 +28 @@
 
     if (PSF_FIT_RADIUS == 0.0) {
-        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
+        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
         if (!status) {
             gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
@@ -37 +37 @@
 
     if (PSF_APERTURE == 0.0) {
-        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
+        float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
         if (!status) {
             gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
@@ -122 +122 @@
 }
 
+static float EXT_FIT_SKY_SIG;
 static float EXT_FIT_NSIGMA;
 static float EXT_FIT_PADDING;
 static float EXT_FIT_MAX_RADIUS;
 
-bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) {
+bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout) {
 
     bool status;
@@ -134 +135 @@
     EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS");
 
+    float skyMean  = psMetadataLookupF32 (&status, readout->analysis, "SKY_MEAN");
+    float skyStdev = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV");
+
+    fprintf (stderr, "sky: %f +/- %f\n", skyMean, skyStdev);
+
+    EXT_FIT_SKY_SIG = skyStdev;
+
     return true;
 }
 
 // call this function whenever you (re)-define the EXT model
-float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
+bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) {
 
     psAssert (source, "source not defined??");
@@ -146 +154 @@
 
     // set the radius based on the footprint:
-    if (!peak->footprint) goto escape;
+    if (!peak->footprint) return false;
     pmFootprint *footprint = peak->footprint;
-    if (!footprint->spans) goto escape;
-    if (footprint->spans->n < 1) goto escape;
+    if (!footprint->spans) return false;
+    if (footprint->spans->n < 1) return false;
 
     // find the max radius
-    float radius = 0.0;
+    float rawRadius = 0.0;
     for (int j = 0; j < footprint->spans->n; j++) {
         pmSpan *span = footprint->spans->data[j];
@@ -160 +168 @@
         float dX1 = span->x1 - peak->xf;
 
-        radius = PS_MAX (radius, hypot(dY, dX0));
-        radius = PS_MAX (radius, hypot(dY, dX1));
-    }
-
-    radius += EXT_FIT_PADDING;
-    if (isnan(radius)) psAbort("error in radius");
-
-    radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS);
+        rawRadius = PS_MAX (rawRadius, hypot(dY, dX0));
+        rawRadius = PS_MAX (rawRadius, hypot(dY, dX1));
+    }
+    if (isnan(rawRadius)) return false;
+    rawRadius = PS_MIN (factor*rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS);
 
     // redefine the pixels if needed
-    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius);
-
-    // set the mask to flag the excluded pixels
-    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal);
-    return radius;
-
-escape:
-    return NAN;
-    // bool result = psphotCheckRadiusEXT (readout, source, model, markVal);
-    // return result;
+    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, rawRadius);
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, rawRadius, "OR", markVal);
+
+    *radius = rawRadius;
+    return true;
 }
 
 // alternative EXT radius based on model guess (for use without footprints)
-bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) {
-
-    psAbort ("do not use this function");
+bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) {
 
     psF32 *PAR = model->params->data.F32;
@@ -193 +193 @@
 
     // set the fit radius based on the object flux limit and the model
-    float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky);
-
-    model->fitRadius = rawRadius + EXT_FIT_PADDING;
-    if (isnan(model->fitRadius)) psAbort("error in radius");
+    float flux = deep ? EXT_FIT_NSIGMA*EXT_FIT_SKY_SIG : 0.1 * model->params->data.F32[PM_PAR_I0];
+
+    float rawRadius = model->modelRadius (model->params, flux);
+    if (isnan(rawRadius)) return false;
+
+    rawRadius = PS_MIN (rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS);
+    model->fitRadius = rawRadius;
 
     // redefine the pixels if needed
-    bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
+    pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius);
 
     // set the mask to flag the excluded pixels
     psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal);
-    return status;
-}
+    return true;
+}