Index: /branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c (revision 5806) +++ /branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c (revision 5807) @@ -1,6 +1,4 @@ # include -# include "psLibUtils.h" # include "pmObjects.h" -# include "psModulesUtils.h" # include "pmPSF.h" # include "pmPSFtry.h" @@ -41,5 +39,5 @@ type = pmModelSetType (modelName); test->psf = pmPSFAlloc (type); - test->sources = psMemCopy(sources); + test->sources = psMemIncrRefCounter(sources); test->modelFLT = psArrayAlloc (sources->n); test->modelPSF = psArrayAlloc (sources->n); @@ -79,129 +77,101 @@ 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_Alt (try - , 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 - ); -} - - -bool pmPSFtryMetric (pmPSFtry *try - , float RADIUS) -{ - - float dBin; - int nKeep, nSkip; - - // the measured (aperture - fit) magnitudes (dA == try->metric) + pmPSFtryMetric (psfTry, RADIUS); + psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", + modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); + + return (psfTry); +} + +bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS) +{ + + // 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,193 +183,28 @@ // 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]); - } - - // 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); - - // build binned versions of rflux, metric - dBin = (stats->max - stats->min) / 10.0; - psVector *rfBin = psVectorCreate(NULL, stats->min, stats->max, dBin, PS_TYPE_F64); - psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64); - psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8); - psFree (stats); - - psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin); - - // group data in daBin bins, measure lower 50% mean - for (int i = 0; i < daBin->n; i++) { - - psVector *tmp = psVectorAlloc (try - ->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 ++; - } - - // is this a valid point? - maskB->data.U8[i] = 0; - if (tmp->n < 2) { - maskB->data.U8[i] = 1; - psFree (tmp); + 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; - } - - // dA values are contaminated with low outliers - // measure statistics only on upper 50% of points - // this would be easier if we could sort in reverse: - // - // psVectorSort (tmp, tmp); - // tmp->n = 0.5*tmp->n; - // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - // psVectorStats (stats, tmp, NULL, NULL, 0); - // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); - - psVectorSort (tmp, tmp); - nKeep = 0.5*tmp->n; - nSkip = tmp->n - nKeep; - - psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64); - for (int j = 0; j < tmp2->n; j++) { - tmp2->data.F64[j] = tmp->data.F64[j + nSkip]; - } - - stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - psVectorStats (stats, tmp2, NULL, NULL, 0); - psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian); - - daBin->data.F64[i] = stats->sampleMedian; - - psFree (stats); - psFree (tmp); - psFree (tmp2); - } - - // linear fit to rfBin, daBin - psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD); - psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin); - - // XXX EAM : this is the intended API (cycle 7? cycle 8?) - // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin); - - // XXX EAM : replace this when the above version is implemented - // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL); - - psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin); - psVector *daResid = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit); - - stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV); - 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)); - - psFree (rflux); - psFree (rfBin); - psFree (daBin); - psFree (maskB); - psFree (daBinFit); - psFree (daResid); - psFree (poly); - psFree (stats); - psFree (fitstat); - - return true; -} - -bool pmPSFtryMetric_Alt (pmPSFtry *try - , float RADIUS) -{ - - // the measured (aperture - fit) magnitudes (dA == try->metric) - // depend on both the true ap-fit (dAo) and the bias in the sky measurement: - // dA = dAo + dsky/flux - // where flux is the flux of the star - // we fit this trend to find the infinite flux aperture correction (dAo), - // the nominal sky bias (dsky), and the error on dAo - // the values of dA are contaminated by stars with close neighbors in the aperture - // we use an outlier rejection to avoid this bias - - // 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]); - } - - // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit + rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); + } + + // use 3hi/1lo sigma clipping on the rflux vs metric fit psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - - // XXX EAM stats->min = 1.0; stats->max = 3.0; stats->clipIter = 3; - // linear clipped fit to rfBin, daBin + // linear clipped fit of ApResid to rflux psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); - poly = psVectorClipFitPolynomial1D (poly, stats, try - ->mask, PSFTRY_MASK_ALL, try - ->metric, NULL, rflux); - fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); - - try - ->psf->ApResid = poly->coeff[0]; - try - ->psf->dApResid = stats->sampleStdev; - try - ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); - - fprintf (stderr, "*******************************************************************************\n"); - - FILE *f; - f = fopen ("apresid.dat", "w"); - if (f == NULL) - psAbort ("pmPSFtry", "can't open output file"); - - for (int i = 0; i < try - ->sources->n; i++) { - fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try - ->fitMag->data.F64[i], rflux->data.F64[i], try - ->metric->data.F64[i], try - ->mask->data.U8[i]); - } - fclose (f); + poly = psVectorClipFitPolynomial1D (poly, stats, psfTry->mask, PSFTRY_MASK_ALL, psfTry->metric, NULL, rflux); + psLogMsg ("pmPSFtryMetric", 4, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + + psfTry->psf->ApTrend->coeff[0][0][0] = poly->coeff[0]; + psfTry->psf->ApTrend->coeff[0][0][1] = poly->coeff[1]; + + psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + psfTry->psf->ApResid = poly->coeff[0]; + psfTry->psf->dApResid = stats->sampleStdev; psFree (rflux); @@ -407,7 +212,4 @@ psFree (stats); - // psFree (daFit); - // psFree (daResid); - return true; }