Index: branches/sc_branches/trunkTest/psphot/doc/notes.20100715.txt =================================================================== --- branches/sc_branches/trunkTest/psphot/doc/notes.20100715.txt (revision 29060) +++ branches/sc_branches/trunkTest/psphot/doc/notes.20100715.txt (revision 29060) @@ -0,0 +1,33 @@ + +2010.08.12 + + speeding up the PCM fitting. + + * first, I need to convert the convolutions to FFTs (double check that this will actually be faster!) + * second, I need to thread the PCM fitting process + + I have three options: + + 1) use the FFT threading in fftw -- each convolution would run threaded, but separate ones would be run in series + * con: the images are generally small, so the threading overhead may be large compared to the operation (test?) + 2) each pass of the minimization has 1 + N convolutions to perform: run these in separate threads + * can only use 1+N threads, where N is the number of free parameters (Io, Sxx, Syy, Sxy, but maybe also Xo, Yo) + 3) thread by source, using the system defined for the other fitting. + * since these source are large, this is a bit trickier to guarantee non-collision. + + infrastructure is easy -- copy psphotBlendFit. + +2010.07.15 + +various psphot issues that need to be addressed: + +* update extNsigma based on kron mags +* convert EXT sigmas to a probability (watch that extNsigma is correctly defined) +* psphot forced fails for negative flux + - do a run, plot a histogram of fluxes as various steps + FIXED: this was caused by minimum limit of 0.0 in linear fitting + +* capability to choose extended source fits based on |b| + +* update psf-convolved fitting to work with sersic (iteration) + DONE: but, needs to be sped up + Index: branches/sc_branches/trunkTest/psphot/doc/psf.load.txt =================================================================== --- branches/sc_branches/trunkTest/psphot/doc/psf.load.txt (revision 29060) +++ branches/sc_branches/trunkTest/psphot/doc/psf.load.txt (revision 29060) @@ -0,0 +1,75 @@ + +We need to have consistent behavior when we determine the PSF from the image or load it externally. + +Here is a list of the metadata info generated during the psphot analysis (these are mostly stored on the 'analysis' metadata, sometimes in chip and sometimes in readout). + + ANGLE : PSF Model position angle + APLOSS : flux lost out of reference aperture + APMIFIT : aperture magnitude - psf fit magnitude (aperture correction + DAPMIFIT : error on APMIFIT + *EFFECTIVE_AREA : effective area of the PSF model + FLIMIT : limiting flux (not currently set correctly) + FSATUR : saturation flux (not currently set correctly) + FWHM_MAJ : PSF Model FWHM (Major-Axis) + FWHM_MIN : PSF Model FWHM (Minor-Axis) + FW_MJ_LQ : PSF Model FWHM (Major-Axis lower quartile) + FW_MJ_SG : PSF Model FWHM (Major-Axis sigma) + FW_MJ_UQ : PSF Model FWHM (Major-Axis upper quartile) + FW_MN_LQ : PSF Model FWHM (Minor-Axis lower quartile) + FW_MN_SG : PSF Model FWHM (Minor-Axis sigma) + FW_MN_UQ : PSF Model FWHM (Minor-Axis upper quartile) + IQ_FW1 : Moments-based FWHM (Major-Axis) + IQ_FW1_E : Moments-based FWHM (Major-Axis Stdev) + IQ_FW2 : Moments-based FWHM (Minor-Axis) + IQ_FW2_E : Moments-based FWHM (Minor-Axis Stdev) + IQ_M2 : Moments-based shape parameter + IQ_M2C : Moments-based shape parameter + IQ_M2C_E : Moments-based shape parameter + IQ_M2C_L : Moments-based shape parameter + IQ_M2C_U : Moments-based shape parameter + IQ_M2S : Moments-based shape parameter + IQ_M2S_E : Moments-based shape parameter + IQ_M2S_L : Moments-based shape parameter + IQ_M2S_U : Moments-based shape parameter + IQ_M2_ER : Moments-based shape parameter + IQ_M2_LQ : Moments-based shape parameter + IQ_M2_UQ : Moments-based shape parameter + IQ_M3 : Moments-based shape parameter + IQ_M3_ER : Moments-based shape parameter + IQ_M3_LQ : Moments-based shape parameter + IQ_M3_UQ : Moments-based shape parameter + IQ_M4 : Moments-based shape parameter + IQ_M4_ER : Moments-based shape parameter + IQ_M4_LQ : Moments-based shape parameter + IQ_M4_UQ : Moments-based shape parameter + IQ_NSTAR : Number of sources used to measure moments-based parameters + MOMENTS_GAUSS_SIGMA : Sigma of Gaussian window used for moments analysis + MOMENTS_SX_MAX : max allowed Mxx (second moment in x) allowed for PSF clump + MOMENTS_SY_MAX : max allowed Myy (second moment in y) allowed for PSF clump + MSKY_MAX : Sky Model max value + MSKY_MIN : Sky Model min value + MSKY_MN : Sky Model mean value + MSKY_NX : Sky Model number of cells (x) + MSKY_NY : Sky Model number of cells (y) + MSKY_SIG : Sky Model sigma + NAPMIFIT : Number of stars used for APMIFIT analysis + NDET_CR : Number of sources classified as 'cosmic ray' + NDET_EXT : Number of sources classified as 'extended' + NPSFSTAR : Number of sources used for PSF model + NSTARS : Total number of sources + PEAK_THRESHOLD : threshold of peak detection analysis (saved but not re-used) + PSF.CLUMP.DX : size of clump in X + PSF.CLUMP.DY : size of clump in Y + PSF.CLUMP.NREGIONS : number of PSF clump regions measured for this image + PSF.CLUMP.X : center of clump in X + PSF.CLUMP.Y : center of clump in Y + PSFMODEL : do we have a valid PSF model? + PSF_APERTURE : aperture used to measure the aperture magnitude of PSF sources + PSF_CLUMP_GRID_SCALE: pixelization size used for PSF clump search + PSF_FIT_RADIUS : size of window used for PSF model fits + PSF_MAX_NSTARS : maximum number of sources used to find the PSF model + PSF_MOMENTS_RADIUS : size of window used for moments analysis (Gaussian window is applied within this radius) + SIGMA_SMOOTH : sigma of kernel used for smoothing to generate the significance + SIGNIFICANCE_SCALE_FACTOR : factor needed to adjust the level of the significance image + SKY_MEAN : mean of background + SKY_STDEV : stdev of full background image Index: branches/sc_branches/trunkTest/psphot/doc/visual.txt =================================================================== --- branches/sc_branches/trunkTest/psphot/doc/visual.txt (revision 29060) +++ branches/sc_branches/trunkTest/psphot/doc/visual.txt (revision 29060) @@ -0,0 +1,47 @@ + +Visualization Options + +image : display image +backgnd : display image +signif : display image + +show image : image 1 or object 1 +show background : image 2 or backgnd 1 +show significance : image 3 or signif 1 +show residual : image 2 or resid 1 + +XXX: these all need to have the image: +show peaks : object 1 or peaks 1 +show moments : object 2 or image.moments 1 +show footprints : object 3 or footprints 1 +show rough class : object 3 or class 3 +show flags : (disabled) +show size class : object 2 or class 2 +show petrosian : object 2 or petro 1 + +plot final moments (2D) : moments 1 +plot test moments (2D) : moments 2 +plot source sizes : class 1 +plot ap resid : apresid 1 +plot chisq : chisq 1 + +show PSF model : psf 1 +show PSF stars : psf 2 +show Sat stars : ? + +show radial profiles : profile + +petrosian visualization: + + psphotPetrosianVisualEllipticalContour + psphotPetrosianVisualStats + psphotPetrosianVisualProfileRadii + psphotPetrosianVisualProfileByAngle + +psphot visualization in psModules: + + pmSourceVisualPSFModelResid + pmSourceVisualShowModelFit + pmSourceVisualShowModelFits + pmSourceVisualPlotPSFMetricSubpix + pmSourceVisualPlotPSFMetric Index: branches/sc_branches/trunkTest/psphot/src/Makefile.am =================================================================== --- branches/sc_branches/trunkTest/psphot/src/Makefile.am (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/Makefile.am (revision 29060) @@ -163,10 +163,8 @@ psphotOutput.c \ psphotFakeSources.c \ - psphotModelWithPSF.c \ psphotExtendedSourceAnalysis.c \ psphotExtendedSourceAnalysisByObject.c \ psphotExtendedSourceFits.c \ psphotKernelFromPSF.c \ - psphotPSFConvModel.c \ psphotFitSet.c \ psphotSourceFreePixels.c \ Index: branches/sc_branches/trunkTest/psphot/src/psphot.h =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphot.h (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphot.h (revision 29060) @@ -12,12 +12,4 @@ #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs - -// pmPCMData : PSF Convolved Model data storage structure -typedef struct { - psImage *model; - psArray *dmodels; - psImage *modelConv; - psArray *dmodelsConv; -} pmPCMData; // top-level psphot functions @@ -114,4 +106,5 @@ bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotExtendedSourceFits_Threaded (psThreadJob *job); bool psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -188,11 +181,11 @@ // functions to set the correct source pixels -bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type); +bool psphotInitRadiusPSF (psMetadata *recipe, pmReadout *readout); +bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout); bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR); -bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type); -bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); -float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal); +bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor); +bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep); bool psphotDumpMoments (psMetadata *recipe, psArray *sources); @@ -210,10 +203,10 @@ // functions to support the source fitting process bool psphotInitLimitsPSF (psMetadata *recipe, pmReadout *readout); -bool psphotInitLimitsEXT (psMetadata *recipe); -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal); -psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); +bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout); +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal); +psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); // functions to support simultaneous multi-source fitting @@ -252,4 +245,5 @@ bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout); bool psphotVisualShowSignificance (psImage *image, float min, float max); +bool psphotVisualShowLogSignificance (psImage *image, float min, float max); bool psphotVisualShowPeaks (pmDetections *detections); bool psphotVisualShowFootprints (pmDetections *detections); @@ -292,33 +286,4 @@ bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise); -// structures & functions to support psf-convolved model fitting - -// psf-convolved model fitting -bool psphotModelWithPSF_LMM ( - psMinimization *min, - psImage *covar, - psVector *params, - psMinConstraint *constraint, - pmSource *source, - const psKernel *psf, - psMinimizeLMChi2Func func); - -psF32 psphotModelWithPSF_SetABX( - psImage *alpha, - psVector *beta, - const psVector *params, - const psVector *paramMask, - pmPCMData *pcm, - const pmSource *source, - const psKernel *psf, - psMinimizeLMChi2Func func); - -pmPCMData *pmPCMDataAlloc ( - const psVector *params, - const psVector *paramMask, - pmSource *source); - -psImage *pmPCMDataSaveImage (pmPCMData *pcm); - int psphotKapaOpen (void); bool psphotKapaClose (void); @@ -364,6 +329,4 @@ bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final); -int pmPhotObjSortBySN (const void **a, const void **b); -int pmPhotObjSortByX (const void **a, const void **b); typedef enum { @@ -462,3 +425,8 @@ bool psphotStackObjectsUnifyPosition (psArray *objects); +bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal); + +bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); +pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize); + #endif Index: branches/sc_branches/trunkTest/psphot/src/psphotAddNoise.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotAddNoise.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotAddNoise.c (revision 29060) @@ -34,7 +34,4 @@ bool status = false; - psEllipseShape oldshape; - psEllipseShape newshape; - psEllipseAxes axes; // find the currently selected readout @@ -86,44 +83,5 @@ if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; // model must be defined - - if (add) { - psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - } else { - psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - } - - psF32 *PAR = model->params->data.F32; - - // save original values - float oldI0 = PAR[PM_PAR_I0]; - oldshape.sx = PAR[PM_PAR_SXX]; - oldshape.sy = PAR[PM_PAR_SYY]; - oldshape.sxy = PAR[PM_PAR_SXY]; - - // XXX can this be done more intelligently? - if (oldI0 == 0.0) continue; - if (!isfinite(oldI0)) continue; - - // increase size and height of source - axes = psEllipseShapeToAxes (oldshape, 20.0); - axes.major *= SIZE; - axes.minor *= SIZE; - newshape = psEllipseAxesToShape (axes); - PAR[PM_PAR_I0] = FACTOR*oldI0; - PAR[PM_PAR_SXX] = newshape.sx; - PAR[PM_PAR_SYY] = newshape.sy; - PAR[PM_PAR_SXY] = newshape.sxy; - - // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected - pmSourceOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, add, maskVal, 0, 0); - - // restore original values - PAR[PM_PAR_I0] = oldI0; - PAR[PM_PAR_SXX] = oldshape.sx; - PAR[PM_PAR_SYY] = oldshape.sy; - PAR[PM_PAR_SXY] = oldshape.sxy; + pmSourceNoiseOp (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, FACTOR, SIZE, add, maskVal, 0, 0); } if (add) { @@ -132,4 +90,5 @@ psLogMsg ("psphot.noise", PS_LOG_INFO, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); } + return true; } Index: branches/sc_branches/trunkTest/psphot/src/psphotApResid.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotApResid.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotApResid.c (revision 29060) @@ -459,5 +459,5 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - bool status = pmSourceMagnitudes (source, psf, photMode, maskVal); + bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit Index: branches/sc_branches/trunkTest/psphot/src/psphotArguments.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotArguments.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotArguments.c (revision 29060) @@ -176,10 +176,4 @@ } - // visual : interactive display mode - if ((N = psArgumentGet (argc, argv, "-visual"))) { - psArgumentRemove (N, &argc, argv); - pmVisualSetVisual(true); - } - // break : used from recipe throughout psphotReadout if ((N = psArgumentGet (argc, argv, "-break"))) { Index: branches/sc_branches/trunkTest/psphot/src/psphotBlendFit.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotBlendFit.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotBlendFit.c (revision 29060) @@ -65,19 +65,38 @@ assert (status && fitIter > 0); - float fitTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_TOL"); // Fit tolerance - assert (status && isfinite(fitTol) && fitTol > 0); + float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance + if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) { + fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance + if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) { + psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive"); + } + } + + float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance + if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) { + fitMaxTol = 1.0; + } bool poisson = psMetadataLookupBool(&status, recipe, "POISSON.ERRORS.PHOT.LMM"); // Poisson errors? assert (status); + float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_CHISQ"); // Fit tolerance + float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG"); assert (status && isfinite(skySig) && skySig > 0); // Define source fitting parameters for extended source fits - pmSourceFitModelInit(fitIter, fitTol, PS_SQR(skySig), poisson); + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->nIter = fitIter; + fitOptions->minTol = fitMinTol; + fitOptions->maxTol = fitMaxTol; + fitOptions->maxChisqDOF = maxChisqDOF; + fitOptions->poissonErrors = poisson; + fitOptions->weight = PS_SQR(skySig); + fitOptions->mode = PM_SOURCE_FIT_PSF; psphotInitLimitsPSF (recipe, readout); - psphotInitLimitsEXT (recipe); - psphotInitRadiusPSF (recipe, readout->analysis, psf->type); + psphotInitLimitsEXT (recipe, readout); + psphotInitRadiusPSF (recipe, readout); // starts the timer, sets up the array of fitSets @@ -88,4 +107,5 @@ if (!sources->n) { psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + psFree (fitOptions); return true; } @@ -112,4 +132,5 @@ psArrayAdd(job->args, 1, psf); psArrayAdd(job->args, 1, newSources); // return for new sources + psArrayAdd(job->args, 1, fitOptions); // default fit options psFree (newSources); @@ -121,31 +142,7 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (fitOptions); return NULL; } - -# if (0) - { - int nfit = 0; - int npsf = 0; - int next = 0; - int nfail = 0; - psArray *newSources = psArrayAllocEmpty(16); - - if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - Nfit += nfit; - Npsf += npsf; - Next += next; - Nfail += nfail; - - // add these back onto sources - for (int k = 0; k < newSources->n; k++) { - psArrayAdd (sources, 16, newSources->data[k]); - } - psFree (newSources); - } -# endif } @@ -153,4 +150,5 @@ if (!psThreadPoolWait (false)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (fitOptions); return NULL; } @@ -163,11 +161,11 @@ } else { psScalar *scalar = NULL; - scalar = job->args->data[5]; + scalar = job->args->data[6]; Nfit += scalar->data.S32; - scalar = job->args->data[6]; + scalar = job->args->data[7]; Npsf += scalar->data.S32; - scalar = job->args->data[7]; + scalar = job->args->data[8]; Next += scalar->data.S32; - scalar = job->args->data[8]; + scalar = job->args->data[9]; Nfail += scalar->data.S32; @@ -186,4 +184,5 @@ psphotSaveImage (NULL, readout->image, "image.v2.fits"); } + psFree (fitOptions); psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit); @@ -204,9 +203,10 @@ psScalar *scalar = NULL; - pmReadout *readout = job->args->data[0]; - psMetadata *recipe = job->args->data[1]; - psArray *sources = job->args->data[2]; - pmPSF *psf = job->args->data[3]; - psArray *newSources = job->args->data[4]; + pmReadout *readout = job->args->data[0]; + psMetadata *recipe = job->args->data[1]; + psArray *sources = job->args->data[2]; + pmPSF *psf = job->args->data[3]; + psArray *newSources = job->args->data[4]; + pmSourceFitOptions *fitOptions = job->args->data[5]; // bit-masks to test for good/bad pixels @@ -269,8 +269,14 @@ Nfit ++; + if (0) { + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + fprintf (stderr, "1: Mxx: %f, Myy: %f\n", Mxx, Myy); + } + // try fitting PSFs or extended sources depending on source->mode // these functions subtract the resulting fitted source if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { + if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) { source->type = PM_SOURCE_TYPE_EXTENDED; psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); @@ -280,5 +286,5 @@ } } else { - if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { + if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) { source->type = PM_SOURCE_TYPE_STAR; psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); @@ -289,4 +295,10 @@ } + if (0) { + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + fprintf (stderr, "2: Mxx: %f, Myy: %f\n", Mxx, Myy); + } + psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); Nfail ++; @@ -298,14 +310,14 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[5]; + scalar = job->args->data[6]; scalar->data.S32 = Nfit; - scalar = job->args->data[6]; + scalar = job->args->data[7]; scalar->data.S32 = Npsf; - scalar = job->args->data[7]; + scalar = job->args->data[8]; scalar->data.S32 = Next; - scalar = job->args->data[8]; + scalar = job->args->data[9]; scalar->data.S32 = Nfail; Index: branches/sc_branches/trunkTest/psphot/src/psphotChoosePSF.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotChoosePSF.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotChoosePSF.c (revision 29060) @@ -74,5 +74,5 @@ // structure to store user options defining the psf - pmPSFOptions *options = pmPSFOptionsAlloc (); + pmPSFOptions *options = pmPSFOptionsAlloc(); // load user options from the recipe. no need to check existence -- they are @@ -138,11 +138,28 @@ return false; } - float fitTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance - if (!status || !isfinite(fitTol) || fitTol <= 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_TOL is not positive"); - return false; - } - pmSourceFitModelInit(fitIter, fitTol, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM); - + float fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MIN_TOL"); // Fit tolerance + if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) { + fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance + if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive"); + return false; + } + } + float fitMaxTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_TOL"); // Fit tolerance + if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) { + fitMaxTol = 1.0; + } + float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance + + // options which modify the behavior of the model fitting + options->fitOptions = pmSourceFitOptionsAlloc(); + options->fitOptions->nIter = fitIter; + options->fitOptions->minTol = fitMinTol; + options->fitOptions->maxTol = fitMaxTol; + options->fitOptions->maxChisqDOF = maxChisqDOF; + options->fitOptions->poissonErrors = options->poissonErrorsPhotLMM; + options->fitOptions->weight = PS_SQR(SKY_SIG); + options->fitOptions->mode = PM_SOURCE_FIT_PSF; + psArray *stars = psArrayAllocEmpty (sources->n); @@ -227,4 +244,5 @@ // try each model option listed in config + // pmPSFtryModel makes a local copy of the sources -- those points are not the same as those for 'sources' for (int i = 0; i < modelNames->n; i++) { char *modelName = modelNames->data[i]; @@ -304,10 +322,42 @@ // unset the PSFSTAR flag for stars not used for PSF model + // XXX a more efficient way of achieving this would be to record a pair of arrays + // of the source index and the source id for the psf stars. but that would require we do + // not re-sort the source list in the meanwhile + int nDrop = 0; for (int i = 0; i < try->sources->n; i++) { pmSource *source = try->sources->data[i]; if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; - } - } + // need to find this source in the original list (these are copies, not pointers) + for (int j = 0; j < sources->n; j++) { + pmSource *realSource = sources->data[j]; + if (realSource->id != source->id) continue; + realSource->mode &= ~PM_SOURCE_MODE_PSFSTAR; + source->mode &= ~PM_SOURCE_MODE_PSFSTAR; + nDrop ++; + break; + } + } + } + // fprintf (stderr, "drop %d stars as PSF stars\n", nDrop); + + // XXX is this working? + // int N1 = 0; + // for (int i = 0; i < try->sources->n; i++) { + // pmSource *source = try->sources->data[i]; + // fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR)); + // if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + // N1 ++; + // } + // } + // int N2 = 0; + // for (int i = 0; i < sources->n; i++) { + // pmSource *source = sources->data[i]; + // fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR)); + // if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + // N2 ++; + // } + // } + // fprintf (stderr, "N1: %d, N2: %d\n", N1, N2); // build a PSF residual image Index: branches/sc_branches/trunkTest/psphot/src/psphotEllipticalContour.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotEllipticalContour.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotEllipticalContour.c (revision 29060) @@ -82,5 +82,5 @@ params->data.F32[PAR_RMIN] = Rmin; - psMinimization *myMin = psMinimizationAlloc (25, 0.001); + psMinimization *myMin = psMinimizationAlloc (25, 0.01, 1.00); psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); Index: branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceAnalysis.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceAnalysis.c (revision 29060) @@ -121,4 +121,6 @@ if (source->peak->y > AnalysisRegion.y1) continue; + // fprintf (stderr, "xsrc: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); + // replace object in image if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { @@ -157,27 +159,4 @@ } -# if (0) - // Isophotal Mags - if (doIsophotal) { - if (!psphotIsophotal (source, recipe, maskVal)) { - psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - } else { - psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Nisophot ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } - } - // Kron Mags - if (doKron) { - if (!psphotKron (source, recipe, maskVal)) { - psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - } else { - psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Nkron ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } - } -# endif - // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); @@ -201,4 +180,6 @@ } + // fprintf (stderr, "xsrc : tried %ld objects\n", sources->n); + return true; } Index: branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceFits.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceFits.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotExtendedSourceFits.c (revision 29060) @@ -31,4 +31,186 @@ // non-linear model fitting for extended sources bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + int Next = 0; + int Nconvolve = 0; + int NconvolvePass = 0; + int Nplain = 0; + int NplainPass = 0; + + psTimerStart ("psphot.extended"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // 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); + + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // select the collection of desired models + psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); + if (!status) { + psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); + return true; + } + if (models->list->n == 0) { + psWarning ("extended source model fits requested but no models are specified\n"); + return true; + } + + // validate the model entries + psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + + if (item->type != PS_DATA_METADATA) { + psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); + // XXX we could cull the bad entries or build a validated model folder + } + + psMetadata *model = (psMetadata *) item->data.md; + + // check on the model type + char *modelName = psMetadataLookupStr (&status, model, "MODEL"); + int modelType = pmModelClassGetType (modelName); + if (modelType < 0) { + psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); + } + psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); + + // check on the SNLIM, set a float value + char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); + if (!status) { + psAbort("SNLIM not defined for extended source model %s\n", item->name); + } + float SNlim = atof (SNword); + psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); + + // check on the PSF-Convolution status + char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); + if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { + psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); + } + bool convolved = !strcasecmp (convolvedWord, "true"); + psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); + } + psFree (iter); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0); + *AnalysisRegion = psRegionForImage(readout->image, psRegionFromString (region)); + if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined"); + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_EXTENDED_FIT"); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, models); + psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer + + PS_ARRAY_ADD_SCALAR(job->args, psfSize, 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 Next + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nconvolve + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NconvolvePass + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[7]; + Next += scalar->data.S32; + scalar = job->args->data[8]; + Nconvolve += scalar->data.S32; + scalar = job->args->data[9]; + NconvolvePass += scalar->data.S32; + scalar = job->args->data[10]; + Nplain += scalar->data.S32; + scalar = job->args->data[11]; + NplainPass += scalar->data.S32; + } + psFree(job); + } + } + psFree (cellGroups); + psFree(AnalysisRegion); + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d convolved models (%d passed)\n", Nconvolve, NconvolvePass); + psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); + return true; +} + +// non-linear model fitting for extended sources +bool psphotExtendedSourceFits_Threaded (psThreadJob *job) { bool status; @@ -39,93 +221,24 @@ int NplainPass = 0; bool savePics = false; - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); - - if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); - return true; - } - - // 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); - - psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels - assert (markVal); - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // select the collection of desired models - psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); - if (!status) { - psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); - return true; - } - if (models->list->n == 0) { - psWarning ("extended source model fits requested but no models are specified\n"); - return true; - } - - // validate the model entries - psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); - psMetadataItem *item = NULL; - while ((item = psMetadataGetAndIncrement (iter)) != NULL) { - - if (item->type != PS_DATA_METADATA) { - psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); - // XXX we could cull the bad entries or build a validated model folder - } - - psMetadata *model = (psMetadata *) item->data.md; - - // check on the model type - char *modelName = psMetadataLookupStr (&status, model, "MODEL"); - int modelType = pmModelClassGetType (modelName); - if (modelType < 0) { - psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); - } - psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); - - // check on the SNLIM, set a float value - char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); - if (!status) { - psAbort("SNLIM not defined for extended source model %s\n", item->name); - } - float SNlim = atof (SNword); - psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); - - // check on the PSF-Convolution status - char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); - if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { - psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); - } - bool convolved = !strcasecmp (convolvedWord, "true"); - psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); - } - psFree (iter); - - // 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"); - - // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) - int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); - assert (status); - - // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + float radius; + psScalar *scalar = NULL; + + // arguments: readout, sources, models, region, psfSize, maskVal, markVal + pmReadout *readout = job->args->data[0]; + psArray *sources = job->args->data[1]; + psMetadata *models = job->args->data[2]; + psRegion *region = job->args->data[3]; + int psfSize = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + + // Define source fitting parameters for extended source fits + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT; + // XXX for now, use the defaults for the rest: + // fitOptions->nIter = fitIter; + // fitOptions->tol = fitTol; + // fitOptions->poissonErrors = poisson; + // fitOptions->weight = PS_SQR(skySig); // choose the sources of interest @@ -140,8 +253,8 @@ // XXX this should use peak? - 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; + if (source->peak->x < region->x0) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->x > region->x1) continue; + if (source->peak->y > region->y1) continue; // if model is NULL, we don't have a starting guess @@ -155,7 +268,36 @@ Next ++; + // set the radius based on the footprint (also sets the mask pixels) + if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) { + fprintf (stderr, "skipping (1) %f, %f\n", source->peak->xf, source->peak->yf); + psFree (fitOptions) + return false; + } + + // 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 + // this uses the footprint to judge both radius and aperture? + // XXX save the psf-based moments for output + if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) { + fprintf (stderr, "skipping (2) %f, %f\n", source->peak->xf, source->peak->yf); + // subtract the best fit from the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + // XXX raise an error of some kind + continue; + } + // save the modelFlux here in case we need to subtract it (for failure) psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); - + if (!modelFluxStart) { + pmSourceCacheModel (source, maskVal); + modelFluxStart = psMemIncrRefCounter (source->modelFlux); + if (!modelFluxStart) { + fprintf (stderr, "skipping (3) %f, %f\n", source->peak->xf, source->peak->yf); + // XXX raise an error of some kind? + continue; + } + } + if (savePics) { psphotSaveImage (NULL, readout->image, "image.xp.fits"); @@ -183,4 +325,6 @@ assert (status); + // fprintf (stderr, "xfit: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); + // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? @@ -201,5 +345,5 @@ pmModel *modelFit = NULL; if (convolved) { - modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize); + modelFit = psphotFitPCM (readout, source, fitOptions, modelType, maskVal, markVal, psfSize); if (!modelFit) { psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My); @@ -215,5 +359,5 @@ psFree (source->modelFlux); source->modelFlux = NULL; - modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal); + modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal); if (!modelFit) { psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); @@ -233,4 +377,5 @@ // test for fit quality / result + modelFit->fitRadius = radius; psArrayAdd (source->modelFits, 4, modelFit); @@ -311,8 +456,24 @@ } } - - psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); - psLogMsg ("psphot", PS_LOG_INFO, " %d convolved models (%d passed)\n", Nconvolve, NconvolvePass); - psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); + psFree (fitOptions); + + // fprintf (stderr, "xfit : tried %ld objects\n", sources->n); + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[7]; + scalar->data.S32 = Next; + + scalar = job->args->data[8]; + scalar->data.S32 = Nconvolve; + + scalar = job->args->data[9]; + scalar->data.S32 = NconvolvePass; + + scalar = job->args->data[10]; + scalar->data.S32 = Nplain; + + scalar = job->args->data[11]; + scalar->data.S32 = NplainPass; + return true; } Index: branches/sc_branches/trunkTest/psphot/src/psphotFindDetections.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotFindDetections.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotFindDetections.c (revision 29060) @@ -90,10 +90,6 @@ psphotVisualShowSignificance (significance, -1.0, PS_SQR(3.0*NSIGMA_PEAK)); - // XXX getting some strange results from significance image - if (0) { - psImage *lsig = (psImage *) psUnaryOp (NULL, significance, "log"); - psphotVisualShowSignificance (lsig, 0.0, 4.0); - psFree (lsig); - } + // display the log significance image + psphotVisualShowLogSignificance (significance, 0.0, 4.5); // detect the peaks in the significance image Index: branches/sc_branches/trunkTest/psphot/src/psphotFitSet.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotFitSet.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotFitSet.c (revision 29060) @@ -24,6 +24,15 @@ } + // Define source fitting parameters for extended source fits + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT; + // XXX for now, use the defaults for the rest: + // fitOptions->nIter = fitIter; + // fitOptions->tol = fitTol; + // fitOptions->poissonErrors = poisson; + // fitOptions->weight = PS_SQR(skySig); + // XXX pmSourceFitSet must cache the modelFlux? - pmSourceFitSet (source, modelSet, mode, maskVal); + pmSourceFitSet (source, modelSet, fitOptions, maskVal); // write out positive object Index: branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinear.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinear.c (revision 29060) @@ -100,4 +100,13 @@ if (!status) { SKY_FIT_LINEAR = false; + } + + float MIN_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MIN_VALID_FLUX"); + if (!status) { + MIN_VALID_FLUX = 1e-8; + } + float MAX_VALID_FLUX = psMetadataLookupF32(&status, recipe, "PSF_FIT_MAX_VALID_FLUX"); + if (!status) { + MAX_VALID_FLUX = 1e+8; } @@ -228,7 +237,7 @@ psSparseConstraint constraint; - constraint.paramMin = 0.0; - constraint.paramMax = 1e8; - constraint.paramDelta = 1e8; + constraint.paramMin = MIN_VALID_FLUX; + constraint.paramMax = MAX_VALID_FLUX; + constraint.paramDelta = 1e7; // solve for normalization terms (need include local sky?) Index: branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinearStack.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotFitSourcesLinearStack.c (revision 29060) @@ -170,17 +170,2 @@ return true; } - -// sort by X (ascending) -int pmPhotObjSortByX (const void **a, const void **b) -{ - pmPhotObj *objA = *(pmPhotObj **)a; - pmPhotObj *objB = *(pmPhotObj **)b; - - psF32 fA = objA->x; - psF32 fB = objB->x; - - psF32 diff = fA - fB; - if (diff > FLT_EPSILON) return (+1); - if (diff < FLT_EPSILON) return (-1); - return (0); -} Index: branches/sc_branches/trunkTest/psphot/src/psphotGuessModels.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotGuessModels.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotGuessModels.c (revision 29060) @@ -80,5 +80,5 @@ // setup the PSF fit radius details - psphotInitRadiusPSF (recipe, readout->analysis, psf->type); + psphotInitRadiusPSF (recipe, readout); // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) Index: branches/sc_branches/trunkTest/psphot/src/psphotLoadSRCTEXT.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotLoadSRCTEXT.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotLoadSRCTEXT.c (revision 29060) @@ -84,5 +84,6 @@ source->peak->yf = PAR[PM_PAR_YPOS]; // but we know the pixel coordinate - source->pixWeight = 1.0; + source->pixWeightNotBad = 1.0; + source->pixWeightNotPoor = 1.0; source->crNsigma = 0.0; source->extNsigma = 0.0; Index: branches/sc_branches/trunkTest/psphot/src/psphotMagnitudes.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotMagnitudes.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotMagnitudes.c (revision 29060) @@ -124,13 +124,4 @@ return false; } - -# if (0) - int nap = 0; - if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; - } - Nap += nap; -# endif } @@ -186,5 +177,5 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); if (status && isfinite(source->apMag)) Nap ++; @@ -268,4 +259,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); if (!psThreadJobAddPending(job)) { @@ -304,4 +296,5 @@ psArray *sources = job->args->data[0]; psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[1],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { @@ -312,12 +305,14 @@ if (model == NULL) { psTrace ("psphot", 3, "fail mag : no valid model"); - source->pixWeight = NAN; + source->pixWeightNotBad = NAN; + source->pixWeightNotPoor = NAN; continue; } - status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal); + status = pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal); if (!status) { psTrace ("psphot", 3, "fail to measure pixel weight"); - source->pixWeight = NAN; + source->pixWeightNotBad = NAN; + source->pixWeightNotPoor = NAN; continue; } Index: branches/sc_branches/trunkTest/psphot/src/psphotModelWithPSF.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotModelWithPSF.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotModelWithPSF.c (revision 29060) @@ -73,7 +73,7 @@ // iterate until the tolerance is reached, or give up - while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) { + while ((min->iter < min->maxIter) && ((min->lastDelta > min->minTol) || !isfinite(min->lastDelta))) { psTrace("psphot", 5, "Iteration number %d. (max iterations is %d).\n", min->iter, min->maxIter); - psTrace("psphot", 5, "Last delta is %f. Min->tol is %f.\n", min->lastDelta, min->tol); + psTrace("psphot", 5, "Last delta is %f. Min->minTol is %f.\n", min->lastDelta, min->minTol); @@ -166,11 +166,11 @@ psFree(pcm); - if (min->iter == min->maxIter) { - psTrace("psphot", 3, "---- end (false) ----\n"); - return(false); - } - - psTrace("psphot", 3, "---- end (true) ----\n"); - return(true); + // if the last improvement was at least as good as maxTol, accept the fit: + if (min->lastDelta <= min->maxTol) { + psTrace("psphot", 6, "---- end (true) ----\n"); + return(true); + } + psTrace("psphot", 6, "---- end (false) ----\n"); + return(false); } Index: branches/sc_branches/trunkTest/psphot/src/psphotPSFConvModel.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotPSFConvModel.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotPSFConvModel.c (revision 29060) @@ -4,10 +4,14 @@ // save as static values so they may be set externally static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; -static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; +static psF32 PM_SOURCE_FIT_MODEL_MIN_TOL = 0.1; +static psF32 PM_SOURCE_FIT_MODEL_MAX_TOL = 2.0; // input source has both modelPSF and modelEXT. on successful exit, we set the // modelConv to contain the fitted parameters, and the modelFlux to contain the // convolved model image. -pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + +// XXX need to generalize this -- number of fitted parameters must be flexible based on the fitOptions + +pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { // maskVal is used to test for rejected pixels, and must include markVal @@ -90,5 +94,5 @@ // set up the minimization process - psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE); + psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_MIN_TOL, PM_SOURCE_FIT_MODEL_MAX_TOL); psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); Index: branches/sc_branches/trunkTest/psphot/src/psphotPetrosianVisual.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotPetrosianVisual.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotPetrosianVisual.c (revision 29060) @@ -1,4 +1,3 @@ # include "psphotInternal.h" -# define FORCE_VISUAL 0 // this function displays representative images as the psphot analysis progresses: @@ -54,6 +53,5 @@ Graphdata graphdata; - // return true; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.byangle", 2)) return true; if (kapa2 == -1) { @@ -101,6 +99,5 @@ Graphdata graphdata; - // return true; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.radii", 2)) return true; if (kapa == -1) { @@ -173,5 +170,5 @@ KapaSection section; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.stats", 2)) return true; if (kapa2 == -1) { @@ -311,5 +308,5 @@ Graphdata graphdata; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.ellipse", 2)) return true; if (kapa == -1) { Index: branches/sc_branches/trunkTest/psphot/src/psphotRadiusChecks.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotRadiusChecks.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotRadiusChecks.c (revision 29060) @@ -8,5 +8,5 @@ // and a per-object radius is calculated) -bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) { +bool psphotInitRadiusPSF(psMetadata *recipe, pmReadout *readout) { bool status = true; @@ -15,10 +15,10 @@ PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING"); - PSF_FIT_RADIUS = psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS"); + PSF_FIT_RADIUS = psMetadataLookupF32(&status, readout->analysis, "PSF_FIT_RADIUS"); if (!status) { PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); } - PSF_APERTURE = psMetadataLookupF32(&status, analysis, "PSF_APERTURE"); + PSF_APERTURE = psMetadataLookupF32(&status, readout->analysis, "PSF_APERTURE"); if (!status) { PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); @@ -28,5 +28,5 @@ if (PSF_FIT_RADIUS == 0.0) { - float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); if (!status) { gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); @@ -37,5 +37,5 @@ if (PSF_APERTURE == 0.0) { - float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + float gaussSigma = psMetadataLookupF32(&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); if (!status) { gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); @@ -122,9 +122,10 @@ } +static float EXT_FIT_SKY_SIG; static float EXT_FIT_NSIGMA; static float EXT_FIT_PADDING; static float EXT_FIT_MAX_RADIUS; -bool psphotInitRadiusEXT (psMetadata *recipe, pmModelType type) { +bool psphotInitRadiusEXT (psMetadata *recipe, pmReadout *readout) { bool status; @@ -134,9 +135,16 @@ EXT_FIT_MAX_RADIUS = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_RADIUS"); + float skyMean = psMetadataLookupF32 (&status, readout->analysis, "SKY_MEAN"); + float skyStdev = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV"); + + fprintf (stderr, "sky: %f +/- %f\n", skyMean, skyStdev); + + EXT_FIT_SKY_SIG = skyStdev; + return true; } // call this function whenever you (re)-define the EXT model -float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) { +bool psphotSetRadiusFootprint (float *radius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float factor) { psAssert (source, "source not defined??"); @@ -146,11 +154,11 @@ // set the radius based on the footprint: - if (!peak->footprint) goto escape; + if (!peak->footprint) return false; pmFootprint *footprint = peak->footprint; - if (!footprint->spans) goto escape; - if (footprint->spans->n < 1) goto escape; + if (!footprint->spans) return false; + if (footprint->spans->n < 1) return false; // find the max radius - float radius = 0.0; + float rawRadius = 0.0; for (int j = 0; j < footprint->spans->n; j++) { pmSpan *span = footprint->spans->data[j]; @@ -160,30 +168,22 @@ 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); + rawRadius = PS_MAX (rawRadius, hypot(dY, dX0)); + rawRadius = PS_MAX (rawRadius, hypot(dY, dX1)); + } + if (isnan(rawRadius)) return false; + rawRadius = PS_MIN (factor*rawRadius + EXT_FIT_PADDING, 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; + pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, rawRadius); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, peak->xf, peak->yf, rawRadius, "OR", markVal); + + *radius = rawRadius; + return true; } // 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"); +bool psphotSetRadiusModel (pmModel *model, pmReadout *readout, pmSource *source, psImageMaskType markVal, bool deep) { psF32 *PAR = model->params->data.F32; @@ -193,14 +193,17 @@ // set the fit radius based on the object flux limit and the model - float rawRadius = model->modelRadius (model->params, EXT_FIT_NSIGMA*moments->dSky); - - model->fitRadius = rawRadius + EXT_FIT_PADDING; - if (isnan(model->fitRadius)) psAbort("error in radius"); + float flux = deep ? EXT_FIT_NSIGMA*EXT_FIT_SKY_SIG : 0.1 * model->params->data.F32[PM_PAR_I0]; + + float rawRadius = model->modelRadius (model->params, flux); + if (isnan(rawRadius)) return false; + + rawRadius = PS_MIN (rawRadius + EXT_FIT_PADDING, EXT_FIT_MAX_RADIUS); + model->fitRadius = rawRadius; // redefine the pixels if needed - bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->fitRadius); + 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->fitRadius, "OR", markVal); - return status; -} + return true; +} Index: branches/sc_branches/trunkTest/psphot/src/psphotReadout.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotReadout.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotReadout.c (revision 29060) @@ -53,6 +53,7 @@ } - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are determined and saved on - // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT + // load the psf model, if suppled. FWHM_MAJ,FWHM_MIN,etc are determined and saved on + // readout->analysis. XXX Note: this function currently only works with a single + // PSPHOT.INPUT if (!psphotLoadPSF (config, view)) { // ??? need to supply 2 ? psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); Index: branches/sc_branches/trunkTest/psphot/src/psphotSetThreads.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotSetThreads.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotSetThreads.c (revision 29060) @@ -15,5 +15,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_PSF_WEIGHTS", 2); + task = psThreadTaskAlloc("PSPHOT_PSF_WEIGHTS", 3); task->function = &psphotPSFWeights_Threaded; psThreadTaskAdd(task); @@ -30,6 +30,11 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 9); + task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10); task->function = &psphotBlendFit_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12); + task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); psFree(task); Index: branches/sc_branches/trunkTest/psphot/src/psphotSignificanceImage.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotSignificanceImage.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotSignificanceImage.c (revision 29060) @@ -97,5 +97,6 @@ smooth_im->data.F32[j][i] = 0.0; } else { - smooth_im->data.F32[j][i] = factor * PS_SQR(value) / smooth_wt->data.F32[j][i]; + float v2 = value + PS_SQR(value/1000.0); + smooth_im->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i]; } } Index: branches/sc_branches/trunkTest/psphot/src/psphotSourceFits.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotSourceFits.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotSourceFits.c (revision 29060) @@ -2,5 +2,4 @@ // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful -// XXX this function does not call pmSourceFitModelInit : fix this? static int NfitPSF = 0; @@ -8,4 +7,5 @@ static int NfitDBL = 0; static int NfitEXT = 0; +static int NfitPCM = 0; bool psphotFitInit (int nThreads) { @@ -18,11 +18,13 @@ psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", - NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); - return true; -} - -bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { + NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); + return true; +} + +bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { float x, y, dR; + + pmSourceFitOptions options = *fitOptions; // maskVal is used to test for rejected pixels, and must include markVal @@ -31,5 +33,5 @@ // if this source is not a possible blend, just fit as PSF if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) { - bool status = psphotFitPSF (readout, source, psf, maskVal, markVal); + bool status = psphotFitPSF (readout, source, psf, fitOptions, maskVal, markVal); return status; } @@ -91,5 +93,6 @@ // fit PSF model - pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); + options.mode = PM_SOURCE_FIT_PSF; + pmSourceFitSet (source, modelSet, &options, maskVal); // clear the circular mask @@ -154,7 +157,8 @@ } -bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { +bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { double chiTrend; + pmSourceFitOptions options = *fitOptions; // maskVal is used to test for rejected pixels, and must include markVal @@ -171,5 +175,6 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF, maskVal); + options.mode = PM_SOURCE_FIT_PSF; + pmSourceFitModel (source, PSF, &options, maskVal); if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); @@ -199,7 +204,8 @@ } +// save a local, static copy of the EXT model type so we don't have to lookup for each object static pmModelType modelTypeEXT; -bool psphotInitLimitsEXT (psMetadata *recipe) { +bool psphotInitLimitsEXT (psMetadata *recipe, pmReadout *readout) { bool status; @@ -208,11 +214,13 @@ char *modelNameEXT = psMetadataLookupStr (&status, recipe, "EXT_MODEL"); modelTypeEXT = pmModelClassGetType (modelNameEXT); - psphotInitRadiusEXT (recipe, modelTypeEXT); - - return true; -} - -bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { - + + psphotInitRadiusEXT (recipe, readout); + + return true; +} + +bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *newSources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { + + float radius; bool okEXT, okDBL; float chiEXT, chiDBL; @@ -222,6 +230,8 @@ pmModel *EXT = NULL; psArray *DBL = NULL; + pmMoments psfMoments; // skip the source if we don't think it is extended + // XXX are these robust, or do we have better info from the source-size analysis?? if (source->type == PM_SOURCE_TYPE_UNKNOWN) return false; if (source->type == PM_SOURCE_TYPE_DEFECT) return false; @@ -229,11 +239,17 @@ // set the radius based on the footprint (also sets the mask pixels) - float radius = psphotSetRadiusEXT (readout, source, markVal); + if (!psphotSetRadiusFootprint(&radius, readout, source, markVal, 1.0)) return false; // XXX note that this changes the source moments that are published... + // XXX all published moments should use the same measurement // 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 // this uses the footprint to judge both radius and aperture? - if (!pmSourceMoments (source, radius, 0.0, 0.0, maskVal)) return false; + // XXX save the psf-based moments for output + psfMoments = *source->moments; + if (!pmSourceMoments (source, radius, 0.0, 0.5, maskVal)) { + *source->moments = psfMoments; + return false; + } psTrace ("psphot", 5, "trying blob...\n"); @@ -244,42 +260,47 @@ // this temporary source is used as a place-holder by the psphotEval functions below tmpSrc = pmSourceAlloc (); - - // XXX need to handle failures better here - EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal); - if (!EXT) goto escape; - if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape; - - okEXT = psphotEvalEXT (tmpSrc, EXT); - chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; - - // DBL will always be defined, but DBL->data[n] might not - DBL = psphotFitDBL (readout, source, maskVal, markVal); - if (!DBL) goto escape; - if (!DBL->n) goto escape; - - okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); - okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); - // XXX should I keep / save the flags set in the eval functions? + { + // DBL will always be defined, but DBL->data[n] might not + DBL = psphotFitDBL (readout, source, fitOptions, maskVal, markVal); + if (!DBL) goto escape; + if (!DBL->n) goto escape; + + okDBL = psphotEvalDBL (tmpSrc, DBL->data[0]); + okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]); + // XXX should I keep / save the flags set in the eval functions? + + // 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; + chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison + ONE->fitRadius = radius; + } + + // correct second model chisqs for flux trend + 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; + ONE->fitRadius = radius; + } + } + + { + // XXX need to handle failures better here + EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal); + if (!EXT) goto escape; + if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape; + + okEXT = psphotEvalEXT (tmpSrc, EXT); + chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; + } // 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; - chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison - } - - // correct second model chisqs for flux trend - 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; - } psFree (tmpSrc); @@ -298,5 +319,8 @@ // both models failed; reject them both - // XXX -- change type flags to psf in this case and keep original moments? + // XXX -- change type flags to psf in this case, and make sure we subtract it? + // reset the psf moments + *source->moments = psfMoments; + psFree (EXT); psFree (DBL); @@ -308,6 +332,6 @@ // save new model + // XXX save the correct radius... source->modelEXT = EXT; - source->modelEXT->fitRadius = radius; source->type = PM_SOURCE_TYPE_EXTENDED; source->mode |= PM_SOURCE_MODE_EXTMODEL; @@ -327,4 +351,6 @@ # endif + // reset the psf moments + *source->moments = psfMoments; return true; @@ -337,10 +363,8 @@ source->modelPSF = psMemIncrRefCounter (DBL->data[0]); source->mode |= PM_SOURCE_MODE_PAIR; - source->modelPSF->fitRadius = radius; // copy most data from the primary source (modelEXT, blends stay NULL) - pmSource *newSrc = pmSourceCopy (source); + pmSource *newSrc = pmSourceCopyData (source); newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); - newSrc->modelPSF->fitRadius = radius; // build cached models and subtract @@ -365,4 +389,8 @@ # endif + // reset the (original) psf moments + *source->moments = psfMoments; + *newSrc->moments = psfMoments; + psArrayAdd (newSources, 100, newSrc); psFree (newSrc); @@ -371,4 +399,6 @@ escape: + // reset the psf moments + *source->moments = psfMoments; psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); psFree (tmpSrc); @@ -379,5 +409,5 @@ // fit a double PSF source to an extended blob -psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) { +psArray *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { float dx, dy; @@ -387,4 +417,6 @@ psEllipseMoments moments; psArray *modelSet; + + pmSourceFitOptions options = *fitOptions; NfitDBL ++; @@ -403,4 +435,8 @@ axes = psEllipseMomentsToAxes (moments, 20.0); + if (isnan(axes.major)) return NULL; + if (isnan(axes.minor)) return NULL; + if (isnan(axes.theta)) return NULL; + // XXX this is really arbitrary: 4 pixel separation? dx = 2 * cos (axes.theta); @@ -426,9 +462,16 @@ // fit PSF model (set/unset the pixel mask) - pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF, maskVal); + options.mode = PM_SOURCE_FIT_PSF; + pmSourceFitSet (source, modelSet, &options, maskVal); return (modelSet); } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) { +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType 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); + } + + pmSourceFitOptions options = *fitOptions; NfitEXT ++; @@ -438,16 +481,224 @@ // use the source moments, etc to guess basic model parameters - pmModel *EXT = pmSourceModelGuess (source, modelType); - if (!EXT) { + pmModel *model = pmSourceModelGuess (source, modelType); + if (!model) { psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); return NULL; } + // for sersic models, use a grid search to choose an index, then float the params there + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + // for the test fits, use a somewhat smaller radius + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + psphotFitSersicIndex (model, readout, source, fitOptions, maskVal, markVal); + } + + if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0); + } + + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + options.mode = PM_SOURCE_FIT_NO_INDEX; + } else { + options.mode = PM_SOURCE_FIT_EXT; + } + + // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); + pmSourceFitModel (source, model, &options, maskVal); + // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + + // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); + return (model); +} + +pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + 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); } + pmSourceFitOptions options = *fitOptions; + + NfitPCM ++; + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // allocate the model + pmModel *model = pmModelAlloc(modelType); + if (!model) { + return NULL; + } + + pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize); + if (!pcm) { + psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy); + model->flags |= PM_MODEL_STATUS_BADARGS; // XXX this is probably already set in pmPCMinit + return model; + } + + // use the source moments, etc to guess basic model parameters + if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) { + psFree (pcm); + model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } + + // for sersic models, use a grid search to choose an index, then float the params there + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + // for the test fits, use a somewhat smaller radius + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + + if (!psphotFitSersicIndexPCM (pcm, readout, source, fitOptions, maskVal, markVal, psfSize)) { + model->flags |= PM_MODEL_STATUS_BADARGS; + return model; + } + } + + if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 1.0); + } + + if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) { + options.mode = PM_SOURCE_FIT_NO_INDEX; + } else { + options.mode = PM_SOURCE_FIT_EXT; + } + // update the pcm elements if we have changed the circumstance (options.mode or source->pixels) + pmPCMupdate(pcm, source, &options, model); + + // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); + // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + + // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0); + psFree (pcm); + + return model; +} + +// note that these should be 1/2n of the standard sersic index +float indexGuess[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083}; +# define N_INDEX_GUESS 6 + +// A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index +// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations +bool psphotFitSersicIndex (pmModel *model, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + // fit EXT (not PSF) model (set/unset the pixel mask) - pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT, maskVal); - return (EXT); -} + options.mode = PM_SOURCE_FIT_NO_INDEX; + options.nIter = 4; + + int iMin = -1; + float xMin = NAN; + float chiSquare[N_INDEX_GUESS]; + + for (int i = 0; i < N_INDEX_GUESS; i++) { + model->params->data.F32[PM_PAR_7] = indexGuess[i]; + + if (!model->modelGuess(model, source)) { + model->flags |= PM_MODEL_STATUS_BADARGS; + return false; + } + + // each time we change the model guess, we need to adjust the radius + // XXX this did not work : we do not need such a large radius -- just uses moments-based radius + if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + } + + pmSourceFitModel (source, model, &options, maskVal); + // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + + chiSquare[i] = model->chisqNorm; + if (i == 0) { + xMin = chiSquare[i]; + iMin = i; + } else { + if (chiSquare[i] < xMin) { + xMin = chiSquare[i]; + iMin = i; + } + } + } + assert (iMin >= 0); + + model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it + model->params->data.F32[PM_PAR_7] = indexGuess[iMin]; + model->modelGuess(model, source); + + // each time we change the model guess, we need to adjust the radius + // if (!psphotSetRadiusModel (model, readout, source, markVal, true)) { + // psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal); + // } + + return true; +} + +// A sersic model is very sensitive to the index. attempt to find the index first by grid search in just the index +// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations +bool psphotFitSersicIndexPCM (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) { + + pmModel *model = pcm->modelConv; + + assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC")); + + pmSourceFitOptions options = *fitOptions; + + // fit EXT (not PSF) model (set/unset the pixel mask) + options.mode = PM_SOURCE_FIT_NO_INDEX; + options.nIter = 4; + + // update the pcm elements if we have changed the circumstance (here, options.mode) + pmPCMupdate(pcm, source, &options, model); + + int iMin = -1; + float xMin = NAN; + float chiSquare[N_INDEX_GUESS]; + + for (int i = 0; i < N_INDEX_GUESS; i++) { + model->params->data.F32[PM_PAR_7] = indexGuess[i]; + + if (!model->modelGuess(model, source)) { + model->flags |= PM_MODEL_STATUS_BADARGS; + return false; + } + + // each time we change the model guess, we need to adjust the radius + // XXX this did not work : we do not need such a large radius -- just uses moments-based radius + if (false && !psphotSetRadiusModel (model, readout, source, markVal, false)) { + psphotSetRadiusFootprint(&model->fitRadius, readout, source, markVal, 0.5); + } + +# if (0) + pmSourceFitModel (source, model, &options, maskVal); +# else + pmSourceModelGuessPCM(pcm, source, maskVal, markVal); + pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); +# endif + // fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix); + + chiSquare[i] = model->chisq; + if (i == 0) { + xMin = chiSquare[i]; + iMin = i; + } else { + if (chiSquare[i] < xMin) { + xMin = chiSquare[i]; + iMin = i; + } + } + } + assert (iMin >= 0); + + model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it + model->params->data.F32[PM_PAR_7] = indexGuess[iMin]; + + pmSourceModelGuessPCM(pcm, source, maskVal, markVal); + + return true; +} Index: branches/sc_branches/trunkTest/psphot/src/psphotSourceSize.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotSourceSize.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotSourceSize.c (revision 29060) @@ -1,4 +1,6 @@ # include "psphotInternal.h" # include + +# define KRON 1 typedef struct { @@ -168,4 +170,5 @@ psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); + psImageMaskType markVal = options->markVal; psImageMaskType maskVal = options->maskVal | options->markVal; @@ -189,5 +192,5 @@ // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit @@ -197,5 +200,10 @@ pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); + // XXX test: switch to kron flux +# if (KRON) + float apMag = -2.5*log10(source->moments->KronFlux); +# else float apMag = -2.5*log10(source->moments->Sum); +# endif float dMag = source->psfMag - apMag; @@ -288,4 +296,5 @@ pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + psImageMaskType markVal = options->markVal; psImageMaskType maskVal = options->maskVal | options->markVal; @@ -335,5 +344,5 @@ // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit @@ -343,16 +352,28 @@ pmSourceSub (source, PM_MODEL_OP_FULL, options->maskVal); +# if (KRON) + float apMag = -2.5*log10(source->moments->KronFlux); +# else float apMag = -2.5*log10(source->moments->Sum); +# endif float dMag = source->psfMag - apMag; // set nSigma to include both systematic and poisson error terms // XXX the 'poisson error' contribution for size is probably wrong... - float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); + // XXX add in a hard floor on the Ap Sys Err (to be a bit generous) + float nSigmaMAG = (dMag - options->ApResid) / hypot(source->errMag, hypot(options->ApSysErr, 0.025)); float nSigmaMXX = (Mxx - psfClump->X) / hypot(psfClump->dX, psfClump->X*psfClump->X*source->errMag); float nSigmaMYY = (Myy - psfClump->Y) / hypot(psfClump->dY, psfClump->Y*psfClump->Y*source->errMag); - // partially-masked sources are more likely to be mis-measured PSFs + // fprintf (stderr, "%f %f : Mxx: %f, Myy: %f, dx: %f, dy: %f, psfMag: %f, apMag: %f, dMag: %f, errMag: %f, nSigmaMag: %f, nSigmaMxx: %f, nSigmaMyy: %f\n", + // source->peak->xf, source->peak->yf, Mxx, Myy, source->peak->xf - source->moments->Mx, source->peak->yf - source->moments->My, + // source->psfMag, apMag, dMag, source->errMag, nSigmaMAG, nSigmaMXX, nSigmaMYY); + + // XXX double check on ths stuff!! partially-masked sources are more likely to be mis-measured PSFs float sizeBias = 1.0; - if (source->pixWeight < 0.9) { + if (source->pixWeightNotBad < 0.9) { + sizeBias = 3.0; + } + if (source->pixWeightNotPoor < 0.9) { sizeBias = 3.0; } @@ -389,5 +410,5 @@ if (isCR) { psTrace("psphotSourceClassRegion.CR",4,"CLASS: %g %g %f\t%g %g %g %g %g %g\t%g %g\t%g CR\t%g %g\n", - source->peak->xf,source->peak->yf,source->pixWeight,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, + source->peak->xf,source->peak->yf,source->pixWeightNotBad,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigmaMAG, options->nSigmaApResid,sizeBias*options->nSigmaMoments); source->mode |= PM_SOURCE_MODE_DEFECT; Index: branches/sc_branches/trunkTest/psphot/src/psphotSourceStats.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotSourceStats.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotSourceStats.c (revision 29060) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +void pmSourceMomentsSetVerbose(bool state); // convert detections to sources and measure their basic properties (moments, local sky, sky @@ -369,4 +370,6 @@ maskVal |= markVal; + // XXX test : pmSourceMomentsSetVerbose(true); + // threaded measurement of the sources moments int Nfail = 0; @@ -408,4 +411,17 @@ // measure basic source moments (no S/N clipping on input pixels) status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); + // XXX moments / aperture test: + if (0) { + // 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, 4.0*source->moments->Mrf, "OR", markVal); + + float apMag = NAN; + pmSourcePhotometryAper (&apMag, NULL, source->pixels, source->maskObj, maskVal); + fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux)); + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + } if (status) { Nmoments ++; @@ -502,5 +518,7 @@ 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, analysis, sources); + if (pmVisualTestLevel("psphot.moments.full", 2)) { + psphotVisualPlotMoments (recipe, analysis, sources); + } #endif Index: branches/sc_branches/trunkTest/psphot/src/psphotStack.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotStack.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotStack.c (revision 29060) @@ -24,4 +24,13 @@ exit (psphotGetExitStatus()); } + // Dump configuration + bool mdok; + psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config + if (dump && !pmConfigDump(config, dump)) { + psError(psErrorCodeLast(), false, "Unable to dump configuration."); + exit (psphotGetExitStatus()); + } + + psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); Index: branches/sc_branches/trunkTest/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotStackMatchPSFsUtils.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotStackMatchPSFsUtils.c (revision 29060) @@ -346,6 +346,6 @@ if (!fake) goto escape; - dumpImage(fake, readoutSrc, index, "fake"); - dumpImage(readoutSrc, readoutSrc, index, "real"); + // dumpImage(fake, readoutSrc, index, "fake"); + // dumpImage(readoutSrc, readoutSrc, index, "real"); if (threads) pmSubtractionThreadsInit(); @@ -380,6 +380,6 @@ } - dumpImage(readoutOut, readoutSrc, index, "conv"); - dumpImageDiff(readoutOut, fake, readoutSrc, index, "diff"); + // dumpImage(readoutOut, readoutSrc, index, "conv"); + // dumpImageDiff(readoutOut, fake, readoutSrc, index, "diff"); psFree(fake); Index: branches/sc_branches/trunkTest/psphot/src/psphotVisual.c =================================================================== --- branches/sc_branches/trunkTest/psphot/src/psphotVisual.c (revision 28957) +++ branches/sc_branches/trunkTest/psphot/src/psphotVisual.c (revision 29060) @@ -66,9 +66,9 @@ int myKapa = psphotKapaChannel (channel); if (!(strcasecmp (overlay, "all"))) { - KiiEraseOverlay (myKapa, "red"); - KiiEraseOverlay (myKapa, "green"); - KiiEraseOverlay (myKapa, "blue"); - KiiEraseOverlay (myKapa, "yellow"); - return true; + KiiEraseOverlay (myKapa, "red"); + KiiEraseOverlay (myKapa, "green"); + KiiEraseOverlay (myKapa, "blue"); + KiiEraseOverlay (myKapa, "yellow"); + return true; } KiiEraseOverlay (myKapa, overlay); @@ -182,5 +182,5 @@ bool psphotVisualShowImage (pmReadout *readout) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image", 1)) return true; int kapa = psphotKapaChannel (1); @@ -199,5 +199,5 @@ pmReadout *backgnd; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image.backgnd", 2)) return true; int kapa = psphotKapaChannel (1); @@ -208,7 +208,7 @@ if (file->mode == PM_FPA_MODE_INTERNAL) { - backgnd = file->readout; + backgnd = file->readout; } else { - backgnd = pmFPAviewThisReadout (view, file->fpa); + backgnd = pmFPAviewThisReadout (view, file->fpa); } @@ -222,5 +222,5 @@ bool psphotVisualShowSignificance (psImage *image, float min, float max) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image.signif", 2)) return true; int kapa = psphotKapaChannel (1); @@ -233,4 +233,20 @@ } +bool psphotVisualShowLogSignificance (psImage *image, float min, float max) { + + if (!pmVisualTestLevel("psphot.image.logsignif", 3)) return true; + + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + psImage *lsig = (psImage *) psUnaryOp (NULL, image, "log"); + psphotVisualRangeImage (kapa, lsig, "log-signif", 2, min, max); + psFree (lsig); + + pmVisualAskUser(NULL); + return true; +} + +// XXX : requires psphotVisualShowImage bool psphotVisualShowPeaks (pmDetections *detections) { @@ -238,5 +254,5 @@ KiiOverlay *overlay; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.peaks", 1)) return true; int kapa = psphotKapaChannel (1); @@ -252,15 +268,15 @@ for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - if (peak == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_BOX; - overlay[Noverlay].x = peak->xf; - overlay[Noverlay].y = peak->yf; - overlay[Noverlay].dx = 2.0; - overlay[Noverlay].dy = 2.0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmPeak *peak = peaks->data[i]; + if (peak == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 2.0; + overlay[Noverlay].dy = 2.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -272,4 +288,5 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowFootprints (pmDetections *detections) { @@ -277,5 +294,5 @@ KiiOverlay *overlay; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.footprints", 3)) return true; int kapa = psphotKapaChannel (1); @@ -292,82 +309,82 @@ for (int i = 0; i < footprints->n; i++) { - pmSpan *span = NULL; - - pmFootprint *footprint = footprints->data[i]; - if (footprint == NULL) continue; - if (footprint->spans == NULL) continue; - if (footprint->spans->n < 1) continue; - - // 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; - overlay[Noverlay].x = span->x0; - overlay[Noverlay].y = span->y; - overlay[Noverlay].dx = span->x1 - span->x0; - overlay[Noverlay].dy = 0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - int ys = span->y; - int x0s = span->x0; - int x1s = span->x1; - - // draw the outer span edges - for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int ye = span1->y; - int x0e = span1->x0; - int x1e = span1->x1; - - // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) - // find all of the spans in this row and generate x0e, x01: - for (int k = j + 1; k < footprint->spans->n; k++) { - pmSpan *span2 = footprint->spans->data[k]; - if (span2->y > span1->y) break; - x0e = PS_MIN (x0e, span2->x0); - x1e = PS_MAX (x1e, span2->x1); - j++; - } - - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = x0s; - overlay[Noverlay].y = ys; - overlay[Noverlay].dx = x0e - x0s; - overlay[Noverlay].dy = ye - ys; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = x1s; - overlay[Noverlay].y = ys; - overlay[Noverlay].dx = x1e - x1s; - overlay[Noverlay].dy = ye - ys; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - ys = ye; - x0s = x0e; - x1s = x1e; - } - - // draw the bottom - span = footprint->spans->data[footprint->spans->n - 1]; - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = span->x0; - overlay[Noverlay].y = span->y; - overlay[Noverlay].dx = span->x1 - span->x0; - overlay[Noverlay].dy = 0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + pmSpan *span = NULL; + + pmFootprint *footprint = footprints->data[i]; + if (footprint == NULL) continue; + if (footprint->spans == NULL) continue; + if (footprint->spans->n < 1) continue; + + // 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; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + int ys = span->y; + int x0s = span->x0; + int x1s = span->x1; + + // draw the outer span edges + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int ye = span1->y; + int x0e = span1->x0; + int x1e = span1->x1; + + // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) + // find all of the spans in this row and generate x0e, x01: + for (int k = j + 1; k < footprint->spans->n; k++) { + pmSpan *span2 = footprint->spans->data[k]; + if (span2->y > span1->y) break; + x0e = PS_MIN (x0e, span2->x0); + x1e = PS_MAX (x1e, span2->x1); + j++; + } + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x0s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x0e - x0s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x1s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x1e - x1s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + ys = ye; + x0s = x0e; + x1s = x1e; + } + + // draw the bottom + span = footprint->spans->data[footprint->spans->n - 1]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); } @@ -379,4 +396,5 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowMoments (psArray *sources) { @@ -387,5 +405,5 @@ psEllipseAxes axes; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.moments", 2)) return true; int kapa = psphotKapaChannel (1); @@ -401,27 +419,27 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) 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.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - - overlay[Noverlay].dx = 2.0*axes.major; - overlay[Noverlay].dy = 2.0*axes.minor; - - overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; - - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) 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.xy = moments->Mxy; + emoments.y2 = moments->Myy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -439,5 +457,5 @@ KapaSection section; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.moments", 1)) return true; int myKapa = psphotKapaChannel (2); @@ -456,26 +474,26 @@ float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, analysis, "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, analysis, 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 - 4.0*psfdX; - float X1 = psfX + 4.0*psfdX; - float Y0 = psfY - 4.0*psfdY; - float Y1 = psfY + 4.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); } - } + int nRegions = psMetadataLookupS32 (&status, analysis, "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, analysis, 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 - 4.0*psfdX; + float X1 = psfX + 4.0*psfdX; + float Y0 = psfY - 4.0*psfdY; + float Y1 = psfY + 4.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); } + } } Xmin = PS_MAX(Xmin, -0.1); @@ -498,21 +516,21 @@ int nF = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->moments == NULL) - continue; - - xFaint->data.F32[nF] = source->moments->Mxx; - yFaint->data.F32[nF] = source->moments->Myy; - mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum); - nF++; - - // XXX make this a user-defined cutoff - if (source->moments->SN < SN_LIM) - continue; - - xBright->data.F32[nB] = source->moments->Mxx; - yBright->data.F32[nB] = source->moments->Myy; - mBright->data.F32[nB] = -2.5*log10(source->moments->Sum); - nB++; + pmSource *source = sources->data[i]; + if (source->moments == NULL) + continue; + + xFaint->data.F32[nF] = source->moments->Mxx; + yFaint->data.F32[nF] = source->moments->Myy; + mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum); + nF++; + + // XXX make this a user-defined cutoff + if (source->moments->SN < SN_LIM) + continue; + + xBright->data.F32[nB] = source->moments->Mxx; + yBright->data.F32[nB] = source->moments->Myy; + mBright->data.F32[nB] = -2.5*log10(source->moments->Sum); + nB++; } xFaint->n = nF; @@ -652,45 +670,45 @@ // 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, analysis, "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, analysis, 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); + 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, analysis, "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, analysis, 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); } @@ -721,28 +739,28 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - if (source->type != type) continue; - if (mode && !(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 = 2.0*axes.major; - overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + if (source->type != type) continue; + if (mode && !(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 = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -753,7 +771,8 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowRoughClass (psArray *sources) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.size", 3)) return true; int myKapa = psphotKapaChannel (1); @@ -776,5 +795,5 @@ bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.model", 1)) return true; int myKapa = psphotKapaChannel (3); @@ -797,28 +816,28 @@ // generate a fake model at each of the 3x3 image grid positions for (int x = -2; x <= +2; x ++) { - for (int y = -2; y <= +2; y ++) { - // use the center of the center pixel of the image - float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5; - float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5; - - // assign the x and y coords to the image center - // create an object with center intensity of 1000 - modelRef->params->data.F32[PM_PAR_SKY] = 0; - modelRef->params->data.F32[PM_PAR_I0] = 1000; - modelRef->params->data.F32[PM_PAR_XPOS] = xc; - modelRef->params->data.F32[PM_PAR_YPOS] = yc; - - // create modelPSF from this model - pmModel *model = pmModelFromPSF (modelRef, psf); - if (!model) continue; - - // place the reference object in the image center - // no need to mask the source here - // XXX should we measure this for the analytical model only or the full model? - pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - psFree (model); - } + for (int y = -2; y <= +2; y ++) { + // use the center of the center pixel of the image + float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5; + float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + modelRef->params->data.F32[PM_PAR_XPOS] = xc; + modelRef->params->data.F32[PM_PAR_YPOS] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + if (!model) continue; + + // place the reference object in the image center + // no need to mask the source here + // XXX should we measure this for the analytical model only or the full model? + pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + psFree (model); + } } @@ -842,5 +861,5 @@ bool status; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.stars", 2)) return true; int myKapa = psphotKapaChannel (3); @@ -871,32 +890,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - bool keep = false; - keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); - if (!keep) continue; - - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows - - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + bool keep = false; + keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); + if (!keep) continue; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } NY += DY; @@ -918,59 +937,59 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - bool keep = false; - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - nPSF ++; - keep = true; - } - if (!keep) continue; - - if (Xo + DX > NX) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - 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); - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - - 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); - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - Xo = DX; - dY = DY; - } - } else { - // extend this row - 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); - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - Xo += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + bool keep = false; + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + nPSF ++; + keep = true; + } + if (!keep) continue; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + 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); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + + 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); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + Xo = DX; + dY = DY; + } + } else { + // extend this row + 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); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + Xo += DX; + dY = PS_MAX (dY, DY); + } } @@ -994,5 +1013,5 @@ bool status; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.sat", 3)) return true; int myKapa = psphotKapaChannel (3); @@ -1023,32 +1042,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // 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? - // DX = source->pixels->numCols - // DY = source->pixels->numRows - - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + // 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? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } NY += DY; @@ -1068,46 +1087,46 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // 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) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - 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); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - - 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); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - Xo = DX; - dY = DY; - } - } else { - // extend this row - 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); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - Xo += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + // 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) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + 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); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + + 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); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + Xo = DX; + dY = DY; + } + } else { + // extend this row + 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); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + Xo += DX; + dY = PS_MAX (dY, DY); + } } @@ -1151,19 +1170,19 @@ float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; for (int iy = 0; iy < source->pixels->numRows; iy++) { - for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { - rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; - // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; - Rb->data.F32[nb] = log10(rb->data.F32[nb]); - fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); - nb++; - } else { - rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; - // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; - Rg->data.F32[ng] = log10(rg->data.F32[ng]); - fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]); - ng++; - } - } + 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] = log10(rb->data.F32[nb]); + fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); + nb++; + } else { + rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; + Rg->data.F32[ng] = log10(rg->data.F32[ng]); + fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]); + ng++; + } + } } @@ -1358,5 +1377,5 @@ KapaSection section; // put the positive profile in one and the residuals in another? - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.profiles", 3)) return true; int myKapa = psphotKapaChannel (2); @@ -1394,22 +1413,22 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - - psphotVisualPlotRadialProfile (myKapa, source, maskVal); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - if (key[0] == 'e') { - KapaClearPlots (myKapa); - } - if (key[0] == 's') { - break; - } + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + psphotVisualPlotRadialProfile (myKapa, source, maskVal); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + if (key[0] == 'e') { + KapaClearPlots (myKapa); + } + if (key[0] == 's') { + break; + } } @@ -1429,5 +1448,5 @@ return true; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.flags", 3)) return true; int myKapa = psphotKapaChannel (1); @@ -1445,62 +1464,62 @@ for (int i = 0; i < sources->n; i++) { - float Xo, Yo, Rmaj, Rmin, cs, sn; - - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - pmMoments *moments = source->moments; - if (0) { - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - Xo = moments->Mx; - Yo = moments->My; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - Rmaj = 2.0*axes.major; - Rmin = 2.0*axes.minor; - cs = cos(axes.theta); - sn = sin(axes.theta); - } else { - Rmaj = Rmin = 5.0; - cs = 1.0; - sn = 0.0; - Xo = source->peak->xf; - Yo = source->peak->yf; - } - - unsigned short int flagMask = 0x01; - for (int j = 0; j < 8; j++) { - if (source->mode & flagMask) { - overlayE[NoverlayE].type = KII_OVERLAY_LINE; - overlayE[NoverlayE].x = Xo; - overlayE[NoverlayE].y = Yo; - - float phi = j*M_PI/4.0; - overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; - overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; - overlayE[NoverlayE].angle = 0; - overlayE[NoverlayE].text = NULL; - NoverlayE ++; - CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100); - } - flagMask <<= 1; - - if (source->mode & flagMask) { - overlayO[NoverlayO].type = KII_OVERLAY_LINE; - overlayO[NoverlayO].x = Xo + 1; - overlayO[NoverlayO].y = Yo; - - float phi = j*M_PI/4.0; - overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; - overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; - overlayO[NoverlayO].angle = 0; - overlayO[NoverlayO].text = NULL; - NoverlayO ++; - CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100); - } - flagMask <<= 1; - } + float Xo, Yo, Rmaj, Rmin, cs, sn; + + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + pmMoments *moments = source->moments; + if (0) { + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + Xo = moments->Mx; + Yo = moments->My; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + Rmaj = 2.0*axes.major; + Rmin = 2.0*axes.minor; + cs = cos(axes.theta); + sn = sin(axes.theta); + } else { + Rmaj = Rmin = 5.0; + cs = 1.0; + sn = 0.0; + Xo = source->peak->xf; + Yo = source->peak->yf; + } + + unsigned short int flagMask = 0x01; + for (int j = 0; j < 8; j++) { + if (source->mode & flagMask) { + overlayE[NoverlayE].type = KII_OVERLAY_LINE; + overlayE[NoverlayE].x = Xo; + overlayE[NoverlayE].y = Yo; + + float phi = j*M_PI/4.0; + overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; + overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; + overlayE[NoverlayE].angle = 0; + overlayE[NoverlayE].text = NULL; + NoverlayE ++; + CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100); + } + flagMask <<= 1; + + if (source->mode & flagMask) { + overlayO[NoverlayO].type = KII_OVERLAY_LINE; + overlayO[NoverlayO].x = Xo + 1; + overlayO[NoverlayO].y = Yo; + + float phi = j*M_PI/4.0; + overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; + overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; + overlayO[NoverlayO].angle = 0; + overlayO[NoverlayO].text = NULL; + NoverlayO ++; + CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100); + } + flagMask <<= 1; + } } @@ -1531,33 +1550,33 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) 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 ++; + pmSource *source = sources->data[i]; + if (source == NULL) 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 ++; } @@ -1570,5 +1589,5 @@ bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.size", 2)) return true; int myKapa = psphotKapaChannel (1); @@ -1597,5 +1616,5 @@ KapaSection section; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.size", 2)) return true; int myKapa = psphotKapaChannel (2); @@ -1612,26 +1631,26 @@ float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, analysis, "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, analysis, 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); } - } + int nRegions = psMetadataLookupS32 (&status, analysis, "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, analysis, 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); } + } } Xmin = PS_MAX(Xmin, -0.1); @@ -1677,53 +1696,53 @@ int nCR = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->moments == NULL) continue; + pmSource *source = sources->data[i]; + if (source->moments == NULL) continue; // only plot the measured sources... - if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue; - - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { - 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->errMag > 0.1) { - xLOW->data.F32[nLOW] = source->moments->Mxx; - yLOW->data.F32[nLOW] = source->moments->Myy; - mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum); - sLOW->data.F32[nLOW] = source->extNsigma; - nLOW++; - continue; - } - 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++; + if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) 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->errMag > 0.1) { + xLOW->data.F32[nLOW] = source->moments->Mxx; + yLOW->data.F32[nLOW] = source->moments->Myy; + mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum); + sLOW->data.F32[nLOW] = source->extNsigma; + nLOW++; + continue; + } + 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++; } @@ -2029,45 +2048,45 @@ // 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, analysis, "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, analysis, 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); + 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, analysis, "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, analysis, 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); } @@ -2108,5 +2127,5 @@ bool psphotVisualShowResidualImage (pmReadout *readout) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image.resid", 2)) return true; int myKapa = psphotKapaChannel (1); @@ -2124,5 +2143,5 @@ float lineX[2], lineY[2]; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.apresid", 1)) return true; int myKapa = psphotKapaChannel (2); @@ -2144,19 +2163,19 @@ 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->apMag)) continue; - if (!isfinite (source->psfMag)) continue; - - 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]); - graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); - graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); - - n++; + pmSource *source = sources->data[i]; + if (!source) continue; + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (!isfinite (source->apMag)) continue; + if (!isfinite (source->psfMag)) continue; + + 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]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + + n++; } x->n = y->n = dy->n = n; @@ -2242,5 +2261,5 @@ Graphdata graphdata; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.chisq", 1)) return true; int myKapa = psphotKapaChannel (2); @@ -2263,22 +2282,22 @@ 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 (!source->moments) continue; - if (!isfinite(source->moments->Sum)) continue; - if (!source->modelPSF) continue; - if (!isfinite(source->modelPSF->chisq)) continue; - - x->data.F32[n] = -2.5*log10(source->moments->Sum); - y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF; - graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); - graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); - graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); - graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); - - fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]); - - n++; + pmSource *source = sources->data[i]; + if (!source) continue; + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (!source->moments) continue; + if (!isfinite(source->moments->Sum)) continue; + if (!source->modelPSF) continue; + if (!isfinite(source->modelPSF->chisq)) continue; + + x->data.F32[n] = -2.5*log10(source->moments->Sum); + y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + + fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]); + + n++; } x->n = y->n = n; @@ -2326,5 +2345,5 @@ KiiOverlay *overlay; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.petro", 2)) return true; int kapa = psphotKapaChannel (1); @@ -2336,32 +2355,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - if (!source) continue; - if (!source->extpars) continue; - if (!source->extpars->petProfile) continue; - - float petrosianRadius = source->extpars->petrosianRadius; + pmSource *source = sources->data[i]; + + if (!source) continue; + if (!source->extpars) continue; + if (!source->extpars->petProfile) continue; + + float petrosianRadius = source->extpars->petrosianRadius; psEllipseAxes *axes = &source->extpars->axes; - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - overlay[Noverlay].dx = 1.0*petrosianRadius; - overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major; - overlay[Noverlay].angle = 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 = 2.0*petrosianRadius; - // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major; - // overlay[Noverlay].angle = 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 = 1.0*petrosianRadius; + overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major; + overlay[Noverlay].angle = 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 = 2.0*petrosianRadius; + // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major; + // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; + // overlay[Noverlay].text = NULL; + // Noverlay ++; + // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); }