Changeset 5704

Timestamp:

Dec 6, 2005, 7:57:54 AM (21 years ago)

Author:

eugene

Message:

adding basic deblending

Location:

trunk/psphot/src

Files:

• psphotBasicDeblend.c (modified) (1 diff)
• psphotEnsemblePSF.c (modified) (1 diff)

trunk/psphot/src/psphotBasicDeblend.c

@@ -6 +6 @@
 { 
 
-    // sources must be sorted by Y
-    // index must be sorted by SN
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortByY);
 
-    // search for sources which fall within current box
+    psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32);
+    for (int i = 0; i < SN->n; i++) {
+      source = sources->data[i];
+      SN->data.F32[i] = source->moments->SN;
+    }
+    psVector *index = psVectorSortIndex (NULL, SN);
+    
+    psArray *overlap = psArrayAlloc (100);
 
-    // generate source contour (1/4 peak counts)
-    
-    // find sources from overlap which fall in contour
-    // mark source
+    for (int i = 0; i < sources->n; i++) {
+      N = index->data.U32[i];
+      source = sources->data[N];
+      overlap->n = 0;
 
-    // select the candidate PSF stars (pointers to original sources)
-    for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
-        pmSource *source = sources->data[i];
-        if (source->type != PM_SOURCE_PSFSTAR) continue;
-        psArrayAdd (stars, 200, source);
-    }
-    psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
+      // search backwards for overlapping sources
+      for (int j = N - 1; j >= 0; j--) {
+        testSource = sources->data[j];
+        if (testSource->moments->x < source->pixels->col0) continue;
+        if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue;
+        if (testSource->moments->y < source->pixels->row0) break;
+        if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) {
+          fprintf (stderr, "warning: invalid condition\n");
+          continue;
+        }
+        psArrayAdd (overlap, 100, testSource);
+      }
 
-    // get the fixed PSF fit radius
-    // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
+      // search forwards for overlapping sources
+      for (int j = N + 1; j < sources->n; j++) {
+        testSource = sources->data[j];
+        if (testSource->moments->x < source->pixels->col0) continue;
+        if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue;
+        if (testSource->moments->y < source->pixels->row0) {
+          fprintf (stderr, "warning: invalid condition\n");
+          continue;
+        }
+        if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) break;
+        psArrayAdd (overlap, 100, testSource);
+      }
 
-    // get the list pointers for the PSF_MODEL entries
-    psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
-    if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
+      // generate source contour (1/4 peak counts)
+      if (overlap->n > 0) {
+        // XXX EAM : make the 1/4 user-defined.
+        threshold = 0.25 * (source->moments->Peak - source->moments->Sky) + source->moments->Sky); 
+        psArray *contour = pmSourceContour_EAM (source->pixels, source->moments->x, source->moments->x, threshold);
 
-    psList *list = (psList *) mdi->data.list;
-    psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 
-
-    // set up an array to store the results
-    psArray *models = psArrayAlloc (list->n);
-
-    // try each model option listed in config
-    for (int i = 0; i < models->n; i++) { 
-
-        item = psListGetAndIncrement (iter);
-        modelName = item->data.V;
-
-        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
-        psFree (modelName);
-        psFree (item);
-    }
-    psFree (iter);
-    // psFree (list); XXX EAM - is list freed with iter?
-    psFree (stars);
-
-    // select the best of the models
-    // here we are using the clippedStdev on the metric as the indicator
-    try = models->data[0];
-    int bestN = 0;
-    float bestM = try->psf->dApResid;
-    for (int i = 1; i < models->n; i++) {
-        try = models->data[i];
-        float M = try->psf->dApResid;
-        if (M < bestM) {
-            bestM = M;
-            bestN = i;
-        }
+        // XXX EAM : the order is wrong: need to check j, j+1 contour points
+        psVector *xv = contour->data[0];
+        psVector *yv = contour->data[1];
+        for (int j = 0; j < xv->n; j++) {
+          for (int k = 0; k < overlap->n; k++) {
+            testSource = overlap->data[k];
+            if (fabs(yv->data.F32[j] - testSource->moments->y) > 0.5) continue;
+            if (xv->data.F32[j] > testSource->moments->x) continue;
+            if (xv->data.F32[j] < testSource->moments->x) continue;
+            // XXX EAM : mark source with flag (blended)
+          }
+        }  
     }
 
-    // keep only the selected model:
-    try = models->data[bestN];
-    pmPSF *psf = psMemCopy(try->psf);
-    psFree (models);                             // keep only the pmPSF resulting from this analysis
-
-    modelName = pmModelGetType (psf->type);
     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 

trunk/psphot/src/psphotEnsemblePSF.c

@@ -96 +96 @@
             pmSource *Mj = models->data[j];
 
+            // XXX double check that this is working!  should not break here
+            check me here;
             if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) break;
             if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;