Index: /branches/eam_branches/ipp-20130419/ppStack/src/ppStackConvolve.c =================================================================== --- /branches/eam_branches/ipp-20130419/ppStack/src/ppStackConvolve.c (revision 35552) +++ /branches/eam_branches/ipp-20130419/ppStack/src/ppStackConvolve.c (revision 35553) @@ -403,5 +403,6 @@ continue; } - if ((options->matchChi2->data.F32[i] > thresh) || ! isfinite(options->matchChi2->data.F32[i])) { + if ((options->matchChi2->data.F32[i] > thresh) || + ! isfinite(options->matchChi2->data.F32[i])) { numRej++; options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2; Index: /branches/eam_branches/ipp-20130419/ppStack/src/ppStackMatch.c =================================================================== --- /branches/eam_branches/ipp-20130419/ppStack/src/ppStackMatch.c (revision 35552) +++ /branches/eam_branches/ipp-20130419/ppStack/src/ppStackMatch.c (revision 35553) @@ -467,7 +467,10 @@ sum += kernels->rms; num++; + // psLogMsg("ppStack", PS_LOG_INFO, "KERNMATCHDUMP %d %d %g %g %g\n", + // index,num,kernels->mean,kernels->rms,sum); } psFree(iter); options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); + } Index: /branches/eam_branches/ipp-20130419/ppStack/src/ppStackPrepare.c =================================================================== --- /branches/eam_branches/ipp-20130419/ppStack/src/ppStackPrepare.c (revision 35552) +++ /branches/eam_branches/ipp-20130419/ppStack/src/ppStackPrepare.c (revision 35553) @@ -257,8 +257,10 @@ bool mdok = false; + bool simpleClip = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE"); 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? @@ -309,54 +311,105 @@ #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]; + // Begin new rejection code. + float limit = INFINITY; // Default to a big value. + if (!isfinite(threshFWHM)) { // If this isn't set, initialize it to zero. No inputs may be clipped if smaller than this + threshFWHM = 0.0; + } + if (!isfinite(maxFWHM)) { // No maxFWHM was set, accept all inputs + maxFWHM = INFINITY; + } + limit = maxFWHM; // Initialize the limit to be the maxFWHM + if (limit < threshFWHM) { + limit = threshFWHM; + } + if (!isfinite(clipFWHMnSig)) { // We cannot do a sigma clip with nSig, so use the maxFWHM to set the limit. + psLogMsg("ppStack",PS_LOG_INFO, + "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)", + limit, + clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM); + } + if (simpleClip) { // Do a sigma clip like the old rejection code. + psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff); + psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev); + + options->clippedMean = fwhmStats->clippedMean; + options->clippedStdev = fwhmStats->clippedStdev; + limit = options->clippedMean + clipFWHMnSig * options->clippedStdev; + if (limit < threshFWHM) { + limit = threshFWHM; } - else { // The smaller FWHM mode is more populated - limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0]; + if (limit > maxFWHM) { + limit = maxFWHM; } - } - if (limit > maxFWHM) { limit = maxFWHM; } // We should not be larger than our max - if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min - + psLogMsg("ppStack",PS_LOG_INFO, + "PSF FWHM distribution: Used simple clip method. Limit set to %f.", + limit); + psFree(fwhmStats); + } + else { // Do the mixture model rejection. + if(!isfinite(asymmetryFWHM)) { // Or not, because the parameters aren't set. + psLogMsg("ppStack",PS_LOG_INFO, + "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)", + limit, + clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM); + } + else { + // 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); + + 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]); + + // 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 + psLogMsg("ppStack",PS_LOG_INFO, + "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f", + limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0], + means->data.F32[1],S->data.F32[1],pi->data.F32[1], + Punimodal); + psFree(means); + psFree(modes); + psFree(S); + psFree(pi); + psFree(P); + } // End mixture model case + } + // Perform rejection using the limit set for (int i = 0; i < num; i++) { if (options->inputSeeing->data.F32[i] > limit) { @@ -367,9 +420,4 @@ } } - psFree(means); - psFree(modes); - psFree(S); - psFree(pi); - psFree(P); // End new rejection code #endif @@ -404,5 +452,5 @@ if (options->convolve) { options->psf = ppStackPSF(config, numCols, numRows, psfs, options); - psFree(psfs); + psFree(psfs); if (!options->psf) { #if 1