Changeset 5718 for trunk/psphot/src/psphotBasicDeblend.c

Timestamp:

Dec 6, 2005, 9:55:31 PM (21 years ago)

Author:

eugene

Message:

finshed deblender, dropped local copies of psModule code

File:

• trunk/psphot/src/psphotBasicDeblend.c (modified) (1 diff)

trunk/psphot/src/psphotBasicDeblend.c

@@ -1 +1 @@
 # include "psphot.h"
 
-// try PSF models and select best option
+bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky) { 
 
-pmPSF *psphotBasicDeblend (psMetadata *config, psArray *sources, psStats *skystats) 
-{ 
+    int N;
+    bool status;
+    float threshold;
+    pmSource *source, *testSource;
 
-    // source analysis is done in S/N order (brightest first)
+    FILE *f = fopen ("deblend.dat", "w");
+    if (f == NULL) psAbort ("psphot", "can't open deblend.dat output file");
+    int Nblend = 0;
+
+    psTimerStart ("psphot");
+
+    // this should be added to the PM_SOURCE flags:
+    int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
+
+    float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION");
+    float NSIGMA   = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA");
+    float minThreshold = NSIGMA*sky->sampleStdev;
+    fprintf (stderr, "min threshold: %f\n", minThreshold);
+
+    // we need sources spatially-sorted to find overlaps
     sources = psArraySort (sources, psphotSortByY);
 
+    // source analysis is done in peak order (brightest first)
+    // we use an index for this so the spatial sorting is kept
     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;
+        source = sources->data[i];
+        SN->data.F32[i] = source->peak->counts;
     }
     psVector *index = psVectorSortIndex (NULL, SN);
-    
+    // this results in an index of increasing SN
+
+    // temporary array for overlapping objects we find
     psArray *overlap = psArrayAlloc (100);
 
-    for (int i = 0; i < sources->n; i++) {
-      N = index->data.U32[i];
-      source = sources->data[N];
-      overlap->n = 0;
+    // XXX make sure this results in decreasing, not increasing, SN
+    for (int i = sources->n - 1; i >= 0; i--) {
+        N = index->data.U32[i];
+        source = sources->data[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;
+        if (source->type == PM_SOURCE_BLEND) continue;
+
+        overlap->n = 0;
+
+        // search backwards for overlapping sources
+        for (int j = N - 1; j >= 0; j--) {
+            testSource = sources->data[j];
+            if (testSource->peak->x <  source->pixels->col0) continue;
+            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+            if (testSource->peak->y <  source->pixels->row0) break;
+            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) {
+                fprintf (stderr, "warning: invalid condition\n");
+                continue;
+            }
+            psArrayAdd (overlap, 100, testSource);
         }
-        psArrayAdd (overlap, 100, testSource);
-      }
 
-      // 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;
+        // search forwards for overlapping sources
+        for (int j = N + 1; j < sources->n; j++) {
+            testSource = sources->data[j];
+            if (testSource->peak->x <  source->pixels->col0) continue;
+            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+            if (testSource->peak->y <  source->pixels->row0) {
+                fprintf (stderr, "warning: invalid condition\n");
+                continue;
+            }
+            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break;
+            psArrayAdd (overlap, 100, testSource);
         }
-        if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) break;
-        psArrayAdd (overlap, 100, testSource);
-      }
 
-      // 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);
+        // generate source contour (1/4 peak counts)
+        if (overlap->n > 0) {
+            // XXX EAM : make the 1/4 user-defined.
+            // XXX keep threshold from dropping too low (N*sky.sigma)
+            threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; 
+            threshold = PS_MAX (threshold, minThreshold);
 
-        // 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)
-          }
-        }  
+            // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now
+            psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold);
+            if (contour == NULL) {
+                fprintf (stderr, "below threshold? invalid peak?\n");
+                continue;
+            }
+
+            // the source contour consists of two vectors, xv and yv.  the contour is 
+            // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]).  both
+            // coordinate pairs have the same yv[] value, with xv[i] corresponding to
+            // the left boundary and xv[i+1] corresponding to the right boundary
+
+            psVector *xv = contour->data[0];
+            psVector *yv = contour->data[1];
+            for (int k = 0; k < overlap->n; k++) {
+                testSource = overlap->data[k];
+                if (testSource->peak->counts > source->peak->counts) continue;
+                for (int j = 0; j < xv->n; j+=2) {
+                    if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;
+                    if (xv->data.F32[j+0] > testSource->peak->x) break;
+                    if (xv->data.F32[j+1] < testSource->peak->x) break;
+
+                    int xp0 = source->moments->x - source->pixels->col0;
+                    int xp1 = source->peak->x - source->pixels->col0;
+                    int xp2 = testSource->moments->x - source->pixels->col0;
+                    int xp3 = testSource->peak->x - source->pixels->col0;
+
+                    int yp0 = source->moments->y - source->pixels->row0;
+                    int yp1 = source->peak->y - source->pixels->row0;
+                    int yp2 = testSource->moments->y - testSource->pixels->row0;
+                    int yp3 = testSource->peak->y - testSource->pixels->row0;
+
+                    fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n", 
+                             source->peak->x, source->peak->y, 
+                             source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1], 
+                             testSource->peak->x, testSource->peak->y, 
+                             testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
+                             );
+
+                    testSource->type = PM_SOURCE_BLEND;
+                    Nblend ++;
+                    j = xv->n;
+                }
+            }  
+        }
     }
+    fclose (f);
+    psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f)\n", Nblend, psTimerMark ("psphot"));
+    return true;
+}
 
-    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
-    return (psf);
-}
+