Index: trunk/psphot/doc/notes.20100131.txt =================================================================== --- trunk/psphot/doc/notes.20100131.txt (revision 27314) +++ trunk/psphot/doc/notes.20100131.txt (revision 27532) @@ -1,2 +1,46 @@ + +2010.03.24 + + I have finished another iteration on source size analysis. I feel + fairly confident of the extended source and CR identifications. They + are not perfect, but they are in pretty good shape. The CR IDs use + the moments to find things smaller in one dimension than a PSF, while + the EXT is set for things larger than a PSF in one dimention. + + There still seems to be some trouble with SATSTARS : these are being + set for sources which are fainter than the SAT limit -- I'm not sure + if this is due to the quality of the PSF fit or if they really are + SAT, but have integrated fluxes lower than expected (underfitted in + PSF). + + Additional Parameters needed from psphot for distinguishing dipoles: + + -- note regarding CfA parameters -- + + Here are our settings. + + Note that Fpos is the flux which is 3-sigma above background, and + Fneg is the flux which is 3-sigma below background. + + Npos is the # of pixels with flux which is 3-sigma above background. + + etc. + + # for positive objects: + # Ngood = Npos, Nbad = Nneg, Fgood = Fpos, Fbad = Fneg + + # for negative objects: + # Ngood = Nneg, Nbad = Npos, Fgood = Fneg, Fbad = Fpos + # FRATIO=Fgood/(Fgood+Fbad) + # NRATIO_BAD=Ngood/(Ngood+Nbad) + # NRATIO_MASK=Ngood/(Ngood+Nmask) + # NRATIO_ALL=Ngood/(Ngood+Nbad+Nmask) + DC_FRATIO_MIN 0.75 + DC_NRATIO_BAD_MIN 0.70 + DC_NRATIO_MASK_MIN 0.60 + DC_NRATIO_ALL_MIN 0.50 + DC_NGOOD_MIN 3 + + *** double-check on the use of fluxScale in pmSourcePhotometry.c:106 20100131 Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 27314) +++ trunk/psphot/src/psphot.h (revision 27532) @@ -92,5 +92,5 @@ bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final); -bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final); +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final); bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize); @@ -119,5 +119,5 @@ bool psphotMagnitudes (pmConfig *config, const pmFPAview *view); -bool psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf); bool psphotMagnitudes_Threaded (psThreadJob *job); Index: trunk/psphot/src/psphotEfficiency.c =================================================================== --- trunk/psphot/src/psphotEfficiency.c (revision 27314) +++ trunk/psphot/src/psphotEfficiency.c (revision 27532) @@ -4,8 +4,5 @@ PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models -//#define TESTING - - -# if 0 +# define TESTING 0 // Calculate the limiting magnitude for an image @@ -19,5 +16,5 @@ float *covarFactor,// Covariance factor const pmReadout *ro, // Readout of interest - const pmPSF *psf, // Point-spread function + pmPSF *psf, // Point-spread function float thresh, // Threshold for source identification float smoothSigma, // Gaussian smoothing sigma @@ -152,6 +149,4 @@ } -# endif - bool psphotEfficiency (pmConfig *config, const pmFPAview *view) { @@ -178,5 +173,4 @@ bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { -# if 0 bool status = true; @@ -206,34 +200,26 @@ // Collect recipe information - float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x - float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y + float smoothNsigma = psMetadataLookupF32(&status, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit + psAssert (status && isfinite(smoothNsigma), "Unable to find PEAKS_SMOOTH_NSIGMA in recipe (or invalid value)"); + + float thresh = psMetadataLookupF32(&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); + psAssert (status && isfinite(thresh), "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe (or invalid value)"); + + psVector *magOffsets = psMetadataLookupVector(&status, recipe, "EFF.MAG"); // Magnitude offsets + psAssert (status, "Unable to find EFF.MAG F32 vector in recipe"); + psAssert (magOffsets->type.type == PS_TYPE_F32, "Unable to find EFF.MAG F32 vector in recipe"); + + int numSources = psMetadataLookupS32(&status, recipe, "EFF.NUM"); // Number of sources for each bin + psAssert (status && (numSources > 0), "Unable to find EFF.NUM in recipe (or invalid value)"); + + float minGauss = psMetadataLookupF32(&status, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel + psAssert (status && isfinite(minGauss), "PEAKS_MIN_GAUSS is not set in recipe (or invalid)"); + + // find the PSF size information (why is this not part of the psf structure?) + float fwhmMajor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MAJ"); // PSF size in x + float fwhmMinor = psMetadataLookupF32(NULL, readout->analysis, "FWHM_MIN"); // PSF size in y if (!isfinite(fwhmMajor) || !isfinite(fwhmMinor) || fwhmMajor == 0.0 || fwhmMinor == 0.0) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in recipe"); - return false; - } - float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA"); // Smoothing limit - if (!isfinite(smoothNsigma)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_SMOOTH_NSIGMA in recipe"); - return false; - } - float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2"); - if (!isfinite(thresh)) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find PEAKS_NSIGMA_LIMIT_2 in recipe"); - return false; - } - psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "EFF.MAG"); // Magnitude offsets - if (!magOffsets || magOffsets->type.type != PS_TYPE_F32) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.MAG F32 vector in recipe"); - return NULL; - } - int numSources = psMetadataLookupS32(NULL, recipe, "EFF.NUM"); // Number of sources for each bin - if (numSources == 0) { - psError(PSPHOT_ERR_CONFIG, false, "Unable to find EFF.NUM in recipe"); - return NULL; - } - float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel - if (!isfinite(minGauss)) { - psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value"); - minGauss = 0.5; + psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_MAJ and FWHM_MIN in readout->analysis"); + return false; } @@ -246,8 +232,6 @@ // remove all sources, adding noise for subtracted sources psphotRemoveAllSources(realSources, recipe); - // psphotAddNoise(readout, realSources, recipe); - - -#if defined(TESTING) && 0 + +#if TESTING { psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); @@ -278,5 +262,5 @@ } -#ifdef TESTING +#if TESTING psphotSaveImage(NULL, readout->image, "orig_image.fits"); psphotSaveImage(NULL, readout->variance, "orig_variance.fits"); @@ -293,5 +277,5 @@ } -#ifdef TESTING +#if TESTING psphotSaveImage(NULL, readout->image, "fake_image.fits"); psphotSaveImage(NULL, readout->variance, "fake_variance.fits"); @@ -408,5 +392,5 @@ psFree(significance); - if (!psphotFitSourcesLinear(readout, fakeSourcesAll, recipe, psf, true)) { + if (!psphotFitSourcesLinearReadout(recipe, readout, fakeSourcesAll, psf, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform linear fit on fake sources."); psFree(fakeSources); @@ -415,19 +399,19 @@ } - // Disable aperture corrections; casting away const! + // Disable aperture corrections (save current value) pmTrend2D *apTrend = psf->ApTrend; // Aperture trend - ((pmPSF*)psf)->ApTrend = NULL; - - if (!psphotMagnitudes(config, readout, view, fakeSourcesAll, psf)) { + psf->ApTrend = NULL; + + if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources."); psFree(fakeSources); psFree(count); - ((pmPSF*)psf)->ApTrend = apTrend; // Casting away const! + psf->ApTrend = apTrend; // Casting away const! return false; } psFree(fakeSourcesAll); - // Re-enable aperture corrections; casting away const! - ((pmPSF*)psf)->ApTrend = apTrend; + // Replace aperture corrections + psf->ApTrend = apTrend; psVector *magDiffMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean difference in magnitude for each bin @@ -443,5 +427,5 @@ psVectorInit(magMask, 0); -#ifdef TESTING +#if TESTING psString name = NULL; psStringAppend(&name, "fake_%d.dat", i); @@ -459,5 +443,5 @@ } -#ifdef TESTING +#if TESTING fprintf(file, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf, source->modelPSF->params->data.F32[PM_PAR_XPOS], @@ -495,5 +479,5 @@ } -#ifdef TESTING +#if TESTING fclose(file); #endif @@ -520,6 +504,4 @@ psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake")); -# endif - return true; } Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 27314) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 27532) @@ -17,4 +17,8 @@ bool status = true; + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); @@ -22,5 +26,22 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotFitSourcesLinearReadout (config, view, "PSPHOT.INPUT", i, final)) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i); return false; @@ -30,5 +51,5 @@ } -bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) { +bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) { bool status; @@ -36,22 +57,4 @@ float y; float f; - // float r; - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); if (!sources->n) { @@ -59,7 +62,4 @@ return true; } - - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - psAssert (sources, "missing psf?"); psTimerStart ("psphot.linear"); Index: trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27314) +++ trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27532) @@ -9,6 +9,4 @@ // these are used to determine the simultaneous linear fit of fluxes. // the analysis is performed wrt the simulated pixel values - -static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal); bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { @@ -194,58 +192,2 @@ return true; } - -// Calculate the weight terms for the sky fit component of the matrix. This function operates -// on the pixels which correspond to all of the sources of interest. These elements fill in -// the border matrix components in the sparse matrix equation. -static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal) { - - // generate the image-wide weight terms - // turn on MARK for all image pixels - psRegion fullArray = psRegionSet (0, 0, 0, 0); - fullArray = psRegionForImage (readout->mask, fullArray); - psImageMaskRegion (readout->mask, fullArray, "OR", markVal); - - // turn off MARK for all object pixels - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; - float x = model->params->data.F32[PM_PAR_XPOS]; - float y = model->params->data.F32[PM_PAR_YPOS]; - psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_IMAGE_MASK(markVal)); - } - - // accumulate the image statistics from the masked regions - psF32 **image = readout->image->data.F32; - psF32 **variance = readout->variance->data.F32; - psImageMaskType **mask = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; - - double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; - w = x = y = x2 = xy = y2 = fo = fx = fy = 0; - - int col0 = readout->image->col0; - int row0 = readout->image->row0; - - for (int j = 0; j < readout->image->numRows; j++) { - for (int i = 0; i < readout->image->numCols; i++) { - if (mask[j][i]) continue; - if (constant_weights) { - wt = 1.0; - } else { - wt = variance[j][i]; - } - f = image[j][i]; - w += 1/wt; - fo += f/wt; - } - } - - // turn off MARK for all image pixels - psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // set the Border T elements - psSparseBorderElementG (border, 0, fo); - psSparseBorderElementT (border, 0, 0, w); - - return true; -} Index: trunk/psphot/src/psphotKernelFromPSF.c =================================================================== --- trunk/psphot/src/psphotKernelFromPSF.c (revision 27314) +++ trunk/psphot/src/psphotKernelFromPSF.c (revision 27532) @@ -4,8 +4,8 @@ assert (source); - assert (source->psfFlux); // XXX build if needed? + assert (source->psfImage); // XXX build if needed? - int x0 = source->peak->xf - source->psfFlux->col0; - int y0 = source->peak->yf - source->psfFlux->row0; + int x0 = source->peak->xf - source->psfImage->col0; + int y0 = source->peak->yf - source->psfImage->row0; // need to decide on the size: dynamically? statically? @@ -17,12 +17,12 @@ // if the realized PSF for this object does not cover the full kernel, give up for now if (x0 + psf->xMin < 0) goto escape; - if (x0 + psf->xMax >= source->psfFlux->numCols) goto escape; + if (x0 + psf->xMax >= source->psfImage->numCols) goto escape; if (y0 + psf->yMin < 0) goto escape; - if (y0 + psf->yMax >= source->psfFlux->numRows) goto escape; + if (y0 + psf->yMax >= source->psfImage->numRows) goto escape; double sum = 0.0; for (int j = psf->yMin; j <= psf->yMax; j++) { for (int i = psf->xMin; i <= psf->xMax; i++) { - double value = source->psfFlux->data.F32[y0 + j][x0 + i]; + double value = source->psfImage->data.F32[y0 + j][x0 + i]; psf->kernel[j][i] = value; sum += value; Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 27314) +++ trunk/psphot/src/psphotMagnitudes.c (revision 27532) @@ -14,5 +14,22 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotMagnitudesReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + + if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i); return false; @@ -22,31 +39,15 @@ } -bool psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) { bool status = false; int Nap = 0; + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes"); + return true; + } + psTimerStart ("psphot.mags"); - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); - - if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); - return true; - } - - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - psAssert (psf, "missing psf?"); // determine the number of allowed threads @@ -79,10 +80,12 @@ } - bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); + bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); + bool DIFF_STATS = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT; if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH; if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP; + if (DIFF_STATS) photMode |= PM_SOURCE_PHOT_DIFFSTATS; // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) Index: trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- trunk/psphot/src/psphotMakeResiduals.c (revision 27314) +++ trunk/psphot/src/psphotMakeResiduals.c (revision 27532) @@ -305,4 +305,7 @@ if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) { resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; + resid->Ro->data.F32[oy][ox] = 0.0; + resid->Rx->data.F32[oy][ox] = 0.0; + resid->Ry->data.F32[oy][ox] = 0.0; } } Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 27314) +++ trunk/psphot/src/psphotReadout.c (revision 27532) @@ -91,11 +91,10 @@ return psphotReadoutCleanup (config, view); } - if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup(config, view); - } - // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) if (!psphotImageQuality (config, view)) { // pass 1 psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); + return psphotReadoutCleanup(config, view); + } + if (!strcasecmp (breakPt, "MOMENTS")) { return psphotReadoutCleanup(config, view); } Index: trunk/psphot/src/psphotSourceMerge.c =================================================================== --- trunk/psphot/src/psphotSourceMerge.c (revision 27532) +++ trunk/psphot/src/psphotSourceMerge.c (revision 27532) @@ -0,0 +1,49 @@ + +// we need a couple of functions to distinguish coincident sources: +// XXX these are similar (identical?) to the goals of pmSourceMatch.c + +# define NEXT1 { i++; continue; } +# define NEXT2 { j++; continue; } + +psphotSourceMerge () { + + float dx, dy; + + // sort the source list by X + pmSourceSortByX (sources1); + pmSourceSortByX (sources2); + + int i, j; + for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) { + + pmSource *src1 = sources1->data[i]; + pmSource *src2 = sources2->data[j]; + + if (!src1) NEXT1; + if (!src1->peak) NEXT1; + if (!finite(src1->peak->xf)) NEXT1; + if (!finite(src1->peak->yf)) NEXT1; + + if (!src2) NEXT2; + if (!src2->peak) NEXT2; + if (!finite(src2->peak->xf)) NEXT2; + if (!finite(src2->peak->yf)) NEXT2; + + dx = src1->peak->xf - src2->peak->xf; + if (dx < -1.02*RADIUS) NEXT1; + if (dx > +1.02*RADIUS) NEXT2; + + // we are within match range, look for matches: + for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { + + dx = src1->peak->xf - src2->peak->xf; + dy = src1->peak->yf - src2->peak->yf; + + dr = dx*dx + dy*dy; + if (dr > RADIUS2) continue; + + // add to group? + } + i++; + } +} Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 27314) +++ trunk/psphot/src/psphotSourceSize.c (revision 27532) @@ -13,4 +13,6 @@ float soft; int grow; + int xtest, ytest; + bool apply; // apply CR mask? } psphotSourceSizeOptions; @@ -22,4 +24,5 @@ bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal); bool psphotMaskCosmicRayFootprintCheck (psArray *sources); +int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh); // we need to call this function after sources have been fitted to the PSF model and @@ -101,4 +104,9 @@ assert (status); + // XXX recipe name is not great + options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST"); + options.ytest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.YTEST"); + assert (status); + options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs if (!status || options.grow < 0) { @@ -111,4 +119,10 @@ psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero."); options.soft = 0.0; + } + + options.apply = psMetadataLookupBool(&status, recipe, "PSPHOT.CRMASK.APPLY"); // Growth size for CRs + if (!status) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CRMASK.APPLY is not defined."); + return false; } @@ -239,4 +253,5 @@ continue; } + // psphotVisualPlotSourceSize (recipe, readout->analysis, sources); } @@ -244,4 +259,5 @@ } +# define SIZE_SN_LIM 10 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { @@ -285,5 +301,5 @@ } - // we are basically classifying by moments; use the default if not found + // we are basically classifying by moments psAssert (source->moments, "why is this source missing moments?"); if (source->mode & noMoments) { @@ -292,4 +308,5 @@ } + // convert to Mmaj, Mmin: psF32 Mxx = source->moments->Mxx; psF32 Myy = source->moments->Myy; @@ -315,18 +332,36 @@ float apMag = -2.5*log10(source->moments->Sum); float dMag = source->psfMag - apMag; - float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); - - source->extNsigma = nSigma; - source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + // set nSigma to include both systematic and poisson error terms + // XXX the 'poisson error' contribution for size is probably wrong... + float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); + float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag); + float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag); + + // partially-masked sources are more likely to be mis-measured PSFs + float sizeBias = 1.0; + if (source->pixWeight < 0.9) { + sizeBias = 3.0; + } + + float minMxx = psfClump->X - sizeBias*options->nSigmaMoments*psfClump->dX; + float minMyy = psfClump->Y - sizeBias*options->nSigmaMoments*psfClump->dY; + + // include MAG, MXX, and MYY? + source->extNsigma = nSigmaMAG; + + // notes to clarify the source size classification rules: + // * a defect should be functionally equivalent to a cosmic ray + // * CR & defect should have a faintess limit (min S/N) + // * SAT stars should not be faint, but defects may? // Anything within this region is a probably PSF-like object. Saturated stars may land // in this region, but are detected elsewhere on the basis of their peak value. - bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) && - (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && - (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY)); + bool isPSF = (fabs(nSigmaMAG) < options->nSigmaApResid) && (fabs(nSigmaMXX) < sizeBias*options->nSigmaMoments) && (fabs(nSigmaMYY) < sizeBias*options->nSigmaMoments); if (isPSF) { - psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g PSF\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, - options->nSigmaApResid,options->nSigmaMoments); + psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g PSF\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, + options->nSigmaApResid,sizeBias*options->nSigmaMoments); + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; Npsf ++; continue; @@ -336,10 +371,13 @@ // Defects may also be marked as SATSTAR -- XXX deactivate this flag? // XXX this rule is not great - if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { - - psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g CR\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, - options->nSigmaApResid,options->nSigmaMoments); + // XXX only accept brightish detections as CRs + // (nSigmaMAG < -options->nSigmaApResid) || + bool isCR = isCR = (source->errMag < 1.0 / SIZE_SN_LIM) && ((Mxx < minMxx) || (Myy < minMyy)); + if (isCR) { + psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g %g %g %g %g\t%g %g\t%g CR\t%g %g\n", + source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, + options->nSigmaApResid,sizeBias*options->nSigmaMoments); source->mode |= PM_SOURCE_MODE_DEFECT; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_CR_CANDIDATE; Ncr ++; continue; @@ -349,8 +387,8 @@ // just large saturated regions. if (source->mode & PM_SOURCE_MODE_SATSTAR) { - psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g SAT\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, - options->nSigmaApResid,options->nSigmaMoments); - + psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g SAT\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, + options->nSigmaApResid,sizeBias*options->nSigmaMoments); + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; Nsat ++; continue; @@ -358,19 +396,22 @@ // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? - bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y)); + bool isEXT = (nSigmaMAG > options->nSigmaApResid) || (Mxx > minMxx) || (Myy > minMyy); if (isEXT) { - psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Ext\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, - options->nSigmaApResid,options->nSigmaMoments); + psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Ext\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, + options->nSigmaApResid,sizeBias*options->nSigmaMoments); source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; Next ++; continue; } - psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Unk\t%g %g\n", - source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma, - options->nSigmaApResid,options->nSigmaMoments); + psTrace("psphotSourceClassRegion.MISS",4,"CLASS: %g %g\t%g %g %g %g %g %g\t%g %g\t%g Unk\t%g %g\n", + source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, + options->nSigmaApResid,sizeBias*options->nSigmaMoments); - psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); + // sources that reach here are probably too faint for a reasonable source size measurement + // psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigmaMAG); + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; Nmiss ++; } @@ -385,4 +426,7 @@ bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { + psTimerStart ("psphot.cr"); + + int nMasked = 0; for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; @@ -393,4 +437,14 @@ continue; } + + // only check candidates marked above + if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_CR_CANDIDATE)) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // skip unless this source is thought to be a cosmic ray. flag the detection and mask the pixels + // XXX this may be degenerate with the above test + if (!(source->mode & PM_SOURCE_MODE_DEFECT)) continue; // Integer position of peak @@ -402,17 +456,27 @@ yPeak < 1 || yPeak > source->pixels->numRows - 2) { psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); - // psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n"); continue; } - // this source is thought to be a cosmic ray. flag the detection and mask the pixels - if (source->mode & PM_SOURCE_MODE_DEFECT) { - // XXX this is running slowly and is too agressive, but it more-or-less works - // psphotMaskCosmicRayCZW(readout, source, options->crMask); - - // XXX these are the old versions which we are not using - // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); - // psphotMaskCosmicRay_Old (source, maskVal, crMask); - } + // XXX for testing, only CRMASK a single source: + if (options->xtest && (fabs(source->peak->xf - options->xtest) > 5)) continue; + if (options->ytest && (fabs(source->peak->yf - options->ytest) > 5)) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); + } + + // XXX this is running slowly and is too agressive, but it more-or-less works + psTrace("psphot", 6, "mask cosmic ray at %f, %f\n", source->peak->xf, source->peak->yf); + if (options->apply) { + psphotMaskCosmicRay(readout, source, options->crMask); + } else { + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + } + nMasked ++; + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); } @@ -430,10 +494,17 @@ } + psLogMsg ("psphot.cr", PS_LOG_INFO, "mask CR: %d masked in %f sec\n", nMasked, psTimerMark ("psphot.cr")); + // XXX test : save the mask image if (0) { psphotSaveImage (NULL, readout->mask, "mask.fits"); } + return true; } + +# define DUMPPICS 0 +# define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); } +# define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); } // Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray @@ -448,32 +519,28 @@ pmFootprint *footprint = peak->footprint; - int xm = footprint->bbox.x0; - int xM = footprint->bbox.x1; - int ym = footprint->bbox.y0; - int yM = footprint->bbox.y1; - - if (xm < 0) { xm = 0; } - if (ym < 0) { ym = 0; } - if (xM > mask->numCols) { xM = mask->numCols; } - if (yM > mask->numRows) { yM = mask->numRows; } - int dx = xM - xm; - int dy = yM - ym; - // Bounding boxes are inclusive of final pixel - // XXX ensure dx,dy do not become too large here - dx++; - dy++; + int xs = footprint->bbox.x0; + int xe = footprint->bbox.x1 + 1; + int ys = footprint->bbox.y0; + int ye = footprint->bbox.y1 + 1; + + LIMIT_XRANGE(xs, mask); + LIMIT_XRANGE(xe, mask); + LIMIT_YRANGE(ys, mask); + LIMIT_YRANGE(ye, mask); + + int dx = xe - xs; + int dy = ye - ys; psImage *image= readout->image; psImage *variance = readout->variance; - int binning = 1; - float sigma_thresh = 2.0; - int iteration = 0; - int max_iter = 5; - float noise_factor = 5.0 / 4.0; // Intrinsic to the Laplacian making noise spikes spikier. + int binning = 2; + float sigma_thresh = 3.0; + int max_iter = 1; // XXX with isophot masking, we only want to do a single pass // Temporary images. psImage *mypix = psImageAlloc(dx,dy,image->type.type); + psImage *myfix = psImageAlloc(dx,dy,image->type.type); psImage *myvar = psImageAlloc(dx,dy,image->type.type); psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type); @@ -482,10 +549,9 @@ psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK); - int x,y; // Load my copy of things. - for (x = 0; x < dx; x++) { - for (y = 0; y < dy; y++) { - psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym)); - psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym)); + for (int y = 0; y < dy; y++) { + for (int x = 0; x < dx; x++) { + mypix->data.F32[y][x] = image->data.F32[y+ys][x+xs]; + myvar->data.F32[y][x] = variance->data.F32[y+ys][x+xs]; mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00; } @@ -495,111 +561,123 @@ pmSpan *sp = footprint->spans->data[i]; for (int j = sp->x0; j <= sp->x1; j++) { - y = sp->y - ym; - x = j - xm; + int y = sp->y - ys; + int x = j - xs; mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01; } } - int CRpix_count = 0; - do { - CRpix_count = 0; - // Zero out my temp images. - for (x = 0; x < binning * dx; x++) { - for (y = 0; y < binning * dy; y++) { - psImageSet(binned,x,y,0.0); - psImageSet(conved,x,y,0.0); - psImageSet(edges,x/binning,y/ binning,0.0); + int nCRpix = 1; // force at least one pass... + for (int iteration = 0; (iteration < max_iter) && (nCRpix > 0); iteration++) { + nCRpix = 0; + psImageInit (binned, 0.0); + psImageInit (conved, 0.0); + psImageInit (edges, 0.0); + + // Make subsampled image. Maybe this should be called "unbinned" or something + for (int y = 0; y < binning * dy; y++) { + int yraw = y / binning; + for (int x = 0; x < binning * dx; x++) { + int xraw = x / binning; + binned->data.F32[y][x] = mypix->data.F32[yraw][xraw]; } - } - // Make subsampled image. Maybe this should be called "unbinned" or something - for (x = 0; x < binning * dx; x++) { - for (y = 0; y < binning * dy; y++) { - psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning)); + } + + // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero + for (int y = 1; y < binning * dy - 1; y++) { + for (int x = 1; x < binning * dx - 1; x++) { + float value = binned->data.F32[y][x] - 0.25 * + (binned->data.F32[y+0][x-1] + binned->data.F32[y+0][x+1] + + binned->data.F32[y-1][x+0] + binned->data.F32[y+1][x+0]); + value = PS_MAX(0.0, value); + + conved->data.F32[y][x] = value; } } - // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero - for (x = 1; x < dx - 1; x++) { - for (y = 1; y < dy - 1; y++) { - psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 * - (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) + - psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1))); - if (psImageGet(conved,x,y) < 0.0) { - psImageSet(conved,x,y,0.0); - } + + // Create an edge map by rebinning + for (int y = 0; y < binning * dy; y++) { + int yraw = y / binning; + for (int x = 0; x < binning * dx; x++) { + int xraw = x / binning; + edges->data.F32[yraw][xraw] += conved->data.F32[y][x]; } } - // Create an edge map by rebinning - for (x = 0; x < binning * dx; x++) { - for (y = 0; y < binning * dy; y++) { - psImageSet(edges,x / binning, y / binning, - psImageGet(edges, x / binning, y / binning) + - psImageGet(conved,x,y)); - } - } - // Modify my mask if we're above the significance threshold - for (x = 0; x < dx; x++) { - for (y = 0; y < dy; y++) { - if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) { - if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) { - mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40; - CRpix_count++; - } - } - } - } - + + // coordinate of peak in subimage pixels: + int xPeak = peak->x - xs; + int yPeak = peak->y - ys; + + // Modify my mask if we're above the significance threshold, but only for connected pixels + nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh); + +# if DUMPPICS + psphotSaveImage (NULL, mypix, "crmask.pix.fits"); +# endif + +// XXX do not repair the pixels in isophot version +# if 0 // "Repair" Masked pixels for the next round. - for (x = 1; x < dx - 1; x++) { - for (y = 1; y < dy - 1; y++) { - if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) { - psImageSet(mypix,x,y,0.25 * - (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) + - psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1))); - } - } - } - -# if 0 - if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) { - psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration); - psphotSaveImage (NULL, mypix, "czw.pix.fits"); - psphotSaveImage (NULL, myvar, "czw.var.fits"); - psphotSaveImage (NULL, binned, "czw.binn.fits"); - psphotSaveImage (NULL, conved, "czw.conv.fits"); - psphotSaveImage (NULL, edges, "czw.edge.fits"); - psphotSaveImage (NULL, mymask, "czw.mask.fits"); - } -# endif - - psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count); - iteration++; - } while ((iteration < max_iter)&&(CRpix_count > 0)); - - // A solitary masked pixel is likely a lie. Remove those. - for (x = 0; x < dx; x++) { - for (y = 0; y < dy; y++) { - if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) { - if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) { + for (int y = 1; y < dy - 1; y++) { + for (int x = 1; x < dx - 1; x++) { + if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) { + myfix->data.F32[y][x] = mypix->data.F32[y][x]; continue; } - if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) { - continue; - } - if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) { - continue; - } - if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) { - continue; - } - mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40; + myfix->data.F32[y][x] = 0.25 * + (mypix->data.F32[y+0][x-1] + mypix->data.F32[y+0][x+1] + + mypix->data.F32[y-1][x+0] + mypix->data.F32[y+1][x+0]); } } - } - - // transfer temporary mask to real mask - for (x = 0; x < dx; x++) { - for (y = 0; y < dy; y++) { + + // "Repair" Masked pixels for the next round. + for (int y = 1; y < dy - 1; y++) { + for (int x = 1; x < dx - 1; x++) { + mypix->data.F32[y][x] = myfix->data.F32[y][x]; + } + } +# endif + +# if DUMPPICS + fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration); + psphotSaveImage (NULL, mypix, "crmask.fix.fits"); + psphotSaveImage (NULL, myvar, "crmask.var.fits"); + psphotSaveImage (NULL, binned, "crmask.binn.fits"); + psphotSaveImage (NULL, conved, "crmask.conv.fits"); + psphotSaveImage (NULL, edges, "crmask.edge.fits"); + psphotSaveImage (NULL, mymask, "crmask.mask.fits"); +# endif + psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration, nCRpix); + } + +# if 0 + // A solitary masked pixel is likely a lie. Remove those + // XXX can't we use nCRpix == 1 to test for these? + for (int x = 0; x < dx; x++) { + for (int y = 0; y < dy; y++) { + if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40)) continue; + if ((x-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) { + continue; + } + if ((y-1 >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) { + continue; + } + if ((x+1 < dx) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) { + continue; + } + if ((y+1 < dy) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) { + continue; + } + mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40; + } + } +# endif + + // transfer temporary mask to real mask & count masked pixels + nCRpix = 0; + for (int x = 0; x < dx; x++) { + for (int y = 0; y < dy; y++) { if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ys+mask->row0][x+xs+mask->col0] |= maskVal; + nCRpix ++; } } @@ -607,7 +685,12 @@ // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit - source->mode |= PM_SOURCE_MODE_CR_LIMIT; + if (nCRpix > 1) { + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + } + // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix); psFree(mypix); + psFree(myfix); psFree(myvar); psFree(binned); @@ -680,4 +763,147 @@ } return true; +} + +# define VERBOSE 0 +int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) { + + int xLo, xRo; + int nCRpix = 0; + + float noise_factor = 5.0 / 4.0; // Intrinsic to the Laplacian making noise spikes spikier. + + // mark the pixels in this row to the left, then the right. stay within footprint + int xL = xPeak; // find the range of valid pixels in this row + int xR = xPeak; + for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) { + float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]); + float value = edges->data.F32[yPeak][ix] / noise; + if (value < sigma_thresh ) break; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40; + xL = ix; + nCRpix ++; + if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR); + } + for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) { + float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]); + float value = edges->data.F32[yPeak][ix] / noise; + if (value < sigma_thresh ) break; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40; + xR = ix; + nCRpix ++; + if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR); + } + // xL and xR mark the first and last valid pixel in the row + + // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR + xLo = PS_MAX(xL - 1, 0); + xRo = PS_MIN(xR + 1, mymask->numCols); + + // first go down: + for (int iy = yPeak - 1; iy >= 0; iy--) { + + int xLn = -1; + int xRn = -1; + int newPix = 0; + + // mark the pixels in the good range + for (int ix = xLo; ix < xRo; ix++) { + if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint + float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]); + float value = edges->data.F32[iy][ix] / noise; + if (value < sigma_thresh ) continue; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40; + if (xLn == -1) xLn = ix; // first valid pixel in this row + xRn = ix; // last valid pixel in this row + nCRpix ++; + newPix ++; + if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn); + } + + // mark the pixels to the left of the good range + for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) { + float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]); + float value = edges->data.F32[iy][ix] / noise; + if (value < sigma_thresh ) break; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40; + if (xRn == -1) xRn = ix; // last valid pixel in this row + xLn = ix; + nCRpix ++; + newPix ++; + if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn); + } + + // mark the pixels to the right of the good range + for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) { + float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]); + float value = edges->data.F32[iy][ix] / noise; + if (value < sigma_thresh ) break; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40; + if (xLn == -1) xLn = ix; // first valid pixel in this row + xRn = ix; + nCRpix ++; + newPix ++; + if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn); + } + if (newPix == 0) break; + xLo = PS_MAX(xLn - 1, 0); + xRo = PS_MIN(xRn + 1, mymask->numCols); + } + + xLo = PS_MAX(xL - 1, 0); + xRo = PS_MIN(xR + 1, mymask->numCols); + + // next go up: + for (int iy = yPeak + 1; iy < mymask->numRows; iy++) { + + int xLn = -1; + int xRn = -1; + int newPix = 0; + + // mark the pixels in the good range + for (int ix = xLo; ix < xRo; ix++) { + if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint + float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]); + float value = edges->data.F32[iy][ix] / noise; + if (value < sigma_thresh ) continue; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40; + if (xLn == -1) xLn = ix; // first valid pixel in this row + xRn = ix; // last valid pixel in this row + nCRpix ++; + newPix ++; + if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn); + } + + // mark the pixels to the left of the good range + for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) { + float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]); + float value = edges->data.F32[iy][ix] / noise; + if (value < sigma_thresh ) break; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40; + if (xRn == -1) xRn = ix; // last valid pixel in this row + xLn = ix; + nCRpix ++; + newPix ++; + if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn); + } + + // mark the pixels to the right of the good range + for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) { + float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]); + float value = edges->data.F32[iy][ix] / noise; + if (value < sigma_thresh ) break; + mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40; + if (xLn == -1) xLn = ix; // first valid pixel in this row + xRn = ix; + nCRpix ++; + newPix ++; + if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn); + } + if (newPix == 0) break; + xLo = PS_MAX(xLn - 1, 0); + xRo = PS_MIN(xRn + 1, mymask->numCols); + } + + return nCRpix; } Index: trunk/psphot/src/psphotStackChisqImage.c =================================================================== --- trunk/psphot/src/psphotStackChisqImage.c (revision 27532) +++ trunk/psphot/src/psphotStackChisqImage.c (revision 27532) @@ -0,0 +1,85 @@ +# include "psphotInternal.h" + +// Chi = sqrt(\sum_i im_i^2/var_i) + +static int npass = 0; + +// XXX supply filename or keep PSPHOT.INPUT fixed? +bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view) +{ + bool status = false; + + psTimerStart ("psphot.chisq.image"); + + // create a holder for the image + pmReadout *chiImage = pmReadoutAlloc(); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotStackChisqImageAddReadout(config, view, chiImage, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); + return false; + } + } + + // need to save the resulting image somewhere (PSPHOT.INPUT?) + + psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image")); + + return true; +} + +bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration + pmFPAview *view, + pmReadout *chiReadout, + char *filename, + int index) +{ + bool status = true; + + // find the currently selected readout + pmFPAfile *input = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (input, "missing file?"); + + pmReadout *inReadout = pmFPAviewThisReadout(view, input->fpa); + psAssert (inReadout, "missing readout?"); + + psImage *inImage = inReadout->image; + psImage *inVariance = inReadout->variance; + psImage *inMask = inReadout->mask; + + psImage *chiImage = chiReadout->image; + psImage *chiVariance = chiReadout->variance; + psImage *chiMask = chiReadout->mask; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // generate the chisq image: + for (int iy = 0; iy < inImage->numRows; iy++) { + for (int ix = 0; ix < inImage->numCols; ix++) { + if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + chiImage->data.F32[iy][ix] += PS_SQR(inImage->data.F32[iy][ix]) / inVariance->data.F32[iy][ix]; + chiVariance->data.F32[iy][ix] = 1.0; // ?? what is the right value? just init to this? + chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x00; // we have valid data so unmask this pixel + } + } + + if (psTraceGetLevel("psphot") > 5) { + char name[256]; + sprintf (name, "chiImage.%02d.fits", npass); + psphotSaveImage (NULL, chiImage, name); + } + npass ++; + + return true; +} + Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 27314) +++ trunk/psphot/src/psphotVisual.c (revision 27532) @@ -1099,4 +1099,16 @@ psFree (outsat); return true; +} + +static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) +{ + float x[2], y[2]; + x[0] = x0; + x[1] = x1; + y[0] = y0; + y[1] = y1; + KapaPrepPlot (myKapa, 2, graphdata); + KapaPlotVector (myKapa, 2, x, "x"); + KapaPlotVector (myKapa, 2, y, "y"); } @@ -1138,4 +1150,64 @@ } + // generate model profiles (major and minor axis): + // create a model with theta = 0.0 so major and minor axes are equiv to x and y: + psEllipseShape rawShape, rotShape; + + rawShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + rawShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0); + + axes.theta = 0.0; + + rotShape = psEllipseAxesToShape (axes); + + psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32); + for (int i = 0; i < source->modelPSF->params->n; i++) { + params->data.F32[i] = source->modelPSF->params->data.F32[i]; + } + params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2; + params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2; + params->data.F32[PM_PAR_SXY] = rotShape.sxy; + params->data.F32[PM_PAR_XPOS] = 0.0; + params->data.F32[PM_PAR_YPOS] = 0.0; + + psVector *rmod = psVectorAlloc(300, PS_TYPE_F32); + psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32); + psVector *fmin = psVectorAlloc(300, PS_TYPE_F32); + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + float r = 0.0; + for (int i = 0; i < rmod->n; i++) { + r = i*0.1; + rmod->data.F32[i] = r; + + coord->data.F32[1] = r; + coord->data.F32[0] = 0.0; + fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + + coord->data.F32[0] = r; + coord->data.F32[1] = 0.0; + fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + } + psFree (coord); + psFree (params); + + float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); + + psEllipseMoments emoments; + emoments.x2 = source->moments->Mxx; + emoments.xy = source->moments->Mxy; + emoments.y2 = source->moments->Myy; + axes = psEllipseMomentsToAxes (emoments, 20.0); + float MOMENTS_MAJOR = 2.355*axes.major; + float MOMENTS_MINOR = 2.355*axes.minor; + + float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + // reset source Add/Sub state to recorded if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); @@ -1174,4 +1246,33 @@ KapaPlotVector (myKapa, nb, fb->data.F32, "y"); + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y"); + plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM); + plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0); + graphdata.ltype = 1; + plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0); + graphdata.ltype = 0; + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y"); + plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0); + graphdata.ltype = 1; + plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0); + graphdata.ltype = 0; + + for (int i = 0; i < rmod->n; i++) { + rmod->data.F32[i] = log10(rmod->data.F32[i]); + } + // ** loglog ** KapaSelectSection (myKapa, "loglog"); @@ -1182,4 +1283,5 @@ graphdata.ymin = -0.05; graphdata.ymax = +5.05; + graphdata.color = KapaColorByName ("black"); KapaSetLimits (myKapa, &graphdata); @@ -1204,4 +1306,24 @@ KapaPlotVector (myKapa, nb, Rb->data.F32, "x"); KapaPlotVector (myKapa, nb, fb->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y"); + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (myKapa, rmod->n, &graphdata); + KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y"); + + psFree (rmod); + psFree (fmin); + psFree (fmaj); psFree (rg); @@ -1392,5 +1514,4 @@ if (source == NULL) continue; - // if (source->type != type) continue; if (mode) { if (keep) { @@ -1516,4 +1637,9 @@ psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *xLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sLOW = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); @@ -1526,8 +1652,12 @@ int nPSF = 0; int nDEF = 0; + int nLOW = 0; int nCR = 0; for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; if (source->moments == NULL) continue; + + // only plot the measured sources... + if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue; if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { @@ -1561,5 +1691,10 @@ continue; } - if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { + if (source->errMag > 0.1) { + xLOW->data.F32[nLOW] = source->moments->Mxx; + yLOW->data.F32[nLOW] = source->moments->Myy; + mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum); + sLOW->data.F32[nLOW] = source->extNsigma; + nLOW++; continue; } @@ -1570,4 +1705,5 @@ nPSF++; } + xSAT->n = nSAT; ySAT->n = nSAT; @@ -1594,4 +1730,9 @@ mDEF->n = nDEF; sDEF->n = nDEF; + + xLOW->n = nLOW; + yLOW->n = nLOW; + mLOW->n = nLOW; + sLOW->n = nLOW; // four sections: MxxMyy, MagMxx, MagMyy, MagSigma @@ -1657,4 +1798,12 @@ KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y"); + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nLOW, &graphdata); + KapaPlotVector (myKapa, nLOW, xLOW->data.F32, "x"); + KapaPlotVector (myKapa, nLOW, yLOW->data.F32, "y"); + // second section: MagMyy section.dx = 0.75; @@ -1668,5 +1817,5 @@ graphdata.color = KapaColorByName ("black"); graphdata.xmin = -17.1; - graphdata.xmax = -7.9; + graphdata.xmax = -6.9; graphdata.ymin = Ymin; graphdata.ymax = Ymax; @@ -1730,5 +1879,5 @@ graphdata.xmin = Xmin; graphdata.xmax = Xmax; - graphdata.ymin = -7.9; + graphdata.ymin = -6.9; graphdata.ymax = -17.1; KapaSetLimits (myKapa, &graphdata); @@ -1789,5 +1938,5 @@ graphdata.color = KapaColorByName ("black"); - graphdata.xmax = -7.9; + graphdata.xmax = -6.9; graphdata.xmin = -17.1; graphdata.ymin = -20.1;