Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 25738) +++ trunk/psphot/src/Makefile.am (revision 25755) @@ -25,5 +25,7 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotTest psphotMomentsStudy +bin_PROGRAMS = psphot psphotTest psphotMomentsStudy +# bin_PROGRAMS = psphotPetrosianStudy +# bin_PROGRAMS = psphot psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -38,4 +40,8 @@ psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotMomentsStudy_LDADD = libpsphot.la + +# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +# psphotPetrosianStudy_LDADD = libpsphot.la psphot_SOURCES = \ @@ -61,4 +67,7 @@ psphotMomentsStudy_SOURCES = \ psphotMomentsStudy.c + +# psphotPetrosianStudy_SOURCES = \ +# psphotPetrosianStudy.c libpsphot_la_SOURCES = \ @@ -104,9 +113,4 @@ psphotExtendedSourceAnalysis.c \ psphotExtendedSourceFits.c \ - psphotRadialProfile.c \ - psphotPetrosian.c \ - psphotIsophotal.c \ - psphotAnnuli.c \ - psphotKron.c \ psphotKernelFromPSF.c \ psphotPSFConvModel.c \ @@ -129,7 +133,25 @@ psphotThreadTools.c \ psphotAddNoise.c \ + psphotRadialProfile.c \ + psphotRadialProfileByAngles.c \ + psphotRadiiFromProfiles.c \ + psphotEllipticalContour.c \ + psphotEllipticalProfile.c \ + psphotPetrosian.c \ + psphotPetrosianRadialBins.c \ + psphotPetrosianStats.c \ psphotEfficiency.c -# dropped? psphotGrowthCurve.c +# re-instate these +# psphotIsophotal.c \ +# psphotAnnuli.c \ +# psphotKron.c \ +# psphotPetrosianVisual.c \ +# + +# test versions +# psphotPetrosianProfile.c \ +# psphotPetrosianAnalysis.c \ +# include_HEADERS = \ Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 25738) +++ trunk/psphot/src/psphot.h (revision 25755) @@ -58,5 +58,5 @@ bool psphotBlendFit_Threaded (psThreadJob *job); -psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections); +psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow); bool psphotSourceStats_Threaded (psThreadJob *job); @@ -71,5 +71,4 @@ bool psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf); -bool psphotApResidMags_Threaded (psThreadJob *job); bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); @@ -87,5 +86,5 @@ psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal); psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); -bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal); +bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); psErrorCode psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints); @@ -94,6 +93,6 @@ // used by ApResid -bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup); -bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag); +pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag); +bool psphotApResidMags_Threaded (psThreadJob *job); // basic support functions @@ -109,4 +108,5 @@ bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type); bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); +float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal); // output functions @@ -159,5 +159,5 @@ bool psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal); bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first); +bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first); bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal); @@ -167,9 +167,11 @@ psKernel *psphotKernelFromPSF (pmSource *source, int nPix); -bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); +// functions related to extended source analysis +bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); +bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax); +float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); +bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax); +bool psphotEllipticalProfile (pmSource *source); +bool psphotEllipticalContour (pmSource *source); // psphotVisual functions @@ -190,6 +192,29 @@ bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); bool psphotVisualShowResidualImage (pmReadout *readout); -bool psphotVisualPlotApResid (psArray *sources); -bool psphotVisualPlotSourceSize (psArray *sources); +bool psphotVisualPlotApResid (psArray *sources, float mean, float error); +bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources); +bool psphotVisualShowPetrosians (psArray *sources); + +// bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); +// bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise); + +bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); +bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise); +bool psphotPetrosianStats (pmSource *source); + +// XXX old versions, currently disabled +// bool psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); +// bool psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); +// bool psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); + +// XXX visualization functions related to radial profiles (disabled) +// bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux); +// bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); +// bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian); +// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, +// psVector *refRadius, psVector *meanSB, +// psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, +// float petRadius, float ratioForRadius, +// float petFlux, float radiusForFlux); bool psphotImageQuality (psMetadata *recipe, psArray *sources); Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 25738) +++ trunk/psphot/src/psphotAddNoise.c (revision 25755) @@ -50,4 +50,5 @@ // skip sources which were not subtracted + // NOTE: this bit is not modified when pmSourceOp applies to noise if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; @@ -79,5 +80,5 @@ axes.minor *= SIZE; newshape = psEllipseAxesToShape (axes); - PAR[PM_PAR_I0] = FACTOR*PS_SQR(oldI0); + PAR[PM_PAR_I0] = FACTOR*oldI0; PAR[PM_PAR_SXX] = newshape.sx; PAR[PM_PAR_SYY] = newshape.sy; Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 25738) +++ trunk/psphot/src/psphotApResid.c (revision 25755) @@ -33,10 +33,8 @@ if (!measureAptrend) { // save nan values since these were not calculated - psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation", NAN); psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); + psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); return true; } @@ -53,8 +51,10 @@ maskVal |= markVal; - // S/N limit to perform full non-linear fits + // clipping for extreme outliers + // XXX this is not currently defined in the recipe float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET"); - // this is the smallest radius allowed: need to at least extend growth curve down to this... + // options for how the photometry is calculated + // XXX are these sensible? bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH"); bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP"); @@ -100,4 +100,5 @@ PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nskip @@ -110,17 +111,4 @@ } psFree(job); - -# if (0) - int nskip = 0; - int nfail = 0; - - if (!psphotApResidMags_Unthreaded (&nskip, &nfail, sources, psf, photMode, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; - } - Nskip += nskip; - Nfail += nfail; -# endif - } @@ -138,7 +126,7 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[4]; + scalar = job->args->data[5]; Nskip += scalar->data.S32; - scalar = job->args->data[5]; + scalar = job->args->data[6]; Nfail += scalar->data.S32; } @@ -150,5 +138,4 @@ // gather the stats to assess the aperture residuals - psVector *mask = psVectorAllocEmpty (300, PS_TYPE_VECTOR_MASK); psVector *mag = psVectorAllocEmpty (300, PS_TYPE_F32); psVector *xPos = psVectorAllocEmpty (300, PS_TYPE_F32); @@ -158,4 +145,9 @@ Npsf = 0; +# ifdef DEBUG + FILE *f = fopen ("apresid.dat", "w"); + psAssert (f, "failed open"); +# endif + for (int i = 0; i < sources->n; i++) { source = sources->data[i]; @@ -170,8 +162,12 @@ if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED"); if (source->mode & PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY"); + if (source->mode & PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT"); if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { continue; } + + // XXX make this user-configurable? + if (source->errMag > 0.01) continue; // aperture residual for this source @@ -185,19 +181,24 @@ } - mag->data.F32[Npsf] = source->psfMag; - apResid->data.F32[Npsf] = dap; - xPos->data.F32[Npsf] = model->params->data.F32[PM_PAR_XPOS]; - yPos->data.F32[Npsf] = model->params->data.F32[PM_PAR_YPOS]; - - mask->data.PS_TYPE_VECTOR_MASK_DATA[Npsf] = 0; - - dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; - - psVectorExtend (mag, 100, 1); - psVectorExtend (mask, 100, 1); - psVectorExtend (xPos, 100, 1); - psVectorExtend (yPos, 100, 1); - psVectorExtend (dMag, 100, 1); - psVectorExtend (apResid, 100, 1); +# ifdef DEBUG + fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", + source->peak->xf, source->peak->yf, + source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], + source->psfMag, source->apMag, source->errMag, + source->modelPSF->params->data.F32[PM_PAR_I0], + source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], + source->modelPSF->params->data.F32[PM_PAR_7]); +# endif + if (!isfinite(source->psfMag)) psAbort ("nan in psfMag"); + if (!isfinite(source->errMag)) psAbort ("nan in errMag"); + if (!isfinite(source->apMag)) psAbort ("nan in apMag"); + if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos"); + if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos"); + + psVectorAppend (mag, source->psfMag); + psVectorAppend (dMag,source->errMag); + psVectorAppend (apResid, dap); + psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]); + psVectorAppend (yPos, model->params->data.F32[PM_PAR_YPOS]); Npsf ++; } @@ -205,4 +206,8 @@ psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n", Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail); + +# ifdef DEBUG + fclose (f); +# endif // XXX choose a better value here? @@ -212,47 +217,149 @@ } - // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux - // XXX is this asymmetric clipping still needed? this analysis should come after neighbors are subtracted... - // 3hi/1lo sigma clipping on the rflux vs metric fit - // systematic error information - float errorScale = 0.0; - float errorFloor = 0.0; - + // XXX set the min number of needed source more carefully + if ((Npsf < 15) && (APTREND_ORDER_MAX >= 4)) APTREND_ORDER_MAX = 3; + if ((Npsf < 11) && (APTREND_ORDER_MAX >= 3)) APTREND_ORDER_MAX = 2; + if ((Npsf < 8) && (APTREND_ORDER_MAX >= 2)) APTREND_ORDER_MAX = 1; + + psFree (psf->ApTrend); + psf->ApTrend = NULL; float errorFloorMin = FLT_MAX; - int entryMin = -1; - - // Fit out the dap vs mag trend, iterate over spatial scale until error Floor increases. - // Stop if Npsf / (Nx * Ny) < 3 + + // as we loop over orders, we need to refer to the initial selection, but we modify the + // option values to match the current guess: save the max values here: + int NX = readout->image->numCols; + int NY = readout->image->numRows; for (int i = 1; i <= APTREND_ORDER_MAX; i++) { - - if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) { - break; - } - - // store the resulting errorFloor values and the scales, redo the best + + int Nx, Ny; + if (NX > NY) { + Nx = i; + Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5)); + } else { + Ny = i; + Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5)); + } + + float errorFloor; + pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag); + if (!apTrend) { + continue; + } + + // apply ApTrend results + // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; + // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; + // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center + // if (!isfinite(ApResid)) psAbort("nan apresid @ center"); + + // store the minimum errorFloor and best ApTrend to keep if (errorFloor < errorFloorMin) { errorFloorMin = errorFloor; - entryMin = i; + psFree (psf->ApTrend); + psf->ApTrend = psMemIncrRefCounter(apTrend); } - } - if (entryMin == -1) { + psFree (apTrend); + } + if (psf->ApTrend == NULL) { psWarning("Failed to find a valid aperture residual value"); goto escape; } - // XXX catch error condition - psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag); + // apply ApTrend results + float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; + float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; + + psf->ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center + psf->dApResid = errorFloorMin; + psf->nApResid = Npsf; + + // save results for later output + psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid); + psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss); + + psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid); + psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid")); + + psFree (xPos); + psFree (yPos); + psFree (apResid); + psFree (mag); + psFree (dMag); + + psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid); + + return true; + +escape: + // save nan values since these were not calculated + psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); + psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); + psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); + psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); + + psFree (xPos); + psFree (yPos); + psFree (apResid); + psFree (mag); + psFree (dMag); + return false; +} + +pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) { + + // the mask marks the values not used to calculate the ApTrend + psVector *mask = psVectorAlloc(xPos->n, PS_TYPE_VECTOR_MASK); + psVectorInit (mask, 0); + + // XXX allow user to set this optionally? + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + // measure Trend2D for the current spatial scale + pmTrend2D *apTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); + + // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate: + // XXX use this or not? probably not, since this is the point of the systematic error analysis + psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32); + for (int i = 0; i < dMag->n; i++) { + dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.005); + } + + // XXX test for errors here + if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) { + psWarning("Failed to fit trend for %d x %d map", Nx, Ny); + psFree (apTrend); + return NULL; + } + if (apTrend->mode == PM_TREND_MAP) { + // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST: + psImageMapRepair (apTrend->map->map); + // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST: + } // construct the fitted values and the residuals - psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos); + psVector *apResidFit = pmTrend2DEvalVector (apTrend, mask, 0xff, xPos, yPos); psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit); - psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32)); - - if (psTraceGetLevel("psphot") >= 2) { - FILE *dumpFile = fopen ("apresid.dat", "w"); + + // measure systematic error + *apResidSysErr = psVectorSystematicError (apResidRes, dMag, 0.10); + if (!isfinite(*apResidSysErr)) { + psWarning("Failed to find systematic error for %d x %d map", Nx, Ny); + psFree (apTrend); + return NULL; + } + + psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny); + psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *apResidSysErr); + + if (psTraceGetLevel("psphot") >= 4) { + char filename[64]; + snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny); + FILE *dumpFile = fopen (filename, "w"); for (int i = 0; i < xPos->n; i++) { - fprintf (dumpFile, "%f %f %f %f %f %f %f %x\n", + fprintf (dumpFile, "%f %f %f %f %f %f %x\n", xPos->data.F32[i], yPos->data.F32[i], - mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i], + dMag->data.F32[i], hypot(dMag->data.F32[i], *apResidSysErr), apResid->data.F32[i], apResidRes->data.F32[i], mask->data.PS_TYPE_VECTOR_MASK_DATA[i]); @@ -261,204 +368,5 @@ } - // apply ApTrend results - float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; - float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; - - psf->ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center - psf->dApResid = errorFloor; - psf->nApResid = Npsf; - - // save results for later output - psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", 0.0); - psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation", 0.0); - psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); - psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); - psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss); - psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid); - - psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid); - psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid")); - - psFree (mag); psFree (mask); - psFree (xPos); - psFree (yPos); - - psFree (apResid); - psFree (apResidFit); - psFree (apResidRes); - - psFree (dMagSys); - psFree (dMag); - - psphotVisualPlotApResid (sources); - - return true; - -escape: - // save nan values since these were not calculated - psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); - - psFree (mag); - psFree (mask); - psFree (xPos); - psFree (yPos); - psFree (apResid); - psFree (dMag); - return false; -} - -/* - (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) - (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) - (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) -*/ - -// XXX this still sucks... need a better way to estimate the error floor... -bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) { - - psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN); - - // measure the trend in bins with 10 values each; if < 10 total, use them all - int nBin = PS_MAX (dMag->n / nGroup, 1); - - // output vectors for ApResid trend - psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32); - psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32); - psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32); - - // use dMag to group the dMag and dap vectors - psVector *index = psVectorSortIndex (NULL, dMag); - - // subset vectors for dMag and dap values within the given range - psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK); - - int n = 0; - for (int i = 0; i < dMo->n; i++) { - int j; - for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) { - int N = index->data.U32[n]; - dMSubset->data.F32[j] = dMag->data.F32[N]; - dASubset->data.F32[j] = dap->data.F32[N]; - mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N]; - } - dMSubset->n = j; - dASubset->n = j; - mkSubset->n = j; - - psStatsInit (statsS); - psStatsInit (statsM); - - if (j > 2) { - if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dSo->data.F32[i] = statsS->robustStdev; - dMo->data.F32[i] = statsM->sampleMean; - dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean; - } else { - dSo->data.F32[i] = NAN; - dMo->data.F32[i] = NAN; - dRo->data.F32[i] = NAN; - } - } - psFree (dMSubset); - psFree (dASubset); - psFree (mkSubset); - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - if (!psVectorStats (stats, dRo, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - - *errorScale = stats->sampleMedian; - for (int i = 0; i < dSo->n; i++) { - *errorFloor = dSo->data.F32[i]; - if (fabs(*errorFloor) <= FLT_EPSILON) continue; - if (isfinite(*errorFloor)) break; - } - - psFree (stats); - psFree (index); - - psFree (dRo); - psFree (dMo); - psFree (dSo); - - psFree (statsS); - psFree (statsM); - - return true; -} - -bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) { - - int Nx, Ny; - - if (readout->image->numCols > readout->image->numRows) { - Nx = scale; - float AR = readout->image->numRows / (float) readout->image->numCols; - Ny = (int) (Nx * AR + 0.5); - Ny = PS_MAX (1, Ny); - } else { - Ny = scale; - float AR = readout->image->numRows / (float) readout->image->numCols; - Nx = (int) (Ny * AR + 0.5); - Nx = PS_MAX (1, Nx); - } - - // require at least 10 stars per spatial bin - if (Npsf < 10*Nx*Ny) { - return false; - } - - // the mask marks the values not used to calculate the ApTrend - psVectorInit (mask, 0); - - // XXX stats structure for use by ApTrend : make parameters user setable - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - stats->min = 2.0; - stats->max = 3.0; - - // measure Trend2D for the current spatial scale - psFree (psf->ApTrend); - psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats); - - // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate: - psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32); - for (int i = 0; i < dMag->n; i++) { - dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01); - } - - // XXX test for errors here - pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft); - - // construct the fitted values and the residuals - psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos); - psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit); - - // measure systematic errorFloor & systematic / photon scale factor - // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin - int nGroup = PS_MAX (3*Nx*Ny, 10); - psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup); - - psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup); - psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale); - psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor); - psFree (stats); psFree (dMagSoft); @@ -466,5 +374,5 @@ psFree (apResidRes); - return true; + return apTrend; } @@ -480,4 +388,5 @@ pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[2],S32); psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { @@ -490,61 +399,42 @@ if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); - if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) { - Nskip ++; - psTrace ("psphot", 3, "skip : bad source mag"); - continue; + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } - if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { - Nfail ++; - psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); - continue; - } - source->mode |= PM_SOURCE_MODE_AP_MAGS; + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); + + bool status = pmSourceMagnitudes (source, psf, photMode, maskVal); + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + if (!status) { + Nskip ++; + psTrace ("psphot", 3, "skip : bad source mag"); + continue; + } + + if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { + Nfail ++; + psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); + continue; + } + source->mode |= PM_SOURCE_MODE_AP_MAGS; } // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[4]; + scalar = job->args->data[5]; scalar->data.S32 = Nskip; - scalar = job->args->data[5]; + scalar = job->args->data[6]; scalar->data.S32 = Nfail; return true; } - -# if (0) -bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) { - - int Nskip = 0; - int Nfail = 0; - - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR"); - if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR"); - if (source->mode & PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND"); - if (source->mode & PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR"); - if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); - - if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) { - Nskip ++; - psTrace ("psphot", 3, "skip : bad source mag"); - continue; - } - - if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { - Nfail ++; - psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); - continue; - } - } - - // change the value of a scalar on the array (wrap this and put it in psArray.h) - *nskip = Nskip; - *nfail = Nfail; - - return true; -} -# endif Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 25738) +++ trunk/psphot/src/psphotBlendFit.c (revision 25755) @@ -254,5 +254,4 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } @@ -365,5 +364,4 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } @@ -374,4 +372,8 @@ *nfail = Nfail; + // moments are modified by the fit; re-display + psphotVisualPlotMoments (recipe, sources); + psphotVisualShowResidualImage (readout); + return true; } Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 25738) +++ trunk/psphot/src/psphotChoosePSF.c (revision 25755) @@ -1,18 +1,3 @@ # include "psphotInternal.h" - -void psphotCountPSFStars (psArray *sources) { - - int nPSF = 0; - - // select the candidate PSF stars (pointers to original sources) - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - nPSF ++; - } - } - - fprintf (stderr, "N PSF: %d\n", nPSF); -} // try PSF models and select best option @@ -73,6 +58,18 @@ // get the fixed PSF fit radius // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS - options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); - assert (status); + // options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + // assert (status); + + // We have calculated a Gaussian window function, use that for both the PSF fit radius and + // the aperture radius (scaling SIGMA) + float gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE"); + float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE"); + options->fitRadius = (int)(fitScale*gaussSigma); + options->apRadius = (int)(apScale*gaussSigma); + + // XXX use the same radii for standard analysis as for the PSF creation + psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius); + psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius); // XXX ROBUST seems to be too agressive given the small numbers. @@ -99,6 +96,4 @@ psArray *stars = psArrayAllocEmpty (sources->n); - - // psphotCountPSFStars (sources); // select the candidate PSF stars (pointers to original sources) @@ -115,11 +110,7 @@ } - // psphotCountPSFStars (sources); - // check that the identified psf stars sufficiently cover the region; if not, extend the // limits somewhat psphotCheckStarDistribution (stars, sources, options); - - // psphotCountPSFStars (sources); psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n); @@ -289,10 +280,8 @@ // XXX test dump of psf star data and psf-subtracted image if (psTraceGetLevel("psphot.psfstars") > 5) { - psphotDumpPSFStars (readout, try, options->radius, maskVal, markVal); + psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal); } // save only the best model; - // XXX we are not saving the fitted sources - // XXX do we want to keep them so we may optionally write them out? pmPSF *psf = psMemIncrRefCounter(try->psf); psFree (models); @@ -305,6 +294,4 @@ } - // psphotCountPSFStars (sources); - char *modelName = pmModelClassGetName (psf->type); psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf")); @@ -341,5 +328,8 @@ // create modelPSF from this model pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0); - if (!modelPSF) continue; + if (!modelPSF) { + fprintf (stderr, "?"); + continue; + } // get the model full-width at half-max @@ -354,5 +344,8 @@ float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); - if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue; + if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) { + fprintf (stderr, "!"); + continue; + } psVectorAppend (fwhmMajor, FWHM_MAJOR); psVectorAppend (fwhmMinor, FWHM_MINOR); Index: trunk/psphot/src/psphotEllipticalContour.c =================================================================== --- trunk/psphot/src/psphotEllipticalContour.c (revision 25755) +++ trunk/psphot/src/psphotEllipticalContour.c (revision 25755) @@ -0,0 +1,174 @@ +# include "psphotInternal.h" + +// model parameters +enum {PAR_PHI, PAR_EPSILON, PAR_RMIN}; +psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord); + +bool psphotEllipticalContour (pmSource *source) { + + pmSourceRadialProfile *profile = source->extpars->profile; + + // use LMM to fit theta vs radius to an ellipse + psVector *theta = profile->theta; + psVector *radius = profile->isophotalRadii; + + // find Rmin and Rmax for the initial guess + float Rmin = radius->data.F32[0]; + float Rmax = radius->data.F32[0]; + + // arrays to hold the data to be fitted + // we fit x and y vs theta in separate passes. + psArray *x = psArrayAllocEmpty(2*radius->n); + psVector *y = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32); + psVector *yErr = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32); + + int n = 0; + for (int i = 0; i < radius->n; i++) { + if (!isfinite(radius->data.F32[i])) continue; + + psVector *coord = NULL; + + // Rx coordinate value + coord = psVectorAlloc (2, PS_TYPE_F32); + coord->data.F32[1] = 0.0; + coord->data.F32[0] = theta->data.F32[i]; + x->data[n] = coord; + y->data.F32[n] = radius->data.F32[i]*cos(theta->data.F32[i]); + yErr->data.F32[n] = 1000.0; + n++; + + // Ry coordinate value + coord = psVectorAlloc (2, PS_TYPE_F32); + coord->data.F32[1] = 1.0; + coord->data.F32[0] = theta->data.F32[i]; + x->data[n] = coord; + y->data.F32[n] = radius->data.F32[i]*sin(theta->data.F32[i]); + yErr->data.F32[n] = 1000.0; + n++; + + // check the radius range + Rmin = MIN (Rmin, radius->data.F32[i]); + Rmax = MAX (Rmax, radius->data.F32[i]); + } + x->n = n; + y->n = n; + yErr->n = n; + + if (n < 4) { + psFree (x); + psFree (y); + psFree (yErr); + return false; + } + + psVector *params = psVectorAlloc (3, PS_TYPE_F32); + + // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 7); + + // create the minimization constraints + psMinConstraint *constraint = psMinConstraintAlloc(); + + // XXX for now, no parameter masks, skip checkLimits + // XXX might help to add a limit to the angle or re-parameterize the ellipse in terms of Rxx, Ryy, Rxy + // constraint->checkLimits = psastroModelBoresiteLimits; + + params->data.F32[PAR_PHI] = 0.0; + params->data.F32[PAR_EPSILON] = Rmin / Rmax; + params->data.F32[PAR_RMIN] = Rmin; + + psMinimization *myMin = psMinimizationAlloc (25, 0.001); + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); + + // XXX skip the weights for now + psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, psphotEllipticalContourFunc); + + /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees) + if (params->data.F32[PAR_EPSILON] < 1.0) { + profile->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + profile->axes.minor = params->data.F32[PAR_RMIN]; + profile->axes.theta = params->data.F32[PAR_PHI]; + } else { + profile->axes.major = params->data.F32[PAR_RMIN]; + profile->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + profile->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI; + } + + psTrace ("psphot", 4, "# fitted values:\n"); + psTrace ("psphot", 4, "Phi: %f\n", profile->axes.theta*PS_DEG_RAD); + psTrace ("psphot", 4, "Rmaj: %f\n", profile->axes.major); + psTrace ("psphot", 4, "Rmin: %f\n", profile->axes.minor); + + // show the results + // psphotPetrosianVisualEllipticalContour (petrosian); + + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree (covar); + psFree (myMin); + psFree (constraint); + + return true; +} + +/** + * the full chisq is built of two associated sums over coordinates: + * chisq = sum ((Rx_obs - Rx_fit(t))^2 + (Ry_obs - Ry_fit(t))^2) + * we use split this into a 2x long vector and use coord[1] to distinguish the X and Y terms: + * coord[0] = measured X or measured Y + * coord[1] = 0 or 1 + */ +psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) { + + static int pass = 0; + + psF32 *par = params->data.F32; + + float alpha = coord->data.F32[0]; + + float cs_alpha = cos(alpha); + float sn_alpha = sin(alpha); + + float cs_phi = cos(alpha - par[PAR_PHI]); + float sn_phi = sin(alpha - par[PAR_PHI]); + + float r = 1.0 / sqrt(SQ(sn_phi) + SQ(par[PAR_EPSILON]*cs_phi)); + float r3 = pow(r, 3.0); + float drdE = -0.5 * r3 * SQ(cs_phi) * 2.0 * par[PAR_EPSILON]; + float drdP = -0.5 * r3 * (SQ(par[PAR_EPSILON]) - 1) * 2.0 * cs_phi * sn_phi; + + // value is X + // if (coord->data.F32[1] == 0) { + if (pass == 0) { + pass = 1; + + float value = par[PAR_RMIN]*cs_alpha*r; + + if (deriv) { + psF32 *dPAR = deriv->data.F32; + dPAR[PAR_RMIN] = r*cs_alpha; + dPAR[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE; + dPAR[PAR_PHI] = 4.0*par[PAR_RMIN]*cs_alpha*drdP; + } + return (value); + } + + // value is Y + // if (coord->data.F32[1] == 1) { + if (pass == 1) { + pass = 0; + + float value = par[PAR_RMIN]*sn_alpha*r; + + if (deriv) { + psF32 *dPAR = deriv->data.F32; + dPAR[PAR_RMIN] = r*sn_alpha; + dPAR[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE; + dPAR[PAR_PHI] = 4.0*par[PAR_RMIN]*sn_alpha*drdP; + } + return (value); + } + + psAbort ("programming error: invalid coordinate"); +} Index: trunk/psphot/src/psphotEllipticalProfile.c =================================================================== --- trunk/psphot/src/psphotEllipticalProfile.c (revision 25755) +++ trunk/psphot/src/psphotEllipticalProfile.c (revision 25755) @@ -0,0 +1,75 @@ +# include "psphotInternal.h" + +bool psphotEllipticalProfile (pmSource *source) { + + pmSourceRadialProfile *profile = source->extpars->profile; + + profile->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); + profile->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); + + psVector *radius = profile->radiusElliptical; + psVector *flux = profile->fluxElliptical; + + // the psEllipse functions use z = 0.5(x/Sxx)^2 + 0.5(y/Syy)^2 + x y Sxy + // which are converted to z = 0.5(x/a)^2 + 0.5(y/b)^2 + // we have major and minor axes of a specific ellipse with r^2 = (x/A)^2 + (y/B)^2 + // a = A / sqrt(2) + + // we have the shape parameters of the elliptical contour at the reference isophote. + // use the axis ratio (major/minor) to rescale the radial profile so that 1 pixel + // along the major axis is 1 pixel, and a smaller amount on the minor axis + + psEllipseAxes axes; + axes.major = M_SQRT1_2; + axes.minor = M_SQRT1_2 * (profile->axes.minor / profile->axes.major); + + // axes.major = 1.0; + // axes.minor = profile->axes.minor / profile->axes.major; + + axes.theta = profile->axes.theta; + psEllipseShape shape = psEllipseAxesToShape (axes); + + float Sxx = shape.sx; + float Sxy = shape.sxy; + float Syy = shape.sy; + + // XXX drop these two vectors? + psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32); + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + + // 0.5 PIX: get radius as a function of pixel coord + float x = ix + 0.5 - source->peak->xf + source->pixels->col0; + float y = iy + 0.5 - source->peak->yf + source->pixels->row0; + + float r2 = 0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy; + + psVectorAppend(radius, sqrt(r2)); + psVectorAppend(flux, source->pixels->data.F32[iy][ix]); + + float Rraw = hypot(x, y); + psVectorAppend(radiusRaw, Rraw); + psVectorAppend(fluxRaw, source->pixels->data.F32[iy][ix]); + } + } + + // psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32); + // psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32); + // for (int i = 0; i < profile->radii->n; i++) { + // psVector *r = profile->radii->data[i]; + // psVector *f = profile->fluxes->data[i]; + // for (int j = 0; j < r->n; j++) { + // psVectorAppend(radiusRaw, r->data.F32[j]); + // psVectorAppend(fluxRaw, f->data.F32[j]); + // } + // } + + // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0); + // psphotPetrosianVisualProfileByAngle (radius, flux); + + psFree (radiusRaw); + psFree (fluxRaw); + return true; +} Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 25738) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 25755) @@ -21,8 +21,15 @@ } - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + // XXX require petrosian analysis for non-linear fits? + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + +# if (0) + // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN + // XXX use this to set skynoise + pmReadout *backModel = psphotSelectBackground (config, view); + pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); +# endif // S/N limit to perform full non-linear fits @@ -38,5 +45,8 @@ sources = psArraySort (sources, pmSourceSortBySN); - // XXX some init functions for the extended source recipe options? + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); // choose the sources of interest @@ -49,4 +59,7 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; // limit selection to some SN limit @@ -68,10 +81,9 @@ // if we request any of these measurements, we require the radial profile if (doPetrosian || doIsophotal || doAnnuli || doKron) { - if (!psphotRadialProfile (source, recipe, maskVal)) { + if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { // all measurements below require the radial profile; skip them all // re-subtract the object, leave local sky psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; continue; } @@ -79,4 +91,16 @@ } + // Petrosian Mags + if (doPetrosian) { + if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + Npetro ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; + } + } + +# if (0) // Isophotal Mags if (doIsophotal) { @@ -89,16 +113,4 @@ } } - - // Petrosian Mags - if (doPetrosian) { - if (!psphotPetrosian (source, recipe, maskVal)) { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - } else { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Npetro ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } - } - // Kron Mags if (doKron) { @@ -111,19 +123,12 @@ } } - - // Radial Annuli - if (doAnnuli) { - if (!psphotAnnuli (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis"); - return false; - } - psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Nannuli ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } +# endif // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + + if (source->extpars) { + pmSourceRadialProfileFreeVectors(source->extpars->profile); + } } @@ -133,4 +138,11 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); + + psphotVisualShowResidualImage (readout); + + if (doPetrosian) { + psphotVisualShowPetrosians (sources); + } + return true; } Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 25738) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 25755) @@ -226,5 +226,4 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; psFree (modelFluxes); @@ -253,5 +252,4 @@ // subtract the best fit from the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; // the initial model flux is no longer needed Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 25738) +++ trunk/psphot/src/psphotFindDetections.c (revision 25755) @@ -52,5 +52,5 @@ // optionally merge peaks into footprints if (useFootprints) { - psphotFindFootprints (detections, significance, readout, recipe, pass, maskVal); + psphotFindFootprints (detections, significance, readout, recipe, threshold, pass, maskVal); } Index: trunk/psphot/src/psphotFindFootprints.c =================================================================== --- trunk/psphot/src/psphotFindFootprints.c (revision 25738) +++ trunk/psphot/src/psphotFindFootprints.c (revision 25755) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psImageMaskType maskVal) { +bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) { bool status; @@ -9,17 +9,4 @@ int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN"); PS_ASSERT (status, false); - - float FOOTPRINT_NSIGMA_LIMIT; - if (pass == 1) { - FOOTPRINT_NSIGMA_LIMIT = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NSIGMA_LIMIT"); - } else { - FOOTPRINT_NSIGMA_LIMIT = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NSIGMA_LIMIT_2"); - } - PS_ASSERT (status, false); - - // XXX do we need to use the same threshold here as for peaks? does it make sense for - // these to be different? - - float threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT); int growRadius = 0; Index: trunk/psphot/src/psphotFindPeaks.c =================================================================== --- trunk/psphot/src/psphotFindPeaks.c (revision 25738) +++ trunk/psphot/src/psphotFindPeaks.c (revision 25755) @@ -31,4 +31,10 @@ peak->SN = sqrt(peak->value); peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; + // if (peak->flux / peak->value > 5.0/12.0) { + // psWarning ("odd peak levels (1)"); + // } + // if (peak->value / peak->flux > 5*12.0) { + // psWarning ("odd peak levels (2)"); + // } if (readout->variance && isfinite (peak->dx)) { peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]); Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 25738) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 25755) @@ -76,5 +76,5 @@ if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; } else { - if (source->mode & PM_SOURCE_MODE_BLEND) continue; + // if (source->mode & PM_SOURCE_MODE_BLEND) continue; } @@ -186,5 +186,5 @@ if (SKY_FIT_LINEAR) { psSparseBorderSolve (&norm, &skyfit, constraint, border, 5); - fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); + psLogMsg ("psphot", PS_LOG_MINUTIA, "skyfit: %f\n", skyfit->data.F32[0]); } else { norm = psSparseSolve (NULL, constraint, sparse, 5); @@ -215,5 +215,4 @@ // subtract object pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -239,5 +238,5 @@ psphotVisualShowResidualImage (readout); - psphotVisualShowFlags (sources); + // psphotVisualShowFlags (sources); return true; @@ -264,5 +263,5 @@ 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)); + psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); } Index: trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 25738) +++ trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 25755) @@ -168,5 +168,4 @@ // subtract object pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 25738) +++ trunk/psphot/src/psphotGuessModels.c (revision 25755) @@ -177,9 +177,9 @@ // set the source PSF model + psAssert (source->modelPSF == NULL, "failed to free one of the models?"); source->modelPSF = modelPSF; source->modelPSF->residuals = psf->residuals; pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) - } Index: trunk/psphot/src/psphotImageLoop.c =================================================================== --- trunk/psphot/src/psphotImageLoop.c (revision 25738) +++ trunk/psphot/src/psphotImageLoop.c (revision 25755) @@ -67,18 +67,19 @@ // Update the header - { - pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); - if (hdu && hdu != lastHDU) { - psphotVersionHeaderFull(hdu->header); - lastHDU = hdu; - } + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; } - psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels - if (!pmReadoutMaskNonfinite(readout, maskSat)) { - psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); - psFree(view); - return false; - } + // if an external mask is supplied, ensure that NAN pixels are also masked + if (readout->mask) { + psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels + if (!pmReadoutMaskNonfinite(readout, maskSat)) { + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); + psFree(view); + return false; + } + } // run the actual photometry analysis on this chip/cell/readout Index: trunk/psphot/src/psphotImageQuality.c =================================================================== --- trunk/psphot/src/psphotImageQuality.c (revision 25738) +++ trunk/psphot/src/psphotImageQuality.c (revision 25755) @@ -279,8 +279,8 @@ #endif - psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) : %f, %f\n", + psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n", M2->n, fwhm_major, fwhm_minor); - psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n", + psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f [pixels^n]\n", vM2, dM2, vM3, dM3, vM4, dM4); Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 25738) +++ trunk/psphot/src/psphotMagnitudes.c (revision 25755) @@ -71,4 +71,5 @@ PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for nAp @@ -102,5 +103,5 @@ fprintf (stderr, "error with job\n"); } else { - psScalar *scalar = job->args->data[7]; + psScalar *scalar = job->args->data[8]; Nap += scalar->data.S32; } @@ -127,9 +128,26 @@ pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32); psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { pmSource *source = (pmSource *) sources->data[i]; - status = pmSourceMagnitudes (source, psf, photMode, maskVal); + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); + + status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal if (status && isfinite(source->apMag)) Nap ++; + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); if (backModel) { @@ -155,5 +173,5 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - psScalar *scalar = job->args->data[7]; + psScalar *scalar = job->args->data[8]; scalar->data.S32 = Nap; Index: trunk/psphot/src/psphotMakeFluxScale.c =================================================================== --- trunk/psphot/src/psphotMakeFluxScale.c (revision 25738) +++ trunk/psphot/src/psphotMakeFluxScale.c (revision 25755) @@ -60,4 +60,9 @@ goto DONE; } + if (trend->mode == PM_TREND_MAP) { + // p_psImagePrint (2, trend->map->map, "FluxScale Before"); // XXX TEST: + psImageMapRepair (trend->map->map); + // p_psImagePrint (2, trend->map->map, "FluxScale After"); // XXX TEST: + } // XXX do something useful to measure residual statistics Index: trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- trunk/psphot/src/psphotMakeResiduals.c (revision 25738) +++ trunk/psphot/src/psphotMakeResiduals.c (revision 25755) @@ -1,3 +1,5 @@ # include "psphotInternal.h" + +# define RESIDUAL_SOFTENING 0.005 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal) { @@ -31,4 +33,7 @@ float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN"); + PS_ASSERT (status, false); + + float radiusMax = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.RADIUS"); PS_ASSERT (status, false); @@ -171,5 +176,4 @@ bool offImage = false; if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) { - // fprintf (stderr, "off image: %f %f : %f %f\n", ix, iy, flux, dflux); // This pixel is off the image offImage = true; @@ -179,6 +183,9 @@ } fluxes->data.F32[i] = flux; - dfluxes->data.F32[i] = dflux; + dfluxes->data.F32[i] = hypot(dflux, RESIDUAL_SOFTENING); if (isnan(flux)) { + fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask; + } + if (isnan(dflux)) { fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask; } @@ -234,4 +241,10 @@ } + float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); + if (radius > radiusMax) { + resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; + continue; + } + resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption); resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; @@ -248,9 +261,13 @@ resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; } - - // fprintf (stderr, "res: %2d %2d : %6.4f %6.4f %6.4f %3d %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); - } else { assert (SPATIAL_ORDER == 1); + + float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); + if (radius > radiusMax) { + resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; + continue; + } + psImageInit(A, 0.0); psVectorInit(B, 0.0); @@ -275,8 +292,7 @@ if (!psMatrixGJSolve(A, B)) { - psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (y,x) = (%d,%d)'s residuals", - oy, ox); - psFree(resid); resid = NULL; - break; + resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; + psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox); + continue; } @@ -286,11 +302,8 @@ float dRo = sqrt(A->data.F32[0][0]); - // fprintf (stderr, "res: %2d %2d : %6.4f %6.4f %6.4f %3d %1d\n", - // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) { resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; } - //resid->variance->data.F32[oy][ox] = XXX; } } Index: trunk/psphot/src/psphotMaskReadout.c =================================================================== --- trunk/psphot/src/psphotMaskReadout.c (revision 25738) +++ trunk/psphot/src/psphotMaskReadout.c (revision 25755) @@ -33,13 +33,17 @@ } + bool softenVariance = psMetadataLookupBool (&status, recipe, "SOFTEN.VARIANCE"); + float softenFraction = psMetadataLookupF32 (&status, recipe, "SOFTEN.VARIANCE.FRACTION"); + // make this an option via the recipe - if (0) { + if (softenVariance) { psImage *im = readout->image; psImage *wt = readout->variance; - psImage *mk = readout->mask; for (int j = 0; j < im->numRows; j++) { for (int i = 0; i < im->numCols; i++) { - if (isfinite(im->data.F32[j][i]) && isfinite(wt->data.F32[j][i])) continue; - mk->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= maskBad; + if (!isfinite(im->data.F32[j][i])) continue; + if (!isfinite(wt->data.F32[j][i])) continue; + float sysError = softenFraction * im->data.F32[j][i]; + wt->data.F32[j][i] += PS_SQR(sysError); } } Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 25738) +++ trunk/psphot/src/psphotOutput.c (revision 25755) @@ -31,4 +31,36 @@ } return background; +} + +// dump source stats for psf stars +bool psphotDumpStats (psArray *sources, char *stage) { + + char filename[64]; + snprintf (filename, 64, "psf.%s.dat", stage); + FILE *f = fopen (filename, "w"); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + pmModel *model = source->modelPSF; + if (!model) continue; + + // int xc = source->peak->x - source->pixels->col0; + // int yc = source->peak->y - source->pixels->row0; + // float mcore = source->modelFlux ? source->modelFlux->data.F32[yc][xc] : NAN; + // float mpeak = model ? model->params->data.F32[PM_PAR_I0] : NAN; + // bool subtracted = source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED; + // fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->flux, source->pixels->data.F32[yc][xc], mcore, mpeak); + + fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", + source->peak->xf, source->peak->yf, + model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], + source->psfMag, source->apMag, source->errMag, + model->params->data.F32[PM_PAR_I0], + model->params->data.F32[PM_PAR_SXX], model->params->data.F32[PM_PAR_SXY], model->params->data.F32[PM_PAR_SYY], + model->params->data.F32[PM_PAR_7]); + } + fclose (f); + return true; } @@ -160,6 +192,4 @@ psMetadataItemSupplement (header, recipe, "DAPMIFIT"); psMetadataItemSupplement (header, recipe, "NAPMIFIT"); - psMetadataItemSupplement (header, recipe, "SKYBIAS"); - psMetadataItemSupplement (header, recipe, "SKYSAT"); // PSF model parameters (shape values for image center) @@ -256,5 +286,5 @@ psImageKeepCircle (source->maskObj, x, y, radius, "OR", markVal); pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal); - psImageKeepCircle (source->maskObj, x, y, radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); } Index: trunk/psphot/src/psphotPSFConvModel.c =================================================================== --- trunk/psphot/src/psphotPSFConvModel.c (revision 25738) +++ trunk/psphot/src/psphotPSFConvModel.c (revision 25755) @@ -37,4 +37,8 @@ } + // adjust the pixels based on the footprint + float radius = psphotSetRadiusEXT (readout, source, markVal); + if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false; + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: psEllipseShape psfShape; @@ -67,6 +71,4 @@ psVector *params = modelConv->params; psVector *dparams = modelConv->dparams; - - psphotCheckRadiusEXT (readout, source, modelConv, markVal); // create the minimization constraints Index: trunk/psphot/src/psphotPetrosian.c =================================================================== --- trunk/psphot/src/psphotPetrosian.c (revision 25738) +++ trunk/psphot/src/psphotPetrosian.c (revision 25755) @@ -1,109 +1,32 @@ # include "psphotInternal.h" -bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) { +bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) { - bool status; + // XXX these need to go into recipe values + float Rmax = 200; - assert (source->extpars); - assert (source->extpars->profile); - assert (source->extpars->profile->radius); - assert (source->extpars->profile->flux); + psAssert (source->extpars, "need to run psphotRadialProfile first"); + psAssert (source->extpars->profile, "need to run psphotRadialProfile first"); - psVector *radius = source->extpars->profile->radius; - psVector *flux = source->extpars->profile->flux; + // integrate the radial profile for radial bins defined for the petrosian measurement: + // SB_i (r_i) where \alpha r_i < r < \beta r_i + if (!psphotPetrosianRadialBins (source, Rmax, skynoise)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; + } + + // use the SB_i from above to calculate the petrosian radius and the flux within that radius + if (!psphotPetrosianStats (source)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; + } + + psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f", + source->peak->xf, source->peak->yf, + source->extpars->petrosian_80->radius, + source->extpars->petrosian_80->flux, + source->extpars->profile->axes.minor/source->extpars->profile->axes.major, + source->extpars->profile->axes.theta*PS_DEG_RAD); - // flux at which to measure isophotal parameters - float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); - float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); - assert (status); - - // first find flux at R0 - int firstAbove = -1; - int lastBelow = -1; - for (int i = 0; i < radius->n; i++) { - if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i; - if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i; - } - // if we don't go out far enough, we have a problem... - if (lastBelow == radius->n - 1) { - psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius..."); - // XXX skip object? raise a flag ? - return false; - } - if (firstAbove < 0) { - psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius"); - // XXX raise a flag ? - return false; - } - - // average flux in this range - float fluxR0 = 0.0; - int fluxRn = 0; - for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) { - fluxR0 += flux->data.F32[i]; - fluxRn ++; - } - fluxR0 /= (float)(fluxRn); - - // target flux for petrosian radius - float fluxRP = fluxR0 * PETROSIAN_RF; - - // find the first bin below the flux level and the last above the level - // XXX can this be done faster with bisection? - // XXX do I need to worry about crazy outliers? - // XXX should i be smoothing or fitting the curve? - int firstBelow = -1; - int lastAbove = -1; - for (int i = 0; i < flux->n; i++) { - if (flux->data.F32[i] > fluxRP) lastAbove = i; - if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i; - } - // if we don't go out far enough, we have a problem... - if (lastAbove == radius->n - 1) { - psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius..."); - // XXX skip object? raise a flag ? - return false; - } - if (firstBelow < 0) { - psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius"); - // XXX raise a flag ? - return false; - } - - // need to examine pixels in this vicinity - float fluxFirst = 0; - float fluxLast = 0; - for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) { - if (i <= firstBelow) { - fluxFirst += flux->data.F32[i]; - } - if (i <= lastAbove) { - fluxLast += flux->data.F32[i]; - } - } - float fluxRPSum = 0.5*(fluxLast + fluxFirst); - float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst); - // XXX need to use the weight appropriately here... - - float rad = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]); - float radErr = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]); - - if (!source->extpars->petrosian) { - source->extpars->petrosian = pmSourcePetrosianValuesAlloc (); - } - - // these are uncalibrated: instrumental mags and pixel units - source->extpars->petrosian->mag = -2.5*log10(fluxRPSum); - source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum; - - source->extpars->petrosian->rad = rad; - source->extpars->petrosian->radErr = radErr; - - psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n", - source->extpars->petrosian->mag, source->extpars->petrosian->magErr, - source->extpars->petrosian->rad, source->extpars->petrosian->radErr, - source->peak->xf, source->peak->yf); - - return true; - + return true; } Index: trunk/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- trunk/psphot/src/psphotPetrosianAnalysis.c (revision 25755) +++ trunk/psphot/src/psphotPetrosianAnalysis.c (revision 25755) @@ -0,0 +1,68 @@ +# include "psphotInternal.h" + +// aperture-like measurements for extended sources +bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + + // 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); + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + +# if (0) + // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN + // XXX use this to set skynoise + pmReadout *backModel = psphotSelectBackground (config, view); + pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); +# endif + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) continue; + + // limit selection by analysis region + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.y1) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + psphotPetrosianProfile (readout, source, skynoise); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + } + + psphotVisualShowResidualImage (readout); + return true; +} Index: trunk/psphot/src/psphotPetrosianProfile.c =================================================================== --- trunk/psphot/src/psphotPetrosianProfile.c (revision 25755) +++ trunk/psphot/src/psphotPetrosianProfile.c (revision 25755) @@ -0,0 +1,80 @@ +# include "psphotInternal.h" + +// generate the Petrosian radius and flux using elliptical contours + +// XXX much of this function is focused on generating the clean contours, which can be used by +// any number of aperture-like measurements. probably will want to rename the pmPetrosian +// structure to something the pmRadialProfile + +bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise) { + + // container to hold results from the radial profile analysis + pmPetrosian *petrosian = pmPetrosianAlloc(); + + // XXX these need to go into recipe values + int Nsec = 24; + float Rmax = 200; + float fluxMin = 0.0; + float fluxMax = source->peak->flux; + + // generate a series of radial profiles at Nsec evenly spaced angles. the profile flux + // is measured by interpolation for small radii; for large radii, the pixels in a box + // are averaged to increase the S/N (XXX not yet done) + if (!psphotRadialProfilesByAngles (source, petrosian, Nsec, Rmax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); + psFree (petrosian); + return false; + } + + // use the radial profiles to determine the radius of a given isophote. this isophote + // is used to determine the elliptical shape of the object, so it has a relatively high + // value (nominally 50% of the peak) + if (!psphotRadiiFromProfiles (source, petrosian, fluxMin, fluxMax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); + psFree (petrosian); + return false; + } + + // convert the isophotal radius vs angle measurements to an elliptical contour + if (!psphotEllipticalContour (source, petrosian)) { + psLogMsg ("psphot", 3, "failed to measure elliptical contour"); + psFree (petrosian); + return false; + } + + // generate a single, normalized radial profile following the elliptical contours. + // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel + if (!psphotEllipticalProfile (source, petrosian)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + psFree (petrosian); + return false; + } + + // integrate the radial profile for radial bins defined for the petrosian measurement: + // SB_i (r_i) where \alpha r_i < r < \beta r_i + if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + psFree (petrosian); + return false; + } + + // use the SB_i from above to calculate the petrosian radius and the flux within that radius + if (!psphotPetrosianStats (source, petrosian)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + psFree (petrosian); + return false; + } + + // XXX this will only work in the psphot context, not the psphotPetrosianStudy... + // XXX add the petrosian to the pmSource structure... + // psphotVisualShowResidualImage (readout); + psphotVisualShowPetrosian (source, petrosian); + + psphotPetrosianFreeVectors(petrosian); + + psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f", + source->peak->xf, source->peak->yf, petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta); + + psFree (petrosian); + return true; +} Index: trunk/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- trunk/psphot/src/psphotPetrosianRadialBins.c (revision 25755) +++ trunk/psphot/src/psphotPetrosianRadialBins.c (revision 25755) @@ -0,0 +1,189 @@ +# include "psphotInternal.h" + +// convert the flux vs elliptical radius to annular bins + +// we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels +// themselves. this function does not attempt to measure the radial profiles accurately +// for radii that are smaller than a minimum (currently 1.0 pixels). + +// for small radii, we are measuring the mean surface brightness in non-overlapping radial +// bins. for large radii (r > 2 pixels), we are measuring the surface brightness for a +// radius range \alpha r_i < i < \beta r_i, but performing this measurement for radii more +// finely spaced than r_{i+1} = r_i * \beta / \alpha. for the integration, we need to +// track the non-overlapping radius values. + +// XXX move the resulting elements from profile to extpars->petrosian? +bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise) { + + pmSourceRadialProfile *profile = source->extpars->profile; + + float skyModelErrorSQ = PS_SQR(skynoise); + + psVector *radius = profile->radiusElliptical; + psVector *flux = profile->fluxElliptical; + + // sort incoming vectors by radius + pmSourceRadialProfileSortPair (radius, flux); + + int nMax = radiusMax; + + // radBin stores the centers of the radial bins, + // radMin, radMax store the bounds + psVector *radMin = psVectorAllocEmpty(nMax, PS_TYPE_F32); + psVector *radMax = psVectorAllocEmpty(nMax, PS_TYPE_F32); + psVector *radAlp = psVectorAllocEmpty(nMax, PS_TYPE_F32); + psVector *radBet = psVectorAllocEmpty(nMax, PS_TYPE_F32); + + psVector *binSB = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin + psVector *binSBstdev = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin + psVector *binRad = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin + psVector *binArea = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping) + + psVectorInit (binSB, 0.0); + psVectorInit (binSBstdev, 0.0); + psVectorInit (binRad, 0.0); + + // generate radial bin bounds + radMin->data.F32[0] = 0.0; + radMax->data.F32[0] = 1.0; + radAlp->data.F32[0] = 0.0; + radBet->data.F32[0] = 1.0; + + radMin->data.F32[1] = 1.0; + radMax->data.F32[1] = 1.5; + radAlp->data.F32[1] = 1.0; + radBet->data.F32[1] = 1.5; + + radMin->data.F32[2] = 1.5; + radMax->data.F32[2] = 2.0; + radAlp->data.F32[2] = 1.5; + radBet->data.F32[2] = 2.0; + +# define PETROSIAN_ALPHA 0.8 +# define PETROSIAN_BETA 1.25 +# define POWER_LAW_SPACING true + + // power-law spacing with overlapping boundaries at the geometric mid-points + float rBeta = sqrt(PETROSIAN_BETA); + for (int i = 3; radBet->data.F32[i-1] < radiusMax; i++) { + if (POWER_LAW_SPACING) { + radMin->data.F32[i] = radMax->data.F32[i-1]; + radMax->data.F32[i] = radMin->data.F32[i] * PETROSIAN_BETA; + radAlp->data.F32[i] = radMin->data.F32[i] / rBeta; + radBet->data.F32[i] = radMax->data.F32[i] * rBeta; + } else { + radMin->data.F32[i] = radMax->data.F32[i-1]; + radMax->data.F32[i] = radMin->data.F32[i] + 1; + float rMid = 0.5*(radMin->data.F32[i] + radMax->data.F32[i]); + radAlp->data.F32[i] = rMid * PETROSIAN_ALPHA; + radBet->data.F32[i] = rMid * PETROSIAN_BETA; + } + radMin->n = radMax->n = radAlp->n = radBet->n = i + 1; + } + + // generate radial area-weighted mean radius & non-overlapping areas + for (int i = 0; i < radMin->n; i++) { + float rMin = radMin->data.F32[i]; + float rMax = radMax->data.F32[i]; + + float rMin2 = rMin*rMin; + float rMin3 = rMin2*rMin; + + float rMax2 = rMax*rMax; + float rMax3 = rMax2*rMax; + + float rBin = 2.0 * (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0; + + // XXX calculate area-weighted radius rather than asserting? + binRad->data.F32[i] = rBin; + binArea->data.F32[i] = M_PI * (rMax2 - rMin2); + + psTrace ("psphot", 6, "%3d %5.1f %5.1f : %5.1f : %5.1f %5.1f\n", + i, radAlp->data.F32[i], radMin->data.F32[i], binRad->data.F32[i], + radMax->data.F32[i], radBet->data.F32[i]); + } + + // storage vector for stats + psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32); + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + // psStats *stats = psStatsAlloc(PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4); + + // integrate flux, radius for each of these bins. since flux is sorted by radius, + // we can do this fairly quickly + + bool done = false; + int nOut = 0; + float Rmin = radAlp->data.F32[nOut]; + float Rmax = radBet->data.F32[nOut]; + float Rnxt = radAlp->data.F32[nOut+1]; // minimum radius for next range + int iNext = 0; + for (int i = 0; !done && (i < radius->n); i++) { + if (radius->data.F32[i] < Rnxt) { + iNext = i; + } + if (radius->data.F32[i] > Rmax) { + // calculate the value for the nOut bin + float value, dvalue; + if (values->n > 0) { + psVectorStats (stats, values, NULL, NULL, 0); + value = stats->robustMedian; + dvalue = stats->robustStdev; + } else { + value = NAN; + dvalue = NAN; + } + // binSB->data.F32[nOut] = stats->sampleMedian; + binSB->data.F32[nOut] = value; + binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ); + // binSB->data.F32[nOut] = stats->fittedMean; + // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ); + + // error in the SB is the stdev per bin / sqrt (number of pixels) + // added in quadrature to a fraction of the local sky (not the + // residual flux, but the sky from the sky model) + + psTrace ("psphot", 5, "%3d %5.1f %5.1f : %5.1f %5.2f\n", + nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]); + + nOut ++; + if (nOut >= radAlp->n) break; + Rmin = radAlp->data.F32[nOut]; + Rmax = radBet->data.F32[nOut]; + Rnxt = (nOut < nMax - 1) ? radAlp->data.F32[nOut+1] : Rmax; // minimum radius for next range + values->n = 0; + psStatsInit(stats); + i = iNext; + } + if (radius->data.F32[i] < Rmin) { + continue; + } + psVectorAppend (values, flux->data.F32[i]); + } + binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut; + // XXX I think this misses the last radial bin -- do we care? + + // save the vectors + profile->radialBins = binRad; + profile->area = binArea; + profile->binSB = binSB; + profile->binSBstdev = binSBstdev; + + // psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0); + + psFree(radMin); + psFree(radMax); + psFree(radAlp); + psFree(radBet); + psFree(values); + psFree(stats); + + return true; +} + +// the area-weighted mean radius is given by: + +// integral r * 2 pi r dr / integral 2 pi r dr + +// = 2/3 pi (r_max^3 - r_min^3) / pi (r_max^2 - r_min^2) +// = 2/3 (r_max^3 - r_min^3) / (r_max^2 - r_min^2) + Index: trunk/psphot/src/psphotPetrosianStats.c =================================================================== --- trunk/psphot/src/psphotPetrosianStats.c (revision 25755) +++ trunk/psphot/src/psphotPetrosianStats.c (revision 25755) @@ -0,0 +1,170 @@ +# include "psphotInternal.h" + +# define PETROSIAN_RATIO 0.2 +# define PETROSIAN_RADII 2.0 + +// generate the Petrosian radius and flux from the mean surface brightness (r_i) + +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); + +bool psphotPetrosianStats (pmSource *source) { + + pmSourceRadialProfile *profile = source->extpars->profile; + + float petRadius = NAN; + float petFlux = NAN; + + psVector *binSB = profile->binSB; + psVector *binSBstdev = profile->binSBstdev; + psVector *binRad = profile->radialBins; + psVector *area = profile->area; + + psVector *fluxSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *fluxSumErr2 = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *refRadius = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *petRatio = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *petRatioErr = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *meanSB = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *areaSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + + bool anyPetro = false; + bool manyPetro = false; + bool above = true; + float Asum = 0.0; + float Fsum = 0.0; + float dFsum2 = 0.0; + + float nSigma = 3.0; + int lowestSignificantRadius = 0; + float lowestSignificantRatio = 1.0; + + int nOut = 0; + for (int i = 0; i < binSB->n; i++) { + // skip nan bins (do not contribute to flux or area) + if (!isfinite(binSB->data.F32[i])) continue; + + float Area = area->data.F32[i]; + Asum += Area; + Fsum += binSB->data.F32[i] * Area; + dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area); + + float areaInner = 0.5 * Area; + float fluxInner = 0.5 * Area * binSB->data.F32[i]; + float fluxInnerErr2 = PS_SQR(binSBstdev->data.F32[i] * 0.5 * Area); + if (nOut > 0) { + areaInner += areaSum->data.F32[nOut-1]; + fluxInner += fluxSum->data.F32[nOut-1]; + fluxInnerErr2 += fluxSumErr2->data.F32[nOut-1]; + } + + // ratio = binSB / meanSB + // meanSB = flux / area + // flux = sum(binSB(i) * area(i) + // fluxErr^2 = sum(binSBerr(i)^2 area(i)^2) + // meanSBerr^2 = fluxErr^2 / area^2 + // (ratioErr/ratio)^2 = (binSBerr/binSB)^2 + (meanSBerr/meanSB)^2 + + psVectorAppend(meanSB, (fluxInner / areaInner)); + + float ratio = binSB->data.F32[i] / meanSB->data.F32[nOut]; + psVectorAppend(petRatio, ratio); + + float meanSBerr = sqrt(fluxInnerErr2) / areaInner; + float ratioErr = fabs(ratio) * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut])); + + psVectorAppend(petRatioErr, ratioErr); + + psVectorAppend(areaSum, Asum); + psVectorAppend(fluxSum, Fsum); + psVectorAppend(fluxSumErr2, dFsum2); + psVectorAppend(refRadius, binRad->data.F32[i]); + + psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", + i, refRadius->data.F32[nOut], + binSB->data.F32[i], binSBstdev->data.F32[i], + meanSB->data.F32[nOut], meanSBerr, + petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], + fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner); + + // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux + // we will keep and report the last (largest radius) value + if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) { + // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) + if (i == 0) { + // assume Fmax @ R = 0.0 + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + } else { + petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + } + above = false; + if (anyPetro) manyPetro = true; + anyPetro = true; + } + + // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux + // we will keep and report the last (largest radius) value + // find the last signficant measurement of the petrosian ratio + if (above && (petRatio->data.F32[nOut] < lowestSignificantRatio) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) { + lowestSignificantRadius = nOut; + lowestSignificantRatio = petRatio->data.F32[nOut]; + } + + // reset on transitions up, but do not re-calculate rad_90, flux_90 + if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) { + above = true; + } + nOut ++; + } + + if (!anyPetro) { + // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) + if (lowestSignificantRadius == 0) { + // assume Fmax @ R = 0.0 + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + } else { + petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + } + } + + // now measure the flux within PETROSIAN_RADII * petRadius + float apRadius = PETROSIAN_RADII * petRadius; + for (int i = 0; i < refRadius->n; i++) { + // XXX use bisection to do this faster: + if (refRadius->data.F32[i] > apRadius) { + if (i == 0) { + psWarning ("does this case make any sense? (refRadius[0] > apRadius)"); + continue; + } else { + petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); + break; + } + } + } + + if (!source->extpars->petrosian_80) { + source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc (); + } + pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + + // XXX save flags (anyPetro, manyPetro) + petrosian->radius = petRadius; + petrosian->flux = petFlux; + + // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); + + psFree(fluxSum); + psFree(fluxSumErr2); + psFree(refRadius); + psFree(petRatio); + psFree(petRatioErr); + psFree(meanSB); + psFree(areaSum); + + return true; +} + +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) { + float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0); + return Y; +} + Index: trunk/psphot/src/psphotPetrosianStudy.c =================================================================== --- trunk/psphot/src/psphotPetrosianStudy.c (revision 25755) +++ trunk/psphot/src/psphotPetrosianStudy.c (revision 25755) @@ -0,0 +1,184 @@ +# include "psphotInternal.h" + +# define DX 512 +# define DY 512 + +// XXX add noise and seeing. +// XXX double check on sersic functional form +// XXX modify ratio if ratio > 1.0 (swap major and minor) + +pmPeak *psphotLocalPeak(pmReadout *readout, int Xo, int Yo); + +int main (int argc, char **argv) { + + pmErrorRegister(); // register psModule's error codes/messages + pmModelClassInit(); + + int N; + float Xo = 0.5*DX; + float Yo = 0.5*DY; + char *image = NULL; + pmSource *source = NULL; + + float peak = 1000.0; + float sigma = 2.0; // major axis size + float ARatio = 1.0; + float angle = 0.0; + float sersic = 0.5; + float skynoise = 0.0; + + if ((N = psArgumentGet (argc, argv, "-peak"))) { + psArgumentRemove (N, &argc, argv); + peak = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-sigma"))) { + psArgumentRemove (N, &argc, argv); + sigma = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-aratio"))) { + psArgumentRemove (N, &argc, argv); + ARatio = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-angle"))) { + psArgumentRemove (N, &argc, argv); + angle = PS_RAD_DEG*atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-sersic"))) { + psArgumentRemove (N, &argc, argv); + sersic = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-skynoise"))) { + psArgumentRemove (N, &argc, argv); + skynoise = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-visual"))) { + psArgumentRemove (N, &argc, argv); + pmVisualSetVisual(true); + } + if ((N = psArgumentGet (argc, argv, "-coords"))) { + psArgumentRemove (N, &argc, argv); + Xo = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + Yo = atof(argv[N]); + psArgumentRemove (N, &argc, argv); + } + if ((N = psArgumentGet (argc, argv, "-image"))) { + psArgumentRemove (N, &argc, argv); + image = psStringCopy (argv[N]); + psArgumentRemove (N, &argc, argv); + } + + if (argc != 1) { + fprintf (stderr, "USAGE: psphotPetrosianStudy\n"); + exit (2); + } + + // create a containing image & associated readout + pmReadout *readout = pmReadoutAlloc(NULL); + if (!image) { + readout->image = psImageAlloc(DX, DY, PS_TYPE_F32); + + // create a dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0 + readout->variance = psImageAlloc(DX, DY, PS_TYPE_F32); + + // create a model & associated source + pmModelType type = pmModelClassGetType("PS_MODEL_SERSIC"); + pmModel *model = pmModelAlloc(type); + + // set the model parameters + model->params->data.F32[PM_PAR_SKY] = 0.0; + model->params->data.F32[PM_PAR_I0] = peak; + model->params->data.F32[PM_PAR_XPOS] = Xo; + model->params->data.F32[PM_PAR_YPOS] = Yo; + + psEllipseAxes axes; + axes.major = sigma; + axes.minor = sigma*ARatio; + axes.theta = angle; + + psEllipseShape shape = psEllipseAxesToShape (axes); + + // XXX set the sigma with user input + model->params->data.F32[PM_PAR_SXX] = shape.sx * M_SQRT2; + model->params->data.F32[PM_PAR_SYY] = shape.sy * M_SQRT2; + model->params->data.F32[PM_PAR_SXY] = shape.sxy; + + if (model->params->n > 7) { + model->params->data.F32[PM_PAR_7] = sersic; + } + + // generate source container & populate image + source = pmSourceFromModel(model, readout, Xo, PM_SOURCE_TYPE_STAR); + + // generate the modelFlux + pmSourceCacheModel(source, 0); + + // instantiate the source + pmSourceAdd(source, PM_MODEL_OP_FUNC, 0); + + // XXX add noise here... + psphotSaveImage(NULL, readout->image, "sersic.fits"); + + } else { + psRegion full = psRegionSet(0,0,0,0); + psFits *fits = psFitsOpen(image, "r"); + readout->image = psFitsReadImage(fits, full, 0); + + source = pmSourceAlloc(); + source->peak = psphotLocalPeak(readout, Xo, Yo); + pmSourceDefinePixels (source, readout, Xo, Yo, 128); + } + + psphotPetrosianProfile (readout, source, skynoise); + + psFree (source); + + exit (0); +} + +// Xo, Yo are in parent coords +pmPeak *psphotLocalPeak(pmReadout *readout, int Xo, int Yo) { + + int Xp = Xo; + int Yp = Yo; + float peakFlux = readout->image->data.F32[Yp][Xp]; + + // find local peak within +/- 3 pix of the given coordinate + for (int iy = Yo - 3; iy <= Yo + 3; iy++) { + for (int ix = Xo - 3; ix <= Xo + 3; ix++) { + if (peakFlux < readout->image->data.F32[iy][ix]) { + Xp = ix; + Yp = iy; + peakFlux = readout->image->data.F32[Yp][Xp]; + } + } + } + + pmPeak *peak = pmPeakAlloc(Xp, Yp, peakFlux, PM_PEAK_LONE); + + // calculate fractional peak position relative to Xp,Yp + psPolynomial2D *bicube = psImageBicubeFit (readout->image, Xp, Yp); + psPlane min = psImageBicubeMin (bicube); + psFree (bicube); + + // if min point is too deviant, use the peak value + if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { + peak->xf = min.x + Xp; + peak->yf = min.y + Yp; + + // xf,yf must land on image with 0 pixel border + peak->xf = PS_MAX (PS_MIN (peak->xf, readout->image->numCols - 1), readout->image->col0); + peak->yf = PS_MAX (PS_MIN (peak->yf, readout->image->numRows - 1), readout->image->row0); + } else { + peak->xf = Xp; + peak->yf = Yp; + } + + return peak; +} Index: trunk/psphot/src/psphotPetrosianVisual.c =================================================================== --- trunk/psphot/src/psphotPetrosianVisual.c (revision 25755) +++ trunk/psphot/src/psphotPetrosianVisual.c (revision 25755) @@ -0,0 +1,405 @@ +# include "psphotInternal.h" + +// this function displays representative images as the psphot analysis progresses: +// 0 : image, 1 : variance +// 0 : backsub, 1 : variance, 2 : backgnd +// 0 : backsub, 1 : variance, 2 : signif +// (overlay peaks on images) +// (overlay footprints on images) +// (overlay moments on images) +// (overlay rough class on images) +// 0 : backsub, 1 : psfpos, 2: psfsub +// 0 : backsub, 1 : lin_resid, 2: psfsub + +# if (HAVE_KAPA) +# include + +// functions used to visualize the analysis as it goes +// these are invoked by the -visual options + +static int kapa = -1; +static int kapa2 = -1; + +// if no valid data is supplied (NULL or n <- 0), leave limits as they were +bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec) { + + if (xVec && xVec->n > 0) { + graphdata->xmin = graphdata->xmax = xVec->data.F32[0]; + for (int i = 1; i < xVec->n; i++) { + if (!isfinite(xVec->data.F32[i])) continue; + graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]); + graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]); + } + float range = graphdata->xmax - graphdata->xmin; + graphdata->xmax += 0.05*range; + graphdata->xmin -= 0.05*range; + } + if (yVec && yVec->n > 0) { + graphdata->ymin = graphdata->ymax = yVec->data.F32[0]; + for (int i = 1; i < yVec->n; i++) { + if (!isfinite(yVec->data.F32[i])) continue; + graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]); + graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]); + } + float range = graphdata->ymax - graphdata->ymin; + graphdata->ymax += 0.05*range; + graphdata->ymin -= 0.05*range; + } + return true; +} + +bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux) { + + Graphdata graphdata; + + // return true; + if (!pmVisualIsVisual()) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearPlots (kapa2); + KapaInitGraph (&graphdata); + KapaSetFont (kapa2, "courier", 14); + + pmVisualLimitsFromVectors (&graphdata, radius, flux); + KapaSetLimits (kapa2, &graphdata); + + KapaBox (kapa2, &graphdata); + KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "flux", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa2, radius->n, &graphdata); + KapaPlotVector (kapa2, radius->n, radius->data.F32, "x"); + KapaPlotVector (kapa2, radius->n, flux->data.F32, "y"); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef) { + + float FluxRef = 500.0; + + Graphdata graphdata; + + // return true; + if (!pmVisualIsVisual()) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearPlots (kapa); + KapaInitGraph (&graphdata); + KapaSetFont (kapa, "courier", 14); + + graphdata.ymax = +1.05*peakFlux; + graphdata.ymin = -0.05*peakFlux; + pmVisualLimitsFromVectors (&graphdata, radius, NULL); + KapaSetLimits (kapa, &graphdata); + + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa, "flux", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa, radius->n, &graphdata); + KapaPlotVector (kapa, radius->n, radius->data.F32, "x"); + KapaPlotVector (kapa, radius->n, flux->data.F32, "y"); + + // do this with log-r, log-flux? + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.size = 2.0; + KapaPrepPlot (kapa, radiusBin->n, &graphdata); + KapaPlotVector (kapa, radiusBin->n, radiusBin->data.F32, "x"); + KapaPlotVector (kapa, radiusBin->n, fluxBin->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 2; + graphdata.ptype = 7; + graphdata.size = 3.0; + KapaPrepPlot (kapa, 1, &graphdata); + KapaPlotVector (kapa, 1, &RadiusRef, "x"); + KapaPlotVector (kapa, 1, &FluxRef, "y"); + + fprintf (stderr, "radius: %f\n", RadiusRef); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, + psVector *refRadius, psVector *meanSB, + psVector *petRatio, psVector *petRatioErr, + psVector *fluxSum, + float petRadius, float ratioForRadius, + float petFlux, float radiusForFlux) +{ + Graphdata graphdata; + KapaSection section; + + if (!pmVisualIsVisual()) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearPlots (kapa); + KapaInitGraph (&graphdata); + KapaSetFont (kapa, "courier", 14); + + // radius vs flux + // radius vs mean SB + // radius vs petRatio + + // *** section 1: radius vs mean SB + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("meanSB"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualLimitsFromVectors (&graphdata, radBin, fluxBin); + KapaSetLimits (kapa, &graphdata); + + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa, "mean SB", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa, radBin->n, &graphdata); + KapaPlotVector (kapa, radBin->n, radBin->data.F32, "x"); + KapaPlotVector (kapa, radBin->n, fluxBin->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 2; + graphdata.ptype = 1; + graphdata.size = 2.0; + KapaPrepPlot (kapa, refRadius->n, &graphdata); + KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); + KapaPlotVector (kapa, refRadius->n, meanSB->data.F32, "y"); + + // *** section 2: radius vs petrosian ratio + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.33; + section.name = psStringCopy ("ratio"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.ymax = +1.05; + graphdata.ymin = -0.05; + pmVisualLimitsFromVectors (&graphdata, radBin, NULL); + KapaSetLimits (kapa, &graphdata); + + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa, "ratio", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + graphdata.etype = 0x01; + KapaPrepPlot (kapa, refRadius->n, &graphdata); + KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); + KapaPlotVector (kapa, refRadius->n, petRatio->data.F32, "y"); + KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dym"); + KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dyp"); + graphdata.etype = 0; + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.size = 2.0; + KapaPrepPlot (kapa, 1, &graphdata); + KapaPlotVector (kapa, 1, &petRadius, "x"); + KapaPlotVector (kapa, 1, &ratioForRadius, "y"); + + // *** section 3: radius vs integrated flux + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.66; + section.name = psStringCopy ("flux"); + KapaSetSection (kapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualLimitsFromVectors (&graphdata, radBin, fluxSum); + KapaSetLimits (kapa, &graphdata); + + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa, "integrated flux", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa, refRadius->n, &graphdata); + KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); + KapaPlotVector (kapa, refRadius->n, fluxSum->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.style = 2; + graphdata.size = 2.0; + KapaPrepPlot (kapa, 1, &graphdata); + KapaPlotVector (kapa, 1, &radiusForFlux, "x"); + KapaPlotVector (kapa, 1, &petFlux, "y"); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian) { + + Graphdata graphdata; + + if (!pmVisualIsVisual()) return true; + + if (kapa == -1) { + kapa = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearPlots (kapa); + KapaInitGraph (&graphdata); + KapaSetFont (kapa, "courier", 14); + + psVector *theta = petrosian->theta; + psVector *radius = petrosian->isophotalRadii; + + // find Rmin and Rmax for the initial guess + float Rmin = radius->data.F32[0]; + float Rmax = radius->data.F32[0]; + + psVector *Rx = psVectorAlloc(radius->n, PS_TYPE_F32); + psVector *Ry = psVectorAlloc(radius->n, PS_TYPE_F32); + + for (int i = 0; i < theta->n; i++) { + Rx->data.F32[i] = radius->data.F32[i]*cos(theta->data.F32[i]); + Ry->data.F32[i] = radius->data.F32[i]*sin(theta->data.F32[i]); + + // check the radius range + Rmin = MIN (Rmin, radius->data.F32[i]); + Rmax = MAX (Rmax, radius->data.F32[i]); + } + + psVector *rx = psVectorAlloc(361, PS_TYPE_F32); + psVector *ry = psVectorAlloc(361, PS_TYPE_F32); + + float epsilon = petrosian->axes.minor / petrosian->axes.major; + + for (int i = 0; i < 361; i++) { + + float alpha = PS_RAD_DEG * i; + + float cs_alpha = cos(alpha); + float sn_alpha = sin(alpha); + + float cs_phi = cos(alpha - petrosian->axes.theta); + float sn_phi = sin(alpha - petrosian->axes.theta); + + float r = 1.0 / sqrt(SQ(sn_phi) + SQ(epsilon*cs_phi)); + + // generate the model fit here + rx->data.F32[i] = petrosian->axes.minor * cs_alpha * r; + ry->data.F32[i] = petrosian->axes.minor * sn_alpha * r; + } + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -1.1*Rmax; + graphdata.ymin = -1.1*Rmax; + graphdata.xmax = +1.1*Rmax; + graphdata.ymax = +1.1*Rmax; + KapaSetLimits (kapa, &graphdata); + + KapaBox (kapa, &graphdata); + KapaSendLabel (kapa, "R_x", KAPA_LABEL_XM); + KapaSendLabel (kapa, "R_y", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.size = 1.0; + KapaPrepPlot (kapa, Rx->n, &graphdata); + KapaPlotVector (kapa, Rx->n, Rx->data.F32, "x"); + KapaPlotVector (kapa, Rx->n, Ry->data.F32, "y"); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 0; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa, rx->n, &graphdata); + KapaPlotVector (kapa, rx->n, rx->data.F32, "x"); + KapaPlotVector (kapa, rx->n, ry->data.F32, "y"); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +# endif Index: trunk/psphot/src/psphotRadialProfile.c =================================================================== --- trunk/psphot/src/psphotRadialProfile.c (revision 25738) +++ trunk/psphot/src/psphotRadialProfile.c (revision 25755) @@ -1,21 +1,5 @@ # include "psphotInternal.h" -# define COMPARE_RADIUS(A,B) (radius->data.F32[A] < radius->data.F32[B]) -# define SWAP_RADIUS(TYPE,A,B) { \ - float tmp; \ - if (A != B) { \ - tmp = radius->data.F32[A]; \ - radius->data.F32[A] = radius->data.F32[B]; \ - radius->data.F32[B] = tmp; \ - tmp = flux->data.F32[A]; \ - flux->data.F32[A] = flux->data.F32[B]; \ - flux->data.F32[B] = tmp; \ - tmp = variance->data.F32[A]; \ - variance->data.F32[A] = variance->data.F32[B]; \ - variance->data.F32[B] = tmp; \ - } \ -} - -bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) { +bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) { // allocate pmSourceExtendedParameters, if not already defined @@ -28,44 +12,39 @@ } - int nPts = source->pixels->numRows * source->pixels->numCols; - source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32); - source->extpars->profile->flux = psVectorAllocEmpty (nPts, PS_TYPE_F32); - source->extpars->profile->variance = psVectorAllocEmpty (nPts, PS_TYPE_F32); + // XXX these need to go into recipe values + int Nsec = 24; + float Rmax = 200; + float fluxMin = 0.0; + float fluxMax = source->peak->flux; - psVector *radius = source->extpars->profile->radius; - psVector *flux = source->extpars->profile->flux; - psVector *variance = source->extpars->profile->variance; + // generate a series of radial profiles at Nsec evenly spaced angles. the profile flux + // is measured by interpolation for small radii; for large radii, the pixels in a box + // are averaged to increase the S/N + if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); + return false; + } - // XXX use the extended source model here for Xo, Yo? - // XXX define a radius scaled to the elliptical contour? + // use the radial profiles to determine the radius of a given isophote. this isophote + // is used to determine the elliptical shape of the object, so it has a relatively high + // value (nominally 50% of the peak) + if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); + return false; + } - int n = 0; - - float Xo = 0.0; - float Yo = 0.0; - - if (source->modelEXT) { - Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; - Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; - } else { - Xo = source->peak->xf - source->pixels->col0; - Yo = source->peak->yf - source->pixels->row0; + // convert the isophotal radius vs angle measurements to an elliptical contour + if (!psphotEllipticalContour (source)) { + psLogMsg ("psphot", 3, "failed to measure elliptical contour"); + return false; } - 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]) continue; - radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ; - flux->data.F32[n] = source->pixels->data.F32[iy][ix]; - variance->data.F32[n] = source->variance->data.F32[iy][ix]; - n++; - } + + // generate a single, normalized radial profile following the elliptical contours. + // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel + if (!psphotEllipticalProfile (source)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; } - radius->n = n; - variance->n = n; - flux->n = n; - - // sort the vector set by the radius - PSSORT (radius->n, COMPARE_RADIUS, SWAP_RADIUS, NONE); - + return true; } Index: trunk/psphot/src/psphotRadialProfileByAngles.c =================================================================== --- trunk/psphot/src/psphotRadialProfileByAngles.c (revision 25755) +++ trunk/psphot/src/psphotRadialProfileByAngles.c (revision 25755) @@ -0,0 +1,239 @@ +# include "psphotInternal.h" + +// Given a source at (x,y), generate a collection of radial profiles at even angular separations + +// These functions are used to calculate the stats in a rectangle at arbitrary orientation. +// XXX Move these elsewhere (psLib?) +float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta); +psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta); +psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW); +psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); + +bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax) { + + // we want to have an even number of sectors so we can do 180 deg symmetrizing + Nsec = (Nsec % 2) ? Nsec + 1 : Nsec; + float dtheta = 2.0*M_PI / Nsec; + + pmSourceRadialProfile *profile = source->extpars->profile; + psFree(profile->radii); + psFree(profile->fluxes); + psFree(profile->theta); + + profile->radii = psArrayAllocEmpty(Nsec); + profile->fluxes = psArrayAllocEmpty(Nsec); + profile->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32); + + + for (int i = 0; i < Nsec; i++) { + + float theta = i*dtheta; + + psVector *radius = psVectorAllocEmpty(Rmax, PS_TYPE_F32); + psVector *flux = psVectorAllocEmpty(Rmax, PS_TYPE_F32); + + // Start at Xo,Yo and find the x,y locations for r_i, theta where r_i initially + // increments by 1 pixel. At large radii (r*dtheta > 2) use stats in a box rather than + // sub-pixel interpolation + + int dR = 1.0; + for (float r = 0; r < Rmax; r += dR) { + + float Xo = source->peak->xf; + float Yo = source->peak->yf; + + // Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0 + float x = r * cos (theta) + Xo; + float y = r * sin (theta) + Yo; + dR = 2*(int)(0.5*r*sin(dtheta)) + 1; + + if (x < 0) goto badvalue; + if (y < 0) goto badvalue; + if (x >= source->pixels->parent->numCols) goto badvalue; + if (y >= source->pixels->parent->numRows) goto badvalue; + + float value = NAN; + if (dR < 2) { + // value is NAN if we run off the image + // 0.5 PIX: this function takes pixel coords; source peak is in pixel coords + value = psImageInterpolatePixelBilinear(x, y, source->pixels); + } else { + // 0.5 PIX: this function takes pixel coords; source peak is in pixel coords + value = psphotMeanSectorValue(source->pixels, x, y, dR, dR, theta); + } + + // keep the all values (even NAN) so all vectors are matched in length + psVectorAppend (radius, r); + psVectorAppend (flux, value); + continue; + + badvalue: + psVectorAppend (radius, r); + psVectorAppend (flux, NAN); + } + + psArrayAdd (profile->radii, 100, radius); + psArrayAdd (profile->fluxes, 100, flux); + psVectorAppend (profile->theta, theta); + + // psphotPetrosianVisualProfileByAngle (radius, flux); + + psFree(radius); + psFree(flux); + } + + for (int i = 0; i < Nsec / 2; i++) { + + psVector *r1 = profile->radii->data[i]; + psVector *r2 = profile->radii->data[i+Nsec/2]; + + psVector *f1 = profile->fluxes->data[i]; + psVector *f2 = profile->fluxes->data[i+Nsec/2]; + + psAssert (r1->n == r2->n, "mis-matched vectors"); + psAssert (f1->n == f2->n, "mis-matched vectors"); + + // we have a pair of vectors i, i+Nsec/2; replace them with the finite minimum of the pair + for (int j = 0; j < r1->n; j++) { + + float flux; + + if (!isfinite(f1->data.F32[j]) && !isfinite(f2->data.F32[j])) { + flux = NAN; + goto setflux; + } + + if (!isfinite(f1->data.F32[j])) { + flux = f2->data.F32[j]; + goto setflux; + } + if (!isfinite(f2->data.F32[j])) { + flux = f1->data.F32[j]; + goto setflux; + } + + flux = PS_MIN(f1->data.F32[j], f2->data.F32[j]); + + setflux: + f1->data.F32[j] = flux; + f2->data.F32[j] = flux; + } + } + return true; +} + +float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta) { + + psVector *values = psphotBoxValues (image, x, y, dL, dW, theta); + if (!values) goto escape; + if (!values->n) goto escape; + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, values, NULL, NULL, 0); + + float value = stats->sampleMedian; + + psFree (stats); + psFree (values); + + return value; + +escape: + psFree(values); + return NAN; +} + +psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta) { + + // extract pixels from a series of lines (from -0.5*dW to +0.5*dW) of length dL, + // centered on x0, y0 in parent pixel coordinates (not pixel indicies) + + float xs = x0 - image->col0 - 0.5*dL*cos(theta); + float ys = y0 - image->row0 - 0.5*dL*sin(theta); + + float xe = xs + 0.5*dL*cos(theta); + float ye = ys + 0.5*dL*sin(theta); + + psVector *values = psphotLineValues (image, xs, ys, xe, ye, (int) dW); + return values; +} + +/** + * identify the quadrant and draw the correct line + */ +psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW) { + + int FlipDirect, FlipCoords; + int X1, Y1, X2, Y2, dX, dY; + + /* rather than draw the line from float positions, we find the closest + integer end-points and draw the line between those pixels */ + + X1 = ROUND(x1); + Y1 = ROUND(y1); + X2 = ROUND(x2); + Y2 = ROUND(y2); + + dX = X2 - X1; + dY = Y2 - Y1; + + FlipCoords = (abs(dX) < abs(dY)); + FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2); + + psVector *values = NULL; + if (!FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X1, Y1, X2, Y2, dW, FALSE); + if ( FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X2, Y2, X1, Y1, dW, FALSE); + if (!FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y1, X1, Y2, X2, dW, TRUE); + if ( FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y2, X2, Y1, X1, dW, TRUE); + + return values; +} + +/** + * use the Bresenham line drawing technique + * integer-only Bresenham line-draw version which is fast + */ +psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) { + + int X, Y, dX, dY; + int e, e2; + + psVector *values = psVectorAllocEmpty(100, PS_TYPE_F32); + + dX = X2 - X1; + dY = Y2 - Y1; + + Y = Y1; + e = 0; + for (X = X1; X <= X2; X++) { + if (X > 0) { + if (swapcoords) { + if (X >= image->numRows) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= image->numCols) continue; + psVectorAppend(values, image->data.F32[X][y]); + } + } else { + if (X >= image->numCols) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= image->numRows) continue; + psVectorAppend(values, image->data.F32[y][X]); + } + } + } + e += dY; + e2 = 2 * e; + if (e2 > dX) { + Y++; + e -= dX; + } + if (e2 < -dX) { + Y--; + e += dX; + } + } + return values; +} + Index: trunk/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- trunk/psphot/src/psphotRadiiFromProfiles.c (revision 25755) +++ trunk/psphot/src/psphotRadiiFromProfiles.c (revision 25755) @@ -0,0 +1,153 @@ +# include "psphotInternal.h" + +// Given the Radial Profiles (radii, fluxes) determine the radius for each profile at the desired isophote + +bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax) { + + pmSourceRadialProfile *profile = source->extpars->profile; + + psFree(profile->isophotalRadii); + profile->isophotalRadii = psVectorAlloc(profile->theta->n, PS_TYPE_F32); + + for (int i = 0; i < profile->theta->n; i++) { + psVector *radii = profile->radii->data[i]; + psVector *fluxes = profile->fluxes->data[i]; + float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax); + + // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius); + + // warn on NAN? + profile->isophotalRadii->data.F32[i] = radius; + } + return true; +} + +float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax) { + + // 'flux' is a noisy sample of the galaxy radial profile at points 'radius' + // rebin flux into samples defined by the isophote Fo = 0.5*(fluxMax + fluxMin). the noisy + // sample is cleaned by rebinning to a well-matched radial binning + + // base selections on fluxes defined by the flux range dF + float fluxRange = fluxMax - fluxMin; + + // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size + // XXX reconsider the fractional isophote value + float Fm = fluxMin + 0.10*fluxRange; + float Fo = fluxMin + 0.25*fluxRange; + float Fp = fluxMin + 0.50*fluxRange; + int Rbin = 1; + + // find the median radius of the points in the flux range Fm - Fp: + { + // storage vector for stats + psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32); + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + + for (int i = 0; i < flux->n; i++) { + if (!isfinite(flux->data.F32[i])) continue; + if (flux->data.F32[i] < Fm) continue; + if (flux->data.F32[i] > Fp) continue; + + psVectorAppend (values, radius->data.F32[i]); + } + if (values->n > 1) { + psVectorStats (stats, values, NULL, NULL, 0); + + // if we have a valid range, rebin with bin size 1/2 of median radius + if (isfinite(stats->sampleMedian)) { + Rbin = MAX(1, 0.5*stats->sampleMedian); + } + } + psFree (values); + psFree (stats); + } + Rbin = 3; + + psVector *fluxBinned = NULL; + psVector *radiusBinned = NULL; + + // do not bother rebinning if the bin size is only 2 or less + if (Rbin <= 2) { + fluxBinned = psMemIncrRefCounter (flux); + radiusBinned = psMemIncrRefCounter (radius); + } else { + // storage vector for stats + psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32); + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + + // rebinned vectors + fluxBinned = psVectorAllocEmpty (flux->n, PS_TYPE_F32); + radiusBinned = psVectorAllocEmpty (flux->n, PS_TYPE_F32); + + // sort the flux by the radius + pmSourceRadialProfileSortPair (radius, flux); + + int nOut = 0; + radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin; + float Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin; + float Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin; + + for (int i = 0; i < flux->n; i++) { + if (radius->data.F32[i] < Rmin) { + // XXX not sure how we can hit this, if there is full coverage of radiusBinned + continue; + } + if (radius->data.F32[i] > Rmax) { + // calculate the value for the nOut bin + // XXX need to fix this as well psStats (stats, values); + float value; + if (values->n > 0) { + psVectorStats (stats, values, NULL, NULL, 0); + value = stats->sampleMedian; + } else { + value = NAN; + } + fluxBinned->data.F32[nOut] = value; + nOut ++; + radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin; + Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin; + Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin; + values->n = 0; + psStatsInit(stats); + } + if (!isfinite(flux->data.F32[i])) continue; + psVectorAppend (values, flux->data.F32[i]); + } + fluxBinned->n = nOut; + radiusBinned->n = nOut; + psFree (values); + psFree(stats); + } + + float Ro = NAN; + bool above = true; + for (int i = 0; i < fluxBinned->n; i++) { + + if (!isfinite(fluxBinned->data.F32[i])) continue; + + // find the largest radius that matches the flux transition + if (above && (fluxBinned->data.F32[i] < Fo)) { + // XXX is there a macro in psLib that does this interpolation? + if (i == 0) { + psTrace ("psphot", 4, "bogus radial profile for source at %f, %f, skipping", source->peak->xf, source->peak->yf); + psFree (fluxBinned); + psFree (radiusBinned); + return NAN; + } + Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]); + above = FALSE; + } + + if (!above && (fluxBinned->data.F32[i] >= Fo)) { + above = TRUE; + } + } + + // show the results + // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, Ro); + + psFree(fluxBinned); + psFree(radiusBinned); + return Ro; +} Index: trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- trunk/psphot/src/psphotRadiusChecks.c (revision 25738) +++ trunk/psphot/src/psphotRadiusChecks.c (revision 25755) @@ -7,4 +7,7 @@ // and a per-object radius is calculated) +static float PSF_APERTURE = 0; // radius to use in PSF aperture mags + + bool psphotInitRadiusPSF(const psMetadata *recipe, const pmModelType type) { @@ -13,5 +16,6 @@ PSF_FIT_NSIGMA = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA"); PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING"); - PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); + PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); + PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); return true; @@ -34,17 +38,21 @@ radiusFit = model->modelRadius(model->params, 1.0); } + model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); + } else { + model->fitRadius = radiusFit; } - model->radiusFit = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); - - if (isnan(model->radiusFit)) psAbort("error in radius"); + if (isnan(model->fitRadius)) psAbort("error in radius"); if (source->mode & PM_SOURCE_MODE_SATSTAR) { - model->radiusFit *= 2; + model->fitRadius *= 2; } - bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit); + // radius used to measure aperture photometry + source->apRadius = PSF_APERTURE; + + bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius); // set the mask to flag the excluded pixels - psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", markVal); + psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal); return status; } @@ -58,15 +66,31 @@ // set the fit radius based on the object flux limit and the model - model->radiusFit = (RADIUS_TYPE) (model->modelRadius (model->params, PSF_FIT_NSIGMA*moments->dSky) + dR + PSF_FIT_PADDING); - if (isnan(model->radiusFit)) psAbort("error in radius"); - + float radiusFit = PSF_FIT_RADIUS; + if (radiusFit <= 0) { // use fixed radius + if (moments == NULL) { + radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + } else { + radiusFit = model->modelRadius(model->params, 1.0); + } + model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); + } else { + model->fitRadius = radiusFit; + } + if (isnan(model->fitRadius)) psAbort("error in radius"); + + // above sets a radius for a single star, bump by blend separation + model->fitRadius += dR; + if (source->mode & PM_SOURCE_MODE_SATSTAR) { - model->radiusFit *= 2; + model->fitRadius *= 2; } - bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit); + // radius used to measure aperture photometry + source->apRadius = PSF_APERTURE; + + bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius); // set the mask to flag the excluded pixels - psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", markVal); + psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal); return status; } @@ -74,4 +98,5 @@ static float EXT_FIT_NSIGMA; static float EXT_FIT_PADDING; +static float EXT_FIT_MAX_RADIUS; bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) { @@ -79,6 +104,7 @@ bool status; - EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA"); - EXT_FIT_PADDING = psMetadataLookupF32 (&status, recipe, "EXT_FIT_PADDING"); + EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA"); + EXT_FIT_PADDING = psMetadataLookupF32 (&status, recipe, "EXT_FIT_PADDING"); + EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS"); return true; @@ -86,5 +112,52 @@ // call this function whenever you (re)-define the EXT model +float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) { + + psAssert (source, "source not defined??"); + psAssert (source->peak, "peak not defined??"); + + pmPeak *peak = source->peak; + + // set the radius based on the footprint: + if (!peak->footprint) goto escape; + pmFootprint *footprint = peak->footprint; + if (!footprint->spans) goto escape; + if (footprint->spans->n < 1) goto escape; + + // find the max radius + float radius = 0.0; + for (int j = 0; j < footprint->spans->n; j++) { + pmSpan *span = footprint->spans->data[j]; + + float dY = span->y - peak->yf; + float dX0 = span->x0 - peak->xf; + float dX1 = span->x1 - peak->xf; + + radius = PS_MAX (radius, hypot(dY, dX0)); + radius = PS_MAX (radius, hypot(dY, dX1)); + } + + radius += EXT_FIT_PADDING; + if (isnan(radius)) psAbort("error in radius"); + + radius = PS_MIN (radius, EXT_FIT_MAX_RADIUS); + + // redefine the pixels if needed + pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal); + return radius; + +escape: + return NAN; + // bool result = psphotCheckRadiusEXT (readout, source, model, markVal); + // return result; +} + +// alternative EXT radius based on model guess (for use without footprints) bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) { + + psAbort ("do not use this function"); psF32 *PAR = model->params->data.F32; @@ -96,12 +169,12 @@ float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky); - model->radiusFit = rawRadius + EXT_FIT_PADDING; - if (isnan(model->radiusFit)) psAbort("error in radius"); + model->fitRadius = rawRadius + EXT_FIT_PADDING; + if (isnan(model->fitRadius)) psAbort("error in radius"); // redefine the pixels if needed - bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit); + bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius); // set the mask to flag the excluded pixels - psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", markVal); + psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius, "OR", markVal); return status; } Index: trunk/psphot/src/psphotReadout.c =================================================================== --- trunk/psphot/src/psphotReadout.c (revision 25738) +++ trunk/psphot/src/psphotReadout.c (revision 25755) @@ -81,5 +81,5 @@ // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections); + psArray *sources = psphotSourceStats (config, readout, detections, true); if (!sources) return false; if (!strcasecmp (breakPt, "PEAKS")) { @@ -126,5 +126,4 @@ return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); } - psphotVisualShowPSFModel (readout, psf); @@ -132,11 +131,8 @@ psphotLoadExtSources (config, view, sources); - // construct an initial model for each object + // construct an initial model for each object, set the radius to fitRadius, set circular fit mask psphotGuessModels (config, readout, sources, psf); - // XXX test output of models - // psphotTestSourceOutput (readout, sources, recipe, psf); - - // linear PSF fit to source peaks + // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); @@ -144,13 +140,14 @@ // psphotGuessModels or fitted until psphotFitSourcesLinear. psphotVisualShowPSFStars (recipe, psf, sources); - psphotVisualShowSatStars (recipe, psf, sources); // identify CRs and extended sources - psphotSourceSize (config, readout, sources, recipe, 0); + psphotSourceSize (config, readout, sources, recipe, psf, 0); if (!strcasecmp (breakPt, "ENSEMBLE")) { goto finish; } + psphotVisualShowSatStars (recipe, psf, sources); // non-linear PSF and EXT fit to brighter sources + // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) psphotBlendFit (config, readout, sources, psf); @@ -158,5 +155,5 @@ psphotReplaceAllSources (sources, recipe); - // linear fit to include all sources + // linear fit to include all sources (subtract again) psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); @@ -165,10 +162,5 @@ goto pass1finish; } - - // XXX for the moment, drop the re-calc of the background (prove this works) - // replace background in residual image - // psphotSkyReplace (config, view); - // re-measure background model (median, smoothed image) - // psphotImageMedian (config, view); + // NOTE: possibly re-measure background model here with objects subtracted // add noise for subtracted objects @@ -182,5 +174,5 @@ // define new sources based on only the new peaks - psArray *newSources = psphotSourceStats (config, readout, detections); + psArray *newSources = psphotSourceStats (config, readout, detections, false); // set source type @@ -190,5 +182,5 @@ } - // create full input models + // create full input models, set the radius to fitRadius, set circular fit mask psphotGuessModels (config, readout, newSources, psf); @@ -206,5 +198,5 @@ // measure source size for the remaining sources - psphotSourceSize (config, readout, sources, recipe, 0); + psphotSourceSize (config, readout, sources, recipe, psf, 0); psphotExtendedSourceAnalysis (readout, sources, recipe); Index: trunk/psphot/src/psphotReadoutFindPSF.c =================================================================== --- trunk/psphot/src/psphotReadoutFindPSF.c (revision 25738) +++ trunk/psphot/src/psphotReadoutFindPSF.c (revision 25755) @@ -41,5 +41,5 @@ // construct sources and measure basic stats (moments, local sky) - psArray *sources = psphotSourceStats(config, readout, detections); + psArray *sources = psphotSourceStats(config, readout, detections, true); if (!sources) return false; Index: trunk/psphot/src/psphotReadoutKnownSources.c =================================================================== --- trunk/psphot/src/psphotReadoutKnownSources.c (revision 25738) +++ trunk/psphot/src/psphotReadoutKnownSources.c (revision 25755) @@ -41,5 +41,5 @@ // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections); + psArray *sources = psphotSourceStats (config, readout, detections, true); if (!sources) return false; Index: trunk/psphot/src/psphotReadoutMinimal.c =================================================================== --- trunk/psphot/src/psphotReadoutMinimal.c (revision 25738) +++ trunk/psphot/src/psphotReadoutMinimal.c (revision 25755) @@ -56,5 +56,5 @@ // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections); + psArray *sources = psphotSourceStats (config, readout, detections, true); if (!sources) return false; Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 25738) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 25755) @@ -17,5 +17,4 @@ replace: pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot.replace"), sources->n); @@ -41,5 +40,4 @@ pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); @@ -64,6 +62,5 @@ if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); } psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); @@ -71,4 +68,5 @@ } +# if (0) // add source, if the source has been subtracted; do not modify state bool psphotAddWithTest (pmSource *source, bool useState, psImageMaskType maskVal) { @@ -108,2 +106,3 @@ return true; } +# endif Index: trunk/psphot/src/psphotRoughClass.c =================================================================== --- trunk/psphot/src/psphotRoughClass.c (revision 25738) +++ trunk/psphot/src/psphotRoughClass.c (revision 25755) @@ -26,11 +26,12 @@ for (int iy = 0; iy < NY; iy ++) { - psRegion region = psRegionSet (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY); - if (!psphotRoughClassRegion (nRegion, ®ion, sources, recipe, havePSF)) { + psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY); + if (!psphotRoughClassRegion (nRegion, region, sources, recipe, havePSF)) { psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", - region.x0, region.y0, region.x1, region.y1); + region->x0, region->y0, region->x1, region->y1); + psFree (region); continue; } - + psFree (region); nRegion ++; } @@ -45,5 +46,5 @@ psphotVisualPlotMoments (recipe, sources); psphotVisualShowRoughClass (sources); - psphotVisualShowFlags (sources); + // XXX better visualization: psphotVisualShowFlags (sources); return true; @@ -63,11 +64,14 @@ psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); if (!regionMD) { + // allocate the region metadata folder and add this region to it. regionMD = psMetadataAlloc(); psMetadataAddMetadata (recipe, PS_LIST_TAIL, regionName, PS_META_REPLACE, "psf clump region", regionMD); psFree (regionMD); } + psMetadataAddPtr (regionMD, PS_LIST_TAIL, "REGION", PS_DATA_REGION | PS_META_REPLACE, "psf clump region", region); if (!havePSF) { // determine the PSF parameters from the source moment values + // XXX why not save the psfClump as a PTR? psfClump = pmSourcePSFClump (region, sources, recipe); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); Index: trunk/psphot/src/psphotSetThreads.c =================================================================== --- trunk/psphot/src/psphotSetThreads.c (revision 25738) +++ trunk/psphot/src/psphotSetThreads.c (revision 25755) @@ -10,5 +10,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_MAGNITUDES", 8); + task = psThreadTaskAlloc("PSPHOT_MAGNITUDES", 9); task->function = &psphotMagnitudes_Threaded; psThreadTaskAdd(task); @@ -20,5 +20,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_APRESID_MAGS", 6); + task = psThreadTaskAlloc("PSPHOT_APRESID_MAGS", 7); task->function = &psphotApResidMags_Threaded; psThreadTaskAdd(task); Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 25738) +++ trunk/psphot/src/psphotSourceFits.c (revision 25755) @@ -90,6 +90,11 @@ psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, dR); - // fit PSF model (set/unset the pixel mask) + // fit PSF model pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); + + // clear the circular mask + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + + if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); // correct model chisq for flux trend @@ -101,4 +106,6 @@ pmSource *blend = sourceSet->data[i]; pmModel *model = modelSet->data[i]; + + if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); // correct model chisq for flux trend @@ -120,5 +127,4 @@ pmSourceCacheModel (blend, maskVal); pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal); - blend->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; blend->mode |= PM_SOURCE_MODE_BLEND_FIT; } @@ -144,5 +150,4 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; source->mode |= PM_SOURCE_MODE_BLEND_FIT; return true; @@ -167,4 +172,9 @@ // fit PSF model (set/unset the pixel mask) pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal); + + if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); + + // clear the circular mask + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // correct model chisq for flux trend @@ -186,10 +196,7 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; - return true; -} - -static float EXT_MIN_SN; -static float EXT_MOMENTS_RAD; + return true; +} + static pmModelType modelTypeEXT; @@ -197,8 +204,4 @@ bool status; - - // extended source model parameters - EXT_MIN_SN = psMetadataLookupF32 (&status, recipe, "EXT_MIN_SN"); - EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, recipe, "EXT_MOMENTS_RADIUS"); // extended source model descriptions @@ -221,9 +224,13 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) return false; if (source->type == PM_SOURCE_TYPE_SATURATED) return false; - if (source->peak->SN < EXT_MIN_SN) return false; - + + // set the radius based on the footprint (also sets the mask pixels) + float radius = psphotSetRadiusEXT (readout, source, markVal); + + // XXX note that this changes the source moments that are published... // recalculate the source moments using the larger extended-source moments radius // at this stage, skip Gaussian windowing, and do not clip pixels by S/N - if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false; + // this uses the footprint to judge both radius and aperture? + if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false; psTrace ("psphot", 5, "trying blob...\n"); @@ -237,4 +244,6 @@ // XXX need to handle failures better here pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal); + if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); + okEXT = psphotEvalEXT (tmpSrc, EXT); chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; @@ -246,8 +255,12 @@ // XXX should I keep / save the flags set in the eval functions? + // clear the circular mask + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + // correct first model chisqs for flux trend chiDBL = NAN; ONE = DBL->data[0]; if (ONE) { + if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); ONE->chisqNorm = ONE->chisq / chiTrend; @@ -258,4 +271,5 @@ ONE = DBL->data[1]; if (ONE) { + if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); ONE->chisqNorm = ONE->chisq / chiTrend; @@ -277,4 +291,5 @@ // both models failed; reject them both + // XXX -- change type flags to psf in this case and keep original moments? psFree (EXT); psFree (DBL); @@ -287,4 +302,5 @@ // save new model source->modelEXT = EXT; + source->modelEXT->fitRadius = radius; source->type = PM_SOURCE_TYPE_EXTENDED; source->mode |= PM_SOURCE_MODE_EXTMODEL; @@ -293,6 +309,15 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + +# if (PS_TRACE_ON) psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); + if (psTraceGetLevel("psphot") >= 6) { + psLogMsg ("psphot", 1, "source 2:\n"); + for (int i = 0; i < source->modelEXT->params->n; i++) { + psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelEXT->params->data.F32[i], source->modelEXT->dparams->data.F32[i]); + } + } +# endif + return true; @@ -304,11 +329,11 @@ psFree (source->modelPSF); source->modelPSF = psMemIncrRefCounter (DBL->data[0]); - source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; source->mode |= PM_SOURCE_MODE_PAIR; + source->modelPSF->fitRadius = radius; // copy most data from the primary source (modelEXT, blends stay NULL) - // XXX use pmSourceCopy? pmSource *newSrc = pmSourceCopy (source); newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); + newSrc->modelPSF->fitRadius = radius; // build cached models and subtract @@ -317,5 +342,19 @@ pmSourceCacheModel (newSrc, maskVal); pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal); + +# if (PS_TRACE_ON) psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]); + if (psTraceGetLevel("psphot") >= 6) { + psLogMsg ("psphot", 1, "source 1:\n"); + for (int i = 0; i < newSrc->modelPSF->params->n; i++) { + psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, newSrc->modelPSF->params->data.F32[i], newSrc->modelPSF->dparams->data.F32[i]); + } + psLogMsg ("psphot", 1, "source 2:\n"); + for (int i = 0; i < source->modelPSF->params->n; i++) { + psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]); + } + psphotVisualShowResidualImage (readout); + } +# endif psArrayAdd (newSources, 100, newSrc); @@ -356,5 +395,4 @@ // save the PSF model from the Ensemble fit PSF = source->modelPSF; - psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0); if (isnan(PSF->params->data.F32[1])) psAbort("nan in dbl fit"); @@ -389,8 +427,4 @@ PS_ASSERT (EXT, NULL); - // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit"); - - psphotCheckRadiusEXT (readout, source, EXT, markVal); - if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) { psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); @@ -401,3 +435,2 @@ return (EXT); } - Index: trunk/psphot/src/psphotSourcePlots.c =================================================================== --- trunk/psphot/src/psphotSourcePlots.c (revision 25738) +++ trunk/psphot/src/psphotSourcePlots.c (revision 25755) @@ -111,12 +111,14 @@ if (Xo == 0) { // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotRadialPlot (&kapa, "radial.plots.ps", source); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + Yo += DY; Xo = 0; @@ -126,11 +128,13 @@ Yo += dY; Xo = 0; - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotRadialPlot (&kapa, "radial.plots.ps", source); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -139,11 +143,12 @@ } else { // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotRadialPlot (&kapa, "radial.plots.ps", source); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 25738) +++ trunk/psphot/src/psphotSourceSize.c (revision 25755) @@ -2,9 +2,23 @@ # include -static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, - psImageMaskType maskVal, const pmModel *model, float Ro); - -bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); -bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); +typedef struct { + psImageMaskType maskVal; + psImageMaskType markVal; + psImageMaskType crMask; + float ApResid; + float ApSysErr; + float nSigmaApResid; + float nSigmaMoments; + float nSigmaCR; + float soft; + int grow; +} psphotSourceSizeOptions; + +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf); +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); +bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); +bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); +bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); // we need to call this function after sources have been fitted to the PSF model and @@ -14,178 +28,384 @@ // deviation from the psf model at the r = FWHM/2 position -bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first) +// XXX use an internal flag to mark sources which have already been measured +bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first) { bool status; + psphotSourceSizeOptions options; psTimerStart ("psphot.size"); // 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); + options.maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (options.maskVal); + + options.markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + assert (options.markVal); // bit to mask the cosmic-ray pixels - psImageMaskType crMask = pmConfigMaskGet("CR", config); // Mask value for cosmic rays - - float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT"); + options.crMask = pmConfigMaskGet("CR", config); // Mask value for cosmic rays + + options.nSigmaCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT"); assert (status); - float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); + // XXX recipe name is not great + options.nSigmaApResid = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); assert (status); - int grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs - if (!status || grow < 0) { + // XXX recipe name is not great + options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS"); + assert (status); + + options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs + if (!status || options.grow < 0) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.GROW is not positive."); return false; } - float soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter - if (!status || !isfinite(soft) || soft < 0.0) { + options.soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter + if (!status || !isfinite(options.soft) || options.soft < 0.0) { psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero."); - soft = 0.0; - } - - // loop over all source - for (int i = first; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip source if it was already measured - if (isfinite(source->crNsigma)) { - psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n"); - continue; + options.soft = 0.0; + } + + // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness + // of and object. We need to model this distribution for the PSF stars before we can test + // the significance for a specific object + // XXX move this to the code that generates the PSF? + // XXX store the results on pmPSF? + psphotSourceSizePSF (&options, sources, psf); + + // classify the sources based on ApResid and Moments (extended sources) + psphotSourceClass(readout, sources, recipe, psf, &options); + + psphotSourceSizeCR (readout, sources, &options); + + psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size")); + + psphotVisualPlotSourceSize (recipe, sources); + psphotVisualShowSourceSize (readout, sources); + psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr); + + return true; +} + +bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + // replace the source flux + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + // flag this as a CR + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + // grab the matching footprint + pmFootprint *footprint = peak->footprint; + if (!footprint) { + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + + if (!footprint->spans) { + // if we have no footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + + // mask all of the pixels covered by the spans of the footprint + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int iy = span1->y; + int xs = span1->x0; + int xe = span1->x1; + + for (int ix = xs; ix < xe; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + psImage *mask = source->maskView; + psImage *pixels = source->pixels; + psImage *variance = source->variance; + + // XXX This should be a recipe variable +# define SN_LIMIT 5.0 + + int xo = peak->x - pixels->col0; + int yo = peak->y - pixels->row0; + + // mark the pixels in this row to the left, then the right + for (int ix = xo; ix >= 0; ix--) { + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } + } + for (int ix = xo + 1; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } + } + + // for each of the neighboring rows, mark the high pixels if they have a marked neighbor + // first go up: + for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + // next go down: + for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +// model the apmifit distribution for the psf stars: +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) { + + // select stats from the psf stars + psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); + + psImageMaskType maskVal = options->maskVal | options->markVal; + + // XXX why PHOT_WEIGHT?? + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); } - // source must have been subtracted - if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + + psVectorAppend (Ap, 100, dMag); + psVectorAppend (ApErr, 100, source->errMag); + } + + // model the distribution as a mean or median value and a systematic error from that value: + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + psVectorStats (stats, Ap, NULL, NULL, 0); + + psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32); + for (int i = 0; i < Ap->n; i++) { + dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian; + } + + options->ApResid = stats->robustMedian; + options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05); + psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr); + + psFree (Ap); + psFree (ApErr); + psFree (stats); + psFree (dAp); + + return true; +} + +// classify sources based on the combination of psf-mag, Mxx, Myy +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { + + bool status; + pmPSFClump psfClump; + char regionName[64]; + + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nmiss", "Nskip"); + + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + for (int i = 0; i < nRegions; i ++) { + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion"); + + psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION"); + psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion"); + + // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); psAssert (status, "missing PSF.CLUMP.X"); + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); + + if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { + psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + + if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, options)) { + psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + } + + return true; +} + +bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { + + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + + int Nsat = 0; + int Next = 0; + int Npsf = 0; + int Ncr = 0; + int Nmiss = 0; + int Nskip = 0; + + pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + psImageMaskType maskVal = options->maskVal | options->markVal; + + for (psS32 i = 0 ; i < sources->n ; i++) { + + pmSource *source = (pmSource *) sources->data[i]; + + // psfClumps are found for image subregions: + // skip sources not in this region + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // source must have been subtracted + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n"); - continue; + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + continue; + } + + // we are basically classifying by moments; use the default if not found + psAssert (source->moments, "why is this source missing moments?"); + if (source->mode & noMoments) { + Nskip ++; + continue; + } + + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, options->maskVal); } - psF32 **resid = source->pixels->data.F32; - psF32 **variance = source->variance->data.F32; - psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; - - // check for extendedness: measure the delta flux significance at the 1 sigma contour - source->extNsigma = psphotModelContour(source->pixels, source->variance, source->maskObj, maskVal, - source->modelPSF, 1.0); - - // XXX prevent a source from being both CR and EXT? - if (source->extNsigma > EXT_NSIGMA_LIMIT) { - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - } - - // Integer position of peak - int xPeak = source->peak->xf - source->pixels->col0 + 0.5; - int yPeak = source->peak->yf - source->pixels->row0 + 0.5; - - // XXX for now, skip sources which are too close to a boundary - // XXX raise a flag? - if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || - yPeak < 1 || yPeak > source->pixels->numRows - 2) { - source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); - continue; - } - - // XXX for now, just skip any sources with masked pixels - // XXX raise a flag? - bool keep = true; - for (int iy = -1; (iy <= +1) && keep; iy++) { - for (int ix = -1; (ix <= +1) && keep; ix++) { - if (mask[yPeak+iy][xPeak+ix] & maskVal) { - keep = false; - } - } - } - if (!keep) { - psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); - source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - continue; - } - - // Compare the central pixel with those on either side, for the four possible lines through it. - - // Soften variances (add systematic error) - float softening = soft * PS_SQR(source->peak->flux); // Softening for variances - - // Across the middle: y = 0 - float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; - float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; - float nX = cX / sqrtf(dcX + softening); - - // Up the centre: x = 0 - float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; - float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; - float nY = cY / sqrtf(dcY + softening); - - // Diagonal: x = y - float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; - float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; - float nL = cL / sqrtf(dcL + softening); - - // Diagonal: x = - y - float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; - float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; - float nR = cR / sqrtf(dcR + softening); - - // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) - // Ndof = 4 ? (four measurements, no free parameters) - // XXX this value is going to be biased low because of systematic errors. - // we need to calibrate it somehow - // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); - - // not strictly accurate: overcounts the chisq contribution from the center pixel (by - // factor of 4); also biases a bit low if any pixels are masked - // XXX I am not sure I want to keep this value... - source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); - - float fCR = 0.0; - int nCR = 0; - if (nX > 0.0) { - fCR += nX; - nCR ++; - } - if (nY > 0.0) { - fCR += nY; - nCR ++; - } - if (nL > 0.0) { - fCR += nL; - nCR ++; - } - if (nR > 0.0) { - fCR += nR; - nCR ++; - } - source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; - if (!isfinite(source->crNsigma)) { - continue; - } - - // this source is thought to be a cosmic ray. flag the detection and mask the pixels - if (source->crNsigma > CR_NSIGMA_LIMIT) { - // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); - psphotMaskCosmicRay_Old (source, maskVal, crMask); - } - } - - // now that we have masked pixels associated with CRs, we can grow the mask - if (grow > 0) { - bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask - psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow); - psImageConvolveSetThreads(oldThreads); - if (!newMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); - return false; - } - psFree(readout->mask); - readout->mask = newMask; - } - - psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", - sources->n - first, psTimerMark ("psphot.size")); - - psphotVisualPlotSourceSize (sources); - psphotVisualShowSourceSize (readout, sources); + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); + + 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; + + // 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); + if (isPSF) { + Npsf ++; + continue; + } + + // Defects may not always match CRs from peak curvature analysis + // Defects may also be marked as SATSTAR -- XXX deactivate this flag? + // XXX this rule is not great + if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + + // saturated star (determined in PSF fit). These may also be saturated galaxies, or + // just large saturated regions. + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + Nsat ++; + continue; + } + + // 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)); + if (isEXT) { + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + Next ++; + continue; + } + + psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); + Nmiss ++; + } + + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4d %4d %4d %4d %4d %4d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip); return true; @@ -194,7 +414,8 @@ // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total // deviation in sigmas. This is measured on the residual image - should we ignore negative -// deviations? -static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, - psImageMaskType maskVal, const pmModel *model, float Ro) +// deviations? NOTE: This function was an early attempt to classify extended objects, and is +// no longer used by psphot. +float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, + psImageMaskType maskVal, const pmModel *model, float Ro) { psF32 *PAR = model->params->data.F32; // Model parameters @@ -211,6 +432,6 @@ float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); if (Q < 0.0) { - // ellipse is imaginary - return NAN; + // ellipse is imaginary + return NAN; } @@ -220,44 +441,44 @@ for (int x = -radius; x <= radius; x++) { - // Polynomial coefficients - // XXX Should we be using the centre of the pixel as x or x+0.5? - float A = PS_SQR (1.0 / syy); - float B = x * sxy; - float C = PS_SQR (x / sxx) - Ro; - float T = PS_SQR(B) - 4*A*C; - if (T < 0.0) { - continue; - } - - // y position in source frame - float yP = (-B + sqrt (T)) / (2.0 * A); - float yM = (-B - sqrt (T)) / (2.0 * A); - - // Get the closest pixel positions (image frame) - int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; - int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - - if (xPix < 0 || xPix >= image->numCols) { - continue; - } - - if (yPixM >= 0 && yPixM < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); - nSigma += dSigma; - nPts++; - } - - if (yPixM == yPixP) { - continue; - } - - if (yPixP >= 0 && yPixP < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); - nSigma += dSigma; - nPts++; - } + // Polynomial coefficients + // XXX Should we be using the centre of the pixel as x or x+0.5? + float A = PS_SQR (1.0 / syy); + float B = x * sxy; + float C = PS_SQR (x / sxx) - Ro; + float T = PS_SQR(B) - 4*A*C; + if (T < 0.0) { + continue; + } + + // y position in source frame + float yP = (-B + sqrt (T)) / (2.0 * A); + float yM = (-B - sqrt (T)) / (2.0 * A); + + // Get the closest pixel positions (image frame) + int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; + int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + + if (xPix < 0 || xPix >= image->numCols) { + continue; + } + + if (yPixM >= 0 && yPixM < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts++; + } + + if (yPixM == yPixP) { + continue; + } + + if (yPixP >= 0 && yPixP < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts++; + } } nSigma /= nPts; @@ -265,107 +486,134 @@ } -bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { - - // replace the source flux - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; - - // flag this as a CR - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - // grab the matching footprint - pmFootprint *footprint = peak->footprint; - if (!footprint) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); - return true; - } - - if (!footprint->spans) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); - return true; - } - - // mask all of the pixels covered by the spans of the footprint - for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int iy = span1->y; - int xs = span1->x0; - int xe = span1->x1; - - for (int ix = xs; ix < xe; ix++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } +bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { + + // classify the sources based on the CR test (place this in a function?) + // XXX use an internal flag to mark sources which have already been measured + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // source must have been subtracted + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + continue; + } + + psF32 **resid = source->pixels->data.F32; + psF32 **variance = source->variance->data.F32; + psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + + // Integer position of peak + int xPeak = source->peak->xf - source->pixels->col0 + 0.5; + int yPeak = source->peak->yf - source->pixels->row0 + 0.5; + + // Skip sources which are too close to a boundary. These are mostly caught as DEFECT + if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || + yPeak < 1 || yPeak > source->pixels->numRows - 2) { + psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); + continue; + } + + // Skip sources with masked pixels. These are mostly caught as DEFECT + bool keep = true; + for (int iy = -1; (iy <= +1) && keep; iy++) { + for (int ix = -1; (ix <= +1) && keep; ix++) { + if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { + keep = false; + } + } + } + if (!keep) { + psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); + continue; + } + + // Compare the central pixel with those on either side, for the four possible lines through it. + + // Soften variances (add systematic error) + float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances + + // Across the middle: y = 0 + float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; + float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; + float nX = cX / sqrtf(dcX + softening); + + // Up the centre: x = 0 + float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; + float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; + float nY = cY / sqrtf(dcY + softening); + + // Diagonal: x = y + float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; + float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; + float nL = cL / sqrtf(dcL + softening); + + // Diagonal: x = - y + float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; + float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; + float nR = cR / sqrtf(dcR + softening); + + // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) + // Ndof = 4 ? (four measurements, no free parameters) + // XXX this value is going to be biased low because of systematic errors. + // we need to calibrate it somehow + // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); + + // not strictly accurate: overcounts the chisq contribution from the center pixel (by + // factor of 4); also biases a bit low if any pixels are masked + // XXX I am not sure I want to keep this value... + source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); + + float fCR = 0.0; + int nCR = 0; + if (nX > 0.0) { + fCR += nX; + nCR ++; + } + if (nY > 0.0) { + fCR += nY; + nCR ++; + } + if (nL > 0.0) { + fCR += nL; + nCR ++; + } + if (nR > 0.0) { + fCR += nR; + nCR ++; + } + source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + if (!isfinite(source->crNsigma)) { + continue; + } + + // this source is thought to be a cosmic ray. flag the detection and mask the pixels + if (source->crNsigma > options->nSigmaCR) { + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); + // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); + } + } + + // now that we have masked pixels associated with CRs, we can grow the mask + if (options->grow > 0) { + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask + psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow); + psImageConvolveSetThreads(oldThreads); + if (!newMask) { + psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); + return false; + } + psFree(readout->mask); + readout->mask = newMask; } return true; } - -bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { - - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - psImage *mask = source->maskView; - psImage *pixels = source->pixels; - psImage *variance = source->variance; - - // XXX This should be a recipe variable -# define SN_LIMIT 5.0 - - int xo = peak->x - pixels->col0; - int yo = peak->y - pixels->row0; - - // mark the pixels in this row to the left, then the right - for (int ix = xo; ix >= 0; ix--) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } - } - for (int ix = xo + 1; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } - } - - // for each of the neighboring rows, mark the high pixels if they have a marked neighbor - // first go up: - for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } - } - // next go down: - for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } - } - return true; -} Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 25738) +++ trunk/psphot/src/psphotSourceStats.c (revision 25755) @@ -1,5 +1,7 @@ # include "psphotInternal.h" -psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) { +bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources); + +psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) { bool status = false; @@ -21,4 +23,7 @@ float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); if (!status) return NULL; + + OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius + char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); if (!status) return NULL; @@ -60,4 +65,11 @@ psphotVisualShowMoments (sources); return sources; + } + + if (setWindow) { + if (!psphotSetMomentsWindow(recipe, sources)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); + return NULL; + } } @@ -144,6 +156,4 @@ float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); if (!status) return false; - float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN"); - if (!status) return false; float SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); if (!status) return false; @@ -194,6 +204,6 @@ } - // measure basic source moments - status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN); + // measure basic source moments (no S/N clipping on input pixels) + status = pmSourceMoments (source, RADIUS, SIGMA, 0.0); if (status) { Nmoments ++; @@ -205,5 +215,5 @@ BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN); + status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0); if (status) { source->mode |= PM_SOURCE_MODE_BIG_RADIUS; @@ -231,84 +241,91 @@ } -# if (0) -bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) { - - bool status = false; - float BIG_RADIUS; - - float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); - if (!status) return false; - float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); - if (!status) return false; - float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); - if (!status) return false; - - // bit-masks to test for good/bad pixels - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); - assert (maskVal); - - // bit-mask to mark pixels not used in analysis - psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); - assert (markVal); - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // threaded measurement of the sources moments - int Nfail = 0; - int Nmoments = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { - continue; - } - - // measure a local sky value - // the local sky is now ignored; kept here for reference only - status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); - if (!status) { - psErrorClear(); // XXX re-consider the errors raised here - Nfail ++; - continue; - } - - // measure the local sky variance (needed if noise is not sqrt(signal)) - // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken - status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); - if (!status) { - Nfail ++; - psErrorClear(); - continue; - } - - // measure basic source moments - status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN); - if (status) { - Nmoments ++; - continue; - } - - // if no valid pixels, or massive swing, likely saturated source, - // try a much larger box - BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); - psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN); - if (status) { - Nmoments ++; - continue; - } - - Nfail ++; - psErrorClear(); - continue; - } - - // change the value of a scalar on the array (wrap this and put it in psArray.h) - *nmoments = Nmoments; - *nfail = Nfail; - +// this function attempts to iteratively determine the best value for sigma of the moments weighting Gaussian +bool psphotSetMomentsWindow (psMetadata *recipe, psArray *sources) { + + bool status; + + float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); + if (!status) return false; + + // XXX move this to a config file? + float sigma[4] = {1.0, 2.0, 3.0, 4.5}; + float Sout[4]; + + // this sorts by peak->SN + sources = psArraySort (sources, pmSourceSortBySN); + + // loop over radii: + for (int i = 0; i < 4; i++) { + + // XXX move max source number to config + for (int j = 0; (j < sources->n) && (j < 400); j++) { + + pmSource *source = sources->data[j]; + psAssert (source->moments, "force moments to exist"); + source->moments->nPixels = 0; + + // skip faint sources for moments measurement + if (source->peak->SN < MIN_SN) { + continue; + } + + // measure basic source moments (no S/N clipping on input pixels) + status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0); + } + + // choose a grid scale that is a fixed fraction of the psf sigma^2 + psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i])); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); + + // determine the PSF parameters from the source moment values + pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe); + psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); + + psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000"); + if (!regionMD) { + regionMD = psMetadataAlloc(); + psMetadataAddMetadata (recipe, PS_LIST_TAIL, "PSF.CLUMP.REGION.000", PS_META_REPLACE, "psf clump region", regionMD); + psFree (regionMD); + } + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); + + // psphotVisualPlotMoments (recipe, sources); + + Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; + } + + // we are looking for sigma for which Sout = 0.65 (or some other value) + + float Sigma = NAN; + float minS = Sout[0]; + float maxS = Sout[0]; + for (int i = 0; i < 4; i++) { + minS = PS_MIN(Sout[i], minS); + maxS = PS_MAX(Sout[i], maxS); + } + if (minS > 0.65) Sigma = sigma[3]; + if (maxS < 0.65) Sigma = sigma[0]; + + for (int i = 0; i < 3; i++) { + if ((Sout[i] > 0.65) && (Sout[i+1] > 0.65)) continue; + if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue; + Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]); + } + psAssert (isfinite(Sigma), "did we miss a case?"); + + // choose a grid scale that is a fixed fraction of the psf sigma^2 + psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma)); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma)); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma)); + psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma); + psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma); + + psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma); return true; } -# endif Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 25738) +++ trunk/psphot/src/psphotVisual.c (revision 25755) @@ -15,11 +15,48 @@ # include +bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec); + // functions used to visualize the analysis as it goes // these are invoked by the -visual options -static int kapa = -1; +static int kapa1 = -1; static int kapa2 = -1; static int kapa3 = -1; +int psphotKapaChannel (int channel) { + + switch (channel) { + case 1: + if (kapa1 == -1) { + kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images"); + if (kapa1 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa1; + case 2: + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa2; + case 3: + if (kapa3 == -1) { + kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps"); + if (kapa3 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + } + } + return kapa3; + default: + psAbort ("unknown kapa channel"); + } + psAbort ("unknown kapa channel"); +} bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) { @@ -131,12 +168,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; // psphotVisualShowMask (kapa, readout->mask, "mask", 2); @@ -144,11 +175,5 @@ psphotVisualScaleImage (kapa, readout->image, "image", 0); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -160,12 +185,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; bool status = false; @@ -181,11 +200,5 @@ psphotVisualScaleImage (kapa, readout->image, "backsub", 0); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -195,23 +208,10 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - // XXX test: image->data.F32[10][10] = 10000; + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + psphotVisualRangeImage (kapa, image, "signif", 2, -1.0, 25.0*25.0); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -224,8 +224,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; psArray *peaks = detections->peaks; @@ -233,5 +231,5 @@ // note: this uses the Ohana allocation tools: // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1); - ALLOCATE (overlay, KiiOverlay, peaks->n); + ALLOCATE (overlay, KiiOverlay, peaks->n + 2); Noverlay = 0; @@ -271,10 +269,10 @@ } -# if (0) +# if (1) overlay[Noverlay].type = KII_OVERLAY_BOX; overlay[Noverlay].x = 10.0; overlay[Noverlay].y = 10.0; - overlay[Noverlay].dx = 0.5; - overlay[Noverlay].dy = 0.5; + overlay[Noverlay].dx = 1.0; + overlay[Noverlay].dy = 1.0; overlay[Noverlay].angle = 0.0; overlay[Noverlay].text = NULL; @@ -285,11 +283,5 @@ FREE (overlay); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -302,8 +294,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; psArray *footprints = detections->footprints; @@ -325,4 +315,5 @@ // draw the top + // XXX need to allow top (and bottom) to have more than one span span = footprint->spans->data[0]; overlay[Noverlay].type = KII_OVERLAY_LINE; @@ -399,11 +390,5 @@ FREE (overlay); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -419,8 +404,9 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + // XXX mark the different source classes with different color/shape dots + // XXX are moments S/N and peak S/N consistent? // note: this uses the Ohana allocation tools: @@ -448,5 +434,7 @@ overlay[Noverlay].dx = 2.0*axes.major; overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; // XXXXXXXX the axes angle is negative to display of object on kapa + + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; Noverlay ++; @@ -456,12 +444,5 @@ FREE (overlay); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - + pmVisualAskUser(NULL); return true; } @@ -475,22 +456,16 @@ if (!pmVisualIsVisual()) return true; - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - KapaClearPlots (kapa3); + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); KapaInitGraph (&graphdata); - KapaSetFont (kapa3, "courier", 14); + KapaSetFont (myKapa, "courier", 14); float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); // select the max psfX,Y values for the plot limits - float Xmin = 0.0, Xmax = 0.0; - float Ymin = 0.0, Ymax = 0.0; + float Xmin = 1000.0, Xmax = 0.0; + float Ymin = 1000.0, Ymax = 0.0; { int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); @@ -511,10 +486,12 @@ float Y1 = psfY + 4.0*psfdY; - if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); } + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } - if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } } } + Xmin = PS_MAX(Xmin, -0.1); + Ymin = PS_MAX(Ymin, -0.1); // storage vectors for data to be plotted @@ -564,5 +541,5 @@ section.y = 0.00; section.name = psStringCopy ("MxxMyy"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -572,9 +549,9 @@ graphdata.xmax = Xmax; graphdata.ymax = Ymax; - KapaSetLimits (kapa3, &graphdata); - - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM); + KapaSetLimits (myKapa, &graphdata); + + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -582,7 +559,7 @@ graphdata.size = 0.3; graphdata.style = 2; - KapaPrepPlot (kapa3, nF, &graphdata); - KapaPlotVector (kapa3, nF, xFaint->data.F32, "x"); - KapaPlotVector (kapa3, nF, yFaint->data.F32, "y"); + KapaPrepPlot (myKapa, nF, &graphdata); + KapaPlotVector (myKapa, nF, xFaint->data.F32, "x"); + KapaPlotVector (myKapa, nF, yFaint->data.F32, "y"); graphdata.color = KapaColorByName ("red"); @@ -590,7 +567,7 @@ graphdata.size = 0.5; graphdata.style = 2; - KapaPrepPlot (kapa3, nB, &graphdata); - KapaPlotVector (kapa3, nB, xBright->data.F32, "x"); - KapaPlotVector (kapa3, nB, yBright->data.F32, "y"); + KapaPrepPlot (myKapa, nB, &graphdata); + KapaPlotVector (myKapa, nB, xBright->data.F32, "x"); + KapaPlotVector (myKapa, nB, yBright->data.F32, "y"); // second section: MagMyy @@ -600,5 +577,5 @@ section.y = 0.80; section.name = psStringCopy ("MagMyy"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -608,10 +585,10 @@ graphdata.ymin = Ymin; graphdata.ymax = Ymax; - KapaSetLimits (kapa3, &graphdata); + KapaSetLimits (myKapa, &graphdata); strcpy (graphdata.labels, "0210"); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP); - KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -619,7 +596,7 @@ graphdata.size = 0.3; graphdata.style = 2; - KapaPrepPlot (kapa3, nF, &graphdata); - KapaPlotVector (kapa3, nF, mFaint->data.F32, "x"); - KapaPlotVector (kapa3, nF, yFaint->data.F32, "y"); + KapaPrepPlot (myKapa, nF, &graphdata); + KapaPlotVector (myKapa, nF, mFaint->data.F32, "x"); + KapaPlotVector (myKapa, nF, yFaint->data.F32, "y"); graphdata.color = KapaColorByName ("red"); @@ -627,7 +604,7 @@ graphdata.size = 0.5; graphdata.style = 2; - KapaPrepPlot (kapa3, nB, &graphdata); - KapaPlotVector (kapa3, nB, mBright->data.F32, "x"); - KapaPlotVector (kapa3, nB, yBright->data.F32, "y"); + KapaPrepPlot (myKapa, nB, &graphdata); + KapaPlotVector (myKapa, nB, mBright->data.F32, "x"); + KapaPlotVector (myKapa, nB, yBright->data.F32, "y"); // third section: MagMxx @@ -637,5 +614,5 @@ section.y = 0.00; section.name = psStringCopy ("MagMxx"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -645,10 +622,10 @@ graphdata.ymin = -7.9; graphdata.ymax = -17.1; - KapaSetLimits (kapa3, &graphdata); + KapaSetLimits (myKapa, &graphdata); strcpy (graphdata.labels, "2001"); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP); graphdata.color = KapaColorByName ("black"); @@ -656,7 +633,7 @@ graphdata.size = 0.3; graphdata.style = 2; - KapaPrepPlot (kapa3, nF, &graphdata); - KapaPlotVector (kapa3, nF, xFaint->data.F32, "x"); - KapaPlotVector (kapa3, nF, mFaint->data.F32, "y"); + KapaPrepPlot (myKapa, nF, &graphdata); + KapaPlotVector (myKapa, nF, xFaint->data.F32, "x"); + KapaPlotVector (myKapa, nF, mFaint->data.F32, "y"); graphdata.color = KapaColorByName ("red"); @@ -664,11 +641,11 @@ graphdata.size = 0.5; graphdata.style = 2; - KapaPrepPlot (kapa3, nB, &graphdata); - KapaPlotVector (kapa3, nB, xBright->data.F32, "x"); - KapaPlotVector (kapa3, nB, mBright->data.F32, "y"); + KapaPrepPlot (myKapa, nB, &graphdata); + KapaPlotVector (myKapa, nB, xBright->data.F32, "x"); + KapaPlotVector (myKapa, nB, mBright->data.F32, "y"); // draw N circles to outline the clumps { - KapaSelectSection (kapa3, "MxxMyy"); + KapaSelectSection (myKapa, "MxxMyy"); // draw a circle centered on psfX,Y with size of the psf limit @@ -686,5 +663,5 @@ graphdata.xmax = Xmax; graphdata.ymax = Ymax; - KapaSetLimits (kapa3, &graphdata); + KapaSetLimits (myKapa, &graphdata); for (int n = 0; n < nRegions; n++) { @@ -705,47 +682,11 @@ yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; } - KapaPrepPlot (kapa3, xLimit->n, &graphdata); - KapaPlotVector (kapa3, xLimit->n, xLimit->data.F32, "x"); - KapaPlotVector (kapa3, yLimit->n, yLimit->data.F32, "y"); + KapaPrepPlot (myKapa, xLimit->n, &graphdata); + KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); } psFree (xLimit); psFree (yLimit); } - -# if (0) - // *** make a histogram of the source counts in the x and y directions - psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0); - psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0); - psVectorHistogram (nX, xFaint, NULL, NULL, 0); - psVectorHistogram (nY, yFaint, NULL, NULL, 0); - psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); - psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); - psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); - psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); - for (int i = 0; i < nX->nums->n; i++) { - dX->data.F32[i] = nX->nums->data.S32[i]; - vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); - } - for (int i = 0; i < nY->nums->n; i++) { - dY->data.F32[i] = nY->nums->data.S32[i]; - vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); - } - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 0; - graphdata.size = 0.0; - graphdata.style = 0; - KapaPrepPlot (kapa3, dX->n, &graphdata); - KapaPlotVector (kapa3, dX->n, dX->data.F32, "x"); - KapaPlotVector (kapa3, vX->n, vX->data.F32, "y"); - - psFree (nX); - psFree (dX); - psFree (vX); - - psFree (nY); - psFree (dY); - psFree (vY); -# endif psFree (xBright); @@ -756,16 +697,10 @@ psFree (mFaint); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } // assumes 'kapa' value is checked and set -bool psphotVisualShowRoughClass_Single (psArray *sources, pmSourceType type, pmSourceMode mode, char *color) { +bool psphotVisualShowRoughClass_Single (int myKapa, psArray *sources, pmSourceType type, pmSourceMode mode, char *color) { int Noverlay; @@ -802,10 +737,10 @@ overlay[Noverlay].dx = 2.0*axes.major; overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = -axes.theta * PS_DEG_RAD; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; overlay[Noverlay].text = NULL; Noverlay ++; } - KiiLoadOverlay (kapa, overlay, Noverlay, color); + KiiLoadOverlay (myKapa, overlay, Noverlay, color); FREE (overlay); @@ -817,26 +752,18 @@ if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } - - KiiEraseOverlay (kapa, "yellow"); // moments - - psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, 0, "red"); - psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue"); - psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_DEFECT, 0, "blue"); - psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_SATURATED, 0, "red"); - psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow"); - psphotVisualShowRoughClass_Single (sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green"); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; + + KiiEraseOverlay (myKapa, "yellow"); // moments + + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, 0, "red"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_EXTENDED, 0, "blue"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_DEFECT, 0, "blue"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_SATURATED, 0, "red"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_PSFSTAR, "yellow"); + psphotVisualShowRoughClass_Single (myKapa, sources, PM_SOURCE_TYPE_STAR, PM_SOURCE_MODE_SATSTAR, "green"); + fprintf (stdout, "red: STAR or SAT AREA; blue: EXTENDED or DEFECT; green: SATSTAR; yellow: PSFSTAR\n"); - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -846,12 +773,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (3); + if (myKapa == -1) return false; int DX = 64; @@ -898,7 +819,7 @@ psImage *psfLogFlux = (psImage *) psUnaryOp (NULL, psfMosaic, "log"); - psphotVisualRangeImage (kapa2, psfLogFlux, "psf_mosaic", 0, -2.0, 3.0); - psphotVisualRangeImage (kapa2, funMosaic, "psf_analytical", 1, -10.0, 100.0); - psphotVisualRangeImage (kapa2, resMosaic, "psf_residual", 2, -10.0, 100.0); + psphotVisualRangeImage (myKapa, psfLogFlux, "psf_mosaic", 0, -2.0, 3.0); + psphotVisualRangeImage (myKapa, funMosaic, "psf_analytical", 1, -10.0, 100.0); + psphotVisualRangeImage (myKapa, resMosaic, "psf_residual", 2, -10.0, 100.0); psFree (psfMosaic); @@ -908,11 +829,5 @@ psFree (modelRef); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + pmVisualAskUser(NULL); return true; } @@ -924,12 +839,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:psfstars"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (3); + if (myKapa == -1) return false; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1017,10 +926,12 @@ if (Xo == 0) { // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // reset source Add/Sub state to recorded + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Yo += DY; @@ -1032,10 +943,12 @@ Xo = 0; - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -1044,10 +957,11 @@ } else { // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - psphotSubWithTest (source, false, maskVal); // remove source (force) + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - psphotSetState (source, false, maskVal); // replace source (has been subtracted) + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; @@ -1056,15 +970,8 @@ } - psphotVisualRangeImage (kapa2, outpos, "psfpos", 0, -0.05, 0.95); - psphotVisualRangeImage (kapa2, outsub, "psfsub", 1, -0.05, 0.95); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - + psphotVisualRangeImage (myKapa, outpos, "psfpos", 0, -0.05, 0.95); + psphotVisualRangeImage (myKapa, outsub, "psfsub", 1, -0.05, 0.95); + + pmVisualAskUser(NULL); psFree (outpos); psFree (outsub); @@ -1084,12 +991,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (3); + if (myKapa == -1) return false; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1119,7 +1020,7 @@ pmSource *source = sources->data[i]; - bool keep = false; - keep |= (source->mode & PM_SOURCE_MODE_SATSTAR); - if (!keep) continue; + // only show "real" saturated stars (not defects) + if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; // how does this subimage get placed into the output image? @@ -1164,10 +1065,8 @@ pmSource *source = sources->data[i]; - bool keep = false; - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - nSAT ++; - keep = true; - } - if (!keep) continue; + // only show "real" saturated stars (not defects) + if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; + nSAT ++; if (Xo + DX > NX) { @@ -1175,7 +1074,8 @@ if (Xo == 0) { // place source alone on this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // reset source Add/Sub state to recorded + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Yo += DY; @@ -1186,7 +1086,9 @@ Yo += dY; Xo = 0; - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // replace source (has been subtracted) + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Xo = DX; @@ -1195,7 +1097,8 @@ } else { // extend this row - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - psphotSetState (source, true, maskVal); // replace source (has been subtracted) + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); Xo += DX; @@ -1204,14 +1107,7 @@ } - psphotVisualScaleImage (kapa2, outsat, "satstar", 2); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - + psphotVisualScaleImage (myKapa, outsat, "satstar", 2); + + pmVisualAskUser(NULL); psFree (outsat); return true; @@ -1222,6 +1118,6 @@ Graphdata graphdata; - bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - psphotAddWithTest (source, true, maskVal); // replace source if subtracted + 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; @@ -1240,12 +1136,12 @@ for (int ix = 0; ix < source->pixels->numCols; ix++) { if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { - // 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] = 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] = 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] = log10(source->pixels->data.F32[iy][ix]); @@ -1256,5 +1152,5 @@ // reset source Add/Sub state to recorded - psphotSetState (source, state, maskVal); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); KapaInitGraph (&graphdata); @@ -1337,12 +1233,6 @@ if (!pmVisualIsVisual()) return true; - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -1351,5 +1241,5 @@ assert (maskVal); - KapaClearPlots (kapa3); + KapaClearPlots (myKapa); // first section : mag vs CR nSigma section.dx = 1.0; @@ -1359,5 +1249,5 @@ section.name = NULL; psStringAppend (§ion.name, "linlog"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -1369,5 +1259,5 @@ section.name = NULL; psStringAppend (§ion.name, "loglog"); - KapaSetSection (kapa3, §ion); + KapaSetSection (myKapa, §ion); psFree (section.name); @@ -1378,5 +1268,5 @@ if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - psphotVisualPlotRadialProfile (kapa3, source, maskVal); + psphotVisualPlotRadialProfile (myKapa, source, maskVal); // pause and wait for user input: @@ -1388,5 +1278,5 @@ } if (key[0] == 'e') { - KapaClearPlots (kapa3); + KapaClearPlots (myKapa); } if (key[0] == 's') { @@ -1407,10 +1297,11 @@ psEllipseAxes axes; + // XXX skip this for now: it is not very clear + return true; + if (!pmVisualIsVisual()) return true; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; // note: this uses the Ohana allocation tools: @@ -1485,40 +1376,28 @@ } - KiiLoadOverlay (kapa, overlayE, NoverlayE, "red"); - KiiLoadOverlay (kapa, overlayO, NoverlayO, "yellow"); + KiiLoadOverlay (myKapa, overlayE, NoverlayE, "red"); + KiiLoadOverlay (myKapa, overlayO, NoverlayO, "yellow"); FREE (overlayE); FREE (overlayO); - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; fprintf (stdout, "even bits (0x0001, 0x0004, ... : red\n"); fprintf (stdout, "odd bits (0x0002, 0x0008, ... : yellow\n"); - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - - return true; -} - -bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) { - - int Noverlay, NOVERLAY; + pmVisualAskUser(NULL); + + return true; +} + +bool psphotVisualShowSourceSize_Single (int myKapa, psArray *sources, pmSourceMode mode, bool keep, float scale, char *color) { + + int Noverlay; KiiOverlay *overlay; - if (!pmVisualIsVisual()) return true; - - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } + psEllipseMoments emoments; + psEllipseAxes axes; // note: this uses the Ohana allocation tools: + ALLOCATE (overlay, KiiOverlay, sources->n); + Noverlay = 0; - NOVERLAY = 100; - ALLOCATE (overlay, KiiOverlay, sources->n); - - // mark CRs with red boxes for (int i = 0; i < sources->n; i++) { @@ -1526,60 +1405,65 @@ if (source == NULL) continue; - if (!(source->mode & PM_SOURCE_MODE_CR_LIMIT)) continue; - - overlay[Noverlay].type = KII_OVERLAY_BOX; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - - overlay[Noverlay].dx = 4; - overlay[Noverlay].dy = 4; - overlay[Noverlay].angle = 0; + // if (source->type != type) continue; + if (mode) { + if (keep) { + if (!(source->mode & mode)) continue; + } else { + if (source->mode & mode) continue; + } + } + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->Mx; + overlay[Noverlay].y = moments->My; + + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = scale*2.0*axes.major; + overlay[Noverlay].dy = scale*2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; overlay[Noverlay].text = NULL; Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - } - KiiLoadOverlay (kapa, overlay, Noverlay, "red"); - - - Noverlay = 0; - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - // mark EXTs with yellow circles - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - - overlay[Noverlay].dx = 10; - overlay[Noverlay].dy = 10; - overlay[Noverlay].angle = 0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - } - - KiiLoadOverlay (kapa, overlay, Noverlay, "blue"); + } + + KiiLoadOverlay (myKapa, overlay, Noverlay, color); FREE (overlay); - psphotVisualShowMask (kapa, readout->mask, "mask", 2); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "CR: 4pix red BOX; EXT: 10pix blue circle\n"); - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - - return true; -} - -bool psphotVisualPlotSourceSize (psArray *sources) { - + return true; +} + +bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) { + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; + + KiiEraseOverlay (myKapa, "red"); + KiiEraseOverlay (myKapa, "blue"); + KiiEraseOverlay (myKapa, "green"); + KiiEraseOverlay (myKapa, "yellow"); + + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_SATSTAR, 0, 1.0, "green"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_EXT_LIMIT, 1, 1.0, "blue"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_CR_LIMIT, 1, 1.0, "red"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_DEFECT, 1, 2.0, "red"); + psphotVisualShowSourceSize_Single (myKapa, sources, PM_SOURCE_MODE_SATSTAR, 1, 1.0, "yellow"); + + fprintf (stdout, "red: CR; blue: EXTENDED; green: PSF-like; yellow: SATSTAR\n"); + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources) { + + bool status; Graphdata graphdata; KapaSection section; @@ -1587,185 +1471,484 @@ if (!pmVisualIsVisual()) return true; - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); + KapaInitGraph (&graphdata); + KapaSetFont (myKapa, "courier", 14); + + // select the max psfX,Y values for the plot limits + float Xmin = 1000.0, Xmax = 0.0; + float Ymin = 1000.0, Ymax = 0.0; + { + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + + float X0 = psfX - 10.0*psfdX; + float X1 = psfX + 10.0*psfdX; + float Y0 = psfY - 10.0*psfdY; + float Y1 = psfY + 10.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } } } - - KapaClearPlots (kapa3); + Xmin = PS_MAX(Xmin, -0.1); + Ymin = PS_MAX(Ymin, -0.1); + + // storage vectors for data to be plotted + psVector *xSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *ySAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sSAT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sPSF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sEXT = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sDEF = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + psVector *xCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *yCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *mCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *sCR = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + + // construct the vectors + int nSAT = 0; + int nEXT = 0; + int nPSF = 0; + int nDEF = 0; + int nCR = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->moments == NULL) continue; + + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { + xCR->data.F32[nCR] = source->moments->Mxx; + yCR->data.F32[nCR] = source->moments->Myy; + mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum); + sCR->data.F32[nCR] = source->extNsigma; + nCR++; + } + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + xSAT->data.F32[nSAT] = source->moments->Mxx; + ySAT->data.F32[nSAT] = source->moments->Myy; + mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum); + sSAT->data.F32[nSAT] = source->extNsigma; + nSAT++; + } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + xEXT->data.F32[nEXT] = source->moments->Mxx; + yEXT->data.F32[nEXT] = source->moments->Myy; + mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum); + sEXT->data.F32[nEXT] = source->extNsigma; + nEXT++; + continue; + } + if (source->mode & PM_SOURCE_MODE_DEFECT) { + xDEF->data.F32[nDEF] = source->moments->Mxx; + yDEF->data.F32[nDEF] = source->moments->Myy; + mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum); + sDEF->data.F32[nDEF] = source->extNsigma; + nDEF++; + continue; + } + if ((source->mode & PM_SOURCE_MODE_CR_LIMIT) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { + continue; + } + xPSF->data.F32[nPSF] = source->moments->Mxx; + yPSF->data.F32[nPSF] = source->moments->Myy; + mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum); + sPSF->data.F32[nPSF] = source->extNsigma; + nPSF++; + } + xSAT->n = nSAT; + ySAT->n = nSAT; + mSAT->n = nSAT; + sSAT->n = nSAT; + + xPSF->n = nPSF; + yPSF->n = nPSF; + mPSF->n = nPSF; + sPSF->n = nPSF; + + xEXT->n = nEXT; + yEXT->n = nEXT; + mEXT->n = nEXT; + sEXT->n = nEXT; + + xCR->n = nCR; + yCR->n = nCR; + mCR->n = nCR; + sCR->n = nCR; + + xDEF->n = nDEF; + yDEF->n = nDEF; + mDEF->n = nDEF; + sDEF->n = nDEF; + + // four sections: MxxMyy, MagMxx, MagMyy, MagSigma + + // first section: MxxMyy + section.dx = 0.75; + section.dy = 0.60; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("MxxMyy"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, xCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, yCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y"); + + // second section: MagMyy + section.dx = 0.75; + section.dy = 0.20; + section.x = 0.00; + section.y = 0.80; + section.name = psStringCopy ("MagMyy"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -17.1; + graphdata.xmax = -7.9; + graphdata.ymin = Ymin; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + strcpy (graphdata.labels, "0210"); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XP); + KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, yPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, yEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, yDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, mCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, yCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, ySAT->data.F32, "y"); + + // third section: MagMxx + section.dx = 0.25; + section.dy = 0.60; + section.x = 0.80; + section.y = 0.00; + section.name = psStringCopy ("MagMxx"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.xmax = Xmax; + graphdata.ymin = -7.9; + graphdata.ymax = -17.1; + KapaSetLimits (myKapa, &graphdata); + + strcpy (graphdata.labels, "2001"); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_YP); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, xPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, xEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, xDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, xCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, mCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, xSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "y"); + + // fourth section: MagSigma + section.dx = 0.75; + section.dy = 0.15; + section.x = 0.00; + section.y = 0.65; + section.name = psStringCopy ("MagSigma"); + KapaSetSection (myKapa, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmax = -7.9; + graphdata.xmin = -17.1; + graphdata.ymin = -20.1; + graphdata.ymax = +20.1; + KapaSetLimits (myKapa, &graphdata); + + strcpy (graphdata.labels, "0100"); + KapaBox (myKapa, &graphdata); + // KapaSendLabel (myKapa, "inst mag", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "EXT&ss&c", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nPSF, &graphdata); + KapaPlotVector (myKapa, nPSF, mPSF->data.F32, "x"); + KapaPlotVector (myKapa, nPSF, sPSF->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nEXT, &graphdata); + KapaPlotVector (myKapa, nEXT, mEXT->data.F32, "x"); + KapaPlotVector (myKapa, nEXT, sEXT->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (myKapa, nDEF, &graphdata); + KapaPlotVector (myKapa, nDEF, mDEF->data.F32, "x"); + KapaPlotVector (myKapa, nDEF, sDEF->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nCR, &graphdata); + KapaPlotVector (myKapa, nCR, mCR->data.F32, "x"); + KapaPlotVector (myKapa, nCR, sCR->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 7; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (myKapa, nSAT, &graphdata); + KapaPlotVector (myKapa, nSAT, mSAT->data.F32, "x"); + KapaPlotVector (myKapa, nSAT, sSAT->data.F32, "y"); + + // draw N circles to outline the clumps + { + KapaSelectSection (myKapa, "MxxMyy"); + + // draw a circle centered on psfX,Y with size of the psf limit + psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); + psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); + + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 0; + + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + float Rx = psfdX * PSF_CLUMP_NSIGMA; + float Ry = psfdY * PSF_CLUMP_NSIGMA; + + for (int i = 0; i < xLimit->n; i++) { + xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; + yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; + } + KapaPrepPlot (myKapa, xLimit->n, &graphdata); + KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); + } + psFree (xLimit); + psFree (yLimit); + } + + psFree (xSAT); + psFree (ySAT); + psFree (mSAT); + + psFree (xEXT); + psFree (yEXT); + psFree (mEXT); + + psFree (xPSF); + psFree (yPSF); + psFree (mPSF); + + psFree (xDEF); + psFree (yDEF); + psFree (mDEF); + + psFree (xCR); + psFree (yCR); + psFree (mCR); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualShowResidualImage (pmReadout *readout) { + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (1); + if (myKapa == -1) return false; + + psphotVisualScaleImage (myKapa, readout->image, "resid", 1); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualPlotApResid (psArray *sources, float mean, float error) { + + Graphdata graphdata; + float lineX[2], lineY[2]; + + if (!pmVisualIsVisual()) return true; + + int myKapa = psphotKapaChannel (2); + if (myKapa == -1) return false; + + KapaClearPlots (myKapa); KapaInitGraph (&graphdata); - - // first section : mag vs CR nSigma - section.dx = 1.0; - section.dy = 0.5; - section.x = 0.0; - section.y = 0.0; - section.name = NULL; - psStringAppend (§ion.name, "a1"); - KapaSetSection (kapa3, §ion); - psFree (section.name); psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32); psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32); - - graphdata.xmin = +32.0; - graphdata.xmax = -32.0; - graphdata.ymin = +32.0; - graphdata.ymax = -32.0; - - // construct the plot vectors - int n = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source) continue; - if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!isfinite (source->crNsigma)) continue; - - x->data.F32[n] = -2.5*log10(source->peak->flux); - y->data.F32[n] = source->crNsigma; - graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); - graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); - graphdata.ymin = -0.5; - graphdata.ymax = 10.0; - - n++; - } - x->n = y->n = n; - - float range; - range = graphdata.xmax - graphdata.xmin; - graphdata.xmax += 0.05*range; - graphdata.xmin -= 0.05*range; - - // XXX set the plot range to match the image - KapaSetLimits (kapa3, &graphdata); - - KapaSetFont (kapa3, "helvetica", 14); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "Peak as Mag", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "CR N Sigma", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 2; - graphdata.size = 0.5; - graphdata.style = 2; - KapaPrepPlot (kapa3, n, &graphdata); - KapaPlotVector (kapa3, n, x->data.F32, "x"); - KapaPlotVector (kapa3, n, y->data.F32, "y"); - - // second section : mag vs EXT nSigma - section.dx = 1.0; - section.dy = 0.5; - section.x = 0.0; - section.y = 0.5; - section.name = NULL; - psStringAppend (§ion.name, "a2"); - KapaSetSection (kapa3, §ion); - psFree (section.name); - - graphdata.xmin = +32.0; - graphdata.xmax = -32.0; - graphdata.ymin = +32.0; - graphdata.ymax = -32.0; - - // construct the plot vectors - n = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source) continue; - if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!isfinite (source->extNsigma)) continue; - - x->data.F32[n] = -2.5*log10(source->peak->flux); - y->data.F32[n] = source->extNsigma; - graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); - graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); - graphdata.ymin = -0.5; - graphdata.ymax = 10.0; - - n++; - } - x->n = y->n = n; - - range = graphdata.xmax - graphdata.xmin; - graphdata.xmax += 0.05*range; - graphdata.xmin -= 0.05*range; - - // XXX set the plot range to match the image - KapaSetLimits (kapa3, &graphdata); - - KapaSetFont (kapa3, "helvetica", 14); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "EXT N Sigma", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 2; - graphdata.size = 0.5; - graphdata.style = 2; - KapaPrepPlot (kapa3, n, &graphdata); - KapaPlotVector (kapa3, n, x->data.F32, "x"); - KapaPlotVector (kapa3, n, y->data.F32, "y"); - - psFree (x); - psFree (y); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - return true; -} - -bool psphotVisualShowResidualImage (pmReadout *readout) { - - if (!pmVisualIsVisual()) return true; - - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - psphotVisualScaleImage (kapa, readout->image, "resid", 1); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - return true; -} - -bool psphotVisualPlotApResid (psArray *sources) { - - Graphdata graphdata; - - if (!pmVisualIsVisual()) return true; - - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - KapaClearPlots (kapa3); - KapaInitGraph (&graphdata); - - psVector *x = psVectorAllocEmpty (sources->n, PS_TYPE_F32); - psVector *y = psVectorAllocEmpty (sources->n, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty (sources->n, PS_TYPE_F32); graphdata.xmin = +32.0; @@ -1785,4 +1968,5 @@ x->data.F32[n] = source->psfMag; y->data.F32[n] = source->apMag - source->psfMag; + dy->data.F32[n] = source->errMag; graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); @@ -1792,5 +1976,5 @@ n++; } - x->n = y->n = n; + x->n = y->n = dy->n = n; float range; @@ -1802,11 +1986,16 @@ graphdata.ymin -= 0.05*range; - // XXX set the plot range to match the image - KapaSetLimits (kapa3, &graphdata); - - KapaSetFont (kapa3, "helvetica", 14); - KapaBox (kapa3, &graphdata); - KapaSendLabel (kapa3, "PSF Mag", KAPA_LABEL_XM); - KapaSendLabel (kapa3, "Ap Mag - PSF Mag", KAPA_LABEL_YM); + // XXX test + graphdata.xmin = -17.0; + graphdata.xmax = -9.0; + graphdata.ymin = -0.31; + graphdata.ymax = +0.31; + + KapaSetLimits (myKapa, &graphdata); + + KapaSetFont (myKapa, "helvetica", 14); + KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "PSF Mag", KAPA_LABEL_XM); + KapaSendLabel (myKapa, "Ap Mag - PSF Mag", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("black"); @@ -1814,18 +2003,105 @@ graphdata.size = 0.5; graphdata.style = 2; - KapaPrepPlot (kapa3, n, &graphdata); - KapaPlotVector (kapa3, n, x->data.F32, "x"); - KapaPlotVector (kapa3, n, y->data.F32, "y"); + graphdata.etype |= 0x01; + KapaPrepPlot (myKapa, n, &graphdata); + KapaPlotVector (myKapa, n, x->data.F32, "x"); + KapaPlotVector (myKapa, n, y->data.F32, "y"); + KapaPlotVector (myKapa, n, dy->data.F32, "dym"); + KapaPlotVector (myKapa, n, dy->data.F32, "dyp"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 0; + graphdata.etype = 0; + lineX[0] = graphdata.xmin; + lineX[1] = graphdata.xmax; + lineY[0] = lineY[1] = mean; + KapaPrepPlot (myKapa, 2, &graphdata); + KapaPlotVector (myKapa, 2, lineX, "x"); + KapaPlotVector (myKapa, 2, lineY, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 0; + graphdata.etype = 0; + lineX[0] = graphdata.xmin; + lineX[1] = graphdata.xmax; + lineY[0] = lineY[1] = mean + error; + KapaPrepPlot (myKapa, 2, &graphdata); + KapaPlotVector (myKapa, 2, lineX, "x"); + KapaPlotVector (myKapa, 2, lineY, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.5; + graphdata.style = 0; + graphdata.etype = 0; + lineX[0] = graphdata.xmin; + lineX[1] = graphdata.xmax; + lineY[0] = lineY[1] = mean - error; + KapaPrepPlot (myKapa, 2, &graphdata); + KapaPlotVector (myKapa, 2, lineX, "x"); + KapaPlotVector (myKapa, 2, lineY, "y"); psFree (x); psFree (y); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } + psFree (dy); + + pmVisualAskUser(NULL); + return true; +} + +bool psphotVisualShowPetrosians (psArray *sources) { + + int Noverlay, NOVERLAY; + KiiOverlay *overlay; + + if (!pmVisualIsVisual()) return true; + + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + Noverlay = 0; + NOVERLAY = 100; + ALLOCATE (overlay, KiiOverlay, NOVERLAY); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + if (!source) continue; + if (!source->extpars) continue; + if (!source->extpars->profile) continue; + if (!source->extpars->petrosian_80) continue; + + pmSourceRadialProfile *profile = source->extpars->profile; + pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 2.0*petrosian->radius; + overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major; + overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 4.0*petrosian->radius; + overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major; + overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + } + + KiiLoadOverlay (kapa, overlay, Noverlay, "red"); + FREE (overlay); + + pmVisualAskUser(NULL); return true; } @@ -1850,2 +2126,40 @@ # endif + +# if (0) +// *** make a histogram of the source counts in the x and y directions +psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0); +psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0); +psVectorHistogram (nX, xFaint, NULL, NULL, 0); +psVectorHistogram (nY, yFaint, NULL, NULL, 0); +psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); +psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); +psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); +psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); +for (int i = 0; i < nX->nums->n; i++) { + dX->data.F32[i] = nX->nums->data.S32[i]; + vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); +} +for (int i = 0; i < nY->nums->n; i++) { + dY->data.F32[i] = nY->nums->data.S32[i]; + vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); +} + +graphdata.color = KapaColorByName ("black"); +graphdata.ptype = 0; +graphdata.size = 0.0; +graphdata.style = 0; +KapaPrepPlot (myKapa, dX->n, &graphdata); +KapaPlotVector (myKapa, dX->n, dX->data.F32, "x"); +KapaPlotVector (myKapa, vX->n, vX->data.F32, "y"); + +psFree (nX); +psFree (dX); +psFree (vX); + +psFree (nY); +psFree (dY); +psFree (vY); + +# endif +