Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 34363) +++ trunk/psphot/src/psphot.h (revision 34404) @@ -281,6 +281,6 @@ bool psphotVisualShowSatStars (psMetadata *recipe, pmPSF *psf, psArray *sources); bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf); -bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal); -bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources); +bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode); +bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode); bool psphotVisualShowFlags (psArray *sources); bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); @@ -295,4 +295,8 @@ bool psphotVisualClose(void); +int psphotKapaChannel (int channel); +void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1); + + bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise); @@ -511,4 +515,11 @@ ); +bool psphotSatstarProfileModel (pmSource *source, psImageMaskType maskVal); +bool psphotSatstarProfileCreate (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logR, psVector *flux, float Rmax); +bool psphotSatstarProfileOp (pmSource *source, psImageMaskType maskVal, float FACTOR, pmModelOpMode mode, bool add); +bool psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal); +bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add); +bool psphotSatstarPhotometry (pmSource *source); + bool psphotSourceMemory(pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index); Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 34363) +++ trunk/psphot/src/psphotAddNoise.c (revision 34404) @@ -76,7 +76,16 @@ } + psphotVisualShowImage (readout); + // loop over all source for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; + + // add or subtract noise for a saturated star. satstars modeled as a radial profile + // need special handling + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) { + psphotSatstarProfileOp (source, maskVal, FACTOR, PM_MODEL_OP_NOISE, add); + continue; + } // skip sources which were not subtracted @@ -92,4 +101,6 @@ } + psphotVisualShowImage (readout); + return true; } Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 34363) +++ trunk/psphot/src/psphotApResid.c (revision 34404) @@ -200,4 +200,7 @@ if (source->mode & PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR"); if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) SKIPSTAR ("SATSTAR"); if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED"); Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 34363) +++ trunk/psphot/src/psphotBlendFit.c (revision 34404) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma); // for now, let's store the detections on the readout->analysis for each readout @@ -92,4 +93,17 @@ float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG"); assert (status && isfinite(skySig) && skySig > 0); + +# 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); + psImageInit (IDimage, 0); + + // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + psphotBlendFitSetSource (source, IDimage, skySig); + } + psphotSaveImage (NULL, IDimage, "idimage.fits"); +# endif // Define source fitting parameters for extended source fits @@ -241,12 +255,13 @@ pmSource *source = sources->data[i]; -# if (0) -# define TEST_X 34 -# define TEST_Y 28 - + int TEST_ON = false; +# if (1) +# 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); + TEST_ON = true; + } # undef TEST_X # undef TEST_Y @@ -258,4 +273,7 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; // skip DBL second sources (ie, added by psphotFitBlob @@ -306,4 +324,7 @@ if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) { + if (TEST_ON) { + psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); + } psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); Next ++; @@ -312,4 +333,7 @@ } else { if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) { + if (TEST_ON) { + psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); + } source->type = PM_SOURCE_TYPE_STAR; psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); @@ -349,2 +373,82 @@ } +bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float skySigma) { + + if (!source) return false; + if (!source->peak) return false; // XXX how can we have a peak-less source? + if (!source->moments) return false; + if (source->type == PM_SOURCE_TYPE_DEFECT) return false; + if (source->type == PM_SOURCE_TYPE_SATURATED) return false; + if (source->mode2 == PM_SOURCE_MODE2_MATCHED) return false; + psAssert(IDimage, "need a window"); + + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) { + // find the radius at which we hit something like 0.1*SATURATION + return false; + } + + if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false; + + // we have a source with moments Mx, My, Mxx, Mxy, Myy. we just need to define a Gaussian that has + // these values (and peak of 1.0) + + int Nx = IDimage->numCols; + int Ny = IDimage->numRows; + + float Xo = source->moments->Mx; + float Yo = source->moments->My; + + psEllipseMoments moments; + moments.x2 = source->moments->Mxx; + moments.y2 = source->moments->Myy; + moments.xy = source->moments->Mxy; + + psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0); + if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) { + fprintf (stderr, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta); + return false; + } + + // Use 1st radial moment as size, not sigma. Why this factor of 0.5 ? + float scale = 0.5 * source->moments->Mrf / axes.major; + axes.major *= scale; + axes.minor *= scale; + + psEllipseShape shape = psEllipseAxesToShape(axes); + if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) { + fprintf (stderr, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy); + return false; + } + + float Smajor = axes.major; + + // we want to go to 2.15 sigma (0.1 Io) + int minX = PS_MIN(PS_MAX(Xo - 2.15*Smajor, 0), Nx - 1); + int maxX = PS_MIN(PS_MAX(Xo + 2.15*Smajor, 0), Nx - 1); + int minY = PS_MIN(PS_MAX(Yo - 2.15*Smajor, 0), Ny - 1); + int maxY = PS_MIN(PS_MAX(Yo + 2.15*Smajor, 0), Ny - 1); + + 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? + + int ID = source->id; + float Io = source->peak->rawFlux; + for (int iy = minY; iy < maxY; iy++) { + for (int ix = minX; ix < maxX; ix++) { + + float dX = (ix + 0.5 - Xo); + float dY = (iy + 0.5 - Yo); + + float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY; + + float f = Io*exp(-z); + + if (f < 2.0*skySigma) continue; + IDimage->data.S32[iy][ix] = ID; + } + } + + return true; +} Index: trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- trunk/psphot/src/psphotDeblendSatstars.c (revision 34363) +++ trunk/psphot/src/psphotDeblendSatstars.c (revision 34404) @@ -1,3 +1,15 @@ # include "psphotInternal.h" + +typedef struct { + float min; + float max; + float lower20; + float upper20; + int Npts; +} QuickStats; + +bool psImageQuickStats (QuickStats *stats, psImage *image, psRegion region); +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); +bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, float factor, int channel); // for now, let's store the detections on the readout->analysis for each readout @@ -16,5 +28,171 @@ } +/** this function does 3 things: + + 1) choose likely saturated stars + - mode | SATSTAR + - if moments->Peak > 0.5*SATURATION, examine region of peak pixels (ensure we go down to 0.1*SAT) for SAT mask + + 2) adjust the window for saturated objects + + 3) measure the radial profile + + 4) subtract the radial profile + + 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) { + + int N; + pmSource *source; + bool status; + + // XXX disable this function for now + return true; + + psTimerStart ("psphot.deblend.sat"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // 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->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend"); + 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) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskSat = psMetadataLookupImageMask(&status, recipe, "MASK.SAT"); // Mask value for bad pixels + assert (maskSat); + + float SATURATION = NAN; + { + // XXX do we need to set this differently from the value used to mark saturated pixels? + pmCell *cell = readout->parent; + + // do not completely trust the values in the header... + float CELL_SATURATION = psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION"); + float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION"); + if (!status || !isfinite(MIN_SATURATION)) { + MIN_SATURATION = 40000.0; + } + if (!isfinite(CELL_SATURATION)) { + SATURATION = MIN_SATURATION; + } else { + SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION); + } + } + + // source analysis is done in peak order (brightest first) + // we use an index for this so the spatial sorting is kept + psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32); + for (int i = 0; i < SN->n; i++) { + source = sources->data[i]; + SN->data.F32[i] = source->moments->SN; + } + psVector *index = psVectorSortIndex (NULL, SN); + // this results in an index of increasing SN + + // psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_SATSTAR); + + int BIG_RADIUS = 250; + int BIG_SIGMA = BIG_RADIUS / 4.0; + + int display = psphotKapaChannel (1); + psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0); + + // examine sources in decreasing SN order + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; + + bool isSat = source->mode & PM_SOURCE_MODE_SATSTAR; + + // for fairly bright stars, check the pixels near the peak again + if (source->moments->Peak > 0.5*SATURATION) { + // pmSourceRoughClass does this analysis, but uses a small (5x5) window. + // here we use a larger window since stacks can have some funny features + psRegion inner; + QuickStats stats; + // grow out the search radius until we have lower20 < 0.1*SATURATION + for (int radius = 2; radius < 30; radius ++) { + inner = psRegionForSquare (source->peak->x, source->peak->y, radius); + inner = psRegionForImage (source->maskView, inner); + psImageQuickStats (&stats, readout->image, inner); + if ((stats.Npts > 1) && (stats.lower20 < 0.1*SATURATION)) break; + } + int Nsatpix = psImageCountPixelMask (source->maskView, inner, maskSat); + // fprintf (stderr, "test object: %d,%d : %d vs %d : %f - %f - %f - %f\n", + // source->peak->x, source->peak->y, Nsatpix, stats.Npts, + // stats.min, stats.lower20, stats.upper20, stats.max); + if (Nsatpix > 1) isSat = true; + } + if (!isSat) continue; + + // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but + // only slightly: pmSourceRedefinePixels uses the readout to pass the pointers to + // the parent image data. I guess the API could be simplified: we could recover + // this from the source in the function + + pmReadout tmpReadout; + tmpReadout.image = (psImage *)source->pixels->parent; + tmpReadout.mask = (psImage *)source->maskView->parent; + tmpReadout.variance = (psImage *)source->variance->parent; + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, &tmpReadout, source->peak->x, source->peak->y, BIG_RADIUS + 2); + + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, 5.0, maskVal); + source->mode |= PM_SOURCE_MODE_BIG_RADIUS; + + if (!psphotSatstarProfileModel (source, maskVal)) continue; + + source->mode |= PM_SOURCE_MODE_SATSTAR; // yes, this source IS saturated + source->mode2 |= PM_SOURCE_MODE2_SATSTAR_PROFILE; // and we have in fact subtracted the profile + + // XXX visualize, model, and subtract + // if (!psphotVisualRadialProfileSatstar (source, maskVal)) { + // break; + // } + + // generate radial profile, store on the source structure + } + // show the image after object have been subtracted + psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1); + + psFree (SN); + psFree (index); + + psLogMsg ("psphot", PS_LOG_INFO, "deblend satstar: %f sec\n", psTimerMark ("psphot.deblend.sat")); + return true; +} + +bool psphotDeblendSatstarsReadoutOld (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { int N; @@ -219,2 +397,690 @@ return true; } + +// create a model for the radial profile of a saturated star (is this actually more generic?) +bool psphotSatstarProfileModel (pmSource *source, psImageMaskType maskVal) { + + // XXX user define somewhere? + float Rmax = 320.0; + + // XXX is this ever the case?? + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + int nPts = source->pixels->numRows * source->pixels->numCols; + psVector *logR = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *flux = psVectorAllocEmpty (nPts, PS_TYPE_F32); + + int ng = 0; + + // choose the best center for this profile + float Xo = NAN; + float Yo = NAN; + + if (source->modelPSF) { + // XXX do we ever have a PSF model for a SATSTAR? + Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + Xo = source->moments->Mx; + Yo = source->moments->My; + } + + float Xc = Xo - source->pixels->col0 - 0.5; + float Yc = Yo - source->pixels->row0 - 0.5; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) continue; + // XXX do this faster by generating R^2 and returning 0.5*log10(R^2)? + logR->data.F32[ng] = log10(hypot (ix - Xc, iy - Yc)); + flux->data.F32[ng] = source->pixels->data.F32[iy][ix]; + ng++; + } + } + logR->n = ng; + flux->n = ng; + + // XXX do something sensible here if there are no pixels + + // measure the radial profile + psVector *logFmodel = NULL; + psVector *logRmodel = NULL; + psphotSatstarProfileCreate (source, &logRmodel, &logFmodel, logR, flux, Rmax); + + // XXX do something sensible here if the profile is crap + + source->satstar = pmSourceSatstarAlloc(); + source->satstar->Xo = Xo; + source->satstar->Yo = Yo; + source->satstar->Rmax = Rmax; + source->satstar->logFmodel = logFmodel; + source->satstar->logRmodel = logRmodel; + + // subtract the profile (false => subtract) + psphotSatstarProfileOp (source, maskVal, 1.0, 0, false); + + psFree (logR); + psFree (flux); + + return true; +} + +// Take logR + flux and generate radial bins logRmodelOut, logFmodelOut +bool psphotSatstarProfileCreate (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logR, psVector *flux, float Rmax) { + + // we have log(radius) & log(flux). find the median flux in log radial bins + + float logRmax = log10(Rmax); + float logRdel = 0.1; // XXX user-define? + + psVector *logRmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32); + psVector *logFmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32); + + pmSourceRadialProfileSortPair (logR, flux); + + psStats *fluxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + psVector *fluxVals = psVectorAllocEmpty (100, PS_TYPE_F32); + + int bin = 0; + for (int i = 0; i < logRmodel->n; i++) { + + // bin (i) has log radius range (i*logRdel : (i+1)*logRdel) + float lRmin = logRdel*(i + 0); + float lRmax = logRdel*(i + 1); + + // reset the flux vector + fluxVals->n = 0; + psStatsInit (fluxStats); + + while (logR->data.F32[bin] < lRmax) { + if (isfinite(flux->data.F32[bin])) { + psVectorAppend (fluxVals, flux->data.F32[bin]); + } + bin ++; + } + + // we have the set of fluxes for this bin; find the median values + + float Rmin = pow(10.0, lRmin); + float Rmax = pow(10.0, lRmax); + + float Rmean = (2.0/3.0) * (pow(Rmax, 3.0) - pow(Rmin, 3.0)) / (PS_SQR(Rmax) - PS_SQR(Rmin)); + logRmodel->data.F32[i] = log10(Rmean); + + float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin)); + if (fluxVals->n < 0.25*Area) { + logFmodel->data.F32[i] = NAN; + continue; + } + + psVectorStats (fluxStats, fluxVals, NULL, NULL, 0); + if (fluxStats->robustMedian > 0.0) { + logFmodel->data.F32[i] = log10(fluxStats->robustMedian); + } else { + logFmodel->data.F32[i] = -3.0; + } + // fprintf (stderr, "R: %f, F: %f +/- %f\n", Rmean, fluxStats->robustMedian, fluxStats->robustStdev); + } + + // now how do i use this to subtract a model?? + *logRmodelOut = logRmodel; + *logFmodelOut = logFmodel; + + // need to free stuff + psFree (fluxStats); + psFree (fluxVals); + + return true; +} + +bool psphotSatstarProfileOp (pmSource *source, psImageMaskType maskVal, float FACTOR, pmModelOpMode mode, bool add) { + + int alt; + float logRdel = 0.1; + + float Xc = source->satstar->Xo - source->pixels->col0 - 0.5; + float Yc = source->satstar->Yo - source->pixels->row0 - 0.5; + psVector *logRmodel = source->satstar->logRmodel; + psVector *logFmodel = source->satstar->logFmodel; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) continue; + + float radius = hypot (ix - Xc, iy - Yc) ; + float logR = log10(radius); + + int bin = (int)(logR / logRdel); + + // we are going to interpolate if possible, or extrapolate if necessary + if (bin >= logRmodel->n - 1) bin = logRmodel->n - 1; + if (bin < 0) bin = 0; + + // interpolate to the current radial position + // XXX BIG HACK : skip nan bins for now + + float logF0 = logFmodel->data.F32[bin]; + float logR0 = logRmodel->data.F32[bin]; + if (!isfinite(logF0)) continue; + + // interpolate between closest two bins if possible, extrapolate on ends + if (logR < logR0) { + alt = (bin > 0) ? bin - 1 : bin + 1; + } else { + alt = (bin < logRmodel->n - 1) ? bin + 1 : bin - 1; + } + + float logF1 = logFmodel->data.F32[alt]; + float logR1 = logRmodel->data.F32[alt]; + if (!isfinite(logF1)) continue; + + // XXX use linear flux, not logFlux + float logF = InterpolateValues (logR0, logF0, logR1, logF1, logR); + float flux = pow (10.0, logF); + + if (mode & PM_MODEL_OP_NOISE) { + if (add) { + source->variance->data.F32[iy][ix] += FACTOR * flux; + } else { + source->variance->data.F32[iy][ix] -= FACTOR * flux; + } + } else { + if (add) { + source->pixels->data.F32[iy][ix] += flux; + } else { + source->pixels->data.F32[iy][ix] -= flux; + } + } + } + } + return true; +} + +// Take logR + flux and generate radial bins logRmodelOut, logFmodelOut +bool psphotSatstarPhotometry (pmSource *source) { + + if (!source->satstar) return false; + + psVector *logRmodel = source->satstar->logRmodel; + psVector *logFmodel = source->satstar->logFmodel; + + float fluxTotal = 0.0; + float logRdel = 0.1; // XXX user-define (or carry in satstar) + + // integrate flux in radial bins + for (int i = 0; i < logRmodel->n; i++) { + // just add up the mean flux in each bin and multiply by the area of the bin + + float logF = logFmodel->data.F32[i]; + if (!isfinite(logF)) continue; + float flux = pow(10.0, logF); // this is the mean flux per pixel (ie, surface brightness) + + // bin (i) has log radius range (i*logRdel : (i+1)*logRdel) + float lRmin = logRdel*(i + 0); + float lRmax = logRdel*(i + 1); + + float Rmin = pow(10.0, lRmin); + float Rmax = pow(10.0, lRmax); + + float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin)); + + fluxTotal += flux * Area; + } + + source->psfMag = -2.5*log10(fluxTotal); + source->psfMagErr = 0.05; + + // XXX I have no idea of a realistic error on the photometry here. the bottom line is that + // the error is totally dominated by the loss of charge due to saturation. I am not + // modeling the profile in any real detail other than to follow the radial bins. + + source->extMag = NAN; + source->apMag = NAN; + source->apMagRaw = NAN; + source->apFlux = fluxTotal; + source->apFluxErr = 0.05*fluxTotal; + + return true; +} + +# if (0) +static bool skipDisplay = false; +bool psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal) { + + int kapaImage = -1; + int kapaGraph = -1; + Graphdata graphdata; + + if (!pmVisualTestLevel("psphot.satstar", 3)) { + skipDisplay = true; + } + + float Rmax = 320.0; + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + int nPts = source->pixels->numRows * source->pixels->numCols; + psVector *rg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *Rg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *fg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *Fg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *rb = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *Rb = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *fb = psVectorAllocEmpty (nPts, PS_TYPE_F32); + + int ng = 0; + int nb = 0; + + float Xo = NAN; + float Yo = NAN; + + if (source->modelPSF) { + Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; + Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + } else { + Xo = source->moments->Mx - source->pixels->col0; + Yo = source->moments->My - source->pixels->row0; + } + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) { + rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; + Rb->data.F32[nb] = log10(rb->data.F32[nb]); + fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); + nb++; + } else { + rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; + Rg->data.F32[ng] = log10(rg->data.F32[ng]); + Fg->data.F32[ng] = source->pixels->data.F32[iy][ix]; + fg->data.F32[ng] = log10(Fg->data.F32[ng]); + ng++; + } + } + } + rg->n = ng; + Rg->n = ng; + fg->n = ng; + Fg->n = ng; + rb->n = nb; + Rb->n = nb; + fb->n = nb; + + if (!skipDisplay) { + kapaImage = psphotKapaChannel (1); + if (kapaImage == -1) return false; + + kapaGraph = psphotKapaChannel (2); + if (kapaGraph == -1) return false; + + KapaSection section; // put the positive profile in one and the residuals in another? + + // first section : mag vs CR nSigma + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.0; + section.bg = -1; + section.name = NULL; + psStringAppend (§ion.name, "linlog"); + KapaSetSection (kapaGraph, §ion); + psFree (section.name); + + // first section : mag vs CR nSigma + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.5; + section.bg = -1; + section.name = NULL; + psStringAppend (§ion.name, "loglog"); + KapaSetSection (kapaGraph, §ion); + psFree (section.name); + + KapaInitGraph (&graphdata); + + // ** linlog ** + KapaSelectSection (kapaGraph, "linlog"); + + // examine sources to set data range + graphdata.xmin = -0.05; + graphdata.xmax = Rmax + 0.05; + graphdata.ymin = -0.05; + graphdata.ymax = +8.05; + KapaSetLimits (kapaGraph, &graphdata); + + KapaSetFont (kapaGraph, "helvetica", 14); + KapaBox (kapaGraph, &graphdata); + KapaSendLabel (kapaGraph, "radius (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapaGraph, "log flux (counts)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, ng, &graphdata); + KapaPlotVector (kapaGraph, ng, rg->data.F32, "x"); + KapaPlotVector (kapaGraph, ng, fg->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, nb, &graphdata); + KapaPlotVector (kapaGraph, nb, rb->data.F32, "x"); + KapaPlotVector (kapaGraph, nb, fb->data.F32, "y"); + + // ** loglog ** + KapaSelectSection (kapaGraph, "loglog"); + + // examine sources to set data range + graphdata.xmin = -1.51; + graphdata.xmax = log10(Rmax) + 0.02; + graphdata.ymin = -5.05; + graphdata.ymax = +8.05; + graphdata.color = KapaColorByName ("black"); + KapaSetLimits (kapaGraph, &graphdata); + + KapaSetFont (kapaGraph, "helvetica", 14); + KapaBox (kapaGraph, &graphdata); + KapaSendLabel (kapaGraph, "log radius (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapaGraph, "log flux (counts)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, ng, &graphdata); + KapaPlotVector (kapaGraph, ng, Rg->data.F32, "x"); + KapaPlotVector (kapaGraph, ng, fg->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, nb, &graphdata); + KapaPlotVector (kapaGraph, nb, Rb->data.F32, "x"); + KapaPlotVector (kapaGraph, nb, fb->data.F32, "y"); + } + + // measure the radial profile + psVector *logFmodel = NULL; + psVector *logRmodel = NULL; + psphotTestRadialModel (source, &logRmodel, &logFmodel, Rg, Fg, Rmax); + + if (!skipDisplay) { + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, logRmodel->n, &graphdata); + KapaPlotVector (kapaGraph, logRmodel->n, logRmodel->data.F32, "x"); + KapaPlotVector (kapaGraph, logRmodel->n, logFmodel->data.F32, "y"); + } + + // subtract the model from the images + psphotTestRadialModelSub (source, logRmodel, logFmodel, Xo, Yo, Rmax, maskVal); + + psFree (logRmodel); + psFree (logFmodel); + + if (!skipDisplay && source->modelPSF) { + // 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(Rmax*10, PS_TYPE_F32); + psVector *fmaj = psVectorAlloc(Rmax*10, PS_TYPE_F32); + psVector *fmin = psVectorAlloc(Rmax*10, 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); + + // ** linlog ** + KapaSelectSection (kapaGraph, "linlog"); + KapaGetGraphData (kapaGraph, &graphdata); + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmin->data.F32, "y"); + plotline (kapaGraph, &graphdata, graphdata.xmin, logHM, graphdata.xmax, logHM); + plotline (kapaGraph, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax); + graphdata.ltype = 1; + plotline (kapaGraph, &graphdata, 0.5*MOMENTS_MINOR, graphdata.ymin, 0.5*MOMENTS_MINOR, graphdata.ymax); + graphdata.ltype = 0; + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmaj->data.F32, "y"); + plotline (kapaGraph, &graphdata, 0.5*FWHM_MAJOR, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax); + graphdata.ltype = 1; + plotline (kapaGraph, &graphdata, 0.5*MOMENTS_MAJOR, graphdata.ymin, 0.5*MOMENTS_MAJOR, graphdata.ymax); + graphdata.ltype = 0; + + for (int i = 0; i < rmod->n; i++) { + rmod->data.F32[i] = log10(rmod->data.F32[i]); + } + + // ** loglog ** + KapaSelectSection (kapaGraph, "loglog"); + KapaGetGraphData (kapaGraph, &graphdata); + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmin->data.F32, "y"); + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmaj->data.F32, "y"); + + psFree (rmod); + psFree (fmin); + psFree (fmaj); + } + + if (!skipDisplay) { + KiiCenter (kapaImage, source->peak->xf, source->peak->yf, 1); + psphotVisualScaleImage (kapaImage, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1); + KiiOverlay overlay; + overlay.x = source->peak->xf; + overlay.y = source->peak->yf; + overlay.dx = 5; + overlay.dy = 5; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (kapaImage, &overlay, 1, "red"); + overlay.x = source->moments->Mx; + overlay.y = source->moments->My; + overlay.dx = 8; + overlay.dy = 8; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (kapaImage, &overlay, 1, "blue"); + } + + psFree (rg); + psFree (Rg); + psFree (fg); + psFree (Fg); + psFree (rb); + psFree (Rb); + psFree (fb); + + if (skipDisplay) return true; + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[q]uit satstar? [e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + if (key[0] == 'e') { + KapaClearPlots (kapaGraph); + } + if (key[0] == 'q') { + return false; + } + if (key[0] == 's') { + skipDisplay = true; + } + return true; +} +# endif + +// only valid for F32 +bool psImageQuickStats (QuickStats *stats, psImage *image, psRegion region) { + + psAssert (image->type.type == PS_TYPE_F32, "unsupported image type"); + + int Npix = (region.y1 - region.y0) * (region.x1 - region.x0); + psVector *tmp = psVectorAllocEmpty (Npix, PS_TYPE_F32); + + int bin = 0; + for (int iy = region.y0; iy < region.y1; iy++) { + for (int ix = region.x0; ix < region.x1; ix++) { + if (!isfinite(image->data.F32[iy][ix])) continue; + tmp->data.F32[bin] = image->data.F32[iy][ix]; + bin ++; + } + } + tmp->n = bin; + + if (bin < 1) { + stats->min = NAN; + stats->lower20 = NAN; + stats->upper20 = NAN; + stats->max = NAN; + stats->Npts = 0; + psFree (tmp); + return true; + } + + + psVectorSortInPlace (tmp); + + int N20 = 0.2*tmp->n; + int N80 = 0.8*tmp->n; + int Nmax = tmp->n - 1; + + stats->min = tmp->data.F32[0]; + stats->lower20 = tmp->data.F32[N20]; + stats->upper20 = tmp->data.F32[N80]; + stats->max = tmp->data.F32[Nmax]; + stats->Npts = bin; + psFree (tmp); + + return true; +} + +bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add) { + + bool status; + + psTimerStart ("psphot.deblend.sat"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // 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 satstar blend"); + return true; + } + + // 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); + + // examine sources in decreasing SN order + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + if (!(source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE)) continue; + + if (!psphotSatstarProfileOp (source, maskVal, 1.0, 0, add)) continue; + } + + psLogMsg ("psphot", PS_LOG_DETAIL, "satstar op: %f sec\n", psTimerMark ("psphot.deblend.sat")); + return true; +} + Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 34363) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 34404) @@ -218,4 +218,7 @@ if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + // optionally allow non-extended objects to get petrosians as well if (!doPetroStars) { Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 34363) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 34404) @@ -317,4 +317,7 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; // XXX this should use peak? Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 34363) +++ trunk/psphot/src/psphotFindDetections.c (revision 34404) @@ -120,8 +120,16 @@ if (useFootprints) { if (replaceSourcesForFootprints) { + // subtract the noise for all sources including satstars psphotAddOrSubNoiseReadout(config, view, filerule, index, recipe, false); psphotReplaceAllSourcesReadout (config, view, filerule, index, recipe, false); + + // add in the satstars + psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, true); + psFree (significance); significance = psphotSignificanceImage (readout, recipe, maskVal); + + // display the significance image + psphotVisualShowSignificance (significance->variance, 0.98*threshold, 1.02*threshold); } @@ -130,4 +138,7 @@ if (replaceSourcesForFootprints) { psphotRemoveAllSourcesReadout (config, view, filerule, index, recipe, false); + + // subtract the satstars + psphotAddOrSubSatstarsReadout (config, view, filerule, index, recipe, false); } } Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 34363) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 34404) @@ -222,4 +222,7 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; // do not include CRs in the full ensemble fit Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 34363) +++ trunk/psphot/src/psphotGuessModels.c (revision 34404) @@ -169,6 +169,7 @@ source->tmpFlags |= PM_SOURCE_TMPF_MODEL_GUESS; - // skip non-astronomical objects (very likely defects) - if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + // skip non-astronomical objects (very likely defects) and satstars with profiles subtracted + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; Index: trunk/psphot/src/psphotKronIterate.c =================================================================== --- trunk/psphot/src/psphotKronIterate.c (revision 34363) +++ trunk/psphot/src/psphotKronIterate.c (revision 34404) @@ -331,4 +331,7 @@ } + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + // replace object in image bool reSubtract = false; Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 34363) +++ trunk/psphot/src/psphotMagnitudes.c (revision 34404) @@ -182,4 +182,10 @@ if (saveTest) { psphotSaveImage(NULL, testImage, "test.image.1.fits"); + } + + // satstars modeled as a radial profile need special handling + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) { + psphotSatstarPhotometry (source); + continue; } Index: trunk/psphot/src/psphotRadialApertures.c =================================================================== --- trunk/psphot/src/psphotRadialApertures.c (revision 34363) +++ trunk/psphot/src/psphotRadialApertures.c (revision 34404) @@ -240,4 +240,7 @@ if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + // XXX measure radial apertures even for saturated stars // if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; Index: trunk/psphot/src/psphotRadialProfileWings.c =================================================================== --- trunk/psphot/src/psphotRadialProfileWings.c (revision 34363) +++ trunk/psphot/src/psphotRadialProfileWings.c (revision 34404) @@ -190,4 +190,15 @@ if (!source->moments) continue; + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + + // skip non-detected sources matched from other images + if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true; // skip matched sources (no signal) + + // XXX unclear if I should run this analysis on both 1st and 2nd pass for 1st pass objects, + // or just use the 1st pass value. I think I should just use the 1st pass value, so skip + // any that have already been assigned + if (isfinite(source->skyRadius)) return true; + // replace object in image bool reSubtract = false; @@ -226,11 +237,4 @@ // psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; - - // XXX unclear if I should run this analysis on both 1st and 2nd pass for 1st pass objects, - // or just use the 1st pass value. I think I should just use the 1st pass value, so skip - // any that have already been assigned - if (isfinite(source->skyRadius)) return true; - - if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true; // skip matched sources (no signal) // radii will be MIN_RADIUS to MAX_RADIUS in NN log steps: Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 34363) +++ trunk/psphot/src/psphotReadout.c (revision 34404) @@ -128,10 +128,4 @@ } - // find blended neighbors of very saturated stars (detections->newSources) - // if (!psphotDeblendSatstars (config, view, filerule)) { - // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - // return psphotReadoutCleanup (config, view, filerule); - // } - // mark blended peaks PS_SOURCE_BLEND (detections->newSources) // XXX I've deactivated this because it was preventing galaxies close to stars from being @@ -149,4 +143,10 @@ } + // find and subtract radial profile models for saturated stars (XXX change name eventually) + if (!psphotDeblendSatstars (config, view, filerule)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + return psphotReadoutCleanup (config, view, filerule); + } + // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) if (!psphotImageQuality (config, view, filerule)) { // pass 1 @@ -221,5 +221,5 @@ // NOTE: this block performs the 2nd pass low-significance PSF detection stage { - // add noise for subtracted objects + // add noise for subtracted objects & subtracted saturated stars psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources) @@ -253,5 +253,4 @@ // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - // XXX check on free of sources... psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources) @@ -294,5 +293,4 @@ // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - // XXX check on free of sources... psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources) Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 34363) +++ trunk/psphot/src/psphotSourceSize.c (revision 34404) @@ -498,4 +498,7 @@ continue; } + + // skip saturated stars modeled with a radial profile + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; // we are classifying by moments and PSF_MAG - KRON_MAG Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 34363) +++ trunk/psphot/src/psphotSourceStats.c (revision 34404) @@ -447,4 +447,7 @@ } + // skip saturated stars modeled with a radial profile (this probably never happens, since it is set after) + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) { // measure basic source moments (no S/N clipping on input pixels) @@ -521,5 +524,5 @@ float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); - + // XXX this will cause an error in the vector length is > 8 # define NSIGMA 8 // moved to config file @@ -538,4 +541,5 @@ sigma[i] = sigmavec->data.F32[i]; } + assert (sigmavec->n <= 8); nsigma = sigmavec->n; } else { Index: trunk/psphot/src/psphotStackMatchPSFs.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFs.c (revision 34363) +++ trunk/psphot/src/psphotStackMatchPSFs.c (revision 34404) @@ -46,5 +46,5 @@ } - // loop over the available readouts (ignore chisq image)psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Source Stats ---"); + // loop over the available readouts (ignore chisq image) for (int i = 0; i < num; i++) { if (!psphotStackMatchPSFsReadout (config, view, options, i)) { Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 34363) +++ trunk/psphot/src/psphotVisual.c (revision 34404) @@ -262,4 +262,7 @@ } +static psImage *posImage = NULL; +static psImage *delImage = NULL; + bool psphotVisualShowImage (pmReadout *readout) { @@ -277,4 +280,8 @@ psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1.0, 1); psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", sqrt(factor), 0); + + if (posImage == NULL) { + posImage = psImageCopy (NULL, readout->image, PS_TYPE_F32); + } pmVisualAskUser(NULL); @@ -1139,5 +1146,5 @@ // after displaying (as an image) the psf stars, we cycle throught them and display their // radial profiles: - psphotVisualPlotRadialProfiles (recipe, sources); + psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_PSFSTAR); return true; @@ -1273,5 +1280,5 @@ } -static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) +void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) { float x[2], y[2]; @@ -1285,7 +1292,9 @@ } -bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal) { +bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode) { Graphdata graphdata; + + float Rmax = (showmode & PM_SOURCE_MODE_SATSTAR) ? 100.0 : 30.0; bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); @@ -1302,9 +1311,19 @@ int ng = 0; int nb = 0; - float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; - float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + + float Xo = NAN; + float Yo = NAN; + + if (source->modelPSF) { + Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; + Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + } else { + Xo = source->moments->Mx - source->pixels->col0; + Yo = source->moments->My - source->pixels->row0; + } + for (int iy = 0; iy < source->pixels->numRows; iy++) { for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) { rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; @@ -1322,67 +1341,4 @@ } - // 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); - KapaInitGraph (&graphdata); @@ -1392,5 +1348,5 @@ // examine sources to set data range graphdata.xmin = -0.05; - graphdata.xmax = +30.05; + graphdata.xmax = Rmax + 0.05; graphdata.ymin = -0.05; graphdata.ymax = +8.05; @@ -1418,33 +1374,4 @@ 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"); @@ -1452,5 +1379,5 @@ // examine sources to set data range graphdata.xmin = -1.51; - graphdata.xmax = +1.51; + graphdata.xmax = log10(Rmax) + 0.02; graphdata.ymin = -0.05; graphdata.ymax = +8.05; @@ -1479,23 +1406,123 @@ 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); + if (source->modelPSF) { + // 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(Rmax*10, PS_TYPE_F32); + psVector *fmaj = psVectorAlloc(Rmax*10, PS_TYPE_F32); + psVector *fmin = psVectorAlloc(Rmax*10, 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); + + // ** linlog ** + KapaSelectSection (myKapa, "linlog"); + KapaGetGraphData (myKapa, &graphdata); + 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, graphdata.xmin, logHM, graphdata.xmax, logHM); + plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax); + graphdata.ltype = 1; + plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, graphdata.ymin, 0.5*MOMENTS_MINOR, graphdata.ymax); + 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, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax); + graphdata.ltype = 1; + plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, graphdata.ymin, 0.5*MOMENTS_MAJOR, graphdata.ymax); + graphdata.ltype = 0; + + for (int i = 0; i < rmod->n; i++) { + rmod->data.F32[i] = log10(rmod->data.F32[i]); + } + + // ** loglog ** + KapaSelectSection (myKapa, "loglog"); + KapaGetGraphData (myKapa, &graphdata); + 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); @@ -1508,5 +1535,5 @@ } -bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) { +bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode) { KapaSection section; // put the positive profile in one and the residuals in another? @@ -1549,7 +1576,27 @@ pmSource *source = sources->data[i]; - if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - - psphotVisualPlotRadialProfile (myKapa, source, maskVal); + if (!(source->mode & showmode)) continue; + + psphotVisualPlotRadialProfile (myKapa, source, maskVal, showmode); + + if (pmVisualTestLevel("psphot.image", 1)) { + int display = psphotKapaChannel (1); + KiiCenter (display, source->peak->xf, source->peak->yf, 1); + KiiOverlay overlay; + overlay.x = source->peak->xf; + overlay.y = source->peak->yf; + overlay.dx = 5; + overlay.dy = 5; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (display, &overlay, 1, "red"); + overlay.x = source->moments->Mx; + overlay.y = source->moments->My; + overlay.dx = 8; + overlay.dy = 8; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (display, &overlay, 1, "blue"); + } // pause and wait for user input: @@ -2540,9 +2587,16 @@ } - if (reshow) { + if (false && reshow) { psphotVisualShowMask (myKapa, readout->mask, "mask", 2); psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1.0, 1); } - psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 0); + + if (posImage) { + delImage = (psImage *) psBinaryOp(delImage, posImage, "-", readout->image); + psphotVisualScaleImage (myKapa, posImage, readout->mask, "posimage", sqrt(factor), 0); + psphotVisualScaleImage (myKapa, delImage, readout->mask, "delimage", sqrt(factor), 2); + } + + psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 1); pmVisualAskUser(NULL);