Index: trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- trunk/psModules/src/objects/pmSourceMatch.c (revision 20949) +++ trunk/psModules/src/objects/pmSourceMatch.c (revision 20953) @@ -63,5 +63,5 @@ pmSource *source = sources->data[i]; // Source of interest if (!source || (source->mode & SOURCE_MASK) || !isfinite(source->psfMag) || - source->psfMag > SOURCE_FAINTEST) { + !isfinite(source->errMag) || source->psfMag > SOURCE_FAINTEST) { continue; } @@ -122,4 +122,5 @@ match->image = psVectorAllocEmpty(num, PS_TYPE_U32); match->index = psVectorAllocEmpty(num, PS_TYPE_U32); + match->mask = psVectorAllocEmpty(num, PS_TYPE_MASK); return match; @@ -142,7 +143,8 @@ match->image->data.S32[num] = image; match->index->data.S32[num] = index; + match->index->data.PS_TYPE_MASK_DATA[num] = 0; match->num++; - match->mag->n = match->magErr->n = match->image->n = match->index->n = match->num; + match->mag->n = match->magErr->n = match->image->n = match->index->n = match->mask->n = match->num; } @@ -287,7 +289,11 @@ psVector *stars, // Star magnitudes const psArray *matches, // Array of matches - const psVector *zp // Zero points for each image (including airmass term) + const psVector *zp, // Zero points for each image (including airmass term) + const psVector *photo // Photometric image? ) { + psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points"); + psAssert(matches, "Need list of matches"); + int numImages = zp->n; // Number of images int numStars = matches->n; // Number of stars @@ -297,72 +303,195 @@ psAssert(stars && stars->type.type == PS_TYPE_F32, "Need star magnitudes"); psAssert(stars->n == numStars, "Not enough stars: %ld\n", stars->n); - psAssert(matches, "Need list of matches"); - psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points"); psAssert(zp->n == numImages, "Not enough ZPs: %ld", zp->n); - - psVector *transErr = psVectorAlloc(numImages, PS_TYPE_F32); // Errors in the transparency - psVector *starsErr = psVectorAlloc(numStars, PS_TYPE_F32); // Errors in the star magnitudes + psAssert(!photo || photo->type.type == PS_TYPE_U8, "Photometric determination is wrong type"); + psAssert(!photo || photo->n == numImages, "Not enough photometric determinations: %ld", photo->n); // Solve the star magnitudes psVectorInit(stars, 0.0); - psVectorInit(starsErr, 0.0); + 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 star = 0.0, starErr = 0.0; // Accumulators for star + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_MASK_DATA[j]) { + continue; + } + numMeasurements++; + int index = match->image->data.U32[j]; // Image index + float mag = match->mag->data.F32[j]; // Measured magnitude + double magErr = match->magErr->data.F32[j]; // Error in measured magnitude + double invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error + float cal = zp->data.F32[index]; // Calibration to apply to image + if (!photo || !photo->data.U8[index]) { + cal -= trans->data.F32[index]; + } + star += (mag + cal) * invErr2; + starErr += invErr2; + } + if (numMeasurements > 1) { + // It's only a good star (contributing to the chi^2) if there's more than 1 measurement + numGoodStars++; + } + stars->data.F32[i] = star / starErr; + } + + // Solve for the transparencies + // We solve even for the "photometric" images since they may jump out of that status upon iteration + psVector *accum = psVectorAlloc(numImages, PS_TYPE_F64); // Transparency accumulator + psVector *accumErr = psVectorAlloc(numImages, PS_TYPE_F64); // Transparency accumulator + psVectorInit(accum, 0.0); + psVectorInit(accumErr, 0.0); 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_MASK_DATA[j]) { + continue; + } + int index = match->image->data.U32[j]; // Image index + float mag = match->mag->data.F32[j]; // Measured magnitude + double magErr = match->magErr->data.F32[j]; // Error in measured magnitude + double invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error + float cal = zp->data.F32[index]; // Calibration to apply to image + accum->data.F64[index] += (mag + cal - stars->data.F32[i]) * invErr2; + accumErr->data.F64[index] += invErr2; + } + } + for (int i = 0; i < numImages; i++) { + trans->data.F32[i] = accum->data.F64[i] / accumErr->data.F64[i]; + + psTrace("psModules.objects", 3, "Transparency for image %d: %f\n", i, trans->data.F32[i]); + } + psFree(accum); + psFree(accumErr); + + // 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 + 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 float mag = match->mag->data.F32[j]; // Measured magnitude float magErr = match->magErr->data.F32[j]; // Error in measured magnitude - float invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error - stars->data.F32[i] += (mag + zp->data.F32[index] - trans->data.F32[index]) * invErr2; - starsErr->data.F32[i] += invErr2; - } - } - for (int i = 0; i < numStars; i++) { - float inverse = 1.0 / starsErr->data.F32[i]; // Inverse error - stars->data.F32[i] *= inverse; - starsErr->data.F32[i] = sqrtf(inverse); - } - - // Solve for the transparencies - psVectorInit(trans, 0.0); - psVectorInit(transErr, 0.0); + float cal = zp->data.F32[index]; // Calibration to apply to image + if (!photo || !photo->data.U8[index]) { + cal -= trans->data.F32[index]; + } + float dev = (mag + cal - stars->data.F32[index]) / magErr; // Deviation from fit + if (isfinite(dev)) { + chi2 += PS_SQR(dev); + dof++; + } + } + } + dof -= numGoodStars + numImages; + chi2 /= dof; + + return chi2; +} + +// Determine which images are photometric, based on estimated transparencies +// Returns number of photometric images, or -1 on error +static int sourceMatchRelphotPhotometric(psVector *photo, // Photometric determination + const psVector *trans, // Estimated transparencies + int transIter, // Iterations for transparency + float transClip, // Clipping level for transparency + float photoLevel // Level below which we declare photometric + ) +{ + psAssert(photo && photo->type.type == PS_TYPE_U8, "Need photometric determination"); + psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies"); + + int numImages = photo->n; // Number of images + + psAssert(trans->n == numImages, "Not enough transparencies: %ld", trans->n); + psAssert(transIter >= 0, "Iterations for transparency must be non-negative: %d", transIter); + psAssert(transClip > 0, "Clipping level for transparency must be positive: %f", transClip); + psAssert(photoLevel > 0, "Photometric level must be positive: %f", photoLevel); + + psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); // Statistics + stats->clipIter = transIter; + stats->clipSigma = transClip; + + if (!psVectorStats(stats, trans, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on transparencies."); + psFree(stats); + return -1; + } + + float thresh = stats->clippedMean + photoLevel * stats->clippedStdev; // Threshold for clouds + psFree(stats); + + int numPhoto = 0; // Number of photometric images + for (int i = 0; i < numImages; i++) { + if (trans->data.F32[i] < thresh) { + photo->data.U8[i] = 0xFF; + numPhoto++; + } else { + photo->data.U8[i] = 0; + } + } + + return numPhoto; +} + + +// Reject star measurements +// Returns the fraction of measurements that were rejected +static float sourceMatchRelphotReject(psVector *trans, // Transparencies + psVector *stars, // Star magnitudes + const psArray *matches, // Array of matches + const psVector *zp, // Zero points for each image + const psVector *photo, // Photometric image? + float starClip, // Clipping for stars + float sysErr // Systematic error + ) +{ + psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points"); + psAssert(matches, "Need list of matches"); + + int numImages = zp->n; // Number of images + int numStars = matches->n; // Number of stars + + psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies"); + psAssert(trans->n == numImages, "Not enough transparencies: %ld\n", trans->n); + psAssert(stars && stars->type.type == PS_TYPE_F32, "Need star magnitudes"); + psAssert(stars->n == numStars, "Not enough stars: %ld\n", stars->n); + psAssert(zp->n == numImages, "Not enough ZPs: %ld", zp->n); + psAssert(!photo || photo->type.type == PS_TYPE_U8, "Photometric determination is wrong type"); + psAssert(!photo || photo->n == numImages, "Not enough photometric determinations: %ld", photo->n); + + starClip = PS_SQR(starClip); + sysErr = PS_SQR(sysErr); + + int numRejected = 0; // Number rejected + int numMeasurements = 0; // Number of measurements 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_MASK_DATA[j]) { + continue; + } + numMeasurements++; int index = match->image->data.U32[j]; // Image index float mag = match->mag->data.F32[j]; // Measured magnitude float magErr = match->magErr->data.F32[j]; // Error in measured magnitude - float invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error - trans->data.F32[index] += (mag - stars->data.F32[i]) * invErr2; - transErr->data.F32[index] += invErr2; - } - } - for (int i = 0; i < numImages; i++) { - float inverse = 1.0 / transErr->data.F32[i]; // Inverse error - trans->data.F32[i] *= inverse; - transErr->data.F32[i] = sqrtf(inverse); - - psTrace("psModules.objects", 3, "Transparency for image %d: %f +/- %f\n", - i, trans->data.F32[i], transErr->data.F32[i]); - } - - // Once more through to evaluate chi^2 - float chi2 = 0.0; // chi^2 for iteration - for (int i = 0; i < numStars; i++) { - pmSourceMatch *match = matches->data[i]; // Matched stars - for (int j = 0; j < match->num; j++) { - int index = match->image->data.U32[j]; // Image index - float mag = match->mag->data.F32[j]; // Measured magnitude - float magErr = match->magErr->data.F32[j]; // Error in measured magnitude - chi2 += PS_SQR(mag - stars->data.F32[index] - trans->data.F32[index] + zp->data.F32[index]) / - PS_SQR(magErr); - } - } - chi2 /= numStars + numImages; - - psFree(transErr); - psFree(starsErr); - - return chi2; + float cal = zp->data.F32[index]; // Calibration to apply to image + if (!photo || !photo->data.U8[index]) { + cal -= trans->data.F32[index]; + } + float dev = mag + cal - stars->data.F32[i]; // Deviation + if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr)) { + numRejected++; + match->mask->data.PS_TYPE_MASK_DATA[j] = 0xFF; + } + } + } + + return (float)numRejected / (float)numMeasurements; } @@ -372,6 +501,10 @@ int maxIter, // Maximum number of iterations float tol, // Relative tolerance for convergence - float cloudClip, // Clipping for clouds - float starClip // Clipping for stars + float rejLimit, // Limit on rejection between iterations + int transIter, // Clipping iterations for transparency + float transClip, // Clipping level for transparency + float photoLevel, // Level at which we declare image is photometric + float starClip, // Clipping for stars + float sysErr // Systematic error in measurements ) { @@ -379,5 +512,5 @@ 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); + PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL); int numImages = zp->n; // Number of images @@ -385,13 +518,39 @@ psVector *trans = psVectorAlloc(numImages, PS_TYPE_F32); // Transparencies for each image, magnitudes psVectorInit(trans, 0.0); + psVector *photo = psVectorAlloc(numImages, PS_TYPE_U8); // Photometric determination for each image + psVectorInit(photo, 0); psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star - float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp); // chi^2 for solution + float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo); // chi^2 for solution psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2); float lastChi2 = INFINITY; // chi^2 on last iteration - for (int i = 0; i < maxIter && (lastChi2 - chi2) / chi2 > tol; i++) { + float fracRej = INFINITY; // Fraction of measurements rejected + for (int i = 0; i < maxIter && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) { lastChi2 = chi2; - chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp); // chi^2 for solution - psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f\n", i, chi2); + + // Identify photometric nights + int numPhoto = sourceMatchRelphotPhotometric(photo, trans, transIter, transClip, + photoLevel); // Number of photometric images + if (numPhoto < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform photometric determination"); + psFree(trans); + psFree(photo); + psFree(stars); + return NULL; + } + psTrace("psModules.objects", 3, "Determined %d/%d are photometric", numPhoto, numImages); + + fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, starClip, sysErr); + psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100); + + chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp, photo); // chi^2 for solution + psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); + } + + psFree(photo); + psFree(stars); + + if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) { + psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej); } Index: trunk/psModules/src/objects/pmSourceMatch.h =================================================================== --- trunk/psModules/src/objects/pmSourceMatch.h (revision 20949) +++ trunk/psModules/src/objects/pmSourceMatch.h (revision 20953) @@ -12,4 +12,5 @@ psVector *mag; // Magnitudes psVector *magErr; // Magnitude errors + psVector *mask; // Mask for measurements } pmSourceMatch; @@ -56,6 +57,10 @@ int maxIter, // Maximum number of iterations float tol, // Relative tolerance for convergence - float cloudClip, // Clipping for clouds - float starClip // Clipping for stars + float rejLimit, // Limit on rejection between iterations + int transIter, // Clipping iterations for transparency + float transClip, // Clipping level for transparency + float photoLevel, // Level at which we declare image is photometric + float starClip, // Clipping for stars + float sysErr // Systematic error in measurements );