Changeset 6481 for trunk/psphot/src/psphotSourceFits.c

Timestamp:

Feb 23, 2006, 6:16:11 PM (20 years ago)

Author:

eugene

Message:

fixed some errors with double sources, added second-pass linear PSF fit, multiple-depths

File:

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

trunk/psphot/src/psphotSourceFits.c

@@ -126 +126 @@
     psTrace ("psphot.blend", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[2], ONE->params->data.F32[3]);
 
-    source->modelPSF = (pmModel *) DBL->data[0];
-    source->mode |=  PM_SOURCE_SUBTRACTED;
-    source->mode &= ~PM_SOURCE_TEMPSUB;
-
+    // drop old model, save new second model...
+    psFree (source->modelPSF);
+    source->modelPSF = psMemCopy (DBL->data[0]);
+    source->mode    |= PM_SOURCE_SUBTRACTED;
+    source->mode    |= PM_SOURCE_PAIR;
+    source->mode    &= ~PM_SOURCE_TEMPSUB;
+
+    // copy most data from the primary source (modelEXT, blends stay NULL)
+    pmSource *newSrc = pmSourceAlloc ();
+    newSrc->peak     = psMemCopy (source->peak);
+    newSrc->pixels   = psMemCopy (source->pixels);
+    newSrc->weight   = psMemCopy (source->weight);
+    newSrc->mask     = psMemCopy (source->mask);
+    newSrc->moments  = psMemCopy (source->moments);
+    newSrc->modelPSF = psMemCopy (DBL->data[1]);
+    newSrc->type     = source->type;
+    newSrc->mode     = source->mode;
+    psArrayAdd (sources, 100, newSrc);
+    psFree (newSrc);
+    psFree (DBL);
     return true;
 }
@@ -154 +170 @@
     // save the PSF model from the Ensemble fit
     PSF = source->modelPSF;
-    psphotCheckRadiusPSF (readout, source, PSF);
+    psphotCheckRadiusPSFBlend (readout, source, PSF, 8.0);
 
     modelSet = psArrayAlloc (2);
@@ -160 +176 @@
     DBL = pmModelCopy (PSF);
     DBL->params->data.F32[1] *= 0.5;
-    DBL->params->data.F32[2] += dx;
-    DBL->params->data.F32[3] += dy;
+    DBL->params->data.F32[2] = source->moments->x + dx;
+    DBL->params->data.F32[3] = source->moments->y + dy;
     modelSet->data[0] = DBL;
 
     DBL = pmModelCopy (PSF);
     DBL->params->data.F32[1] *= 0.5;
-    DBL->params->data.F32[2] -= dx;
-    DBL->params->data.F32[3] -= dy;
+    DBL->params->data.F32[2] = source->moments->x - dx;
+    DBL->params->data.F32[3] = source->moments->y - dy;
     modelSet->data[1] = DBL;
 
-    x = PSF->params->data.F32[2];
-    y = PSF->params->data.F32[3];
+    x = source->moments->x;
+    y = source->moments->y;
 
     // fit EXT (not PSF) model (set/unset the pixel mask)
@@ -214 +230 @@
     // save the PSF model from the Ensemble fit
     pmModel *PSF = pmModelCopy (source->modelPSF);
-
-    // extend source radius as needed
-    psphotCheckRadiusPSF (readout, source, PSF);
 
     x = PSF->params->data.F32[2];
@@ -229 +242 @@
     psArrayAdd (sourceSet, 16, source);
 
-    // counter to track the blend elements used in the fit
-    psVector *index  = psVectorAlloc (source->blends->n + 1, PS_TYPE_S32);
-    index->data.S32[0] = -1; // first element corresponds to the primary 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];
 
-        // is this object close enough to include in the fit
-        // XXX this test is bogus: need to think about this
-        dR = hypot (blend->peak->x - x, blend->peak->y - y);
-        if (dR > PSF->radius - 3) continue;
+        // 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
@@ -256 +265 @@
 
         // add this blend to the list
-        index->data.S32[modelSet->n] = i;
         psArrayAdd (modelSet, 16, model);
         psArrayAdd (sourceSet, 16, blend);
     }
+
+    // extend source radius as needed
+    psphotCheckRadiusPSFBlend (readout, source, PSF, dR);
 
     // fit PSF model (set/unset the pixel mask)
@@ -273 +284 @@
         // if we skip this one, free the corresponding blend entry model
         if (!psphotEvalPSF (src, model)) {
-            int n = index->data.S32[i];
-            pmSource *blend = source->blends->data[n];
+            pmSource *blend = source->blends->data[i - 1];
             psFree (blend->modelPSF);
             blend->modelPSF = NULL;
@@ -285 +295 @@
         src->mode &= ~PM_SOURCE_TEMPSUB;
     }
-    psFree (index);
     psFree (modelSet);
     psFree (sourceSet);