Index: trunk/psModules/src/objects/pmPSFtry.c =================================================================== --- trunk/psModules/src/objects/pmPSFtry.c (revision 5255) +++ trunk/psModules/src/objects/pmPSFtry.c (revision 5762) @@ -79,123 +79,102 @@ int Npsf = 0; - pmPSFtry *try - = pmPSFtryAlloc (sources, modelName); + pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName); // stage 1: fit an independent model (freeModel) to all sources psTimerStart ("fit"); - for (int i = 0; i < try - ->sources->n; i++) { - - pmSource *source = try - ->sources->data[i]; - pmModel *model = pmSourceModelGuess (source, try - ->psf->type); - x = source->peak->x; - y = source->peak->y; - - // set temporary object mask and fit object - // fit model as FLT, not PSF - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - // exclude the poor fits - if (!status) { - try - ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; - psFree (model); - continue; - } - try - ->modelFLT->data[i] = model; - Nflt ++; - } + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + pmModel *model = pmSourceModelGuess (source, psfTry->psf->type); + x = source->peak->x; + y = source->peak->y; + + // set temporary object mask and fit object + // fit model as FLT, not PSF + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, model, false); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + // exclude the poor fits + if (!status) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; + psFree (model); + continue; + } + psfTry->modelFLT->data[i] = model; + Nflt ++; + } psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); // make this optional? - // DumpModelFits (try->modelFLT, "modelsFLT.dat"); + // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat"); // stage 2: construct a psf (pmPSF) from this collection of model fits - pmPSFFromModels (try - ->psf, try - ->modelFLT, try - ->mask); + pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask); // stage 3: refit with fixed shape parameters psTimerStart ("fit"); - 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]; - pmModel *modelFLT = try - ->modelFLT->data[i]; - - // set shape for this model based on PSF - pmModel *modelPSF = pmModelFromPSF (modelFLT, try - ->psf); - x = source->peak->x; - y = source->peak->y; - - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, modelPSF, true); - - // skip poor fits - if (!status) { - try - ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; - psFree (modelPSF); - goto next_source; - } - - // otherwise, save the resulting model - try - ->modelPSF->data[i] = modelPSF; - - // XXX : use a different aperture radius from the fit radius? - // XXX : use a different estimator for the local sky? - // XXX : pass 'source' as input? - if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { - try - ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; - goto next_source; - } - - try - ->metric->data.F64[i] = obsMag - fitMag; - try - ->fitMag->data.F64[i] = fitMag; - Npsf ++; + for (int i = 0; i < psfTry->sources->n; i++) { + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = psfTry->sources->data[i]; + pmModel *modelFLT = psfTry->modelFLT->data[i]; + + // set shape for this model based on PSF + pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf); + x = source->peak->x; + y = source->peak->y; + + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, modelPSF, true); + + // skip poor fits + if (!status) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; + psFree (modelPSF); + goto next_source; + } + + // otherwise, save the resulting model + psfTry->modelPSF->data[i] = modelPSF; + + // XXX : use a different aperture radius from the fit radius? + // XXX : use a different estimator for the local sky? + // XXX : pass 'source' as input? + if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; + goto next_source; + } + + psfTry->metric->data.F64[i] = obsMag - fitMag; + psfTry->fitMag->data.F64[i] = fitMag; + Npsf ++; next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - } + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + } psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); // make this optional - // DumpModelFits (try->modelPSF, "modelsPSF.dat"); + // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat"); // XXX this function wants aperture radius for pmSourcePhotometry - pmPSFtryMetric (try - , RADIUS); + pmPSFtryMetric (psfTry, RADIUS); psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", - modelName, try - ->psf->ApResid, try - ->psf->dApResid, try - ->psf->skyBias); - - return (try - ); + modelName, + psfTry->psf->ApResid, + psfTry->psf->dApResid, + psfTry->psf->skyBias); + + return (psfTry); } -bool pmPSFtryMetric (pmPSFtry *try - , float RADIUS) +bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS) { @@ -203,5 +182,5 @@ int nKeep, nSkip; - // the measured (aperture - fit) magnitudes (dA == try->metric) + // the measured (aperture - fit) magnitudes (dA == psfTry->metric) // depend on both the true ap-fit (dAo) and the bias in the sky measurement: // dA = dAo + dsky/flux @@ -213,18 +192,14 @@ // rflux = ten(0.4*fitMag); - psVector *rflux = psVectorAlloc (try - ->sources->n, PS_TYPE_F64); - for (int i = 0; i < try - ->sources->n; i++) { - if (try - ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - rflux->data.F64[i] = pow(10.0, 0.4*try - ->fitMag->data.F64[i]); - } + psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); + } // find min and max of (1/flux): psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - psVectorStats (stats, rflux, NULL, try - ->mask, PSFTRY_MASK_ALL); + psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL); // build binned versions of rflux, metric @@ -240,29 +215,26 @@ for (int i = 0; i < daBin->n; i++) { - psVector *tmp = psVectorAlloc (try - ->sources->n, PS_TYPE_F64); + psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); tmp->n = 0; // accumulate data within bin range - for (int j = 0; j < try - ->sources->n; j++) { - // masked for: bad model fit, outlier in parameters - if (try - ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue; - - // skip points with extreme dA values - if (fabs(try - ->metric->data.F64[j]) > 0.5) continue; - - // skip points outside of this bin - if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) - continue; - if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) - continue; - - tmp->data.F64[tmp->n] = try - ->metric->data.F64[j]; - tmp->n ++; - } + for (int j = 0; j < psfTry->sources->n; j++) { + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL) + continue; + + // skip points with extreme dA values + if (fabs(psfTry->metric->data.F64[j]) > 0.5) + continue; + + // skip points outside of this bin + if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) + continue; + if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) + continue; + + tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j]; + tmp->n ++; + } // is this a valid point? @@ -319,10 +291,7 @@ stats = psVectorStats (stats, daResid, NULL, maskB, 1); - try - ->psf->ApResid = poly->coeff[0]; - try - ->psf->dApResid = stats->clippedStdev; - try - ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + psfTry->psf->ApResid = poly->coeff[0]; + psfTry->psf->dApResid = stats->clippedStdev; + psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); psFree (rflux);