Index: branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotAddNoise.c (revision 34415) @@ -1,3 +1,5 @@ # include "psphotInternal.h" + +bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal); bool psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule) { @@ -30,4 +32,6 @@ } +static int Nmasked = 0; + bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) { @@ -56,4 +60,7 @@ psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels psAssert (maskVal, "missing mask value?"); + + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + psAssert (markVal, "missing mask value?"); // increase variance by factor*(object noise): @@ -94,4 +101,6 @@ pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0); + + psphotMaskSource (source, add, markVal); } if (add) { @@ -100,4 +109,5 @@ psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); } + fprintf (stderr, "masked %d objects\n", Nmasked); psphotVisualShowImage (readout); @@ -105,2 +115,34 @@ return true; } + +bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal) { + + if (!source) return false; + if (!source->peak) return false; // XXX how can we have a peak-less source? + if (source->type == PM_SOURCE_TYPE_DEFECT) return false; + if (source->type == PM_SOURCE_TYPE_SATURATED) return false; + + float Xc = source->peak->xf - source->pixels->col0 - 0.5; + float Yc = source->peak->yf - source->pixels->row0 - 0.5; + + psImageMaskType notMaskVal = ~maskVal; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + + float radius = hypot (ix - Xc, iy - Yc) ; + + if (radius > 4) continue; + + if (add) { + source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= maskVal; + } else { + source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] &= notMaskVal; + } + } + } + Nmasked ++; + + return true; +} + Index: branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotBlendFit.c (revision 34415) @@ -96,5 +96,5 @@ assert (status && isfinite(skySig) && skySig > 0); -# if (1) +# if (0) // **** test block : generate an ID image where pixels are set based the source models (flux > 0.1 peak && flux > 2.0 skySig) psImage *IDimage = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32); @@ -258,10 +258,10 @@ int TEST_ON = false; -# if (1) +# if (0) # define TEST_X 653 # define TEST_Y 466 if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { fprintf (stderr, "test object\n"); - // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); + psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); TEST_ON = true; } @@ -379,4 +379,15 @@ if (!source) return false; if (!source->peak) return false; // XXX how can we have a peak-less source? + +# if (0) +# define TEST_X 653 +# define TEST_Y 466 + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + fprintf (stderr, "test object\n"); + } +# undef TEST_X +# undef TEST_Y +# endif + if (!source->moments) return false; if (source->type == PM_SOURCE_TYPE_DEFECT) return false; @@ -473,9 +484,10 @@ // what is the radius for a specific peak fraction? - for (int i = 0; i < logFmodel->n; i++) { + for (int i = 1; i < logFmodel->n; i++) { float logF = logFmodel->data.F32[i]; + if (!isfinite(logF)) continue; + float flux = pow(10.0, logF); if (flux > threshold) continue; - if (i == 0) continue; float logF0 = logFmodel->data.F32[i - 1]; @@ -512,4 +524,6 @@ int ID = source->id; + fprintf (stderr, "Xo,Yo: %f,%f; radius: %f\n", Xo, Yo, radius); + for (int iy = minY; iy < maxY; iy++) { for (int ix = minX; ix < maxX; ix++) { Index: branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotDeblendSatstars.c (revision 34415) @@ -121,6 +121,6 @@ int BIG_SIGMA = BIG_RADIUS / 4.0; - int display = psphotKapaChannel (1); - psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0); + // int display = psphotKapaChannel (1); + // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0); // examine sources in decreasing SN order @@ -182,5 +182,5 @@ } // show the image after object have been subtracted - psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1); + // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1); psFree (SN); Index: branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotFindDetections.c (revision 34415) @@ -45,4 +45,10 @@ psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels psAssert (maskVal, "missing mask value?"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + psAssert (markVal, "missing mark value?"); + + maskVal |= markVal; // Use the new pmFootprints approach? Index: branches/eam_branches/ipp-20120905/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotFitSourcesLinear.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotFitSourcesLinear.c (revision 34415) @@ -439,4 +439,9 @@ model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; + + if (norm->data.F32[i] < MIN_VALID_FLUX) { + fprintf (stderr, "fit out of bounds for %f,%f : %f\n", source->peak->xf, source->peak->yf, norm->data.F32[i]); + model->params->data.F32[PM_PAR_I0] = MIN_VALID_FLUX; + } // clear the 'mark' pixels so the subtraction covers the full window Index: branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotKronIterate.c (revision 34415) @@ -321,4 +321,14 @@ if (!source->peak) continue; // XXX how can we have a peak-less source? +# if (0) +# define TEST_X 653 +# define TEST_Y 466 + if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) { + fprintf (stderr, "test object\n"); + } +# undef TEST_X +# undef TEST_Y +# endif + // check status of this source's moments if (!source->moments) continue; @@ -541,5 +551,5 @@ if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) { // We did not get a good measurement - source->moments->Mrf = NAN; + source->moments->Mrf = NAN; source->moments->KronFlux = NAN; source->moments->KronFluxErr = NAN; Index: branches/eam_branches/ipp-20120905/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotReadout.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotReadout.c (revision 34415) @@ -185,5 +185,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (config, view, filerule, false, false); // pass 1 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, false, true); // pass 1 (detections->allSources) // measure the radial profiles to the sky @@ -212,5 +212,5 @@ // linear fit to include all sources (subtract again) // NOTE : apply to ALL sources (extended + psf) - psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, true, true); // pass 2 (detections->allSources) // if we only do one pass, skip to extended source analysis @@ -260,5 +260,5 @@ // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, true, true); // pass 3 (detections->allSources) } @@ -300,5 +300,5 @@ // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, filerule, true, false); // pass 3 (detections->allSources) + psphotFitSourcesLinear (config, view, filerule, true, true); // pass 3 (detections->allSources) } Index: branches/eam_branches/ipp-20120905/psphot/src/psphotSourceSize.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotSourceSize.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotSourceSize.c (revision 34415) @@ -675,6 +675,11 @@ return false; } - psFree(readout->mask); - readout->mask = newMask; + // Copy the new mask pixel values to the old mask array. NOTE: we cannot replace + // the mask pointer because the source->maskView objects point to the old data + // area + psImage *oldMask = readout->mask; + readout->mask = psImageCopy (readout->mask, newMask, PS_TYPE_IMAGE_MASK); + psAssert (oldMask == readout->mask, "should have copied the values in-situ"); + psFree(newMask); } Index: branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c =================================================================== --- branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c (revision 34409) +++ branches/eam_branches/ipp-20120905/psphot/src/psphotStackReadout.c (revision 34415) @@ -332,12 +332,7 @@ psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects); - if (splitLinearFit) { - // NOTE: apply to Matched sources. Since the sources that we fit above are subtracted, they will - // not be included in this fit. - psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources) - } else { - // Fit all sources together - psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 3 (detections->allSources) - } + // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit + // the unsubtracted (matched) sources (the sources that we fit above are subtracted) + psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources) // measure the radial profiles to the sky (only measures new objects)