Index: branches/eam_branches/20090715/psModules/src/extras/pmVisual.c =================================================================== --- branches/eam_branches/20090715/psModules/src/extras/pmVisual.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/extras/pmVisual.c (revision 25496) @@ -24,4 +24,5 @@ #include "pmPSF.h" #include "pmPSFtry.h" +#include "pmSource.h" #include "pmFPAExtent.h" Index: branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/imcombine/pmPSFEnvelope.c (revision 25496) @@ -320,5 +320,6 @@ options->poissonErrorsParams = true; options->stats = psStatsAlloc(PSF_STATS); - options->radius = maxRadius; + options->fitRadius = maxRadius; + options->apRadius = maxRadius; // XXX need to decide if aperture mags need a different radius options->psfTrendMode = PM_TREND_MAP; options->psfTrendNx = xOrder; Index: branches/eam_branches/20090715/psModules/src/objects/pmModel.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmModel.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmModel.c (revision 25496) @@ -62,5 +62,5 @@ tmp->nPix = 0; tmp->nIter = 0; - tmp->radiusFit = 0; + tmp->fitRadius = 0; tmp->flags = PM_MODEL_STATUS_NONE; tmp->residuals = NULL; // XXX should the model own this memory? @@ -109,5 +109,5 @@ new->nIter = model->nIter; new->flags = model->flags; - new->radiusFit = model->radiusFit; + new->fitRadius = model->fitRadius; for (int i = 0; i < new->params->n; i++) { Index: branches/eam_branches/20090715/psModules/src/objects/pmModel.h =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmModel.h (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmModel.h (revision 25496) @@ -96,5 +96,5 @@ int nIter; ///< number of iterations to reach min pmModelStatus flags; ///< model status flags - float radiusFit; ///< fit radius actually used + float fitRadius; ///< fit radius actually used pmResiduals *residuals; ///< normalized PSF residuals Index: branches/eam_branches/20090715/psModules/src/objects/pmPSF.h =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPSF.h (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmPSF.h (revision 25496) @@ -79,5 +79,6 @@ bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting bool poissonErrorsParams; ///< use poission errors for model parameter fitting - float radius; + float fitRadius; + float apRadius; bool chiFluxTrend; // Fit a trend in Chi2 as a function of flux? } pmPSFOptions; Index: branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitEXT.c (revision 25496) @@ -70,5 +70,5 @@ // set object mask to define valid pixels // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); // fit model as EXT, not PSF @@ -76,5 +76,5 @@ // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); // exclude the poor fits Index: branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmPSFtryFitPSF.c (revision 25496) @@ -67,9 +67,10 @@ continue; } - source->modelPSF->radiusFit = options->radius; - // XXXX use a different radius for the aperture magnitude than for the PSF fit? + // PSF fit and aperture mags use different radii + source->modelPSF->fitRadius = options->fitRadius; + source->apRadius = options->apRadius; - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); + // set object mask to define valid pixels for PSF model fit + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal); // fit the PSF model to the source @@ -78,5 +79,5 @@ // skip poor fits if (!status) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); @@ -84,8 +85,19 @@ } - // XXX this function calculates the aperture magnitude, but we now use the moments->Sum as the flux + // set object mask to define valid pixels for APERTURE magnitude + if (options->fitRadius < options->apRadius) { + // this is not really a recommended situation... + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal); + } + if (options->fitRadius > options->apRadius) { + // ensure the aperture is defined + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal); + } + + // This function calculates the psf and aperture magnitudes status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); if (!status || isnan(source->apMag)) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); @@ -94,8 +106,9 @@ // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); psfTry->fitMag->data.F32[i] = source->psfMag; - psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag; + // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag; + psfTry->metric->data.F32[i] = source->apMag - source->psfMag; psfTry->metricErr->data.F32[i] = source->errMag; @@ -104,4 +117,7 @@ } psfTry->psf->nPSFstars = Npsf; + + // XXX test code: + pmSourceVisualShowModelFits (psfTry->psf, psfTry->sources, maskVal); psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, psfTry->sources->n); Index: branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmPSFtryMakePSF.c (revision 25496) @@ -184,5 +184,9 @@ // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. // This way, the parameters masked by one of the fits will be applied to the others - for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { + // NOTE : trend->stats (below) points to the same data as psfTrendStats; set the value to 1 + // to limit the iteration within the loop + int nIter = psf->psfTrendStats->clipIter; + psf->psfTrendStats->clipIter = 1; + for (int i = 0; i < nIter; i++) { // XXX we are using the same stats structure on each pass: do we need to re-init it? // XXX we hardwire this to SAMPLE stats above (psphotChoosePSF.c), hardwire here instead? @@ -203,5 +207,5 @@ psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); - // pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); + pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); trend = psf->params->data[PM_PAR_E1]; @@ -212,5 +216,5 @@ psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); - // pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); + pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); trend = psf->params->data[PM_PAR_E2]; @@ -221,5 +225,5 @@ psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); if (psf->psfTrendMode == PM_TREND_MAP) psImageMapCleanup (trend->map); - // pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); + pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); if (!status) { @@ -228,4 +232,5 @@ } } + psf->psfTrendStats->clipIter = nIter; // test dump of the psf parameters Index: branches/eam_branches/20090715/psModules/src/objects/pmSource.h =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSource.h (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSource.h (revision 25496) @@ -82,4 +82,5 @@ float extNsigma; ///< Nsigma deviation from PSF to EXT float sky, skyErr; ///< The sky and its error at the center of the object + float apRadius; psRegion region; ///< area on image covered by selected pixels pmSourceExtendedPars *extpars; ///< extended source parameters Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V1.c (revision 25496) @@ -128,5 +128,5 @@ nDOF = model->nDOF; nPix = model->nPix; - apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + apRadius = source->apRadius; errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { @@ -308,4 +308,5 @@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -313,5 +314,4 @@ model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); source->moments = pmMomentsAlloc (); Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25496) @@ -128,5 +128,5 @@ nDOF = model->nDOF; nPix = model->nPix; - apRadius = model->radiusFit; // XXX should we really use the fitRadius for aperture Radius? + apRadius = source->apRadius; errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0]; } else { @@ -308,4 +308,5 @@ source->crNsigma = psMetadataLookupF32 (&status, row, "CR_NSIGMA"); source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA"); + source->apRadius = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); // note that some older versions used PSF_PROBABILITY: this was not well defined. @@ -313,5 +314,4 @@ model->nDOF = psMetadataLookupS32 (&status, row, "PSF_NDOF"); model->nPix = psMetadataLookupS32 (&status, row, "PSF_NPIX"); - model->radiusFit = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS"); source->moments = pmMomentsAlloc (); Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_RAW.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_RAW.c (revision 25496) @@ -119,5 +119,5 @@ logChi, logChiNorm, source[0].peak->SN, - model[0].radiusFit, + source[0].apRadius, source[0].pixWeight, model[0].nDOF, @@ -181,5 +181,5 @@ log10(model[0].chisqNorm/model[0].nDOF), source[0].peak->SN, - model[0].radiusFit, + source[0].apRadius, source[0].pixWeight, model[0].nDOF, Index: branches/eam_branches/20090715/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourcePhotometry.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourcePhotometry.c (revision 25496) @@ -201,9 +201,10 @@ if (isfinite (source->apMag) && isPSF && psf) { if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) { - source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusFit); + source->apMag += pmGrowthCurveCorrect (psf->growth, source->apRadius); } if (mode & PM_SOURCE_PHOT_APCORR) { + // XXX this should be removed -- we no longer fit for the 'sky bias' rflux = pow (10.0, 0.4*source->psfMag); - source->apMag -= PS_SQR(model->radiusFit)*rflux * psf->skyBias + psf->skySat / rflux; + source->apMag -= PS_SQR(source->apRadius)*rflux * psf->skyBias + psf->skySat / rflux; } } Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.c (revision 25496) @@ -7,4 +7,5 @@ #include "pmPSF.h" #include "pmPSFtry.h" +#include "pmSource.h" #include "pmSourceVisual.h" @@ -17,6 +18,6 @@ static int kapa1 = -1; +static int kapa2 = -1; static bool plotPSF = true; -// static int kapa2 = -1; // static int kapa3 = -1; @@ -33,5 +34,5 @@ Graphdata graphdata; - if (!pmVisualIsVisual()) return true; + // if (!pmVisualIsVisual()) return true; if (kapa1 == -1) { @@ -44,5 +45,5 @@ } - KapaClearPlots (kapa1); + KapaClearSections (kapa1); KapaInitGraph (&graphdata); @@ -112,4 +113,141 @@ } +// to see the structure of the psf model, place the sources in a fake image 1/10th the size +// at their appropriate relative location. later sources stomp on earlier sources +bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal) { + + // if (!pmVisualIsVisual()) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + // create images 1/10 scale: + psImage *image = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32); + psImage *model = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32); + psImage *resid = psImageAlloc (0.1*psf->fieldNx, 0.1*psf->fieldNy, PS_TYPE_F32); + psImageInit (image, 0.0); + psImageInit (model, 0.0); + psImageInit (resid, 0.0); + + FILE *f = fopen ("stats.dat", "w"); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) continue; + if (!source->pixels) continue; + + pmSourceCacheModel (source, maskVal); + if (!source->modelFlux) continue; + + pmModel *srcModel = pmSourceGetModel (NULL, source); + if (!model) continue; + + float norm = srcModel->params->data.F32[PM_PAR_I0]; + + int Xo = 0.1*srcModel->params->data.F32[PM_PAR_XPOS]; + int Yo = 0.1*srcModel->params->data.F32[PM_PAR_YPOS]; + + // insert source pixels in the image at 1/10th offset (XXX ignore centering for now) + // XXX calculate the pixel scatter values, chisq + for (int iy = 0; iy < source->pixels->numRows; iy++) { + int jy = iy + Yo; + if (jy >= image->numRows) continue; + for (int ix = 0; ix < source->pixels->numCols; ix++) { + int jx = ix + Xo; + if (jx >= image->numCols) continue; + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue; + if (source->modelFlux->data.F32[iy][ix] < 0.001) continue; + image->data.F32[jy][jx] = source->pixels->data.F32[iy][ix]; + model->data.F32[jy][jx] = source->modelFlux->data.F32[iy][ix]; + resid->data.F32[jy][jx] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix]; + } + } + + fprintf (f, "%d, %d -> %f %f %f %f\n", Xo, Yo, norm, source->psfMag, source->apMag, -2.5*log10(source->moments->Sum)); + } + fclose (f); + + // KapaClearSections (kapa2); + pmVisualScaleImage (kapa2, image, "image", 0, true); + pmVisualScaleImage (kapa2, model, "model", 1, true); + pmVisualScaleImage (kapa2, resid, "resid", 2, true); + + { + psFits *fits = psFitsOpen ("image.fits", "w"); + psFitsWriteImage (fits, NULL, image, 0, NULL); + psFitsClose (fits); + fits = psFitsOpen ("model.fits", "w"); + psFitsWriteImage (fits, NULL, model, 0, NULL); + psFitsClose (fits); + fits = psFitsOpen ("resid.fits", "w"); + psFitsWriteImage (fits, NULL, resid, 0, NULL); + psFitsClose (fits); + } + + psFree (image); + psFree (model); + psFree (resid); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + +bool pmSourceVisualShowModelFit (pmSource *source) { + + // if (!pmVisualIsVisual()) return true; + if (!source->pixels) return false; + if (!source->modelFlux) return false; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "pmSource:images"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + // KapaClearSections (kapa2); + pmVisualScaleImage (kapa2, source->pixels, "source", 0, false); + pmVisualScaleImage (kapa2, source->modelFlux, "model", 1, false); + + pmModel *model = pmSourceGetModel (NULL, source); + float norm = model->params->data.F32[PM_PAR_I0]; + + psImage *resid = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32); + 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]) { + resid->data.F32[iy][ix] = NAN; + continue; + } + resid->data.F32[iy][ix] = source->pixels->data.F32[iy][ix] - norm*source->modelFlux->data.F32[iy][ix]; + } + } + pmVisualScaleImage (kapa2, resid, "resid", 2, false); + + psFree (resid); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + return true; +} + bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) { @@ -121,5 +259,5 @@ // XXX skip for now - return true; + // return true; if (kapa1 == -1) { @@ -135,38 +273,57 @@ KapaInitGraph (&graphdata); - float min = +1e32; - float max = -1e32; - float Min = +1e32; - float Max = -1e32; + // float min = +1e32; + // float max = -1e32; + // float Min = +1e32; + // float Max = -1e32; + + float Xmin = +1e32; + float Xmax = -1e32; + float Ymin = +1e32; + float Ymax = -1e32; + float Fmin = +1e32; + float Fmax = -1e32; psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32); psVector *model = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *xm = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *ym = psVectorAlloc (x->n, PS_TYPE_F32); + psVector *Fm = psVectorAlloc (x->n, PS_TYPE_F32); + + int n = 0; for (int i = 0; i < x->n; i++) { model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]); resid->data.F32[i] = param->data.F32[i] - model->data.F32[i]; if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - min = PS_MIN (min, resid->data.F32[i]); - max = PS_MAX (max, resid->data.F32[i]); - Min = PS_MIN (min, param->data.F32[i]); - Max = PS_MAX (max, param->data.F32[i]); - } - - psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32); - psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32); - for (int i = 0; i < x->n; i++) { - xn->data.F32[i] = x->data.F32[i] / 5000.0; - yn->data.F32[i] = y->data.F32[i] / 5000.0; - zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min); - Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min); - Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min); - } + Xmin = PS_MIN (Xmin, x->data.F32[i]); + Xmax = PS_MAX (Xmax, x->data.F32[i]); + Ymin = PS_MIN (Ymin, y->data.F32[i]); + Ymax = PS_MAX (Ymax, y->data.F32[i]); + Fmin = PS_MIN (Fmin, param->data.F32[i]); + Fmax = PS_MAX (Fmax, param->data.F32[i]); + xm->data.F32[n] = x->data.F32[i]; + ym->data.F32[n] = y->data.F32[i]; + Fm->data.F32[n] = param->data.F32[i]; + n++; + } + xm->n = ym->n = Fm->n = n; + + // psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32); + // psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32); + // psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32); + // psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32); + // psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32); + // for (int i = 0; i < x->n; i++) { + // xn->data.F32[i] = x->data.F32[i] / 5000.0; + // yn->data.F32[i] = y->data.F32[i] / 5000.0; + // zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min); + // Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min); + // Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min); + // } // view 1 on resid - section.dx = 0.5; - section.dy = 0.33; + section.dx = 1.0; + section.dy = 0.5; section.x = 0.0; section.y = 0.0; @@ -175,17 +332,97 @@ KapaSetSection (kapa1, §ion); psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Xmin; + graphdata.xmax = Xmax; + graphdata.ymin = Fmin; + graphdata.ymax = Fmax; + + { + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + } + + KapaSetLimits (kapa1, &graphdata); + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "X (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM); + + graphdata.ptype = 2; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (kapa1, x->n, &graphdata); + KapaPlotVector (kapa1, x->n, x->data.F32, "x"); + KapaPlotVector (kapa1, x->n, param->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, xm->n, &graphdata); + KapaPlotVector (kapa1, xm->n, xm->data.F32, "x"); + KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, x->n, &graphdata); + KapaPlotVector (kapa1, x->n, x->data.F32, "x"); + KapaPlotVector (kapa1, x->n, model->data.F32, "y"); // view 2 on resid - section.dx = 0.5; - section.dy = 0.33; - section.x = 0.5; - section.y = 0.0; + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.5; section.name = NULL; psStringAppend (§ion.name, "a2"); KapaSetSection (kapa1, §ion); psFree (section.name); - pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); - + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = Ymin; + graphdata.xmax = Ymax; + graphdata.ymin = Fmin; + graphdata.ymax = Fmax; + { + float range; + range = graphdata.xmax - graphdata.xmin; + graphdata.xmax += 0.05*range; + graphdata.xmin -= 0.05*range; + range = graphdata.ymax - graphdata.ymin; + graphdata.ymax += 0.05*range; + graphdata.ymin -= 0.05*range; + } + + KapaSetLimits (kapa1, &graphdata); + KapaSetFont (kapa1, "helvetica", 14); + KapaBox (kapa1, &graphdata); + KapaSendLabel (kapa1, "Y (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapa1, "Model Param", KAPA_LABEL_YM); + + graphdata.ptype = 2; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (kapa1, y->n, &graphdata); + KapaPlotVector (kapa1, y->n, y->data.F32, "x"); + KapaPlotVector (kapa1, y->n, param->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, xm->n, &graphdata); + KapaPlotVector (kapa1, xm->n, ym->data.F32, "x"); + KapaPlotVector (kapa1, xm->n, Fm->data.F32, "y"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 1; + KapaPrepPlot (kapa1, y->n, &graphdata); + KapaPlotVector (kapa1, y->n, y->data.F32, "x"); + KapaPlotVector (kapa1, y->n, model->data.F32, "y"); + +# if (0) // view 3 on resid section.dx = 0.5; @@ -231,11 +468,8 @@ psFree (section.name); pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG); +# endif psFree (resid); - - psFree (xn); - psFree (yn); - psFree (zn); - psFree (Zn); + psFree (model); // pause and wait for user input: Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h (revision 25476) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceVisual.h (revision 25496) @@ -19,4 +19,6 @@ bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask); bool pmSourceVisualPlotPSFMetric (pmPSFtry *try); +bool pmSourceVisualShowModelFit (pmSource *source); +bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal); /// @} Index: branches/eam_branches/20090715/psModules/test/objects/tap_pmGrowthCurve.c =================================================================== --- branches/eam_branches/20090715/psModules/test/objects/tap_pmGrowthCurve.c (revision 25476) +++ branches/eam_branches/20090715/psModules/test/objects/tap_pmGrowthCurve.c (revision 25496) @@ -131,30 +131,30 @@ source->mode = PM_SOURCE_MODE_PSFSTAR; - source->modelPSF->radiusFit = 15.0; + source->apRadius = 15.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); double refMag = source->apMag; - source->modelPSF->radiusFit = 10.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 8.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 6.0; + source->apRadius = 10.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 8.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 6.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0003, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 4.0; + source->apRadius = 4.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0020, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 3.0; + source->apRadius = 3.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 2.0; + source->apRadius = 2.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0075, 0.0001, "growth offset is is %f", refMag - source->apMag); @@ -234,29 +234,29 @@ source->mode = PM_SOURCE_MODE_PSFSTAR; - source->modelPSF->radiusFit = 15.0; + source->apRadius = 15.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); double refMag = source->apMag; - source->modelPSF->radiusFit = 10.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 8.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 6.0; + source->apRadius = 10.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 8.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 6.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0004, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 4.0; + source->apRadius = 4.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0026, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 3.0; + source->apRadius = 3.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0001, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 2.0; + source->apRadius = 2.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0103, 0.0001, "growth offset is is %f", refMag - source->apMag); @@ -336,29 +336,29 @@ source->mode = PM_SOURCE_MODE_PSFSTAR; - source->modelPSF->radiusFit = 15.0; + source->apRadius = 15.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); double refMag = source->apMag; - source->modelPSF->radiusFit = 10.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 8.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 6.0; + source->apRadius = 10.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 8.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 6.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0006, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 4.0; + source->apRadius = 4.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, +0.0038, 0.0001, "growth offset is is %f", refMag - source->apMag); - source->modelPSF->radiusFit = 3.0; - pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); - ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); - - source->modelPSF->radiusFit = 2.0; + source->apRadius = 3.0; + pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); + ok_float_tol(refMag - source->apMag, +0.0000, 0.0001, "growth offset is is %f", refMag - source->apMag); + + source->apRadius = 2.0; pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH | PM_SOURCE_PHOT_INTERP, 0); ok_float_tol(refMag - source->apMag, -0.0164, 0.0001, "growth offset is is %f", refMag - source->apMag); Index: branches/eam_branches/20090715/psModules/test/objects/tap_pmModel.c =================================================================== --- branches/eam_branches/20090715/psModules/test/objects/tap_pmModel.c (revision 25476) +++ branches/eam_branches/20090715/psModules/test/objects/tap_pmModel.c (revision 25496) @@ -77,5 +77,5 @@ ok(model->nDOF == 0, "pmModelAlloc() set pmModel->nDOF correctly"); ok(model->nIter == 0, "pmModelAlloc() set pmModel->nIter correctly"); - ok(model->radiusFit == 0, "pmModelAlloc() set pmModel->radiusFit correctly"); + ok(model->fitRadius == 0, "pmModelAlloc() set pmModel->fitRadius correctly"); ok(model->flags == PM_MODEL_STATUS_NONE, "pmModelAlloc() set pmModel->flags correctly"); ok(model->residuals == NULL, "pmModelAlloc() set pmModel->residuals correctly"); @@ -132,5 +132,5 @@ modelSrc->nIter = 3; modelSrc->flags = PM_MODEL_STATUS_NONE; - modelSrc->radiusFit = 4; + modelSrc->fitRadius = 4; pmModel *modelDst = pmModelCopy(modelSrc); ok(modelDst != NULL && psMemCheckModel(modelDst), "pmModelCopy() returned a non-NULL pmModel"); @@ -139,5 +139,5 @@ ok(modelDst->nIter == 3, "pmModelCopy() set the pmModel->nIter member correctly"); ok(modelDst->flags == PM_MODEL_STATUS_NONE, "pmModelCopy() set the pmModel->flags member correctly"); - ok(modelDst->radiusFit == 4, "pmModelCopy() set the pmModel->radiusFit member correctly"); + ok(modelDst->fitRadius == 4, "pmModelCopy() set the pmModel->fitRadius member correctly"); psFree(modelSrc); Index: branches/eam_branches/20090715/psModules/test/objects/tap_pmModelUtils.c =================================================================== --- branches/eam_branches/20090715/psModules/test/objects/tap_pmModelUtils.c (revision 25476) +++ branches/eam_branches/20090715/psModules/test/objects/tap_pmModelUtils.c (revision 25496) @@ -81,5 +81,5 @@ ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly"); ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly"); - ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly"); + ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly"); ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly"); ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly"); @@ -140,5 +140,5 @@ ok(tmpModel->nIter == testModelPSF->nIter, "pmModelFromPSF() set the model->nIter correctly"); ok(tmpModel->flags == testModelPSF->flags, "pmModelFromPSF() set the model->flags correctly"); - ok(tmpModel->radiusFit == testModelPSF->radiusFit, "pmModelFromPSF() set the model->radiusFit correctly"); + ok(tmpModel->fitRadius == testModelPSF->fitRadius, "pmModelFromPSF() set the model->fitRadius correctly"); ok(tmpModel->modelFunc == testModelPSF->modelFunc, "pmModelFromPSF() set the model->modelFunc correctly"); ok(tmpModel->modelFlux == testModelPSF->modelFlux, "pmModelFromPSF() set the model->modelFlux correctly"); Index: branches/eam_branches/20090715/psModules/test/objects/tap_pmSourcePhotometry.c =================================================================== --- branches/eam_branches/20090715/psModules/test/objects/tap_pmSourcePhotometry.c (revision 25476) +++ branches/eam_branches/20090715/psModules/test/objects/tap_pmSourcePhotometry.c (revision 25496) @@ -96,5 +96,5 @@ source->modelPSF = pmModelFromPSF (modelRef, psf); source->modelPSF->dparams->data.F32[PM_PAR_I0] = 1; - source->modelPSF->radiusFit = radius; + source->apRadius = radius; // measure photometry for centered source (fractional pix : 0.5,0.5)