Index: trunk/ppStack/src/ppStackSources.c =================================================================== --- trunk/ppStack/src/ppStackSources.c (revision 20944) +++ trunk/ppStack/src/ppStackSources.c (revision 20995) @@ -1,3 +1,5 @@ #include +#include +#include #include #include @@ -5,496 +7,163 @@ #include "ppStack.h" -#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \ - PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ - PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources -#define SOURCE_FAINTEST 50.0 // Faintest magnitude to consider -#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 - psVector *indices; // Indices from tree to sources - psArray *duplicates; // Duplicate sources (from overlaps) - psTree *tree; // kd tree with sources -} ppStackSourceList; +//#define TESTING // Enable debugging output -// Extract coordinates from a source -static inline void coordsFromSource(float *x, float *y, // Coordinates to return - pmSource *source // Source of interest - ) +bool ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config) { - 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; -} + PS_ASSERT_ARRAY_NON_NULL(sourceLists, false); + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(config, false); -// Deallocator -static void stackSourceListFree(ppStackSourceList *list) -{ - psFree(list->bounds); - psFree(list->sources); - psFree(list->indices); - psFree(list->duplicates); - 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, // Coordinate vectors to return - psVector **indices, // Source indices 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(indices, "Must be given indices"); - psAssert(sources, "Must be given sources"); - - long numSources = sources->n; // Number of sources - *x = psVectorRecycle(*x, numSources, PS_TYPE_F32); - *y = psVectorRecycle(*y, numSources, PS_TYPE_F32); - *indices = psVectorRecycle(*indices, numSources, PS_TYPE_U32); - float xMin = INFINITY, xMax = -INFINITY, yMin = INFINITY, yMax = -INFINITY; // Bounds of sources - long num = 0; // Number of valid sources - for (long i = 0; i < numSources; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) || - source->psfMag > SOURCE_FAINTEST) { - continue; - } - - 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; - - (*x)->data.F32[num] = xSrc; - (*y)->data.F32[num] = ySrc; - (*indices)->data.U32[num] = i; - num++; - } - (*x)->n = num; - (*y)->n = num; - (*indices)->n = num; - - if (*bounds) { - (*bounds)->x0 = xMin; - (*bounds)->x1 = xMax; - (*bounds)->y0 = yMin; - (*bounds)->y1 = yMax; - } else { - *bounds = psRegionAlloc(xMin, xMax, yMin, yMax); - } - - psTrace("ppStack.sources", 8, "%ld sources: bounds are [%.2f:%.2f,%.2f:%.2f]\n", - num, xMin, xMax, yMin, yMax); - - return num; -} - - -// 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); - - // Copy the sources onto a new array, so there's no confusion between the inputs and outputs - long num = sources->n; // Number of sources - list->sources = psArrayAlloc(num); - for (long i = 0; i < num; i++) { - list->sources->data[i] = psMemIncrRefCounter(sources->data[i]); - } - - list->bounds = NULL; - list->indices = NULL; - list->duplicates = psArrayAllocEmpty(0); - - psVector *x = NULL, *y = NULL; // Source coordinates - stackSourcesParse(&list->bounds, &x, &y, &list->indices, 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 - const psArray *newDups // Duplicate sources to add - ) -{ - psAssert(list, "Must be given a list"); - psAssert(newSources, "Must be given sources"); - +#ifdef TESTING { - 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]); + // Deliberately induce a major transparency difference + psArray *sources = sourceLists->data[1]; // Sources to correct + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + source->psfMag += 1.0; } } - - 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 - psVector *indices = NULL; // Indices for sources - psRegion *bounds = NULL; // Bounding box for sources - long numSources = stackSourcesParse(&bounds, &x, &y, &indices, sources); // Number of (good) sources - - if (bounds->x0 > list->bounds->x1 || bounds->x1 < list->bounds->x0 || - bounds->y0 > list->bounds->y1 || bounds->y1 < list->bounds->y0) { - // Bounds don't overlap, so the sources don't either - psTrace("ppStack.sources", 7, "Bounds don't overlap\n"); - psFree(x); - psFree(y); - psFree(indices); - psFree(bounds); - return 0; - } - psFree(bounds); - - long numIn = 0; // Number of sources in list - long numOut = 0; // Number of sources outside list - -#ifdef TESTING - static int fileNum = 0; // Number of file - psString filename = NULL; // Name of file - psStringAppend(&filename, "source_match_%d.txt", fileNum++); - FILE *file = fopen(filename, "w"); // File to write - psFree(filename); - fprintf(file, "# x1 y1 x2 y2 mag1 mag1err mag2 mag2err flag1 flag2\n"); #endif - - 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 - psArray *inside = psArrayAllocEmpty(numSources); // Sources that do correlate - for (long i = 0; i < numSources; i++) { - 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++; - psArrayAdd(outside, 1, sources->data[indices->data.U32[i]]); - continue; - } - - pmSource *listSource = list->sources->data[list->indices->data.U32[match]]; // Source in list - pmSource *source = sources->data[indices->data.U32[i]]; // Source of interest - float listMag = listSource->psfMag; // Magnitude of source in list - float sourceMag = source->psfMag; // Magnitude of source - if (!(listSource->mode & SOURCE_MASK) && !(source->mode & SOURCE_MASK) && - isfinite(listMag) && isfinite(sourceMag)) { -#ifdef TESTING - float xList, yList; // Coordinates from source on list - coordsFromSource(&xList, &yList, listSource); - fprintf(file, "%f %f %f %f %f %f %f %f %.4x %.4x\n", x->data.F32[i], y->data.F32[i], xList, yList, - listSource->psfMag, listSource->errMag, source->psfMag, source->errMag, - listSource->mode, source->mode); -#endif - magDiff->data.F32[numIn] = listSource->psfMag - source->psfMag; - psArrayAdd(inside, 1, source); - numIn++; - } - } - magDiff->n = numIn; - -#ifdef TESTING - fclose(file); -#endif - - psFree(coords); - psFree(x); - psFree(y); - psFree(indices); - - psTrace("ppStack.sources", 7, "%ld sources (vs %d) overlap list %d\n", numIn, minOverlap, listIndex); - - if (numIn > minOverlap) { - // There's sufficient overlap --- calculate the magnitude difference - psStats *stats = psStatsAlloc(iter == 0 ? (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV) : - (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV)); - 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 = iter == 0 ? stats->sampleMean : stats->clippedMean; // Mean magnitude difference - float stdev = iter == 0 ? stats->sampleStdev : stats->clippedStdev; // Standard deviation - long numStats = iter == 0 ? numIn : stats->clippedNvalues; // Number used for statistics - psLogMsg("ppStack.sources", PS_LOG_INFO, "Magnitude difference from %ld overlaps is %f +/- %f\n", - numStats, meanDiff, stdev); - 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; - } - 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: Merging lists."); - - ppStackSourceListExtend(*merge, list->sources, NULL); - 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; - } - 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, "Overlap: Extending list"); - ppStackSourceListExtend(list, outside, inside); - *merge = list; - } - - // We only need to consider sources outside the list, since the lists are (mostly) disjoint - psFree(*sourcesPtr); - *sourcesPtr = outside; - } else { - psFree(outside); - } - psFree(inside); - psFree(magDiff); - - return numIn; -} - - -psArray *ppStackSourceListAdd(psArray *lists, psArray *sources, const pmConfig *config) -{ - PS_ASSERT_PTR_NON_NULL(sources, NULL); - PS_ASSERT_PTR_NON_NULL(config, 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 - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe - psAssert(recipe, "We've thrown an error on this before."); - float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Exclusion radius - int minOverlap = psMetadataLookupS32(NULL, recipe, "SOURCE.MIN"); // Minimum number - int iter = psMetadataLookupS32(NULL, recipe, "SOURCE.ITER"); // Clipping iterations - float rej = psMetadataLookupF32(NULL, recipe, "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 - psTrace("ppStack.sources", 4, "Disjoint source: Creating 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); psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe psAssert(recipe, "We've thrown an error on this before."); - float radius = psMetadataLookupF32(NULL, recipe, "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; - } - } + float radius = psMetadataLookupF32(NULL, recipe, "ZP.RADIUS"); // Radius (pixels) for matching sources + int iter = psMetadataLookupS32(NULL, recipe, "ZP.ITER"); // Maximum iterations + float tol = psMetadataLookupF32(NULL, recipe, "ZP.TOL"); // Tolerance for zero point iterations + int transIter = psMetadataLookupS32(NULL, recipe, "ZP.TRANS.ITER"); // Iterations for transparency + float transRej = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.REJ");// Rejection threshold for transparency + float transThresh = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.THRESH"); // Threshold for transparency + float starRej = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ"); // Rejection threshold for stars + float starLimit = psMetadataLookupF32(NULL, recipe, "ZP.STAR.LIMIT"); // Limit on star rejection fraction + float starSys = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS"); // Estimated systematic error + + psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms + if (!airmassZP) { + psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe."); + return false; } - if (!firstList || numLists == 0) { - psError(PS_ERR_UNKNOWN, true, "No lists found."); - return NULL; + int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n); + + psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Zero points for each image + const char *filter = NULL; // Filter name + float airmassTerm = NAN; // Airmass term + float sumExpTime = 0.0; // Sum of the exposure time + for (int i = 0; i < num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest + pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest + + float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time + float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass + const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name + if (!isfinite(exptime) || exptime == 0 || !isfinite(airmass) || airmass == 0 || + !expFilter || strlen(expFilter) == 0) { + psError(PS_ERR_UNEXPECTED_NULL, false, + "Unable to find exposure time (%f), airmass (%f) or filter (%s)", + exptime, airmass, expFilter); + psFree(zp); + return false; + } + + if (!filter) { + filter = expFilter; + airmassTerm = psMetadataLookupF32(NULL, airmassZP, filter); + if (!isfinite(airmassTerm)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter); + psFree(zp); + return false; + } + } else if (strcmp(filter, expFilter) != 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter); + psFree(zp); + return false; + } + + zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime); + sumExpTime += exptime; } -#if 0 - if (numLists == 1) { - // Trivial case - return psMemIncrRefCounter(firstList->sources); + psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches + if (!matches) { + psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); + psFree(zp); + return false; + } + psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, + transThresh, starRej, starSys); // Transparencies for each image + +#ifdef TESTING + { + // Dump the corrected magnitudes + FILE *outMatches = fopen("source_match.dat", "w"); // Output matches + psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star + for (int i = 0; i < matches->n; i++) { + pmSourceMatch *match = matches->data[i]; // Match of interest + psVectorInit(mag, NAN); + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_MASK_DATA[j]) { + continue; + } + int index = match->image->data.U32[j]; // Image index + mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index]; + } + for (int j = 0; j < num; j++) { + fprintf(outMatches, "%f ", mag->data.F32[j]); + } + fprintf(outMatches, "\n"); + } + psFree(mag); + fclose(outMatches); } #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) { + psFree(matches); + if (!trans) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); + return false; + } + + // M = m + c0 + c1 * airmass + 2.5log(t) + transparency + // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 + // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 + // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0) + // We don't need to know the magnitude zero point for the filter, since it cancels out + for (int i = 0; i < num; i++) { + psArray *sources = sourceLists->data[i]; // Sources of interest + float magCorr = airmassTerm + 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i]; + psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr); + + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source) { 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. - psTrace("ppStack.sources", 4, "Disjoint list: Merging into first"); - ppStackSourceListExtend(firstList, list->sources, NULL); - } + source->psfMag += magCorr; } - psFree(sources); } - - return psMemIncrRefCounter(firstList->sources); -} + psFree(trans); - -void ppStackSourcesPrint(const psArray *sources) -{ - static int num = 0; // Number of source array - psString filename = NULL; // Name of file - psStringAppend(&filename, "sources_%d.reg", num++); - - const char *goodColour = "green"; // Colour for good sources - const char *badColour = "red"; // Colour for bad sources - float radius = 5; // Radius for circle - - FILE *file = fopen(filename, "w"); // File to which to write - psFree(filename); - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; +#ifdef TESTING + // Double check: all transparencies should be zero + { + psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches + if (!matches) { + psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); + psFree(zp); + return false; } - - float x, y; // Source coordinates - coordsFromSource(&x, &y, source); - - bool bad = (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) || - source->psfMag > SOURCE_FAINTEST; // Is this a bad source? - - fprintf(file, "image; circle(%f,%f,%f) # color = %s\n", x + 1.0, y + 1.0, radius, - bad ? badColour : goodColour); + psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, + transThresh, starRej, starSys); + psFree(trans); } - fclose(file); - - return; +#endif + return true; }