Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 35382) +++ /trunk/ppStack/src/ppStackMatch.c (revision 35383) @@ -8,4 +8,6 @@ // #define TESTING // Enable debugging output // #define TESTING_REUSE_CONV // Enable debugging for re-using convolved outputs from previous run + +// #define TESTING_CZW #ifdef TESTING_REUSE_CONV @@ -264,13 +266,13 @@ // Add the background into the target image psImage *bgImage = stackBackgroundModel(readout, config); // Image of background - psBinaryOp(fake->image, fake->image, "+", bgImage); + psBinaryOp(fake->image, fake->image, "+", bgImage); psFree(bgImage); #endif -#ifdef TESTING +#ifdef TESTING_CZW { pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "fake_%03d.fits", index); + psStringAppend(&name, "/local/tmp/fake_%03d.fits", index); pmStackVisualPlotTestImage(fake->image, name); psFits *fits = psFitsOpen(name, "w"); @@ -282,9 +284,19 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "real_%03d.fits", index); + psStringAppend(&name, "/local/tmp/real_%03d.fits", index); pmStackVisualPlotTestImage(readout->image, name); psFits *fits = psFitsOpen(name, "w"); psFree(name); psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); + psFitsClose(fits); + } + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "/local/tmp/real_var_%03d.fits", index); + pmStackVisualPlotTestImage(readout->variance, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, readout->variance, 0, NULL); psFitsClose(fits); } @@ -355,9 +367,9 @@ -#ifdef TESTING +#ifdef TESTING_CZW { pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "conv_%03d.fits", index); + psStringAppend(&name, "/local/tmp/conv_%03d.fits", index); pmStackVisualPlotTestImage(conv->image, name); psFits *fits = psFitsOpen(name, "w"); @@ -369,5 +381,15 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "diff_%03d.fits", index); + psStringAppend(&name, "/local/tmp/conv_var_%03d.fits", index); + pmStackVisualPlotTestImage(conv->variance, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, conv->variance, 0, NULL); + psFitsClose(fits); + } + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "/local/tmp/diff_%03d.fits", index); pmStackVisualPlotTestImage(fake->image, name); psFits *fits = psFitsOpen(name, "w"); @@ -510,10 +532,6 @@ } - if (!stackRenormaliseReadout(config, readout)) { - psFree(rng); - psFree(bg); - return false; - } - +#ifdef TESTING_CZW + // CZW begin check of weight before renormalize step // Measure the variance level for the weighting if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) { @@ -529,6 +547,31 @@ options->weightings->data.F32[index] = 1.0; } - psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", - index, options->weightings->data.F32[index]); + psLogMsg("ppStack", PS_LOG_INFO, "Prenorm Weighting for image %d is %g covar %g median %g\n", + index, options->weightings->data.F32[index],psImageCovarianceFactor(readout->covariance), + psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN)); + // CZW end +#endif + if (!stackRenormaliseReadout(config, readout)) { + psFree(rng); + psFree(bg); + return false; + } + + // Measure the variance level for the weighting + if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) { + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image."); + psFree(rng); + psFree(bg); + return false; + } + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * + psImageCovarianceFactor(readout->covariance)); + } else { + options->weightings->data.F32[index] = 1.0; + } + psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %g covar %g median %g\n", + index, options->weightings->data.F32[index],psImageCovarianceFactor(readout->covariance), + psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN)); psFree(rng); Index: /trunk/ppStack/src/ppStackSources.c =================================================================== --- /trunk/ppStack/src/ppStackSources.c (revision 35382) +++ /trunk/ppStack/src/ppStackSources.c (revision 35383) @@ -11,5 +11,7 @@ #endif -#ifdef TESTING +// #define TESTING_CZW + +#ifdef TESTING_CZW // Dump matches to a file static void dumpMatches(const char *filename, // File to which to dump @@ -28,7 +30,9 @@ psVectorInit(magErr, NAN); for (int j = 0; j < match->num; j++) { + if (match->mask) { if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { - continue; - } + continue; + } + } int index = match->image->data.U32[j]; // Image index mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index]; @@ -39,5 +43,11 @@ } for (int j = 0; j < num; j++) { - fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]); + fprintf(outMatches, "%f %f %f %f %f %f\t", + match->mag->data.F32[j], + zp->data.F32[j], + trans ? trans->data.F32[j] : NAN, + magErr->data.F32[j], + match->x->data.F32[j], + match->y->data.F32[j]); } fprintf(outMatches, "\n"); @@ -271,11 +281,14 @@ return false; } - -#ifdef TESTING + options->sources = psMemIncrRefCounter(matches); + psArray *matchedSources = psArrayAlloc(0); + +#ifdef TESTING_CZW dumpMatches("source_match.dat", num, matches, zp, NULL); #endif if (options->matchZPs) { - psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, + + psVector *trans = pmSourceMatchRelphot(matches, matchedSources, zp, tol, iter1, starRej1, starSys1, iter2, starRej2, starSys2, starLimit, transIter, transRej, transThresh); // Transparencies per image @@ -324,11 +337,13 @@ // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i - +#ifdef TESTING_CZW + dumpMatches("source_match2.dat", num, matches, zp, trans); +#endif if (options->matchZPs) { options->norm = psVectorAlloc(num, PS_TYPE_F32); for (int i = 0; i < num; i++) { - if (inputMask->data.U8[i] || !isfinite(trans->data.F32[i])) { - continue; - } + if (!isfinite(trans->data.F32[i])) { + continue; + } psArray *sources = sourceLists->data[i]; // Sources of interest float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure) - airmassTerm * airmassTarget; @@ -341,18 +356,25 @@ options->norm->data.F32[i] = magCorr; psLogMsg("ppStack", PS_LOG_INFO, - "Applying scale correction to image %d: %f mag (%f)\n", - i, magCorr, trans->data.F32[i]); + "Applying scale correction to image %d: %f mag (%f) (with inputMask value %d)\n", + i, magCorr, trans->data.F32[i],inputMask->data.U8[i]); for (int j = 0; j < sources->n; j++) { pmSource *source = sources->data[j]; // Source of interest + source->psfMag += magCorr; + source->apMag += magCorr; + if (!source) { continue; } // XXX need to apply to apMag as well - source->psfMag += magCorr; - source->apMag += magCorr; + psTrace("ppStack",5,"Source corrections: %d %d %f %f %f %f %f %f %f %f\n",i,j, + source->peak->xf,source->peak->yf, + source->psfMag,source->psfMagErr, + airmassTerm,airmassTarget,2.5 * log10(options->sumExposure), + source->psfMag + magCorr); } } } + matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches if (zpExpNum == numGoodImages) { @@ -366,13 +388,11 @@ -#ifdef TESTING - dumpMatches("source_mags.dat", num, matches, zp, trans); -#endif psFree(trans); -#ifdef TESTING +#ifdef TESTING_CZW // Double check: all transparencies should be zero { psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches + psArray *matchedSources = psArrayAlloc(matches->n); if (!matches) { psError(PPSTACK_ERR_DATA, false, "Unable to match sources"); @@ -380,5 +400,5 @@ return false; } - psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, + psVector *trans = pmSourceMatchRelphot(matches, matchedSources, zp, tol, iter1, starRej1, starSys1, iter2, starRej2, starSys2, starLimit, transIter, transRej, transThresh); // Transparencies @@ -386,14 +406,17 @@ fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); } + psFree(trans); psFree(matches); - } -#endif - } - - psFree(zp); + psFree(matchedSources); + } +#endif + } + + psFree(zpExp); #ifdef ASTROMETRY + // CZW: This is off by default. // Position offsets { @@ -436,4 +459,5 @@ #if (defined TESTING && defined ASTROMETRY) + // CZW: This is off by default. // Double check: all offsets should be zero { @@ -458,40 +482,59 @@ if (options->photometry) { // Save best matches for future photometry - psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n); - + psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n); // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible //int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required //MEH - if many poor images, will hit oddity with num fraction. use numGoodImages instead - int minMatches = PS_MAX(2, fracMatch * numGoodImages);// Minimum number of matches required - - int good0 = 0; - for (int i = 0; i < num; i++) { - if (inputMask->data.U8[i]) { - continue; - } - good0 = i; - break; + int minMatches = PS_MAX(2, fracMatch * numGoodImages);// Minimum number of matches required + + int good0 = 0; + for (int i = 0; i < num; i++) { + if (inputMask->data.U8[i]) { + continue; } - - for (int i = 0; i < matches->n; i++) { - pmSourceMatch *match = matches->data[i]; // Match of interest - if (match->num < minMatches) { - continue; - } - - // We need to grab a single instance of this source: just take the first available - // MEH - if first available is a rejected image, then problem. so take first non-rejected instance above - //int image = match->image->data.S32[0]; // Index of image - //int index = match->index->data.S32[0]; // Index of source within image - int image = match->image->data.S32[good0]; // Index of image - int index = match->index->data.S32[good0]; // Index of source within image - psArray *sources = sourceLists->data[image]; // Sources for image - pmSource *source = sources->data[index]; // Source of interest - - psArrayAdd(sourcesBest, sourcesBest->n, source); - } - psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); - } - + good0 = i; + break; + } + + for (int i = 0; i < matches->n; i++) { + pmSourceMatch *match = matches->data[i]; // Match of interest + if (match->num < minMatches) { + continue; + } + + // We need to grab a single instance of this source: just take the first available + // MEH - if first available is a rejected image, then problem. so take first non-rejected instance above + // int image = match->image->data.S32[0]; // Index of image + // int index = match->index->data.S32[0]; // Index of source within image + int image = match->image->data.S32[good0]; // Index of image + int index = match->index->data.S32[good0]; // Index of source within image + psArray *sources = sourceLists->data[image]; // Sources for image + pmSource *source = sources->data[index]; // Source of interest + + // CZW: Accept the first source for all properties, except for the magnitudes. + // We have average magnitudes from relphot, so we should use those + pmSource *photsource = matchedSources->data[i]; + if ((photsource)&&(isfinite(photsource->psfMag))) { + // CZW: photsource is the average relphot mag, which is equal to + // (primed are averaged values for this discussion) + // Therefore, the correct magnitude for psfMag' is: + // psfMag' = (photsource - zpIN_j + magCorr) + // apMag - psfMag should be retained, so: + // apMag - psfMag = apMag' - psfMag' + // apMag' = apMag + psfMag' - psfMag + source->apMag = source->apMag + (photsource->psfMag + options->norm->data.F32[image] - zp->data.F32[image] - source->psfMag); + source->psfMag = photsource->psfMag + options->norm->data.F32[image] - zp->data.F32[image]; + + psTrace("ppStack", 5, "Best Source %d/%ld @ %f %f %f %f %f %d",i,sourcesBest->n,source->peak->xf,source->peak->yf, + source->psfMag, + photsource->psfMag + options->norm->data.F32[image] - zp->data.F32[image] , + source->psfMagErr,-1); + psArrayAdd(sourcesBest, sourcesBest->n, source); + } + } + psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); + + } + psFree(zp); psFree(matches); Index: /trunk/ppStack/src/ppStackTarget.c =================================================================== --- /trunk/ppStack/src/ppStackTarget.c (revision 35382) +++ /trunk/ppStack/src/ppStackTarget.c (revision 35383) @@ -60,5 +60,6 @@ } float xSource, ySource; // Coordinates of source - coordsFromSource(&xSource, &ySource, NULL, source); + float Smag; // magnitude of source + coordsFromSource(&xSource, &ySource, &Smag, source); coords->data.F64[0] = xSource; @@ -66,6 +67,6 @@ psVector *indices = psTreeAllWithin(tree, coords, exclusion); // Number within exclusion zone - psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", - coords->data.F64[0], coords->data.F64[1], indices->n); + psTrace("ppStack", 5, "Filtering: Source at %.0lf %.0lf with mag %f has %ld sources in exclusion zone", + coords->data.F64[0], coords->data.F64[1], Smag, indices->n); if (indices->n == 1) { // Only itself inside the exclusion zone Index: /trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.c (revision 35382) +++ /trunk/psModules/src/objects/pmSourceMatch.c (revision 35383) @@ -630,4 +630,5 @@ psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches + psArray *matchedSources, // Array of averaged sources const psVector *zp, // Zero points for each image (including airmass term) float tol, // Relative tolerance for convergence @@ -648,5 +649,6 @@ PS_ASSERT_VECTOR_TYPE(zp, PS_TYPE_F32, NULL); PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL); - + PS_ASSERT_ARRAY_NON_NULL(matchedSources, NULL); + sys1 *= sys1; sys2 *= sys2; @@ -742,4 +744,10 @@ } + for (int i = 0; i < matches->n; i++) { + pmSource *source = pmSourceAlloc(); + source->psfMag = stars->data.F32[i]; + psArrayAdd(matchedSources,matchedSources->n,source); + } + psFree(photo); psFree(badImage); Index: /trunk/psModules/src/objects/pmSourceMatch.h =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.h (revision 35382) +++ /trunk/psModules/src/objects/pmSourceMatch.h (revision 35383) @@ -74,4 +74,5 @@ /// Perform relative photometry to calibrate images psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches + psArray *matchedSources, // Array of average sources const psVector *zp, // Zero points for each image (including airmass term) float tol, // Relative tolerance for convergence