Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 18356) +++ trunk/ppStack/src/ppStack.h (revision 18382) @@ -70,13 +70,13 @@ /// 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 Index: trunk/ppStack/src/ppStackLoop.c =================================================================== --- trunk/ppStack/src/ppStackLoop.c (revision 18356) +++ trunk/ppStack/src/ppStackLoop.c (revision 18382) @@ -213,5 +213,6 @@ // 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,5 +241,6 @@ psError(PS_ERR_UNKNOWN, false, "Unable to find PSF."); psFree(view); - psFree(sources); + psFree(globalSources); + psFree(indSources); psFree(sourceLists); psFree(fileIter); @@ -257,5 +259,6 @@ 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,9 +273,16 @@ 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,5 +298,6 @@ if (!targetPSF) { psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); - psFree(sources); + psFree(globalSources); + psFree(indSources); return false; } @@ -294,6 +305,7 @@ 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,7 +313,7 @@ } } 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,5 +355,6 @@ pmFPAview *view = filesIterateDown(config); if (!view) { - psFree(sources); + psFree(globalSources); + psFree(indSources); psFree(targetPSF); psFree(rng); @@ -354,4 +367,5 @@ // 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, ®ions, &kernels, sources, targetPSF, rng, config)) { psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); @@ -374,5 +388,6 @@ numGood++; } - psFree(sources); + psFree(globalSources); + psFree(indSources); psFree(targetPSF); psFree(rng); Index: trunk/ppStack/src/ppStackSources.c =================================================================== --- trunk/ppStack/src/ppStackSources.c (revision 18356) +++ trunk/ppStack/src/ppStackSources.c (revision 18382) @@ -17,4 +17,5 @@ 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,4 +46,5 @@ psFree(list->sources); psFree(list->indices); + psFree(list->duplicates); psFree(list->tree); } @@ -130,5 +132,6 @@ // 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,19 +139,32 @@ 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,4 +217,5 @@ 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,7 +224,5 @@ 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,4 +235,5 @@ numIn++; magDiff->data.F32[numIn] = listSource->psfMag - source->psfMag; + psArrayAdd(inside, 1, source); } } @@ -256,9 +272,13 @@ 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,7 +289,11 @@ 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,4 +305,5 @@ psFree(outside); } + psFree(inside); return numIn; @@ -286,14 +311,8 @@ -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,5 +418,5 @@ // 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); } }