Index: trunk/ippconfig/recipes/psphot.config =================================================================== --- trunk/ippconfig/recipes/psphot.config (revision 34406) +++ trunk/ippconfig/recipes/psphot.config (revision 34418) @@ -75,4 +75,7 @@ DEBLEND_SKY_NSIGMA F32 10.0 DEBLEND_MIN_SATURATION F32 40000.0 + +# subtract radial profiles for extremely bright stars? +SUBTRACT_SATSTAR_PROFILE BOOL F # parameters to control the selection of the peak in the Sx,Sy plane @@ -414,4 +417,7 @@ PETROSIAN_FOR_STARS BOOL FALSE + # subtract radial profiles for extremely bright stars? + SUBTRACT_SATSTAR_PROFILE BOOL T + CONSTANT_PHOTOMETRIC_WEIGHTS BOOL TRUE POISSON.ERRORS.PHOT.LMM BOOL TRUE Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 34406) +++ trunk/psphot/src/psphot.h (revision 34418) @@ -81,5 +81,5 @@ bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); +bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view, const char *filerule); Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 34406) +++ trunk/psphot/src/psphotAddNoise.c (revision 34418) @@ -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: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 34406) +++ trunk/psphot/src/psphotBlendFit.c (revision 34418) @@ -1,4 +1,6 @@ # include "psphotInternal.h" bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma); +bool psphotBlendSetSourceSatstar (psImage *image, pmSource *source); +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); // for now, let's store the detections on the readout->analysis for each readout @@ -256,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; } @@ -377,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; @@ -385,5 +398,6 @@ if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) { // find the radius at which we hit something like 0.1*SATURATION - return false; + psphotBlendSetSourceSatstar (IDimage, source); + return true; } @@ -431,6 +445,6 @@ float rMxx = 0.5 / PS_SQR(shape.sx); float rMyy = 0.5 / PS_SQR(shape.sy); - float Sxy = -1. * shape.sxy; // factor of -1 is included to match the previous window function - // implementation. XXX: Is this correct? + float Sxy = shape.sxy; // factor of -1 is included to match the previous window function + // implementation. XXX: Is this correct? int ID = source->id; @@ -443,4 +457,5 @@ float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY; + if (z > 2.311) continue; float f = Io*exp(-z); @@ -453,2 +468,74 @@ return true; } + +bool psphotBlendSetSourceSatstar (psImage *IDimage, pmSource *source) { + + float Xo = source->satstar->Xo; + float Yo = source->satstar->Yo; + psVector *logRmodel = source->satstar->logRmodel; + psVector *logFmodel = source->satstar->logFmodel; + + float lPeak = logFmodel->data.F32[0]; + float fPeak = pow(10.0, lPeak); + float threshold = 0.1*fPeak; + float logThresh = log10(threshold); + + float radius = NAN; + + // what is the radius for a specific peak fraction? + 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; + + float logF0 = logFmodel->data.F32[i - 1]; + float logR0 = logRmodel->data.F32[i - 1]; + + float logF1 = logFmodel->data.F32[i]; + float logR1 = logRmodel->data.F32[i]; + + float logR = InterpolateValues (logF0, logR0, logF1, logR1, logThresh); + radius = pow(10.0, logR); + break; + } + + if (!isfinite(radius)) { + for (int i = logFmodel->n - 1; i >= 0; i--) { + if (!isfinite(logFmodel->data.F32[i])) continue; + float logR = logRmodel->data.F32[i]; + radius = pow(10.0, logR); + break; + } + } + + if (!isfinite(radius)) return false; + + int Nx = IDimage->numCols; + int Ny = IDimage->numRows; + + // region to mask + int minX = PS_MIN(PS_MAX(Xo - radius, 0), Nx - 1); + int maxX = PS_MIN(PS_MAX(Xo + radius, 0), Nx - 1); + int minY = PS_MIN(PS_MAX(Yo - radius, 0), Ny - 1); + int maxY = PS_MIN(PS_MAX(Yo + radius, 0), Ny - 1); + + 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++) { + + float dX = (ix - Xo); + float dY = (iy - Yo); + float R = hypot (dX, dY) ; + if (R > radius) continue; + + IDimage->data.S32[iy][ix] = ID; + } + } + return true; +} + Index: trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- trunk/psphot/src/psphotDeblendSatstars.c (revision 34406) +++ trunk/psphot/src/psphotDeblendSatstars.c (revision 34418) @@ -16,9 +16,21 @@ bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule) { + bool status = false; + int num = psphotFileruleCount(config, filerule); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full extended source non-linear fits? + if (!psMetadataLookupBool (&status, recipe, "SUBTRACT_SATSTAR_PROFILE")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n"); + return true; + } // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotDeblendSatstarsReadout (config, view, filerule, i)) { + if (!psphotDeblendSatstarsReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for %s entry %d", filerule, i); return false; @@ -40,15 +52,12 @@ 4) subtract the radial profile + There is also support code for calculation of magnitudes and add/subtract the profile (a la pmSourceOp) + TBD: - * function to replace the radial profile for a source * recenter the profile (based on cross-correlation / convolution with 1D profile) - - * raise a bit somewhere for these super saturated stars (if they were so modeled) - * consider the subtraction bit : when to raise it? - **/ -bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { +bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe) { int N; @@ -56,7 +65,4 @@ bool status; - // XXX disable this function for now - return true; - psTimerStart ("psphot.deblend.sat"); @@ -78,8 +84,4 @@ return true; } - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -124,6 +126,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 @@ -185,5 +187,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); @@ -451,4 +453,6 @@ // XXX do something sensible here if the profile is crap + // replace an existing profile + psFree (source->satstar); source->satstar = pmSourceSatstarAlloc(); source->satstar->Xo = Xo; Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 34406) +++ trunk/psphot/src/psphotFindDetections.c (revision 34418) @@ -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: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 34406) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 34418) @@ -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: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 34406) +++ trunk/psphot/src/psphotKronIterate.c (revision 34418) @@ -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: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 34406) +++ trunk/psphot/src/psphotReadout.c (revision 34418) @@ -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: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 34406) +++ trunk/psphot/src/psphotSourceSize.c (revision 34418) @@ -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: trunk/psphot/src/psphotStackReadout.c =================================================================== --- trunk/psphot/src/psphotStackReadout.c (revision 34406) +++ trunk/psphot/src/psphotStackReadout.c (revision 34418) @@ -169,4 +169,22 @@ } + // If DET and SRC are different images, subtract radial profiles for the convolved + // image first. The profiles found on the convolved image will be replaced by those + // found for the unconvolved image. Downstream, in psphotFindDetections, we will + // replace the the profiles (and re-subtract them) for the detection image, so we want + // to keep those versions of the profiles on the sources. + if (strcmp(STACK_SRC, STACK_DET)) { + // find and subtract radial profile models for saturated stars (XXX change name eventually) + if (!psphotDeblendSatstars (config, view, STACK_SRC)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + } + // find and subtract radial profile models for saturated stars (XXX change name eventually) + if (!psphotDeblendSatstars (config, view, STACK_DET)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) // only run this on detections from the input images, not chisq image @@ -206,8 +224,8 @@ // window of size PSF_MOMENTS_RADIUS (same window used to measure the psf-scale moments) // but iterates to an appropriately larger size - logMemStats("before.kron.1"); + logMemStats("before.kron.1"); psphotKronIterate(config, view, STACK_SRC, 1); - logMemStats("after.kron.1"); - + logMemStats("after.kron.1"); + // identify CRs and extended sources psphotSourceSize (config, view, STACK_SRC, true); @@ -332,12 +350,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)