Index: /trunk/ppStack/src/Makefile.am =================================================================== --- /trunk/ppStack/src/Makefile.am (revision 18345) +++ /trunk/ppStack/src/Makefile.am (revision 18346) @@ -13,5 +13,6 @@ ppStackPhotometry.c \ ppStackVersion.c \ - ppStackMatch.c + ppStackMatch.c \ + ppStackSources.c noinst_HEADERS = \ Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 18345) +++ /trunk/ppStack/src/ppStack.h (revision 18346) @@ -61,5 +61,5 @@ psArray **regions, // Array of regions used in each PSF matching, returned psArray **kernels, // Array of kernels used in each PSF matching, returned - const pmReadout *sourcesRO, ///< Readout with sources + const psArray *sources, ///< Array of sources const pmPSF *psf, ///< Target PSF psRandom *rng, ///< Random number generator @@ -67,3 +67,20 @@ ); + +/// Add sources to source lists, removing duplicates and solving for magnitude differences +/// +/// Returns the source lists. +psArray *ppStackSourceListAdd(psArray *lists, ///< List to which to add, or NULL + const pmReadout *inRO, ///< Readout containing the sources + const pmConfig *config ///< Configuration + ); + +/// Combine source lists to yield a unique set of sources +/// +/// Returns the sources +psArray *ppStackSourceListCombine(psArray *lists, ///< Source lists + const pmConfig *config ///< Configuration + ); + + #endif Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 18345) +++ /trunk/ppStack/src/ppStackArguments.c (revision 18346) @@ -24,6 +24,7 @@ "\tIMAGE(STR): Image filename\n" "\tMASK(STR): Mask filename\n" - "\tWEIGHT(STR) Weight map filename\n" - "\tPSF(STR) PSF filename\n" + "\tWEIGHT(STR): Weight map filename\n" + "\tPSF(STR): PSF filename\n" + "\tSOURCES(STR): Sources filename\n" "\tWEIGHTING(F32): Relative weighting to be applied\n", program); @@ -155,4 +156,8 @@ psMetadataAddBool(arguments, PS_LIST_TAIL, "-variance", 0, "Use variance for rejection?", false); psMetadataAddBool(arguments, PS_LIST_TAIL, "-safe", 0, "Play safe with small numbers of pixels to combine?", false); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-source-radius", 0, "Source exclusion radius", NAN); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-source-iter", 0, "Source clipping iterations", 0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-source-rej", 0, "Source clipping rejection level", NAN); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-source-min", 0, "Source minimum overlap", 0); psMetadataAddBool(arguments, PS_LIST_TAIL, "-renorm", 0, "Renormalise variance maps?", false); psMetadataAddStr(arguments, PS_LIST_TAIL, "-renorm-mean", 0, "Statistic for mean in renormalisation", NULL); @@ -203,4 +208,9 @@ VALUE_ARG_RECIPE_INT("-rows", "ROWS", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-source-radius", "SOURCE.RADIUS", F32); + VALUE_ARG_RECIPE_INT("-source-iter", "SOURCE.ITER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-source-rej", "SOURCE.REJ", F32); + VALUE_ARG_RECIPE_INT("-source-min", "SOURCE.MIN", S32, 0); + VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0); VALUE_ARG_RECIPE_FLOAT("-psf-radius", "PSF.RADIUS", F32); Index: /trunk/ppStack/src/ppStackCamera.c =================================================================== --- /trunk/ppStack/src/ppStackCamera.c (revision 18345) +++ /trunk/ppStack/src/ppStackCamera.c (revision 18346) @@ -72,5 +72,6 @@ { bool haveWeights = false; // Do we have weight maps? - bool havePSFs = false; // Do we have PSFs? + bool havePSFs = false; // Do we have PSFs? + bool haveSources = (bool)psMetadataLookup(config->arguments, "PPSTACK.SOURCES"); // Have global sources? psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info @@ -99,4 +100,5 @@ psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF + psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources float weighting = psMetadataLookupF32(&mdok, input, "WEIGHTING"); // Relative weighting @@ -198,4 +200,38 @@ } + // Add the sources file + if (!haveSources) { + if (!sources || strlen(sources) == 0) { + if (havePSFs) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SOURCES %d", i); + return false; + } + } else { + if (!havePSFs) { + psWarning("SOURCES provided without PSF --- ignoring."); + } else { + psArray *sourcesFiles = psArrayAlloc(1); // Array of filenames for this FPA + sourcesFiles->data[0] = psMemIncrRefCounter(sources); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SOURCES.FILENAMES", PS_META_REPLACE, + "Filenames of SOURCES files", sourcesFiles); + psFree(sourcesFiles); + + bool status; + pmFPAfile *sourcesFile = pmFPAfileBindFromArgs(&status, imageFile, config, + "PPSTACK.INPUT.SOURCES", + "SOURCES.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)", + i, sources); + return false; + } + if (sourcesFile->type != PM_FPA_FILE_CMF) { + psError(PS_ERR_IO, true, "PPSTACK.INPUT.SOURCES is not of type CMF"); + return false; + } + } + } + } + #if 0 // Output convolved files @@ -238,7 +274,29 @@ psMetadataRemoveKey(config->arguments, "PSF.FILENAMES"); } + if (psMetadataLookup(config->arguments, "SOURCES.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "SOURCES.FILENAMES"); + } + + if (haveSources) { + // Global list of sources for use as stamps + bool status = false; // Found the file? + if (havePSFs) { + pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PPSTACK.INPUT.SOURCES", + "PPSTACK.SOURCES"); + if (!status) { + psError(PS_ERR_IO, false, "Failed to load file definition PPSTACK.INPUT.SOURCES"); + return false; + } + if (sources && sources->type != PM_FPA_FILE_CMF) { + psError(PS_ERR_IO, true, "PPSTACK.SOURCES is not of type CMF"); + return false; + } + } + } psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", i); psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.PSF", 0, "Have PSFs available?", havePSFs); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.SOURCES", 0, "Have global sources?", + haveSources); // Output image @@ -326,18 +384,4 @@ } - // Sources for use as stamps - bool status = false; // Found the file? - if (havePSFs) { - pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PPSTACK.SOURCES", "PPSTACK.SOURCES"); - if (!status) { - psError(PS_ERR_IO, false, "Failed to load file definition PPSTACK.SOURCES"); - return false; - } - if (sources && sources->type != PM_FPA_FILE_CMF) { - psError(PS_ERR_IO, true, "PPSTACK.SOURCES is not of type CMF"); - return false; - } - } - return true; } Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 18345) +++ /trunk/ppStack/src/ppStackLoop.c (revision 18346) @@ -17,5 +17,5 @@ // Files required in preparation for convolution -static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.SOURCES", NULL }; +static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", NULL }; // Files required for the convolution @@ -213,6 +213,7 @@ // Preparation iteration: Load the sources, and get a target PSF model psTrace("ppStack", 1, "Determining target PSF....\n"); - pmReadout *sources = NULL; // Readout with sources to use for PSF matching + psArray *sources = NULL; // Sources to use for PSF matching pmPSF *targetPSF = NULL; // Target PSF + bool haveSources = psMetadataLookupBool(NULL, config->arguments, "HAVE.SOURCES"); // Global sources? if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { pmFPAfileActivate(config->files, false, NULL); @@ -223,12 +224,4 @@ } - // We want to hang on to the 'sources' even when its host FPA is blown away - sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES")); - if (!sources) { - psError(PS_ERR_UNKNOWN, true, "Unable to find sources."); - psFree(view); - return false; - } - // Generate list of PSFs psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, @@ -236,4 +229,5 @@ psMetadataItem *fileItem; // Item from iteration psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope + psArray *sourceLists = NULL; // Source lists for merging sources from multiple readouts int numCols = 0, numRows = 0; // Size of image int index = 0; // Index for file @@ -247,4 +241,5 @@ psFree(view); psFree(sources); + psFree(sourceLists); psFree(fileIter); psFree(psfs); @@ -263,4 +258,5 @@ psFree(view); psFree(sources); + psFree(sourceLists); psFree(fileIter); psFree(psfs); @@ -270,4 +266,17 @@ numCols = naxis1; numRows = naxis2; + } + + if (!haveSources) { + pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources + sourceLists = ppStackSourceListAdd(sourceLists, ro, config); + if (!sourceLists) { + psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list."); + psFree(view); + psFree(sources); + psFree(fileIter); + psFree(psfs); + return false; + } } } @@ -281,4 +290,23 @@ psFree(sources); return false; + } + + 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) { + psError(PS_ERR_UNKNOWN, true, "Unable to find sources."); + psFree(view); + return false; + } + } else { + sources = ppStackSourceListCombine(sourceLists, config); + psFree(sourceLists); + if (!sources) { + psError(PS_ERR_UNKNOWN, false, "Unable to add sources to list."); + psFree(view); + psFree(psfs); + return false; + } } Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 18345) +++ /trunk/ppStack/src/ppStackMatch.c (revision 18346) @@ -16,5 +16,5 @@ bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, - const pmReadout *sourcesRO, const pmPSF *psf, psRandom *rng, const pmConfig *config) + const psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config) { assert(readout); @@ -46,5 +46,5 @@ if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) { assert(psf); - assert(sourcesRO); + assert(sources); int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order @@ -80,9 +80,4 @@ } - psArray *sources = NULL; // Sources in image - if (sourcesRO) { - sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); - } - // Generate a fake image to match to float maxMag = -INFINITY; // Maximum magnitude of sources Index: /trunk/ppStack/src/ppStackSources.c =================================================================== --- /trunk/ppStack/src/ppStackSources.c (revision 18346) +++ /trunk/ppStack/src/ppStackSources.c (revision 18346) @@ -0,0 +1,390 @@ +#include +#include +#include + +#include "ppStack.h" + +#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ + PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources +#define SOURCES_MAX_LEAF 2 // Maximum number of points on a tree leaf + +#define TESTING // Enable test output + +// Stack of sources +typedef struct { + psRegion *bounds; // Bounding box for sources + psArray *sources; // Sources within region + psTree *tree; // kd tree with sources +} ppStackSourceList; + + +// Extract coordinates from a source +static inline void coordsFromSource(float *x, float *y, // Coordinates to return + pmSource *source // Source of interest + ) +{ + if (!source) { + *x = NAN; + *y = NAN; + } else if (source->modelPSF) { + *x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + *y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + *x = source->peak->xf; + *y = source->peak->yf; + } + return; +} + +// Deallocator +static void stackSourceListFree(ppStackSourceList *list) +{ + psFree(list->bounds); + psFree(list->sources); + psFree(list->tree); +} + +// Parse the sources into vectors for each coordinate, and a bounding box +// Returns number of valid sources +static long stackSourcesParse(psRegion **bounds, // Region to update with bounding box + psVector **x, psVector **y, // Vectors to return + const psArray *sources // Input sources + ) +{ + psAssert(bounds, "Must be given a region for bounding box"); + psAssert(x && y, "Must be given vectors"); + psAssert(sources, "Must be given sources"); + + long num = sources->n; // Number of sources + *x = psVectorRecycle(*x, num, PS_TYPE_F32); + *y = psVectorRecycle(*y, num, PS_TYPE_F32); + float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources + long numValid = 0; // Number of valid sources + for (long i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + + numValid++; + float xSrc, ySrc; // Coordinates of source + coordsFromSource(&xSrc, &ySrc, source); + if (xSrc < xMin) xMin = xSrc; + if (xSrc > xMax) xMax = xSrc; + if (ySrc < yMin) yMin = ySrc; + if (ySrc > yMax) yMax = ySrc; + } + + if (*bounds) { + (*bounds)->x0 = xMin; + (*bounds)->x1 = xMax; + (*bounds)->y0 = yMin; + (*bounds)->y1 = yMax; + } else { + *bounds = psRegionAlloc(xMin, xMax, yMin, yMax); + } + + return numValid; +} + + +// Allocator for ppStackSourceList +ppStackSourceList *ppStackSourceListAlloc(psArray *sources // Sources to add + ) +{ + psAssert(sources, "Must be given sources"); + + ppStackSourceList *list = psAlloc(sizeof(ppStackSourceList)); + psMemSetDeallocator(list, (psFreeFunc)stackSourceListFree); + + list->sources = psMemIncrRefCounter(sources); + list->bounds = NULL; + + psVector *x = NULL, *y = NULL; // Source coordinates + stackSourcesParse(&list->bounds, &x, &y, sources); + list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources + psFree(x); + psFree(y); + + return list; +} + +// Extend a list of sources +bool ppStackSourceListExtend(ppStackSourceList *list, // List to extend + const psArray *newSources // Sources to add + ) +{ + psAssert(list, "Must be given a list"); + 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); + psFree(list->tree); + + psVector *x = NULL, *y = NULL; // Source coordinates + stackSourcesParse(&list->bounds, &x, &y, list->sources); + list->tree = psTreePlant(2, SOURCES_MAX_LEAF, x, y); // kd Tree with sources + psFree(x); + psFree(y); + + return list; +} + +// Compare an array of sources with a list +// Return number inside the list +static long stackSourceListCompare(ppStackSourceList **merge, // Merge target, modified + psArray *lists, // Array of source lists + int listIndex, // Index of list to compare with sources + psArray **sourcesPtr, // Sources to compare with list, modified + float radius, // Matching radius + int minOverlap, // Minimum number of sources + int iter, // Clipping iterations + float rej // Clipping rejection level + ) +{ + psAssert(merge, "Expected a source list pointer"); + psAssert(radius > 0, "Silly radius: %f", radius); + psAssert(lists, "Expected an array of source lists"); + psAssert(sourcesPtr && *sourcesPtr, "Expected a source array"); + + ppStackSourceList *list = lists->data[listIndex]; // List of interest + if (!list) { + return 0; + } + + psArray *sources = *sourcesPtr; // Sources to compare with list + psVector *x = NULL, *y = NULL; // Coordinates of sources + psRegion *bounds = NULL; // Bounding box for sources + long numSources = sources->n; // Number of sources + stackSourcesParse(&bounds, &x, &y, sources); + + if (bounds->x0 > list->bounds->x1 || bounds->x1 < list->bounds->x0 || + bounds->x0 > list->bounds->x1 || bounds->x1 < list->bounds->x0) { + // Bounds don't overlap, so the sources don't either + psFree(x); + psFree(y); + psFree(bounds); + return 0; + } + psFree(bounds); + + long numIn = 0; // Number of sources in list + long numOut = 0; // Number of sources outside list + + psVector *magDiff = psVectorAlloc(numSources, PS_TYPE_F32); // Magnitude differences + psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates of source + psArray *outside = psArrayAllocEmpty(numSources); // Sources that don't correlate + for (long i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // New source + if (!source) { + continue; + } + coords->data.F32[0] = x->data.F32[i]; + coords->data.F32[1] = y->data.F32[i]; + long match = psTreeNearestWithin(list->tree, coords, radius); // Index of match + if (match < 0) { + numOut++; + if (outside) { + psArrayAdd(outside, 1, source); + } + continue; + } + + pmSource *listSource = list->sources->data[match]; // Source in list + float listMag = listSource->psfMag; // Magnitude of source in list + float sourceMag = source->psfMag; // Magnitude of source + if (isfinite(listMag) && isfinite(sourceMag)) { + numIn++; + magDiff->data.F32[numIn] = listSource->psfMag - source->psfMag; + } + } + magDiff->n = numIn; + + psFree(coords); + psFree(x); + psFree(y); + + if (numIn > minOverlap) { + // There's sufficient overlap --- calculate the magnitude difference + psStats *stats = psStatsAlloc(iter == 0 ? PS_STAT_SAMPLE_MEAN : PS_STAT_CLIPPED_MEAN); + stats->clipSigma = rej; + stats->clipIter = iter; + if (!psVectorStats(stats, magDiff, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine mean magnitude difference"); + psFree(magDiff); + psFree(stats); + return -1; + } + float meanDiff = stats->clippedMean; // Mean magnitude difference + psFree(stats); + + // Merge our sources into the list + + if (*merge) { + // These sources act as a bridge from one list to another + // Correct magnitudes to be on the system defined by the first source list + // Correct the list to 'sources' (because 'sources' has been corrected to the merge target) + long num = list->sources->n; // Number of sources + for (long j = 0; j < num; j++) { + pmSource *source = list->sources->data[j]; + source->psfMag -= meanDiff; + } + + // Suck sources from this list into the first one we found. + ppStackSourceListExtend(*merge, list->sources); + psFree(list); + lists->data[listIndex] = NULL; + } else { + // Correct 'sources' magnitudes to match the list magnitudes + for (long j = 0; j < numOut; j++) { + pmSource *source = outside->data[j]; // New source + source->psfMag += meanDiff; + } + // Put the sources that didn't correlate into the list + ppStackSourceListExtend(list, outside); + *merge = list; + } + + // We only need to consider sources outside the list, since the lists are (mostly) disjoint + psFree(*sourcesPtr); + *sourcesPtr = outside; + } else { + psFree(outside); + } + + return numIn; +} + + +psArray *ppStackSourceListAdd(psArray *lists, const pmReadout *inRO, const pmConfig *config) +{ + PS_ASSERT_PTR_NON_NULL(inRO, 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) { + lists = psArrayAllocEmpty(1); + } + if (lists->n == 0) { + ppStackSourceList *list = ppStackSourceListAlloc(sources); + psArrayAdd(lists, 1, list); + psFree(list); // Drop reference + return lists; + } + + // Read recipe values + float radius = psMetadataLookupF32(NULL, config->arguments, "SOURCE.RADIUS"); // Exclusion radius + int minOverlap = psMetadataLookupS32(NULL, config->arguments, "SOURCE.MIN"); // Minimum number + int iter = psMetadataLookupS32(NULL, config->arguments, "SOURCE.ITER"); // Clipping iterations + float rej = psMetadataLookupF32(NULL, config->arguments, "SOURCE.REJ"); // Clipping rejection + + // Look for overlaps with existing lists + int numLists = lists->n; // Number of source lists on the stack + ppStackSourceList *merge = NULL; // Merged source list + + psMemIncrRefCounter(sources); // So we can play with the pointer + for (int i = 0; i < numLists; i++) { + long numIn = stackSourceListCompare(&merge, lists, i, &sources, radius, minOverlap, + iter, rej); // Number of overlapping sources + if (numIn == -1) { + psError(PS_ERR_UNKNOWN, false, "Unable to compare sources."); + psFree(sources); + psFree(lists); + return NULL; + } + } + psFree(sources); + + if (!merge) { + // The source list is disjoint, so throw them into a new list + ppStackSourceList *list = ppStackSourceListAlloc(sources); + psArrayAdd(lists, 1, list); + psFree(list); // Drop reference + } + + return lists; +} + + + +psArray *ppStackSourceListCombine(psArray *lists, const pmConfig *config) +{ + PS_ASSERT_ARRAY_NON_NULL(lists, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + float radius = psMetadataLookupF32(NULL, config->arguments, "SOURCE.RADIUS"); // Exclusion radius + + int numLists = 0; // Number of lists + ppStackSourceList *firstList = NULL; // Last list + for (int i = 0; i < lists->n; i++) { + ppStackSourceList *list = lists->data[i]; // List of interest + if (list) { + numLists++; + if (!firstList) { + firstList = list; + } + } + } + + if (!firstList || numLists == 0) { + psError(PS_ERR_UNKNOWN, true, "No lists found."); + return NULL; + } + +#if 0 + if (numLists == 1) { + // Trivial case + return psMemIncrRefCounter(firstList->sources); + } +#endif + + for (int i = lists->n - 1; lists->data[i] != firstList; i--) { + // There may be some small overlap (we know it's less than SOURCE.MIN), so attempt to correct + // magnitudes of these first + ppStackSourceList *lastList = lists->data[i]; // Last source list in the array + if (!lastList) { + continue; + } + psArray *sources = psMemIncrRefCounter(lastList->sources); // Sources of interest + ppStackSourceList *merge = NULL; // Merge target + for (int j = 0; j < i; j++) { + ppStackSourceList *list = lists->data[i]; // List to compare + if (!list) { + continue; + } + // Note: setting iterations = 0, so no clipping (figure not enough values to do clipping) + long numIn = stackSourceListCompare(&merge, lists, j, &sources, radius, 0, + 1, NAN); // Number of overlapping sources + if (numIn < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to compare sources."); + psFree(sources); + return NULL; + } + if (numIn == 0) { + // It doesn't match anything at all. Merge it in to the first list. + ppStackSourceListExtend(firstList, list->sources); + } + } + psFree(sources); + } + +#ifdef TESTING + FILE *srcFile = fopen("sources.txt", "w"); + for (long i = 0; i < firstList->sources->n; i++) { + float x, y; + coordsFromSource(&x, &y, firstList->sources->data[i]); + fprintf(srcFile, "%f\t%f\n", x, y); + } + fclose(srcFile); +#endif + + return psMemIncrRefCounter(firstList->sources); +}