Changeset 5718 for trunk/psphot/src/psphotEnsemblePSF.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/psphotEnsemblePSF.c (modified) (6 diffs)

trunk/psphot/src/psphotEnsemblePSF.c

@@ -13 +13 @@
     sources = psArraySort (sources, psphotSortByY);
 
+    // this should be added to the PM_SOURCE flags:
+    int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
+
     float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
     float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
@@ -23 +26 @@
     // pre-calculate all model pixels
     psArray  *models = psArrayAlloc (sources->n);
+    psVector *index  = psVectorAlloc (sources->n, PS_TYPE_U32);
+    models->n = 0;
+    index->n = 0;
+
     for (int i = 0; i < sources->n; i++) {
         pmSource *inSource = sources->data[i];
+
+        // skip non-astronomical objects (very likely defects)
+        // XXX EAM : should we try these anyway?
+        if (inSource->type == PM_SOURCE_BLEND) continue;
+        if (inSource->type == PM_SOURCE_DEFECT) continue; 
+        if (inSource->type == PM_SOURCE_SATURATED) continue;
+
         pmSource *otSource = pmSourceAlloc ();
 
@@ -72 +86 @@
         // save source in list
         otSource->modelPSF = model;
-        models->data[i] = otSource;
+        index->data.U32[models->n] = i;
+        psArrayAdd (models, 100, otSource);
     }
     psLogMsg ("psphot.emsemble", 4, "built models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n);
@@ -79 +94 @@
 
     // fill out the sparse matrix
-    psSparse *sparse = psSparseAlloc (sources->n, 100);
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *Fi = sources->data[i];
+    psSparse *sparse = psSparseAlloc (models->n, 100);
+    for (int i = 0; i < models->n; i++) {
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
         pmSource *Mi = models->data[i];
 
@@ -93 +109 @@
 
         // loop over all other stars following this one
-        for (int j = i + 1; j < sources->n; j++) {
+        for (int j = i + 1; j < models->n; j++) {
             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;
+            // skip over disjoint source images, break after last possible overlap
+            if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break;
+            if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
+            if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue;
             if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;
-            if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) continue;
-            if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
             
             // got an overlap; calculate cross-product and add to output array
@@ -115 +130 @@
 
     // adjust models, set sources and subtract
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *Fi = sources->data[i];
+    for (int i = 0; i < models->n; i++) {
+        int N = index->data.U32[i];
+        pmSource *Fi = sources->data[N];
         pmSource *Mi = models->data[i];