Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5979) +++ /trunk/psphot/Makefile (revision 5980) @@ -58,4 +58,5 @@ $(SRC)/psphotBlendFit.$(ARCH).o \ $(SRC)/psphotApResid.$(ARCH).o \ +$(SRC)/psphotAssessPSF.$(ARCH).o \ $(SRC)/psBicube.$(ARCH).o @@ -82,9 +83,5 @@ $(SRC)/psphotTest.$(ARCH).o \ $(SRC)/psphotTestArguments.$(ARCH).o \ -$(SRC)/pmSourceContour.$(ARCH).o \ -$(SRC)/psImageData.$(ARCH).o \ -$(SRC)/psphotSetup.$(ARCH).o \ -$(SRC)/psModulesUtils.$(ARCH).o \ -$(SRC)/psSparse.$(ARCH).o +$(SRC)/psphotApResid.$(ARCH).o psphot: $(BIN)/psphot.$(ARCH) @@ -102,5 +99,5 @@ $(TEST): $(SRC)/psphot.h -INSTALL = psphot psphotTest psphotModelTest +INSTALL = psphot # dependancy rules for binary code ######################### Index: /trunk/psphot/doc/notes.txt =================================================================== --- /trunk/psphot/doc/notes.txt (revision 5979) +++ /trunk/psphot/doc/notes.txt (revision 5980) @@ -1,4 +1,2 @@ - -Notes on psphot 2005.12.23 Index: /trunk/psphot/doc/speed.txt =================================================================== --- /trunk/psphot/doc/speed.txt (revision 5980) +++ /trunk/psphot/doc/speed.txt (revision 5980) @@ -0,0 +1,129 @@ + +psphot -g +psModule -g +psLib -O2 +on alala: + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 0.983568 sec + System random seed value used seed = 25904F7C7454AE25 hex + background: 1229.280029 +/- 27.700180 + image stats: 0.063085 sec + System random seed value used seed = 411E5A0600A49AF3 hex + fit background model: 0.580152 sec + smooth: 1.842415 sec + threshold: 83.100540 DN + 4356 peaks: 0.571587 sec + 4230 moments: 3.965099 sec + SN range: 3.355829 - 1043.401001 + identified 107 blended objects (0.121737 sec) + selected candidate 300 PSF objects + fit flt: 3.682850 sec for 300 sources + fit psf: 1.704141 sec for 300 sources + fit stats: 0.000884 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717999 + select psf model: 5.402742 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.991345 (4230 objects) + built matrix: 1.472869 (27998 elements) + measure ensemble of PSFs: 4.334556 + fit PSF models: 12.205738 sec for 4230 objects + replace unfitted models: 0.021473 sec (4230 objects) + measure aperture residuals : 0.451886 sec + measure full-frame aperture residual: 0.454032 sec + aperture residual: -0.002842 +/- 0.002729 : 0.166072 bias + wrote output: 4.648821 sec + complete psphot run: 35.334080 sec + + +psphot -g +psModule -O2 +psLib -O2 +on alala: + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 1.001502 sec + System random seed value used seed = 5F43A63F532D4B3F hex + background: 1229.280029 +/- 27.700180 + image stats: 0.063169 sec + System random seed value used seed = F5F5FF99A6578F61 hex + fit background model: 0.530384 sec + smooth: 1.856414 sec + threshold: 83.100540 DN + 4365 peaks: 0.237292 sec + 4239 moments: 4.087611 sec + SN range: 3.356469 - 1043.400879 + identified 107 blended objects (0.123372 sec) + selected candidate 300 PSF objects + fit flt: 3.741444 sec for 300 sources + fit psf: 1.676478 sec for 300 sources + fit stats: 0.000883 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717985 + select psf model: 5.458156 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.791613 (4239 objects) + built matrix: 1.267909 (28101 elements) + measure ensemble of PSFs: 3.551364 + fit PSF models: 11.240197 sec for 4239 objects + replace unfitted models: 0.018051 sec (4239 objects) + measure aperture residuals : 0.411856 sec + measure full-frame aperture residual: 0.413843 sec + aperture residual: -0.002831 +/- 0.002717 : 0.164772 bias + wrote output: 4.355429 sec + complete psphot run: 33.126283 sec + +after reducing indirection in pmSourceMoments: + 4268 moments: 3.762452 sec + +psphot -O2 +psModule -O2 +psLib -O2 + +alala: psphot -vv starfield.fits starfield starfield.cnf -resid starfield.resid.3.fits + load data: 0.484324 sec + System random seed value used seed = 980050732A392EC7 hex + background: 1229.280029 +/- 27.700180 + image stats: 0.058759 sec + System random seed value used seed = 44E1A6EE77B0E9CE hex + fit background model: 0.395125 sec + smooth: 1.861776 sec + threshold: 83.100540 DN + 4351 peaks: 0.234874 sec + 4264 moments: 3.719435 sec + SN range: 2.023345 - 1043.401123 + identified 124 blended objects (0.115230 sec) + selected candidate 300 PSF objects + fit flt: 3.525501 sec for 300 sources + fit psf: 1.595031 sec for 300 sources + fit stats: 0.000884 +/- 0.004499 + try model PS_MODEL_QGAUSS, ap-fit: -0.007532 +/- 0.004499 : sky bias: -0.717976 + select psf model: 5.160442 sec + selected psf model PS_MODEL_QGAUSS, ApResid: -0.007532 +/- 0.004499 + fitting pixels with at least 27.700180 object counts + built models: 0.760805 (4264 objects) + built matrix: 0.908674 (23703 elements) + measure ensemble of PSFs: 3.115425 + fit PSF models: 10.162397 sec for 4264 objects + replace unfitted models: 0.017491 sec (4264 objects) + measure aperture residuals : 0.408242 sec + measure full-frame aperture residual: 0.410207 sec + aperture residual: -0.002908 +/- 0.002734 : 0.170692 bias + wrote output: 4.225560 sec + complete psphot run: 30.149281 sec + + + fit PSF models: 10.162397 sec for 4264 objects + select psf model: 5.160442 sec + wrote output: 4.225560 sec + 4264 moments: 3.719435 sec + measure ensemble of PSFs: 3.115425 sec + + other : 3.5 sec + + fixed per objects : 4 msec / object (17 of 30 sec) + + +** still can re-compile with -DPS_NO_TRACE +** ASSERTS cannot be compiled out at the moment! Index: /trunk/psphot/src/psLine.c =================================================================== --- /trunk/psphot/src/psLine.c (revision 5979) +++ /trunk/psphot/src/psLine.c (revision 5980) @@ -41,4 +41,5 @@ Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap); line->Nline += PS_MIN (nMax - 1, Nchar); + va_end (ap); if (Nchar >= nMax) return (false); Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5979) +++ /trunk/psphot/src/psphot.c (revision 5980) @@ -78,5 +78,5 @@ psphotFullFit (imdata, config, sources, psf, sky); psphotReplaceUnfit (sources); - psphotApResid (sources, config, psf); + psphotApResid (imdata, sources, config, psf); break; @@ -85,5 +85,5 @@ psphotBlendFit (imdata, config, sources, psf, sky); psphotReplaceUnfit (sources); - psphotApResid (sources, config, psf); + psphotApResid (imdata, sources, config, psf); break; Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5979) +++ /trunk/psphot/src/psphot.h (revision 5980) @@ -4,4 +4,5 @@ # include # include +# include # include # include @@ -37,4 +38,5 @@ int psphotSaveImage (psMetadata *header, psImage *image, char *filename); bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); +bool psphotRedefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); void psphotModelGroupInit (void); @@ -84,5 +86,7 @@ bool psphotReplaceUnfit (psArray *sources); bool psphotDumpMoments (psMetadata *config, psArray *sources); -bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf); +bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf); +bool psphotAssessPSF (eamReadout *imdata, psMetadata *config, pmPSF *psf); + bool psphotWritePSF (pmPSF *psf, char *filename); pmPSF *psphotReadPSF (char *filename); @@ -90,6 +94,8 @@ bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...); bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...); +bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...); psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder); psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder); +psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder); bool psphotModelTest (eamReadout *imdata, psMetadata *config); @@ -112,2 +118,21 @@ psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); psPlane psImageBicubeMin (psPolynomial2D *poly); +void psphotTestArguments (int *argc, char **argv); + +bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf); + +psPolynomial4D *psVectorChiClipFitPolynomial4D( + psPolynomial4D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t); + +// XXX deprecated +// bool psphotApResidFullFit (psVector *x, psVector *y, psVector *radius, psVector *rflux, psVector *apresid, psVector *dmag); +// bool psphotApResidPolyFit (psVector *x, psVector *y, psVector *radius, psVector *rflux, psVector *apresid, psVector *dmag); Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 5979) +++ /trunk/psphot/src/psphotApResid.c (revision 5980) @@ -1,5 +1,114 @@ # include "psphot.h" -bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf) { +psPolynomial4D *psVectorChiClipFitPolynomial4D( + psPolynomial4D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_PTR_NON_NULL(stats, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL); + PS_ASSERT_VECTOR_NON_NULL(fErr, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + + // clipping range defined by min and max and/or clipSigma + float minClipSigma; + float maxClipSigma; + if (isfinite(stats->max)) { + maxClipSigma = +fabs(stats->max); + } else { + maxClipSigma = +fabs(stats->clipSigma); + } + if (isfinite(stats->min)) { + minClipSigma = -fabs(stats->min); + } else { + minClipSigma = -fabs(stats->clipSigma); + } + psVector *fit = NULL; + psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); + + // eventual expansion: user supplies one of various stats option pairs, + // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to + // evaluate the clipping sigma + // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used + stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + for (int N = 0; N < stats->clipIter; N++) { + int Nkeep = 0; + + poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); + fit = psPolynomial4DEvalVector (poly, x, y, z, t); + resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); + + stats = psVectorStats (stats, resid, NULL, mask, maskValue); + psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); + + // set mask if pts are not valid + // we are masking out any point which is out of range + // recovery is not allowed with this scheme + for (int i = 0; i < resid->n; i++) { + if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { + continue; + } + float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev); + if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + Nkeep ++; + } + + psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", + Nkeep, x->n); + + stats->clippedNvalues = Nkeep; + psFree (fit); + } + // Free local temporary variables + psFree (resid); + + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); + return(NULL); + } + return(poly); +} + +// measure the aperture residual statistics +bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf) { int Npsf; @@ -8,22 +117,23 @@ pmSource *source; - // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); - psTimerStart ("psphot"); - // XXX EAM : drop this if we know the list is sorted - sources = psArraySort (sources, psphotSortBySN); - + // measure the aperture loss as a function of radius for PSF + float REF_RADIUS = psMetadataLookupF32 (&status, config, "PSF_REF_RADIUS"); + psf->growth = pmGrowthCurveAlloc (3.0, REF_RADIUS, 0.1); + psphotGrowthCurve (imdata, psf); + psVector *mask = psVectorAlloc (300, PS_TYPE_U8); psVector *xPos = psVectorAlloc (300, PS_TYPE_F64); psVector *yPos = psVectorAlloc (300, PS_TYPE_F64); - psVector *rflux = psVectorAlloc (300, PS_TYPE_F64); + psVector *flux = psVectorAlloc (300, PS_TYPE_F64); + psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64); psVector *apResid = psVectorAlloc (300, PS_TYPE_F64); - mask->n = xPos->n = yPos->n = rflux->n = apResid->n = 0; + psVector *dMag = psVectorAlloc (300, PS_TYPE_F64); + mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0; Npsf = 0; - // select the NNN brightest, non-saturated sources, or just select PSFSTARs? - for (int i = 0; (i < sources->n) && (Npsf < 300); i++) { + // select all good PM_SOURCE_STAR entries + for (int i = 0; i < sources->n; i++) { source = sources->data[i]; @@ -34,6 +144,6 @@ if (source->mode & PM_SOURCE_POOR) continue; - // get magnitudes, uncorrected for (x, y, rflux) - model = pmSourceMagnitudes (source, NULL, RADIUS); + // get uncorrected magnitudes in scaled apertures + model = pmSourceMagnitudes (source, NULL, 0); if (model == NULL) continue; @@ -41,45 +151,199 @@ xPos->data.F64[Npsf] = model->params->data.F32[2]; yPos->data.F64[Npsf] = model->params->data.F32[3]; - rflux->data.F64[Npsf] = pow(10.0, 0.4*source->fitMag); - apResid->data.F64[Npsf] = source->apMag - source->fitMag; - - psVectorExtend (mask, 100, 1); - psVectorExtend (xPos, 100, 1); - psVectorExtend (yPos, 100, 1); - psVectorExtend (rflux, 100, 1); + + flux->data.F64[Npsf] = pow(10.0, -0.4*source->fitMag); + r2rflux->data.F64[Npsf] = PS_SQR(model->radius) / flux->data.F64[Npsf]; + + apResid->data.F64[Npsf] = source->apMag + pmGrowthCurveCorrect (psf->growth, model->radius) - source->fitMag ; + + // XXX sanity clip? + // XXX need to see if all data were tossed? + // XXX need to subtract median? + if (fabs(apResid->data.F64[Npsf]) > 0.2) continue; + + dMag->data.F64[Npsf] = model->dparams->data.F32[1] / model->params->data.F32[1]; + + psVectorExtend (mask, 100, 1); + psVectorExtend (xPos, 100, 1); + psVectorExtend (yPos, 100, 1); + psVectorExtend (flux, 100, 1); + psVectorExtend (r2rflux, 100, 1); + psVectorExtend (dMag, 100, 1); psVectorExtend (apResid, 100, 1); Npsf ++; } - psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects\n", psTimerMark ("psphot"), Npsf); + + // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL + char *ApTrendOption = psMetadataLookupPtr (&status, config, "APTREND"); + if (!status) ApTrendOption = psStringCopy ("SKYBIAS"); // 3hi/1lo sigma clipping on the rflux vs metric fit psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->min = 1.0; + stats->min = 3.0; stats->max = 3.0; - stats->clipIter = 3; - - // first clip out objects which are too far from the median - stats->clipIter = 1; - maskToConstant (psf->ApTrend); - psf->ApTrend = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux); - - // next, fit just SkyBias and clip out objects which are too far from the median - stats->clipIter = 2; - maskToSkyBias (psf->ApTrend); - psf->ApTrend = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux); - - // finally, fit x, y, SkyBias and clip out objects which are too far from the median - stats->clipIter = 2; - maskToDefault (psf->ApTrend); - psf->ApTrend = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux); - - // linear clipped fit of apResid to rflux, xPos, yPos - # if (0) - psf->ApTrend = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux); - psf->skyBias = psf->ApTrend->coeff[0][0][1] / (M_PI * PS_SQR(RADIUS)); - psf->ApResid = psf->ApTrend->coeff[0][0][0]; - psf->dApResid = stats->sampleStdev; - psf->ApTrend->coeff[0][0][1] = 0; - # endif + + // constant only + if (!strcasecmp (ApTrendOption, "CONSTANT")) { + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and skybias only + if (!strcasecmp (ApTrendOption, "SKYBIAS")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + if (!strcasecmp (ApTrendOption, "SKYSAT")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and linear X,Y only + if (!strcasecmp (ApTrendOption, "XY_LIN")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_XY_LIN); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and quadratic X,Y only + if (!strcasecmp (ApTrendOption, "XY_QUAD")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_XY_QUAD); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and sky, linear X,Y only + if (!strcasecmp (ApTrendOption, "SKY_XY_LIN")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKY_XY_LIN); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + // constant and sky, linear X,Y only + if (!strcasecmp (ApTrendOption, "SKYSAT_XY_LIN")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // apply the fit + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT_XY_LIN); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + + if (!strcasecmp (ApTrendOption, "ALL")) { + // first clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT); + psf->ApTrend = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // fit just SkyBias and clip out objects which are too far from the median + stats->clipIter = 2; + pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + + // finally, fit x, y, SkyBias and clip out objects which are too far from the median + stats->clipIter = 3; + pmPSF_MaskApTrend (psf, PM_PSF_ALL); + psf->ApTrend = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux); + } + +# if (1) + psPolynomial4D *poly = psf->ApTrend; + for (int nt = 0; nt <= poly->nT; nt++) { + for (int nz = 0; nz <= poly->nZ; nz++) { + for (int ny = 0; ny <= poly->nY; ny++) { + for (int nx = 0; nx <= poly->nX; nx++) { + if (poly->mask[nx][ny][nz][nt]) continue; + fprintf (stderr, "%d %d %d %d : %22.15g\n", nx, ny, nz, nt, poly->coeff[nx][ny][nz][nt]); + } + } + } + } +# endif + + // construct the fitted values and the residuals + psVector *fit = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux); + psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit); + +# if (0) + FILE *fout = fopen ("resid.dat", "w"); + for (int i = 0; i < resid->n; i++) { + fprintf (fout, "%d %f %f %f %f %f %f %f %f %d\n", + i, + (float) psVectorGet(xPos, i), (float) psVectorGet(yPos, i), (float) psVectorGet(r2rflux, i), (float) psVectorGet(flux, i), + (float) psVectorGet(apResid, i), (float) psVectorGet(fit, i), (float) psVectorGet(resid, i), (float) psVectorGet(dMag, i), + mask->data.U8[i]); + } + fclose (fout); +# endif + + // measure scatter for sources with dMag < 0.01 (S/N = 100) + psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); + for (int i = 0; i < dMag->n; i++) { + if (dMag->data.F64[i] > 0.01) { + mask->data.U8[i] |= 0x02; + } + } + residStats = psVectorStats (residStats, resid, NULL, mask, 0x03); + + // apply ApTrend results + psf->skyBias = psf->ApTrend->coeff[0][0][1][0]; + psf->skySat = psf->ApTrend->coeff[0][0][0][1]; + psf->ApResid = psf->ApTrend->coeff[0][0][0][0]; + psf->dApResid = residStats->sampleStdev; + psf->ApTrend->coeff[0][0][1][0] = 0; + psf->ApTrend->coeff[0][0][0][1] = 0; + psf->nApResid = residStats->clippedNvalues; /* @@ -89,23 +353,19 @@ */ - # if (0) - psPolynomial3D *poly = psf->ApTrend; - for (int nz = 0; nz <= poly->nZ; nz++) { - for (int ny = 0; ny <= poly->nY; ny++) { - for (int nx = 0; nx <= poly->nX; nx++) { - fprintf (stderr, "%d %d %d : %22.15g\n", nx, ny, nz, poly->coeff[nx][ny][nz]); - } - } - } - # endif - psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot")); psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias); - - psFree (stats); - psFree (mask); - psFree (rflux); - psFree (apResid); + psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", + 1e3*psf->ApTrend->coeff[1][0][0][0], + 1e6*psf->ApTrend->coeff[2][0][0][0], + 1e6*psf->ApTrend->coeff[1][1][0][0], + 1e3*psf->ApTrend->coeff[0][1][0][0], + 1e6*psf->ApTrend->coeff[0][2][0][0]); + + // psFree (stats); + // psFree (mask); + // psFree (rflux); + // psFree (apResid); return true; } + Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 5979) +++ /trunk/psphot/src/psphotArguments.c (revision 5980) @@ -7,5 +7,6 @@ int N; unsigned int Nfail; - int mode = PS_DATA_STRING | PS_META_REPLACE; + + psMetadata *options = psMetadataAlloc (); // basic pslib options @@ -13,81 +14,100 @@ psArgumentVerbosity (argc, argv); - // optional mask image - add to config - char *mask = NULL; + // command-line options -- place on options MD + // mask image if ((N = psArgumentGet (*argc, argv, "-mask"))) { psArgumentRemove (N, argc, argv); - mask = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "MASK_IMAGE", 0, "", psStringCopy(argv[N])); psArgumentRemove (N, argc, argv); } - // optional weight image - add to config - char *weight = NULL; + // weight image if ((N = psArgumentGet (*argc, argv, "-weight"))) { psArgumentRemove (N, argc, argv); - weight = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "WEIGHT_IMAGE", 0, "", psStringCopy(argv[N])); psArgumentRemove (N, argc, argv); } - // optional output residual image - add to config - char *resid = NULL; + // output residual image if ((N = psArgumentGet (*argc, argv, "-resid"))) { psArgumentRemove (N, argc, argv); - resid = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "RESID_IMAGE", 0, "", psStringCopy(argv[N])); psArgumentRemove (N, argc, argv); } - // optional analysis region - add to config - char *region = NULL; + // analysis region if ((N = psArgumentGet (*argc, argv, "-region"))) { psArgumentRemove (N, argc, argv); - region = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", psStringCopy(argv[N])); psArgumentRemove (N, argc, argv); } - // optional output residual image - add to config - char *photcode = NULL; + // output residual image + psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", 0, "", psStringCopy("NONE")); if ((N = psArgumentGet (*argc, argv, "-photcode"))) { psArgumentRemove (N, argc, argv); - photcode = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", psStringCopy(argv[N])); psArgumentRemove (N, argc, argv); } - char *psffile = NULL; + // input PSF file if ((N = psArgumentGet (*argc, argv, "-psf"))) { psArgumentRemove (N, argc, argv); - psffile = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "PSF_INPUT_FILE", 0, "", psStringCopy (argv[N])); psArgumentRemove (N, argc, argv); } - bool ModelTest = false; - float ModelTest_X = 0; - float ModelTest_Y = 0; - char *model = NULL; - char *fitset = NULL; - char *fitmode = NULL; + // run the 0ltest model (requires X,Y coordinate) if ((N = psArgumentGet (*argc, argv, "-modeltest"))) { - ModelTest = true; + psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT", 0, "", true); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); + psArgumentRemove (N, argc, argv); - ModelTest_X = atof (argv[N]); psArgumentRemove (N, argc, argv); - ModelTest_Y = atof (argv[N]); psArgumentRemove (N, argc, argv); + // specify the modeltest model if ((N = psArgumentGet (*argc, argv, "-model"))) { psArgumentRemove (N, argc, argv); - model = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", psStringCopy (argv[N])); psArgumentRemove (N, argc, argv); } + // specify the test fit mode if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { psArgumentRemove (N, argc, argv); - fitmode = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", psStringCopy (argv[N])); psArgumentRemove (N, argc, argv); } if ((N = psArgumentGet (*argc, argv, "-fitset"))) { psArgumentRemove (N, argc, argv); - fitset = psStringCopy (argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", psStringCopy (argv[N])); psArgumentRemove (N, argc, argv); } + } + + // other arbitrary options: -D key value (all added as string) + while ((N = psArgumentGet (*argc, argv, "-D"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + // other arbitrary options: -Df key value (all added as float) + while ((N = psArgumentGet (*argc, argv, "-Df"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + // other arbitrary options: -Di key value (all added as int) + while ((N = psArgumentGet (*argc, argv, "-Di"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); } @@ -96,52 +116,21 @@ // load config information psMetadata *config = psMetadataAlloc (); - psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); + + psMetadataAdd (config, PS_LIST_TAIL, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE); + + psMetadataItem *item = NULL; + psMetadataIterator *iter = psMetadataIteratorAlloc (options, PS_LIST_HEAD, NULL); + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + psMetadataAddItem (config, item, PS_LIST_TAIL, PS_META_REPLACE); + // psFree (item); + } + // psFree (iter); // identify input image & optional weight & mask images // command-line entries override config-file entries - psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); - psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]); + psMetadataAddStr (config, PS_LIST_TAIL, "IMAGE", PS_META_REPLACE, "", argv[1]); + psMetadataAddStr (config, PS_LIST_TAIL, "OUTPUT_FILE", PS_META_REPLACE, "", argv[2]); - if (mask != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask); - } - if (weight != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "WEIGHT_IMAGE", mode, "", weight); - } - if (resid != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid); - } - if (region != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "ANALYSIS_REGION", mode, "", region); - } - if (photcode != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", photcode); - } else { - psMetadataAdd (config, PS_LIST_HEAD, "PHOTCODE", mode, "", "NONE"); - } - if (psffile != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "PSF_INPUT_FILE", mode, "", psffile); - } - - // model related options - if (ModelTest) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT", PS_DATA_BOOL, "", true); - - int fmode = PS_DATA_F32 | PS_META_REPLACE; - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", fmode, "", ModelTest_X); - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", fmode, "", ModelTest_Y); - - if (model != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODEL", mode, "", model); - } - - if (fitmode != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_MODE", mode, "", fitmode); - } - if (fitset != NULL) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_SET", mode, "", fitset); - } - } return (config); } Index: /trunk/psphot/src/psphotDefinePixels.c =================================================================== --- /trunk/psphot/src/psphotDefinePixels.c (revision 5979) +++ /trunk/psphot/src/psphotDefinePixels.c (revision 5980) @@ -16,6 +16,52 @@ mySource->weight = psImageSubset(imdata->weight, srcRegion); mySource->mask = psImageSubset(imdata->mask, srcRegion); + mySource->region = srcRegion; - return(mySource); + return true; } +bool psphotRedefinePixels(pmSource *mySource, + const eamReadout *imdata, + psF32 x, + psF32 y, + psF32 Radius) +{ + bool extend; + psRegion newRegion; + + if (Radius == 0) return false; + + // check to see if new region is completely contained within old region + newRegion = psRegionForSquare (x, y, Radius); + newRegion = psRegionForImage (imdata->image, newRegion); + + extend = false; + extend |= (int)(newRegion.x0) < (int)(mySource->region.x0); + extend |= (int)(newRegion.x1) > (int)(mySource->region.x1); + extend |= (int)(newRegion.y0) < (int)(mySource->region.y0); + extend |= (int)(newRegion.y1) > (int)(mySource->region.y1); + + extend |= (mySource->pixels == NULL); + extend |= (mySource->weight == NULL); + extend |= (mySource->mask == NULL); + + // extend = true; + if (extend) { + // Grab a new subimage + // psFree (mySource->pixels); + // psFree (mySource->weight); + // psFree (mySource->mask); + + // fprintf (stderr, "extend %f,%f\n", x, y); + mySource->pixels = psImageSubset(imdata->image, newRegion); + mySource->weight = psImageSubset(imdata->weight, newRegion); + mySource->mask = psImageSubset(imdata->mask, newRegion); + mySource->region = newRegion; + } + + return extend; +} + + +//**** this function was ALWAYS restricting to model->radius +// should have left it larger (PSF_FIT_RAD) for the faint sources Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 5979) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 5980) @@ -58,5 +58,4 @@ sources = psArraySort (sources, psphotSortByY); - float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); float INNER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS"); float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); @@ -104,5 +103,4 @@ // really saturated stars should be re-measured for a better centroid // XXX EAM : move this to a 'clear satstar function' - // XXX EAM : extend size of fit box around SATSTAR if (inSource->mode & PM_SOURCE_SATSTAR) { status = pmSourceMoments (inSource, INNER_RADIUS); @@ -143,5 +141,5 @@ y = model->params->data.F32[3]; - // get function which specifies the radius at this the model hits the given flux + // get function which specifies the radius at which the model hits the given flux pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); @@ -151,7 +149,5 @@ // if needed, ask for more object pixels - if (model->radius > OUTER_RADIUS) { - psphotDefinePixels (inSource, imdata, x, y, model->radius); - } + psphotRedefinePixels (inSource, imdata, x, y, model->radius); // make temporary copies of the image pixels and mask Index: /trunk/psphot/src/psphotFixedPSF.c =================================================================== --- /trunk/psphot/src/psphotFixedPSF.c (revision 5979) +++ /trunk/psphot/src/psphotFixedPSF.c (revision 5980) @@ -13,6 +13,4 @@ // we may set this differently here from the value used to mark likely saturated stars - float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); - float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); @@ -52,10 +50,5 @@ } - // check if we need to redefine the pixels - // XXX EAM : a better test would examine the source pixels - if (model->radius > OUTER_RADIUS) { - // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) - psphotDefinePixels (source, imdata, x, y, model->radius); - } + psphotRedefinePixels (source, imdata, x, y, model->radius); // fit PSF model (set/unset the pixel mask) Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 5979) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 5980) @@ -17,5 +17,5 @@ model = source->modelPSF; if (model == NULL) return NULL; - radius = apRadius; + radius = (apRadius > 0) ? apRadius : model->radius; isPSF = true; break; @@ -46,14 +46,18 @@ // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS if (isPSF && (psf != NULL)) { + if (psf->growth != NULL) { + source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius); + } + rflux = pow (10.0, 0.4*source->fitMag); - source->apMag -= rflux * psf->skyBias * (M_PI * PS_SQR(apRadius)); - source->fitMag += psPolynomial3DEval (psf->ApTrend, x, y, 0.0); + source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux; + source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); } + + // unmask aperture + psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); // subtract object, leave local sky pmSourceSubModel (source->pixels, source->mask, model, false, false); - - // unmask aperture - psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); if (!status) return NULL; @@ -62,6 +66,7 @@ /* -aprMag' - fitMag = rflux*skyBias + ApTrend(x,y) -(aprMag - rflux*skyBias) - fitMAg = ApTrend(x,y) -(aprMag - rflux*skyBias) = fitMAg + ApTrend(x,y) +aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) + */ Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5979) +++ /trunk/psphot/src/psphotOutput.c (revision 5980) @@ -65,5 +65,5 @@ } - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); // write sources with models @@ -116,5 +116,5 @@ } - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); // write sources with models @@ -162,5 +162,5 @@ // find config information for output header - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); if (!status) ZERO_POINT = 25.0; @@ -171,6 +171,9 @@ psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT", PS_DATA_F32 | PS_META_REPLACE, "zero point", ZERO_POINT); psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nPSFstars); psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skySat); psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code", PHOTCODE); psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); @@ -250,5 +253,5 @@ // find config information for output header - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE"); @@ -367,5 +370,5 @@ bool status; - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); f = fopen (filename, "w"); @@ -400,6 +403,6 @@ fprintf (f, "%9.6f ", dPAR[j]); } - fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", - source[0].type, source[0].mode, + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + source[0].apMag, source[0].type, source[0].mode, log10(model[0].chisq/model[0].nDOF), source[0].moments->SN, @@ -453,5 +456,6 @@ fprintf (f, "%9.6f ", dPAR[j]); } - fprintf (f, ": %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + source->apMag, source[0].type, source[0].mode, log10(model[0].chisq/model[0].nDOF), @@ -694,4 +698,35 @@ } +psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder) { + + bool status; + char keyword[80]; + + // get polynomial orders + // XXX add status failures tests + int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X"); + int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y"); + int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z"); + int nTorder = psMetadataLookupS32 (&status, folder, "NORDER_T"); + + psPolynomial4D *poly = psPolynomial4DAlloc (nXorder, nYorder, nZorder, nTorder, PS_POLYNOMIAL_ORD); + + for (int nx = 0; nx < poly->nX + 1; nx++) { + for (int ny = 0; ny < poly->nY + 1; ny++) { + for (int nz = 0; nz < poly->nZ + 1; nz++) { + for (int nt = 0; nt < poly->nT + 1; nt++) { + sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt); + poly->coeff[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword); + if (!status) poly->mask[nx][ny][nz][nt] = 1; + + sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt); + poly->coeffErr[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword); + } + } + } + } + return (poly); +} + // XXX : these may need F64, or %g format for output bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) { @@ -774,4 +809,48 @@ } +bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...) { + + int Nbyte; + char tmp; + char *root; + va_list argp; + + va_start (argp, format); + Nbyte = vsnprintf (&tmp, 0, format, argp); + va_end (argp); + + if (!Nbyte) return false; + + va_start (argp, format); + root = (char *) psAlloc (Nbyte + 1); + memset (root, 0, Nbyte + 1); + vsnprintf (root, Nbyte + 1, format, argp); + va_end (argp); + + psMetadata *folder = psMetadataAlloc (); + psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 4D polynomial", folder); + psFree (root); + + // specify the polynomial orders + psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX); + psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY); + psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ); + psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_T", PS_DATA_S32, "number of z orders", poly->nT); + + // place polynomial entries on folder + for (int nx = 0; nx < poly->nX + 1; nx++) { + for (int ny = 0; ny < poly->nY + 1; ny++) { + for (int nz = 0; nz < poly->nZ + 1; nz++) { + for (int nt = 0; nt < poly->nT + 1; nt++) { + if (poly->mask[nx][ny][nz][nt]) continue; + psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz][nt], nx, ny, nz, nt); + psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz][nt], nx, ny, nz, nt); + } + } + } + } + return true; +} + bool psphotWritePSF (pmPSF *psf, char *filename) { @@ -788,5 +867,5 @@ psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i); } - psPolynomial3DtoMD (psfdata, psf->ApTrend, "APTREND"); + psPolynomial4DtoMD (psfdata, psf->ApTrend, "APTREND"); psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); @@ -828,5 +907,5 @@ sprintf (keyword, "APTREND"); psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword); - psPolynomial3D *poly = psPolynomial3DfromMD (folder); + psPolynomial4D *poly = psPolynomial4DfromMD (folder); psFree (psf->ApTrend); psf->ApTrend = poly; Index: /trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- /trunk/psphot/src/psphotRadiusChecks.c (revision 5979) +++ /trunk/psphot/src/psphotRadiusChecks.c (revision 5980) @@ -31,13 +31,19 @@ } - // check if we need to redefine the pixels - // XXX EAM : a better test would examine the source pixels +# if (1) + bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); + return status; + +# else + if (model->radius > PSF_OUTER_RADIUS) { // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); return true; + } else { + fprintf (stderr, "skipping %f,%f\n", model->params->data.F32[2], model->params->data.F32[3]); } - return false; +# endif } @@ -68,4 +74,11 @@ if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); +# if (1) + // redefine the pixels if needed + bool status = psphotRedefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); + return status; + +# else + // check if we need to redefine the pixels if (model->radius > GAL_OUTER_RADIUS) { @@ -76,3 +89,4 @@ } return false; +# endif } Index: /trunk/psphot/src/psphotTest.c =================================================================== --- /trunk/psphot/src/psphotTest.c (revision 5979) +++ /trunk/psphot/src/psphotTest.c (revision 5980) @@ -1,35 +1,65 @@ # include "psphot.h" +bool apResidSetPower (float value); + +void psExit (int status, char *process, char *format, ...) { + + va_list ap; + + va_start (ap, format); + fprintf (stderr, "exiting %s\n", process); + vfprintf (stderr, format, ap); + va_end (ap); + + exit (status); +} int main (int argc, char **argv) { - char line[80]; - psImage *image; + // read in a file consisting of: + // x y radius fitmag apmag dmag - image = psImageAlloc (512, 512, PS_DATA_F32); - for (int j = 0; j < 512; j++) { - for (int i = 0; i < 512; i++) { - image->data.F32[j][i] = i; - } + double fit; + double ap; + + psphotTestArguments (&argc, argv); + + FILE *f = fopen (argv[1], "r"); + if (f == NULL) psExit (2, "test", "no file\n"); + + float power = atof (argv[2]); + apResidSetPower (power); + + psVector *x = psVectorAlloc (100, PS_TYPE_F64); + psVector *y = psVectorAlloc (100, PS_TYPE_F64); + psVector *rad = psVectorAlloc (100, PS_TYPE_F64); + psVector *rflux = psVectorAlloc (100, PS_TYPE_F64); + psVector *apres = psVectorAlloc (100, PS_TYPE_F64); + psVector *dmag = psVectorAlloc (100, PS_TYPE_F64); + + int Npt = 0; + while (fscanf (f, "%lf %lf %lf %lf %lf %lf", + &x->data.F64[Npt], + &y->data.F64[Npt], + &rad->data.F64[Npt], + &fit, &ap, + &dmag->data.F64[Npt]) != EOF) { + + if (fabs(ap-fit) > 0.1) continue; + if (fit > -10) continue; + + apres->data.F64[Npt] = ap-fit; + rflux->data.F64[Npt] = pow(10.0, 0.4*fit); + + psVectorExtend (x, 100, 1); + psVectorExtend (y, 100, 1); + psVectorExtend (rad, 100, 1); + psVectorExtend (rflux, 100, 1); + psVectorExtend (apres, 100, 1); + psVectorExtend (dmag, 100, 1); + Npt ++; } + x->n = y->n = rad->n = rflux->n = apres->n = dmag->n = Npt; - psMetadata *row = NULL; - - psArray *table = psArrayAlloc (100); - table->n = 0; - - for (int i = 0; i < 100; i++) { - row = psMetadataAlloc (); - sprintf (line, "line %03d", i); - psMetadataAdd (row, PS_LIST_HEAD, "FLOAT_DATA", PS_DATA_F32, "", 0.2*i); - psMetadataAdd (row, PS_LIST_HEAD, "INT_DATA", PS_DATA_S32, "", i); - psMetadataAdd (row, PS_LIST_HEAD, "CHAR_DATA", PS_DATA_STRING, "", line); - psArrayAdd (table, 100, row); - // psFree (row); - } - - psFits *fits = psFitsOpen ("test.fits", "w"); - psFitsWriteImage (fits, NULL, image, 0); - psFitsWriteTable (fits, NULL, table); - psFitsClose (fits); + psphotApResidPolyFit (x, y, rad, rflux, apres, dmag); exit (0); Index: /trunk/psphot/src/psphotTestArguments.c =================================================================== --- /trunk/psphot/src/psphotTestArguments.c (revision 5979) +++ /trunk/psphot/src/psphotTestArguments.c (revision 5980) @@ -2,10 +2,5 @@ static int usage (); -psMetadata *psphotTestArguments (int *argc, char **argv) { - - int N; - unsigned int Nfail; - int mode = PS_DATA_STRING | PS_META_REPLACE; - psF32 UseCoords_X, UseCoords_Y; +void psphotTestArguments (int *argc, char **argv) { // basic pslib options @@ -13,58 +8,12 @@ psArgumentVerbosity (argc, argv); - // identify options in args list - these go on config - bool UseCoords = false; - if ((N = psArgumentGet (*argc, argv, "-coords"))) { - UseCoords = true; - psArgumentRemove (N, argc, argv); - UseCoords_X = atof (argv[N]); - psArgumentRemove (N, argc, argv); - UseCoords_Y = atof (argv[N]); - psArgumentRemove (N, argc, argv); - } - - float threshold = -1; - if ((N = psArgumentGet (*argc, argv, "-threshold"))) { - psArgumentRemove (N, argc, argv); - threshold = atof (argv[N]); - psArgumentRemove (N, argc, argv); - } - - float fraction = 0.5; - if ((N = psArgumentGet (*argc, argv, "-frac"))) { - psArgumentRemove (N, argc, argv); - fraction = atof (argv[N]); - psArgumentRemove (N, argc, argv); - } - if (*argc != 3) usage (); - // load config information - psMetadata *config = psMetadataAlloc (); - psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); - config = psMetadataConfigParse (config, &Nfail, argv[2], FALSE); - - // identify input image - psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); - - mode = PS_DATA_F32 | PS_META_REPLACE; - psMetadataAdd (config, PS_LIST_HEAD, "TEST_THRESHOLD", mode, "", threshold); - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FRACTION", mode, "", fraction); - - if (UseCoords) { - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", mode, "", UseCoords_X); - psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", mode, "", UseCoords_Y); - } - - return (config); + return; } static int usage () { - fprintf (stderr, "USAGE: psphotModelTest (image.fits) (config)\n"); - fprintf (stderr, "options: \n"); - fprintf (stderr, " -coords x y : specify object center\n"); - fprintf (stderr, " -threshold (threshold) : (defaults to FRAC of peak\n"); - fprintf (stderr, " -frac (FRAC) : specify how to determine threshold\n"); + fprintf (stderr, "USAGE: psphotModelTest (apresid.dat) (power)\n"); exit (2); }