Index: branches/eam_branches/ipp-20130307/ppStack/src/ppStackMatch.c =================================================================== --- branches/eam_branches/ipp-20130307/ppStack/src/ppStackMatch.c (revision 35324) +++ branches/eam_branches/ipp-20130307/ppStack/src/ppStackMatch.c (revision 35413) @@ -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: branches/eam_branches/ipp-20130307/ppStack/src/ppStackPrepare.c =================================================================== --- branches/eam_branches/ipp-20130307/ppStack/src/ppStackPrepare.c (revision 35324) +++ branches/eam_branches/ipp-20130307/ppStack/src/ppStackPrepare.c (revision 35413) @@ -1,3 +1,9 @@ # include "ppStack.h" + + +// Option to use the old FWHM rejection code. +# define NEW_REJECTION +// #define OLD_REJECTION + # define RE_PHOTOMETER 0 @@ -253,5 +259,6 @@ float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs - + float threshFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.THRESH"); // FWHM size that we refuse to clip below + float asymmetryFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.ASYMMETRY"); // max bimodal asymmetry psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message bool havePSFs = false; // Do we have any PSFs? @@ -286,5 +293,5 @@ options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM; } - +#ifdef OLD_REJECTION // reject any input images which exceed the specified max FWHM if (isfinite(maxFWHM) && (options->inputSeeing->data.F32[i] > maxFWHM)) { @@ -292,5 +299,6 @@ psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f maxFWHM) --- rejected.", i, options->inputSeeing->data.F32[i], maxFWHM); } - + // End old rejection fixed limit. +#endif psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]); } @@ -300,4 +308,72 @@ psFree(log); +#ifdef NEW_REJECTION + // Do GMM test with two modes to decide where to put the break. + double Punimodal; + int m = 2; + psVector *modes = psVectorAlloc(num,PS_TYPE_F32); + psVector *means = psVectorAlloc(m,PS_TYPE_F32); + psVector *S = psVectorAlloc(m,PS_TYPE_F32); + psVector *pi = psVectorAlloc(m,PS_TYPE_F32); + psImage *P = psImageAlloc(m,num,PS_TYPE_F32); + float limit; + if (!psMM1DClass(options->inputSeeing, + options->inputSeeing->n, + modes,means, + S,pi,P, + 2, + &Punimodal)) { + // Handle error here + psFree(modes); + psFree(means); + psFree(S); + psFree(pi); + psFree(P); + return(false); + } + // fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]); + // fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]); + + // fprintf(stderr,"pi: %g %g\n",pi->data.F32[0],pi->data.F32[1]); + + // unwrittenGMMfunction(options->inputSeeing, options->inputMask, 0xff, + // &Punimodal, + // &m1,&s1,&pi1, + // &m2,&s2,&pi2); + + // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit + if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode + (num <= 4)) { // Or we have a small number of inputs, making this statistic poor. + limit = maxFWHM; + } + else { // This is a bimodal distribution + if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)|| // However, both modes are equally populated + (pi->data.F32[1] > pi->data.F32[0])) { // Or the larger FWHM mode is more populated + limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1]; + } + else { // The smaller FWHM mode is more populated + limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0]; + } + } + if (limit > maxFWHM) { limit = maxFWHM; } // We should not be larger than our max + if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min + + for (int i = 0; i < num; i++) { + if (options->inputSeeing->data.F32[i] > limit) { + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF; + psLogMsg("ppStack", PS_LOG_INFO, + "PSF FWHM for image %d is too large (%f vs %f fwhmlimit) --- rejected", + i, options->inputSeeing->data.F32[i],limit); + } + } + psFree(means); + psFree(modes); + psFree(S); + psFree(pi); + psFree(P); + // End new rejection code +#endif + +#ifdef OLD_REJECTION // We should have the ability to filter the input list based on the seeing: // * reject above some max value and/or min value @@ -323,5 +399,6 @@ } } - + // End old rejection sigma clip +#endif // Generate target PSF if (options->convolve) { Index: branches/eam_branches/ipp-20130307/ppStack/src/ppStackSources.c =================================================================== --- branches/eam_branches/ipp-20130307/ppStack/src/ppStackSources.c (revision 35324) +++ branches/eam_branches/ipp-20130307/ppStack/src/ppStackSources.c (revision 35413) @@ -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: branches/eam_branches/ipp-20130307/ppStack/src/ppStackTarget.c =================================================================== --- branches/eam_branches/ipp-20130307/ppStack/src/ppStackTarget.c (revision 35324) +++ branches/eam_branches/ipp-20130307/ppStack/src/ppStackTarget.c (revision 35413) @@ -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