Index: /trunk/psModules/src/objects/Makefile.am =================================================================== --- /trunk/psModules/src/objects/Makefile.am (revision 20945) +++ /trunk/psModules/src/objects/Makefile.am (revision 20946) @@ -4,83 +4,85 @@ libpsmodulesobjects_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmodulesobjects_la_SOURCES = \ - pmDetections.c \ - pmSpan.c \ - pmFootprint.c \ - pmFootprintArrayGrow.c \ - pmFootprintArraysMerge.c \ - pmFootprintAssignPeaks.c \ - pmFootprintFind.c \ - pmFootprintFindAtPoint.c \ - pmFootprintIDs.c \ - pmPeaks.c \ - pmMoments.c \ - pmModel.c \ - pmModelClass.c \ - pmModelUtils.c \ - pmSource.c \ - pmSourceMoments.c \ - pmSourceExtendedPars.c \ - pmSourceUtils.c \ - pmSourceSky.c \ - pmSourceContour.c \ - pmSourceFitModel.c \ - pmSourceFitSet.c \ - pmSourcePhotometry.c \ - pmSourceIO.c \ - pmSourceIO_RAW.c \ - pmSourceIO_OBJ.c \ - pmSourceIO_SX.c \ - pmSourceIO_CMP.c \ - pmSourceIO_SMPDATA.c \ - pmSourceIO_PS1_DEV_0.c \ - pmSourceIO_PS1_DEV_1.c \ - pmSourceIO_PS1_CAL_0.c \ - pmSourcePlots.c \ - pmSourcePlotPSFModel.c \ - pmSourcePlotMoments.c \ - pmSourcePlotApResid.c \ - pmSourceVisual.c \ - pmResiduals.c \ - pmPSF.c \ - pmPSF_IO.c \ - pmPSFtry.c \ - pmTrend2D.c \ - pmGrowthCurveGenerate.c \ - pmGrowthCurve.c + pmDetections.c \ + pmSpan.c \ + pmFootprint.c \ + pmFootprintArrayGrow.c \ + pmFootprintArraysMerge.c \ + pmFootprintAssignPeaks.c \ + pmFootprintFind.c \ + pmFootprintFindAtPoint.c \ + pmFootprintIDs.c \ + pmPeaks.c \ + pmMoments.c \ + pmModel.c \ + pmModelClass.c \ + pmModelUtils.c \ + pmSource.c \ + pmSourceMoments.c \ + pmSourceExtendedPars.c \ + pmSourceUtils.c \ + pmSourceSky.c \ + pmSourceContour.c \ + pmSourceFitModel.c \ + pmSourceFitSet.c \ + pmSourcePhotometry.c \ + pmSourceIO.c \ + pmSourceIO_RAW.c \ + pmSourceIO_OBJ.c \ + pmSourceIO_SX.c \ + pmSourceIO_CMP.c \ + pmSourceIO_SMPDATA.c \ + pmSourceIO_PS1_DEV_0.c \ + pmSourceIO_PS1_DEV_1.c \ + pmSourceIO_PS1_CAL_0.c \ + pmSourcePlots.c \ + pmSourcePlotPSFModel.c \ + pmSourcePlotMoments.c \ + pmSourcePlotApResid.c \ + pmSourceVisual.c \ + pmResiduals.c \ + pmPSF.c \ + pmPSF_IO.c \ + pmPSFtry.c \ + pmTrend2D.c \ + pmGrowthCurveGenerate.c \ + pmGrowthCurve.c \ + pmSourceMatch.c EXTRA_DIST = \ - models/pmModel_GAUSS.c \ - models/pmModel_PGAUSS.c \ - models/pmModel_QGAUSS.c \ - models/pmModel_SGAUSS.c \ - models/pmModel_RGAUSS.c \ - models/pmModel_SERSIC.c + models/pmModel_GAUSS.c \ + models/pmModel_PGAUSS.c \ + models/pmModel_QGAUSS.c \ + models/pmModel_SGAUSS.c \ + models/pmModel_RGAUSS.c \ + models/pmModel_SERSIC.c pkginclude_HEADERS = \ - pmDetections.h \ - pmSpan.h \ - pmFootprint.h \ - pmPeaks.h \ - pmMoments.h \ - pmModel.h \ - pmModelClass.h \ - pmModelUtils.h \ - pmSource.h \ - pmSourceExtendedPars.h \ - pmSourceUtils.h \ - pmSourceSky.h \ - pmSourceContour.h \ - pmSourceFitModel.h \ - pmSourceFitSet.h \ - pmSourcePhotometry.h \ - pmSourceIO.h \ - pmSourcePlots.h \ - pmSourceVisual.h \ - pmResiduals.h \ - pmPSF.h \ - pmPSF_IO.h \ - pmPSFtry.h \ - pmTrend2D.h \ - pmGrowthCurve.h + pmDetections.h \ + pmSpan.h \ + pmFootprint.h \ + pmPeaks.h \ + pmMoments.h \ + pmModel.h \ + pmModelClass.h \ + pmModelUtils.h \ + pmSource.h \ + pmSourceExtendedPars.h \ + pmSourceUtils.h \ + pmSourceSky.h \ + pmSourceContour.h \ + pmSourceFitModel.h \ + pmSourceFitSet.h \ + pmSourcePhotometry.h \ + pmSourceIO.h \ + pmSourcePlots.h \ + pmSourceVisual.h \ + pmResiduals.h \ + pmPSF.h \ + pmPSF_IO.h \ + pmPSFtry.h \ + pmTrend2D.h \ + pmGrowthCurve.h \ + pmSourceMatch.h CLEANFILES = *~ Index: /trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.c (revision 20946) +++ /trunk/psModules/src/objects/pmSourceMatch.c (revision 20946) @@ -0,0 +1,292 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmSource.h" + +#include "pmSourceMatch.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 + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// 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; +} + +// Parse the sources into vectors for each coordinate, and a bounding box +// Returns number of valid sources +static long sourcesParse(psRegion **bounds, // Region to update with bounding box + psVector **x, psVector **y, // Coordinate vectors to return + psVector **mag, psVector **magErr, // Magnitude and error vectors to return + const psArray *sources // Input sources + ) +{ + psAssert(bounds, "Must be given a region for bounding box"); + psAssert(x && y, "Must be given position vectors"); + psAssert(mag && magErr, "Must be given magnitude vectors"); + 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); + 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; + (*mag)->data.F32[num] = source->psfMag; + (*magErr)->data.F32[num] = source->errMag; + num++; + } + (*x)->n = num; + (*y)->n = num; + (*mag)->n = num; + (*magErr)->n = num; + + if (*bounds) { + (*bounds)->x0 = xMin; + (*bounds)->x1 = xMax; + (*bounds)->y0 = yMin; + (*bounds)->y1 = yMax; + } else { + *bounds = psRegionAlloc(xMin, xMax, yMin, yMax); + } + + psTrace("psModules.objects", 8, "%ld sources: bounds are [%.2f:%.2f,%.2f:%.2f]\n", + num, xMin, xMax, yMin, yMax); + + return num; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// pmSourceMatch operations +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +static void sourceMatchFree(pmSourceMatch *match) +{ + psFree(match->mag); + psFree(match->magErr); + psFree(match->image); + psFree(match->index); +} + +pmSourceMatch *pmSourceMatchAlloc(int num // Maximum number of images + ) +{ + pmSourceMatch *match = psAlloc(sizeof(pmSourceMatch)); // Match data + psMemSetDeallocator(match, (psFreeFunc)sourceMatchFree); + + match->num = 0; + match->mag = psVectorAllocEmpty(num, PS_TYPE_F32); + match->magErr = psVectorAllocEmpty(num, PS_TYPE_F32); + match->image = psVectorAllocEmpty(num, PS_TYPE_U32); + match->index = psVectorAllocEmpty(num, PS_TYPE_U32); + + return match; +} + +void pmSourceMatchAdd(pmSourceMatch *match, // Match data + float mag, float magErr, // Magnitude and error + int image, // Image index + int index // Source index + ) +{ + int num = match->num; // Number of matches + psAssert(num <= match->mag->nalloc, "Too many matches."); + psAssert(num <= match->magErr->nalloc, "Too many matches."); + psAssert(num <= match->image->nalloc, "Too many matches."); + psAssert(num <= match->index->nalloc, "Too many matches."); + + match->mag->data.F32[num] = mag; + match->magErr->data.F32[num] = magErr; + match->image->data.S32[num] = image; + match->index->data.S32[num] = index; + match->num++; + + match->mag->n = match->magErr->n = match->image->n = match->index->n = match->num; +} + + +psArray *pmSourceMatchSources(const psArray *sourceArrays, float radius) +{ + PS_ASSERT_ARRAY_NON_NULL(sourceArrays, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, NULL); + + int numImages = sourceArrays->n; // Number of images + + // First, merge the source lists, pulling out matches + + psRegion *boundsMaster = NULL; // Bounds of source list + psVector *xMaster = NULL, *yMaster = NULL; // Coordinates of sources + long numMaster = 0; // Number in master list + psArray *matches = NULL; // Source matches (potential and actual matches) + + for (int i = 0; i < numImages; i++) { + psArray *sources = sourceArrays->data[i]; // Sources in image + if (!sources) { + continue; + } + psRegion *boundsImage = NULL; // Bounds of source list + psVector *xImage = NULL, *yImage = NULL; // Coordinates of sources + psVector *magImage = NULL, *magErrImage = NULL; // Magnitude and mag + + int numSources = sourcesParse(&boundsImage, &xImage, &yImage, &magImage, &magErrImage, + sources); // Number of sources + + if (!boundsMaster) { + // First run through --- can just copy + boundsMaster = boundsImage; + xMaster = xImage; + yMaster = yImage; + matches = psArrayAlloc(numSources); + for (int j = 0; j < numSources; j++) { + pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + matches->data[j] = match; + } + numMaster += numSources; + goto match_image_done; + } + + if (boundsImage->x0 > boundsMaster->x1 || boundsImage->x1 < boundsMaster->x0 || + boundsImage->y0 > boundsMaster->y1 || boundsImage->y1 < boundsMaster->y0) { + psTrace("psModules.objects", 7, "Bounds don't overlap\n"); + // Add everything in to the master list + long size = numMaster + numSources; // New size + xMaster = psVectorRealloc(xMaster, size); + yMaster = psVectorRealloc(yMaster, size); + matches = psArrayRealloc(matches, size); + + memcpy(&xMaster->data.F32[numMaster], xImage->data.F32, + numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + memcpy(&yMaster->data.F32[numMaster], yImage->data.F32, + numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + for (int j = 0, k = numMaster; j < numSources; j++, k++) { + pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + matches->data[k] = match; + } + + xMaster->n = size; + yMaster->n = size; + matches->n = size; + + goto match_image_done; + } + + // Match with the master list + + psTree *tree = psTreePlant(2, SOURCES_MAX_LEAF, xMaster, yMaster); // kd Tree with sources + long numMatch = 0; // Number of matches + + long size = numMaster + numSources; // New size + xMaster = psVectorRealloc(xMaster, size); + yMaster = psVectorRealloc(yMaster, size); + matches = psArrayRealloc(matches, size); + + psVector *coords = psVectorAlloc(2, PS_TYPE_F32); // Coordinates for tree lookup + for (int j = 0; j < numSources; j++) { + coords->data.F32[0] = xImage->data.F32[j]; + coords->data.F32[1] = yImage->data.F32[j]; + long index = psTreeNearestWithin(tree, coords, radius); // Match index + if (index >= 0) { + // Record the match + pmSourceMatch *match = matches->data[index]; // Match data + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + numMatch++; + } else { + // Add to the master list + pmSourceMatch *match = pmSourceMatchAlloc(numImages); // Match data + pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], i, j); + xMaster->data.F32[numMaster] = xImage->data.F32[j]; + yMaster->data.F32[numMaster] = yImage->data.F32[j]; + matches->data[numMaster] = match; + numMaster++; + xMaster->n = yMaster->n = matches->n = numMaster; + } + + } + psFree(coords); + psFree(tree); + + + match_image_done: + psFree(boundsImage); + psFree(xImage); + psFree(yImage); + psFree(magImage); + psFree(magErrImage); + } + + // Now cull the matches that contain only a single star + int numGood = 0; // Number of good matches + for (int i = 0; i < matches->n; i++) { + pmSourceMatch *match = matches->data[i]; // Match of interest + if (match->num > 1) { + if (i != numGood) { + psFree(matches->data[numGood]); + matches->data[numGood] = match; + psFree(match); + } + numGood++; + } + } + matches->n = numGood; + + return matches; +} + + + + + +psVector *pmSourceRelphot(const psArray *matches, // Array of matches + const psVector *zp, // Zero points for each image (including airmass term) + float cloudClip // Clipping for clouds + ) +{ + PS_ASSERT_ARRAY_NON_NULL(matches, NULL); + PS_ASSERT_VECTOR_NON_NULL(zp, NULL); + PS_ASSERT_VECTOR_TYPE(zp, PS_TYPE_F32, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(cloudClip, 0.0, NULL); + + return NULL; +} Index: /trunk/psModules/src/objects/pmSourceMatch.h =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.h (revision 20946) +++ /trunk/psModules/src/objects/pmSourceMatch.h (revision 20946) @@ -0,0 +1,54 @@ +#ifndef PM_SOURCE_MATCH_H +#define PM_SOURCE_MATCH_H + +/// Matched sources +/// +/// A match between sources on different images. The index back to the image and source are provided. +/// The magnitudes are pulled out, for convenience. +typedef struct { + int num; // Number of matches + psVector *image; // Image index + psVector *index; // Source index for image + psVector *mag; // Magnitudes + psVector *magErr; // Magnitude errors +} pmSourceMatch; + + +// Assert that the source match and its components are non-NULL +#define PM_ASSERT_SOURCE_MATCH_NON_NULL(MATCH, RVAL) { \ + if (!MATCH) { \ + psError(PS_ERR_UNEXPECTED_NULL, true, "Source match %s is NULL.", #MATCH); \ + return RVAL; \ + } \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->image, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->index, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->mag, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->magErr, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->image, (MATCH)->num, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->index, (MATCH)->num, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->mag, (MATCH)->num, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->magErr, (MATCH)->num, RVAL); \ +} + +/// Allocator for pmSourceMatch +pmSourceMatch *pmSourceMatchAlloc(int num // Maximum number of images + ); + +/// Add a source to a match +void pmSourceMatchAdd(pmSourceMatch *match, // Match data + float mag, float magErr, // Magnitude and error + int image, // Image index + int index // Source index + ); + + + +/// Match sources between different images +/// +/// Returns an array of psSourceMatch +psArray *pmSourceMatchSources(const psArray *sourceArrays, // Array of arrays of sources on each image + float radius // Correlation radius + ); + + +#endif