Changeset 7437

Timestamp:

Jun 8, 2006, 10:57:13 AM (20 years ago)

Author:

eugene

Message:

moved functions to correspond to order used

File:

• trunk/psphot/src/psphotSourceFits.c (modified) (4 diffs)

trunk/psphot/src/psphotSourceFits.c

@@ -2 +2 @@
 
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
+
+bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 
+
+    float x, y, dR;
+
+    // if this source is not a possible blend, just fit as PSF
+    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
+        bool status = psphotFitPSF (readout, source, psf);
+        return status;
+    }
+    psTrace ("psphot.blend", 5, "trying blend...\n");
+
+    // save the PSF model from the Ensemble fit
+    pmModel *PSF = pmModelCopy (source->modelPSF);
+    if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in blend fit primary");
+
+    x = PSF->params->data.F32[2];
+    y = PSF->params->data.F32[3];
+
+    psArray *modelSet = psArrayAlloc (source->blends->n + 1);
+    psArrayAdd (modelSet, 16, PSF);
+    
+    psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
+    psArrayAdd (sourceSet, 16, source);
+
+    // we need to include all blends in the fit (unless primary is saturated?)
+    dR = 0;
+    for (int i = 0; i < source->blends->n; i++) {
+        pmSource *blend = source->blends->data[i];
+
+        // find the blend which is furthest from source
+        dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
+
+        // create the model and guess parameters for this blend
+        pmModel *model = pmModelAlloc (PSF->type);
+        for (int j = 0; j < model->params->n; j++) {
+            model->params->data.F32[j] = PSF->params->data.F32[j];
+            model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
+        }
+
+        // XXX assume local sky is 0.0?
+        model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
+        if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
+        model->params->data.F32[2] = blend->peak->x;
+        model->params->data.F32[3] = blend->peak->y;
+
+        // add this blend to the list
+        psArrayAdd (modelSet, 16, model);
+        psArrayAdd (sourceSet, 16, blend);
+
+        // free to avoid double counting model 
+        psFree (model);
+    }
+
+    // extend source radius as needed
+    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
+
+    // fit PSF model (set/unset the pixel mask)
+    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
+    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
+    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
+
+    // correct model chisq for flux trend
+    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
+    PSF->chisqNorm = PSF->chisq / chiTrend;
+
+    // evaluate the blend objects, subtract if good, free otherwise
+    for (int i = 1; i < modelSet->n; i++) {
+        pmSource *blend = sourceSet->data[i];
+        pmModel *model  = modelSet->data[i];
+
+        // correct model chisq for flux trend
+        chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
+        model->chisqNorm = model->chisq / chiTrend;
+
+        // if this one failed, skip it
+        if (!psphotEvalPSF (blend, model)) {
+            psTrace ("psphot.blend", 5, "failed on blend %d of %d\n", i, modelSet->n);
+            continue;
+        }
+
+        // otherwise, supply the resulting model to the corresponding blend
+        blend->modelPSF = psMemIncrRefCounter (model);
+        psTrace ("psphot.blend", 5, "fitted blend as PSF\n");
+        pmModelSub (source->pixels, source->mask, model, false, false);
+        blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
+        blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
+    }
+    psFree (modelSet);
+    psFree (sourceSet);
+
+    // evaluate the primary object
+    if (!psphotEvalPSF (source, PSF)) {
+        psTrace ("psphot.blend", 5, "failed on blend 0 of %d\n", modelSet->n);
+        psFree (PSF);
+        return false;
+    }
+
+    psTrace ("psphot.blend", 5, "fitted primary as PSF\n");
+    pmModelSub (source->pixels, source->mask, PSF, false, false);
+    psFree (source->modelPSF);
+    source->modelPSF = PSF;
+    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
+    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
+    return true;
+}
 
 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { 
@@ -26 +132 @@
     chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
     PSF->chisqNorm = PSF->chisq / chiTrend;
-    // PSF->chisqNorm = PSF->chisq;
 
     // does the PSF model succeed?
@@ -43 +148 @@
     source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
     source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
-    return true;
-}
-
+
+    return true;
+}
 
 static float EXT_MIN_SN;
@@ -244 +349 @@
 }
 
-bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 
-
-    float x, y, dR;
-
-    // if this source is not a possible blend, just fit as PSF
-    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
-        bool status = psphotFitPSF (readout, source, psf);
-        return status;
-    }
-    psTrace ("psphot.blend", 5, "trying blend...\n");
-
-    // save the PSF model from the Ensemble fit
-    pmModel *PSF = pmModelCopy (source->modelPSF);
-    if (isnan(PSF->params->data.F32[1])) psAbort ("psphot", "nan in blend fit primary");
-
-    x = PSF->params->data.F32[2];
-    y = PSF->params->data.F32[3];
-
-    psArray *modelSet = psArrayAlloc (source->blends->n + 1);
-    psArrayAdd (modelSet, 16, PSF);
-    
-    psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
-    psArrayAdd (sourceSet, 16, source);
-
-    // we need to include all blends in the fit (unless primary is saturated?)
-    dR = 0;
-    for (int i = 0; i < source->blends->n; i++) {
-        pmSource *blend = source->blends->data[i];
-
-        // find the blend which is furthest from source
-        dR = PS_MAX (dR, hypot (blend->peak->x - x, blend->peak->y - y));
-
-        // create the model and guess parameters for this blend
-        pmModel *model = pmModelAlloc (PSF->type);
-        for (int j = 0; j < model->params->n; j++) {
-            model->params->data.F32[j] = PSF->params->data.F32[j];
-            model->dparams->data.F32[j] = PSF->dparams->data.F32[j];
-        }
-
-        // XXX assume local sky is 0.0?
-        model->params->data.F32[1] = blend->moments->Peak - blend->moments->Sky;
-        if (isnan(model->params->data.F32[1])) psAbort ("psphot", "nan in blend fit");
-        model->params->data.F32[2] = blend->peak->x;
-        model->params->data.F32[3] = blend->peak->y;
-
-        // add this blend to the list
-        psArrayAdd (modelSet, 16, model);
-        psArrayAdd (sourceSet, 16, blend);
-
-        // free to avoid double counting model 
-        psFree (model);
-    }
-
-    // extend source radius as needed
-    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
-
-    // fit PSF model (set/unset the pixel mask)
-    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_MASK_MARK);
-    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
-    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", NOT_U8(PM_MASK_MARK));
-
-    // correct model chisq for flux trend
-    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
-    PSF->chisqNorm = PSF->chisq / chiTrend;
-
-    // evaluate the blend objects, subtract if good, free otherwise
-    for (int i = 1; i < modelSet->n; i++) {
-        pmSource *blend = sourceSet->data[i];
-        pmModel *model  = modelSet->data[i];
-
-        // correct model chisq for flux trend
-        chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
-        model->chisqNorm = model->chisq / chiTrend;
-
-        // if this one failed, skip it
-        if (!psphotEvalPSF (blend, model)) {
-            psTrace ("psphot.blend", 5, "failed on blend %d of %d\n", i, modelSet->n);
-            continue;
-        }
-
-        // otherwise, supply the resulting model to the corresponding blend
-        blend->modelPSF = psMemIncrRefCounter (model);
-        psTrace ("psphot.blend", 5, "fitted blend as PSF\n");
-        pmModelSub (source->pixels, source->mask, model, false, false);
-        blend->mode |=  PM_SOURCE_MODE_SUBTRACTED;
-        blend->mode &= ~PM_SOURCE_MODE_TEMPSUB;
-    }
-    psFree (modelSet);
-    psFree (sourceSet);
-
-    // evaluate the primary object
-    if (!psphotEvalPSF (source, PSF)) {
-        psTrace ("psphot.blend", 5, "failed on blend 0 of %d\n", modelSet->n);
-        psFree (PSF);
-        return false;
-    }
-
-    psTrace ("psphot.blend", 5, "fitted primary as PSF\n");
-    pmModelSub (source->pixels, source->mask, PSF, false, false);
-    psFree (source->modelPSF);
-    source->modelPSF = PSF;
-    source->mode |=  PM_SOURCE_MODE_SUBTRACTED;
-    source->mode &= ~PM_SOURCE_MODE_TEMPSUB;
-    return true;
-}