Index: trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- trunk/psModules/src/objects/pmSourceMatch.c (revision 26070) +++ trunk/psModules/src/objects/pmSourceMatch.c (revision 26076) @@ -17,4 +17,5 @@ #define SOURCES_MAX_LEAF 2 // Maximum number of points on a tree leaf #define ARRAY_BUFFER 16 // Buffer for array + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -124,4 +125,6 @@ match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); + match->x = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); + match->y = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32); match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32); @@ -133,4 +136,5 @@ void pmSourceMatchAdd(pmSourceMatch *match, // Match data float mag, float magErr, // Magnitude and error + float x, float y, // Position int image, // Image index int index // Source index @@ -141,4 +145,6 @@ match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1); match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1); + match->x = psVectorExtend(match->x, match->x->nalloc, 1); + match->y = psVectorExtend(match->y, match->y->nalloc, 1); match->image = psVectorExtend(match->image, match->image->nalloc, 1); match->index = psVectorExtend(match->index, match->index->nalloc, 1); @@ -147,4 +153,6 @@ match->mag->data.F32[num] = mag; match->magErr->data.F32[num] = magErr; + match->x->data.F32[num] = x; + match->y->data.F32[num] = y; match->image->data.S32[num] = image; match->index->data.S32[num] = index; @@ -192,5 +200,5 @@ pmSourceMatch *match = pmSourceMatchAlloc(); // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); matches->data[j] = match; } @@ -212,5 +220,5 @@ pmSourceMatch *match = pmSourceMatchAlloc(); // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); matches->data[k] = match; } @@ -238,5 +246,5 @@ pmSourceMatch *match = matches->data[index]; // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); numMatch++; } else { @@ -244,5 +252,5 @@ pmSourceMatch *match = pmSourceMatchAlloc(); // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); xMaster->data.F32[numMaster] = xImage->data.F32[j]; yMaster->data.F32[numMaster] = yImage->data.F32[j]; @@ -304,5 +312,6 @@ pmSourceMatch *match = matches->data[i]; // Match of interest for (int j = 0; j < match->num && !source; j++) { - if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) { + if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j]) || + !isfinite(match->x->data.F32[j]) || !isfinite(match->y->data.F32[j])) { continue; } @@ -365,5 +374,5 @@ double star = 0.0, starErr = 0.0; // Accumulators for star for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -396,5 +405,5 @@ pmSourceMatch *match = matches->data[i]; // Matched stars for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -424,5 +433,5 @@ pmSourceMatch *match = matches->data[i]; // Matched stars for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -543,5 +552,5 @@ pmSourceMatch *match = matches->data[i]; // Matched stars for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -564,5 +573,5 @@ if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) { numRejected++; - match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; + match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_PHOT; } } @@ -627,8 +636,8 @@ return NULL; } - psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages); + psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric\n", numPhoto, numImages); fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1); - psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100); + psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100); chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1); @@ -649,8 +658,8 @@ return NULL; } - psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages); + psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric\n", numPhoto, numImages); fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2); - psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100); + psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100); chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2); @@ -668,2 +677,220 @@ return trans; } + + +// Iterate on the star positions and image shifts +// Returns the solution chi^2 +static float sourceMatchRelastroIterate(psVector *xShift, psVector *yShift, // Shift for image + psVector *xStar, psVector *yStar, // Position for star + const psArray *matches // Array of matches + ) +{ + psAssert(matches, "Need list of matches"); + + int numImages = xShift->n; // Number of images + int numStars = matches->n; // Number of stars + + psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32, + "Need shifts"); + psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n); + psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32, + "Need star positions"); + psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n"); + + // Solve the star positions + psVectorInit(xStar, NAN); + psVectorInit(yStar, NAN); + psVector *starMask = psVectorAlloc(numStars, PS_TYPE_U8); // Mask for stars + psVectorInit(starMask, 0xFF); + int numGoodStars = 0; // Number of stars with good measurements + for (int i = 0; i < numStars; i++) { + pmSourceMatch *match = matches->data[i]; // Matched stars + int numMeasurements = 0; // Number of unmasked measurements for star + double xSum = 0.0, ySum = 0.0; // Accumulators for star + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) { + continue; + } + numMeasurements++; + int index = match->image->data.U32[j]; // Image index + + xSum += match->x->data.F32[j] - xShift->data.F32[index]; + ySum += match->y->data.F32[j] - yShift->data.F32[index]; + } + if (numMeasurements > 1) { + // It's only a good star (contributing to the chi^2) if there's more than 1 measurement + numGoodStars++; + xStar->data.F32[i] = xSum / numMeasurements; + yStar->data.F32[i] = ySum / numMeasurements; + starMask->data.U8[i] = 0; + } + } + + // Solve for the shifts + psVectorInit(xShift, 0.0); + psVectorInit(yShift, 0.0); + psVector *num = psVectorAlloc(numImages, PS_TYPE_S32); // Number of stars + psVectorInit(num, 0); + for (int i = 0; i < numStars; i++) { + if (starMask->data.U8[i]) { + continue; + } + pmSourceMatch *match = matches->data[i]; // Matched stars + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) { + continue; + } + int index = match->image->data.U32[j]; // Image index + + xShift->data.F32[index] += match->x->data.F32[j] - xStar->data.F32[i]; + yShift->data.F32[index] += match->y->data.F32[j] - yStar->data.F32[i]; + num->data.S32[index]++; + } + } + for (int i = 0; i < numImages; i++) { + xShift->data.F32[i] /= num->data.S32[i]; + yShift->data.F32[i] /= num->data.S32[i]; + psTrace("psModules.objects", 3, "Shift for image %d: %f,%f\n", + i, xShift->data.F32[i], yShift->data.F32[i]); + } + psFree(num); + + // Once more through to evaluate chi^2 + float chi2 = 0.0; // chi^2 for iteration + int dof = 0; // Degrees of freedom + for (int i = 0; i < numStars; i++) { + pmSourceMatch *match = matches->data[i]; // Matched stars + if (starMask->data.U8[i]) { + continue; + } + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + continue; + } + + int index = match->image->data.U32[j]; // Image index + float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i]; + float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i]; + + chi2 += PS_SQR(dx) + PS_SQR(dy); + dof++; + } + } + dof -= numGoodStars + numImages; + chi2 /= dof; + + return chi2; +} + +// Reject star measurements +// Returns the fraction of measurements that were rejected +static float sourceMatchRelastroReject(const psVector *xShift, const psVector *yShift, // Shifts for each image + const psVector *xStar, const psVector *yStar, // Positions for each star + const psArray *matches, // Array of matches + float chi2, // chi^2 from fit + float rej // Rejection threshold + ) +{ + psAssert(matches, "Need list of matches"); + + int numImages = xShift->n; // Number of images + int numStars = matches->n; // Number of stars + + psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32, + "Need shifts"); + psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n); + psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32, + "Need star positions"); + psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n"); + + int numRejected = 0; // Number rejected + int numMeasurements = 0; // Number of measurements + + float thresh = PS_SQR(rej) * chi2; // Threshold for rejection + + for (int i = 0; i < numStars; i++) { + pmSourceMatch *match = matches->data[i]; // Matched stars + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) { + continue; + } + numMeasurements++; + int index = match->image->data.U32[j]; // Image index + + float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i]; + float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i]; + + if (PS_SQR(dx) + PS_SQR(dy) > thresh) { + numRejected++; + match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_ASTRO; + } + } + } + + return (float)numRejected / (float)numMeasurements; +} + +psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches + int numImages, // Number of images + float tol, // Relative tolerance for convergence + int iter1, // Number of iterations for pass 1 + float rej1, // Limit on rejection between iterations for pass 1 + int iter2, // Number of iterations for pass 2 + float rej2, // Limit on rejection between iterations for pass 2 + float rejLimit // Limit on rejection between iterations + ) +{ + PS_ASSERT_ARRAY_NON_NULL(matches, NULL); + + int numStars = matches->n; // Number of stars + psVector *xShift = psVectorAlloc(numImages, PS_TYPE_F32); // x shift for each image + psVector *yShift = psVectorAlloc(numImages, PS_TYPE_F32); // y shift for each image + psVectorInit(xShift, 0.0); + psVectorInit(yShift, 0.0); + psVector *xStar = psVectorAlloc(numStars, PS_TYPE_F32); // x position for each star + psVector *yStar = psVectorAlloc(numStars, PS_TYPE_F32); // y position for each star + + float chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); // chi^2 for solution + psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2); + float lastChi2 = INFINITY; // chi^2 on last iteration + float fracRej = INFINITY; // Fraction of measurements rejected + + // In the first passes, the shifts are not well deteremined: use high systematic error and + // rejection thresholds + for (int i = 0; i < iter1; i++) { + fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej1); + psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100); + + chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); + psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); + } + + for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) { + lastChi2 = chi2; + + fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej2); + psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100); + + chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); + psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); + } + + psFree(xStar); + psFree(yStar); + + if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) { + psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej); + } + + psArray *results = psArrayAlloc(numImages); // Array of results + for (int i = 0; i < numImages; i++) { + psVector *offset = results->data[i] = psVectorAlloc(2, PS_TYPE_F32); // Offset for image + offset->data.F32[0] = xShift->data.F32[i]; + offset->data.F32[1] = yShift->data.F32[i]; + } + psFree(xShift); + psFree(yShift); + + return results; +} + Index: trunk/psModules/src/objects/pmSourceMatch.h =================================================================== --- trunk/psModules/src/objects/pmSourceMatch.h (revision 26070) +++ trunk/psModules/src/objects/pmSourceMatch.h (revision 26076) @@ -1,4 +1,12 @@ #ifndef PM_SOURCE_MATCH_H #define PM_SOURCE_MATCH_H + +#include + +/// Mask values for matched sources +typedef enum { + PM_SOURCE_MATCH_MASK_PHOT = 0x01, // Source was rejected from photometry fit + PM_SOURCE_MATCH_MASK_ASTRO = 0x02, // Source was rejected from astrometry fit +} pmSourceMatchMask; /// Matched sources @@ -12,4 +20,5 @@ psVector *mag; // Magnitudes psVector *magErr; // Magnitude errors + psVector *x, *y; // Positions psVector *mask; // Mask for measurements } pmSourceMatch; @@ -26,8 +35,12 @@ PS_ASSERT_VECTOR_NON_NULL((MATCH)->mag, RVAL); \ PS_ASSERT_VECTOR_NON_NULL((MATCH)->magErr, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->x, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->y, 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); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->x, (MATCH)->num, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->y, (MATCH)->num, RVAL); \ } @@ -38,4 +51,5 @@ void pmSourceMatchAdd(pmSourceMatch *match, // Match data float mag, float magErr, // Magnitude and error + float x, float y, // Position int image, // Image index int index // Source index @@ -75,3 +89,14 @@ ); +/// Perform relative astrometry to calibrate images +psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches + int numImages, // Number of images + float tol, // Relative tolerance for convergence + int iter1, // Number of iterations for pass 1 + float rej1, // Limit on rejection between iterations for pass 1 + int iter2, // Number of iterations for pass 2 + float rej2, // Limit on rejection between iterations for pass 2 + float rejLimit // Limit on rejection between iterations + ); + #endif