Index: /branches/eam_branches/ipp-20110213/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphot.h (revision 30902) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphot.h (revision 30903) @@ -173,8 +173,8 @@ // used by psphotFindDetections -psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal); +pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal); psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); -bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); -psErrorCode psphotCullPeaks(const pmReadout *readout, const psMetadata *recipe, psArray *footprints); +bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); +psErrorCode psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints); // in psphotApResid.c: Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c (revision 30902) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c (revision 30903) @@ -6,4 +6,5 @@ psErrorCode psphotCullPeaks(const pmReadout *readout, // the image wherein lives the footprint + const pmReadout *signifR, // smoothed image and significance const psMetadata *recipe, psArray *footprints) { // array of pmFootprints @@ -26,5 +27,9 @@ psAssert (status, "cannot find SKY_STDEV in readout->analysis"); - const float MIN_THRESHOLD = nsigma_min*skyStdev; + // the sky has been smoothed, so we need to scale down the raw sky stdev by this amount: + float scaleFactor = psMetadataLookupF32(&status, readout->analysis, "SIGNIFICANCE_SCALE_FACTOR"); + if (!status) scaleFactor = 1.0; + + const float MIN_THRESHOLD = nsigma_min*skyStdev / sqrt(scaleFactor); // for saturated stars, we should be somewhat more agressive about culling: instead of @@ -36,6 +41,7 @@ // fStdev = 0.005, stdev_pad = 0.011*40k = 400 // threshold = 40k - 4*400 = 38400 - - // in practice, we should make the threshold in such a case more like 0.5 * saturation... + // this gives too tight of a tolerance on the bright stars + // in practice, we should make the threshold much lower. + // below I am using 0.05 * saturation (eg, 2000 DN above sky in a GPC1 image) float SATURATION = NAN; @@ -72,7 +78,21 @@ } - if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold) != PS_ERR_NONE) { + // if we cull using the significance image, then the thresholds are different +# if (0) + if (pmFootprintCullPeaks(signifR->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) { return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); } +# endif +# if (1) + if (pmFootprintCullPeaks(signifR->image, signifR->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, false) != PS_ERR_NONE) { + return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); + } +# endif +# if (0) + if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) { + return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); + } +# endif + float dtime = psTimerMark ("psphot.cull.footprints"); if (dtime > 1.0) { Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotFindDetections.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotFindDetections.c (revision 30902) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotFindDetections.c (revision 30903) @@ -84,14 +84,14 @@ // generate the smoothed significance image - psImage *significance = psphotSignificanceImage (readout, recipe, maskVal); + pmReadout *significance = psphotSignificanceImage (readout, recipe, maskVal); // display the log significance image - psphotVisualShowLogSignificance (significance, 0.0, 4.5); + psphotVisualShowLogSignificance (significance->variance, 0.0, 4.5); // display the significance image - psphotVisualShowSignificance (significance, 0.98*threshold, 1.02*threshold); + psphotVisualShowSignificance (significance->variance, 0.98*threshold, 1.02*threshold); // detect the peaks in the significance image - detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX); + detections->peaks = psphotFindPeaks (significance->variance, readout, recipe, threshold, NMAX); psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); if (!detections->peaks) { Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotFindFootprints.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotFindFootprints.c (revision 30902) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotFindFootprints.c (revision 30903) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) { +bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) { bool status; @@ -18,6 +18,6 @@ PS_ASSERT (status, false); - // find the raw footprints & assign the peaks to those footprints - psArray *footprints = pmFootprintsFind (significance, threshold, npixMin); + // find the raw footprints in the smoothed significance image & assign the peaks to those footprints + psArray *footprints = pmFootprintsFind (significance->variance, threshold, npixMin); if (pmFootprintsAssignPeaks(footprints, detections->peaks) != PS_ERR_NONE) { @@ -56,5 +56,5 @@ } - psphotCullPeaks(readout, recipe, detections->footprints); + psphotCullPeaks(readout, significance, recipe, detections->footprints); detections->peaks = pmFootprintArrayToPeaks(detections->footprints); psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints")); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotSignificanceImage.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotSignificanceImage.c (revision 30902) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotSignificanceImage.c (revision 30903) @@ -4,5 +4,5 @@ // (S/N)^2. If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the // smoothing kernel. -psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) { +pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) { float SIGMA_SMTH, NSIGMA_SMTH; @@ -94,14 +94,18 @@ psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); - // XXX multithread this if needed + // we are going to return both the image and the weight here: the image contains the signal + // while the 'weight' will contain the significance (NOTE the deviation from the usual + // definition) + + // save the smoothed significance image in the weight array for (int j = 0; j < smooth_im->numRows; j++) { for (int i = 0; i < smooth_im->numCols; i++) { float value = smooth_im->data.F32[j][i]; if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)) { - smooth_im->data.F32[j][i] = 0.0; + smooth_wt->data.F32[j][i] = 0.0; } else { // XXX the value of 100 here (or 1000 before) must depend on the FWHM of the smoothing kernel, right?? float v2 = value + PS_SQR(value/100.0); - smooth_im->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i]; + smooth_wt->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i]; } } @@ -115,13 +119,14 @@ static int pass = 0; sprintf (name, "snsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_im, name); + psphotSaveImage (NULL, smooth_wt, name); pass ++; } - - psFree(smooth_wt); - psImageConvolveSetThreads(oldThreads); - return smooth_im; + pmReadout *significanceRO = pmReadoutAlloc(NULL); + significanceRO->variance = smooth_wt; + significanceRO->image = smooth_im; + + return significanceRO; }