Changeset 18382 for trunk/ppStack

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.

Location:

trunk/ppStack/src

Files:

• ppStack.h (modified) (1 diff)
• ppStackLoop.c (modified) (10 diffs)
• ppStackSources.c (modified) (12 diffs)

trunk/ppStack/src/ppStack.h

@@ -70 +70 @@
 /// Add sources to source lists, removing duplicates and solving for magnitude differences
 ///
-/// Returns the source lists.
+/// Corrects the sources to have consistent magnitudes.  Returns the source lists.
 psArray *ppStackSourceListAdd(psArray *lists, ///< List to which to add, or NULL
-                              const pmReadout *inRO, ///< Readout containing the sources
+                              psArray *sources, ///< Sources to add
                               const pmConfig *config ///< Configuration
     );
 
-/// Combine source lists to yield a unique set of sources
+/// Combine source lists to yield a set of unique sources.
 ///
-/// Returns the sources
+/// Corrects the sources to have consistent magnitudes where possible.  Returns the sources
 psArray *ppStackSourceListCombine(psArray *lists, ///< Source lists
                                   const pmConfig *config ///< Configuration

trunk/ppStack/src/ppStackLoop.c

@@ -213 +213 @@
     // Preparation iteration: Load the sources, and get a target PSF model
     psTrace("ppStack", 1, "Determining target PSF....\n");
-    psArray *sources = NULL;            // Sources to use for PSF matching
+    psArray *globalSources = NULL;      // Global list of sources for matching (haveSources = TRUE)
+    psArray *indSources = psArrayAlloc(num); // Individual lists of sources for matching (haveSources = FALSE)
     pmPSF *targetPSF = NULL;            // Target PSF
     bool haveSources = psMetadataLookupBool(NULL, config->arguments, "HAVE.SOURCES"); // Global sources?
@@ -240 +241 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
                 psFree(view);
-                psFree(sources);
+                psFree(globalSources);
+                psFree(indSources);
                 psFree(sourceLists);
                 psFree(fileIter);
@@ -257 +259 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
                 psFree(view);
-                psFree(sources);
+                psFree(globalSources);
+                psFree(indSources);
                 psFree(sourceLists);
                 psFree(fileIter);
@@ -270 +273 @@
             if (!haveSources) {
                 pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
-                sourceLists = ppStackSourceListAdd(sourceLists, ro, config);
+                psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
+                if (!sources) {
+                    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
+                    return NULL;
+                }
+                sources->data[index] = psMemIncrRefCounter(sources);
+                sourceLists = ppStackSourceListAdd(sourceLists, sources, config);
                 if (!sourceLists) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list.");
                     psFree(view);
-                    psFree(sources);
+                    psFree(globalSources);
+                    psFree(indSources);
                     psFree(fileIter);
                     psFree(psfs);
@@ -288 +298 @@
         if (!targetPSF) {
             psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
-            psFree(sources);
+            psFree(globalSources);
+            psFree(indSources);
             return false;
         }
@@ -294 +305 @@
         if (haveSources) {
             // We want to hang on to the 'sources' even when its host FPA is blown away
-            sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.INPUT.SOURCES"));
-            if (!sources) {
+            globalSources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view,
+                                                                     "PPSTACK.INPUT.SOURCES"));
+            if (!globalSources) {
                 psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
                 psFree(view);
@@ -301 +313 @@
             }
         } else {
-            sources = ppStackSourceListCombine(sourceLists, config);
+            globalSources = ppStackSourceListCombine(sourceLists, config);
             psFree(sourceLists);
-            if (!sources) {
+            if (!globalSources) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list.");
                 psFree(view);
@@ -343 +355 @@
         pmFPAview *view = filesIterateDown(config);
         if (!view) {
-            psFree(sources);
+            psFree(globalSources);
+            psFree(indSources);
             psFree(targetPSF);
             psFree(rng);
@@ -354 +367 @@
         // Background subtraction, scaling and normalisation is performed automatically by the image matching
         psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
+        psArray *sources = haveSources ? globalSources : indSources->data[i]; // Sources for matching
         if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, rng, config)) {
             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
@@ -374 +388 @@
         numGood++;
     }
-    psFree(sources);
+    psFree(globalSources);
+    psFree(indSources);
     psFree(targetPSF);
     psFree(rng);

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);
             }
         }