Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 13834) +++ trunk/psphot/src/psphotChoosePSF.c (revision 13900) @@ -2,5 +2,5 @@ // try PSF models and select best option -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { +pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal, psMaskType mark) { bool status; @@ -43,6 +43,6 @@ psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md); if (!psfTrendMask) { - psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); - return NULL; + psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe"); + return NULL; } @@ -53,11 +53,11 @@ pmSource *source = sources->data[i]; if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - // keep NSTARS PSF stars, unmark the rest - if (stars->n < NSTARS) { - psArrayAdd (stars, 200, source); - } else { - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; - } - } + // keep NSTARS PSF stars, unmark the rest + if (stars->n < NSTARS) { + psArrayAdd (stars, 200, source); + } else { + source->mode &= ~PM_SOURCE_MODE_PSFSTAR; + } + } } psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n); @@ -65,6 +65,6 @@ if (stars->n == 0) { psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates"); - psFree (stars); - psFree (psfTrendMask); + psFree (stars); + psFree (psfTrendMask); return NULL; } @@ -91,5 +91,5 @@ psMetadataItem *item = psListGetAndIncrement (iter); char *modelName = item->data.V; - models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS); + models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark); } @@ -127,13 +127,13 @@ // print/dump psf parameters if (psTraceGetLevel("psphot") >= 5) { - for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) { - psPolynomial2D *poly = try->psf->params_NEW->data[i]; - for (int nx = 0; nx <= poly->nX; nx++) { - for (int ny = 0; ny <= poly->nY; ny++) { - if (poly->mask[nx][ny]) continue; - fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny); - } - } - } + for (int i = PM_PAR_SXX; i < try->psf->params_NEW->n; i++) { + psPolynomial2D *poly = try->psf->params_NEW->data[i]; + for (int nx = 0; nx <= poly->nX; nx++) { + for (int ny = 0; ny <= poly->nY; ny++) { + if (poly->mask[nx][ny]) continue; + fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny); + } + } + } } @@ -147,19 +147,19 @@ psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32); for (int i = 0; i < try->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) - continue; - - pmSource *source = try->sources->data[i]; - Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX]; - Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY]; - dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]; - dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX]; - dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY]; - Sx->n ++; - Sy->n ++; - dSN->n ++; - dSx->n ++; - dSy->n ++; + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = try->sources->data[i]; + Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX]; + Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY]; + dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]; + dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX]; + dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY]; + Sx->n ++; + Sy->n ++; + dSN->n ++; + dSx->n ++; + dSy->n ++; } psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); @@ -174,6 +174,6 @@ for (int i = 0; i < dSN->n; i++) { - dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]); - dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]); + dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]); + dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]); } @@ -206,79 +206,79 @@ // build a PSF residual image - if (!psphotMakeResiduals (try->sources, recipe, try->psf)) { - psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF"); - psFree (models); - return NULL; + if (!psphotMakeResiduals (try->sources, recipe, try->psf, maskVal)) { + psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF"); + psFree (models); + return NULL; } // XXX test dump of psf star data and psf-subtracted image - if (psTraceGetLevel("psphot.psfstars") > 5) { - psphotSaveImage (NULL, readout->image, "rawstars.fits"); - - for (int i = 0; i < try->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) - continue; - - pmSource *source = try->sources->data[i]; - float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - - // set the mask and subtract the PSF model - // XXX should we be using maskObj? should we be unsetting the mask? - // use pmModelSub because modelFlux has not been generated - assert (source->maskObj); - psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK); - pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL); - psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); - } - - FILE *f = fopen ("shapes.dat", "w"); - for (int i = 0; i < try->sources->n; i++) { - psF32 inPar[10]; - - // masked for: bad model fit, outlier in parameters - if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - - pmSource *source = try->sources->data[i]; - psF32 *outPar = source->modelEXT->params->data.F32; - - psEllipseShape shape; - - shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; - shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; - shape.sxy = outPar[PM_PAR_SXY]; - - psEllipsePol pol = pmPSF_ModelToFit (outPar); - inPar[PM_PAR_E0] = pol.e0; - inPar[PM_PAR_E1] = pol.e1; - inPar[PM_PAR_E2] = pol.e2; - pmPSF_FitToModel (inPar, 0.1); - - psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); - psEllipseAxes axes2; - (void)psEllipsePolToAxes(pol, 0.1, &axes2); - psEllipsePol pol2 = psEllipseAxesToPol (axes1); - - fprintf (f, "%3d %7.2f %7.2f %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n", - i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS], - outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], - pol.e0, pol.e1, pol.e2, - pol2.e0, pol2.e1, pol2.e2, - inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY], - axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD, - axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD - ); - } - fclose (f); - - psphotSaveImage (NULL, readout->image, "psfstars.fits"); - pmSourcesWritePSFs (try->sources, "psfstars.dat"); - pmSourcesWriteEXTs (try->sources, "extstars.dat", false); - psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); - psMetadataConfigWrite (psfData, "psfmodel.dat"); - psFree (psfData); - psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n"); - exit (0); + if (psTraceGetLevel("psphot.psfstars") > 5) { + psphotSaveImage (NULL, readout->image, "rawstars.fits"); + + for (int i = 0; i < try->sources->n; i++) { + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = try->sources->data[i]; + float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + + // set the mask and subtract the PSF model + // XXX should we be using maskObj? should we be unsetting the mask? + // use pmModelSub because modelFlux has not been generated + assert (source->maskObj); + psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK); + pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal); + psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + } + + FILE *f = fopen ("shapes.dat", "w"); + for (int i = 0; i < try->sources->n; i++) { + psF32 inPar[10]; + + // masked for: bad model fit, outlier in parameters + if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; + + pmSource *source = try->sources->data[i]; + psF32 *outPar = source->modelEXT->params->data.F32; + + psEllipseShape shape; + + shape.sx = outPar[PM_PAR_SXX] / M_SQRT2; + shape.sy = outPar[PM_PAR_SYY] / M_SQRT2; + shape.sxy = outPar[PM_PAR_SXY]; + + psEllipsePol pol = pmPSF_ModelToFit (outPar); + inPar[PM_PAR_E0] = pol.e0; + inPar[PM_PAR_E1] = pol.e1; + inPar[PM_PAR_E2] = pol.e2; + pmPSF_FitToModel (inPar, 0.1); + + psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0); + psEllipseAxes axes2; + (void)psEllipsePolToAxes(pol, 0.1, &axes2); + psEllipsePol pol2 = psEllipseAxesToPol (axes1); + + fprintf (f, "%3d %7.2f %7.2f %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %7.4f -- %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n", + i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS], + outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], + pol.e0, pol.e1, pol.e2, + pol2.e0, pol2.e1, pol2.e2, + inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY], + axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD, + axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD + ); + } + fclose (f); + + psphotSaveImage (NULL, readout->image, "psfstars.fits"); + pmSourcesWritePSFs (try->sources, "psfstars.dat"); + pmSourcesWriteEXTs (try->sources, "extstars.dat", false); + psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf); + psMetadataConfigWrite (psfData, "psfmodel.dat"); + psFree (psfData); + psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n"); + exit (0); } @@ -290,7 +290,7 @@ if (!psphotPSFstats (readout, recipe, psf)) { - psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); - psFree(psf); - return NULL; + psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); + psFree(psf); + return NULL; } @@ -327,7 +327,7 @@ pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre"); - psFree(modelEXT); - return false; + psError(PSPHOT_ERR_PSF, false, "Failed to estimate PSF model at image centre"); + psFree(modelEXT); + return false; } @@ -346,7 +346,7 @@ psF64 FWHM_Y = FWHM_X * (axes.minor / axes.major); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X", PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y", PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_X", PS_META_REPLACE, "PSF FWHM Major axis", FWHM_X); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_Y", PS_META_REPLACE, "PSF FWHM Minor axis", FWHM_Y); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars); psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true); @@ -386,5 +386,5 @@ moments.xy = source->moments->Sxy; - // limit axis ratio < 20.0 + // limit axis ratio < 20.0 axes = psEllipseMomentsToAxes (moments, 20.0);