Changeset 18382 for trunk/ppStack/src/ppStackSources.c

Timestamp:

Jun 30, 2008, 12:09:14 PM (18 years ago)

Author:

Paul Price

Message:

Figured that we don't want a set of unique sources, since that includes all the cosmic rays from all the inputs. Rather, we want for each image the list of sources in that image, but corrected to be on a common magnitude system (matching zero points across the images) and use that for our target PSF image.

File:

• trunk/ppStack/src/ppStackSources.c (modified) (12 diffs)

trunk/ppStack/src/ppStackSources.c

@@ -17 +17 @@
     psArray *sources;                   // Sources within region
     psVector *indices;                  // Indices from tree to sources
+    psArray *duplicates;                // Duplicate sources (from overlaps)
     psTree *tree;                       // kd tree with sources
 } ppStackSourceList;
@@ -45 +46 @@
     psFree(list->sources);
     psFree(list->indices);
+    psFree(list->duplicates);
     psFree(list->tree);
 }
@@ -130 +132 @@
 // Extend a list of sources
 bool ppStackSourceListExtend(ppStackSourceList *list, // List to extend
-                             const psArray *newSources // Sources to add
+                             const psArray *newSources, // Sources to add
+                             const psArray *newDups // Duplicate sources to add
     )
 {
@@ -136 +139 @@
     psAssert(newSources, "Must be given sources");
 
-    psArray *oldSources = list->sources;// Old list of sources
-    long numOld = oldSources->n;        // Number of old sources
-    long numNew = newSources->n;        // Number of new sources
-
-    list->sources = oldSources = psArrayRealloc(oldSources, numOld + numNew);
-    for (long i = 0; i < numNew; i++) {
-        psArrayAdd(oldSources, 1, newSources->data[i]);
-    }
-    psFree(list->tree);
-
-    psVector *x = NULL, *y = NULL;      // Source coordinates
-    stackSourcesParse(&list->bounds, &x, &y, &list->indices, list->sources);
-    list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources
-    psFree(x);
-    psFree(y);
+    {
+        psArray *oldSources = list->sources;// Old list of sources
+        long numOld = oldSources->n;        // Number of old sources
+        long numNew = newSources->n;        // Number of new sources
+
+        list->sources = oldSources = psArrayRealloc(oldSources, numOld + numNew);
+        for (long i = 0; i < numNew; i++) {
+            psArrayAdd(oldSources, 1, newSources->data[i]);
+        }
+        psFree(list->tree);
+
+        psVector *x = NULL, *y = NULL;      // Source coordinates
+        stackSourcesParse(&list->bounds, &x, &y, &list->indices, list->sources);
+        list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources
+        psFree(x);
+        psFree(y);
+    }
+
+    if (newDups) {
+        psArray *oldDups = list->duplicates; // Old list of duplicates
+        long numOld = oldDups ? oldDups->n : 0; // Number of old duplicates
+        long numNew = newDups->n;       // Number of new duplicates
+
+        list->duplicates = oldDups = psArrayRealloc(oldDups, numOld + numNew);
+        for (long i = 0; i < numNew; i++) {
+            psArrayAdd(oldDups, 1, newDups->data[i]);
+        }
+    }
 
     return list;
@@ -201 +217 @@
     psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates of source
     psArray *outside = psArrayAllocEmpty(numSources); // Sources that don't correlate
+    psArray *inside = psArrayAllocEmpty(numSources); // Sources that do correlate
     for (long i = 0; i < numSources; i++) {
         coords->data.F32[0] = x->data.F32[i];
@@ -207 +224 @@
         if (match < 0) {
             numOut++;
-            if (outside) {
-                psArrayAdd(outside, 1, sources->data[indices->data.U32[i]]);
-            }
+            psArrayAdd(outside, 1, sources->data[indices->data.U32[i]]);
             continue;
         }
@@ -220 +235 @@
             numIn++;
             magDiff->data.F32[numIn] = listSource->psfMag - source->psfMag;
+            psArrayAdd(inside, 1, source);
         }
     }
@@ -256 +272 @@
                 source->psfMag -= meanDiff;
             }
+            for (long j = 0; j < list->duplicates->n; j++) {
+                pmSource *source = list->duplicates->data[j];
+                source->psfMag -= meanDiff;
+            }
 
             // Suck sources from this list into the first one we found.
-            psTrace("ppStack.sources", 4, "Bridge found: Merging lists.");
-
-            ppStackSourceListExtend(*merge, list->sources);
+            psTrace("ppStack.sources", 4, "Bridge: Merging lists.");
+
+            ppStackSourceListExtend(*merge, list->sources, NULL);
             psFree(list);
             lists->data[listIndex] = NULL;
@@ -269 +289 @@
                 source->psfMag += meanDiff;
             }
+            for (long j = 0; j < numIn; j++) {
+                pmSource *source = inside->data[j]; // New source
+                source->psfMag += meanDiff;
+            }
             // Put the sources that didn't correlate into the list
-            psTrace("ppStack.sources", 4, "Extending list");
-            ppStackSourceListExtend(list, outside);
+            psTrace("ppStack.sources", 4, "Overlap: Extending list");
+            ppStackSourceListExtend(list, outside, inside);
             *merge = list;
         }
@@ -281 +305 @@
         psFree(outside);
     }
+    psFree(inside);
 
     return numIn;
@@ -286 +311 @@
 
 
-psArray *ppStackSourceListAdd(psArray *lists, const pmReadout *inRO, const pmConfig *config)
-{
-    PS_ASSERT_PTR_NON_NULL(inRO, NULL);
+psArray *ppStackSourceListAdd(psArray *lists, psArray *sources, const pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(sources, NULL);
     PS_ASSERT_PTR_NON_NULL(config, NULL);
-
-    psArray *sources = psMetadataLookupPtr(NULL, inRO->analysis, "PSPHOT.SOURCES");
-    if (!sources) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
-        return NULL;
-    }
 
     if (!lists) {
@@ -399 +418 @@
                 // It doesn't match anything at all.  Merge it in to the first list.
                 psTrace("ppStack.sources", 4, "Disjoint list: Merging into first");
-                ppStackSourceListExtend(firstList, list->sources);
+                ppStackSourceListExtend(firstList, list->sources, NULL);
             }
         }