Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5992) +++ /trunk/psphot/Makefile (revision 5993) @@ -24,76 +24,55 @@ $(SRC)/psphotArguments.$(ARCH).o \ $(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psphotModelTest.$(ARCH).o \ $(SRC)/psphotImageStats.$(ARCH).o \ +$(SRC)/psphotImageBackground.$(ARCH).o \ +$(SRC)/psphotFindPeaks.$(ARCH).o \ $(SRC)/psphotSourceStats.$(ARCH).o \ +$(SRC)/psphotRoughClass.$(ARCH).o \ +$(SRC)/psphotBasicDeblend.$(ARCH).o \ $(SRC)/psphotChoosePSF.$(ARCH).o \ +$(SRC)/psphotEnsemblePSF.$(ARCH).o \ +$(SRC)/psphotFullFit.$(ARCH).o \ +$(SRC)/psphotBlendFit.$(ARCH).o \ $(SRC)/psphotApplyPSF.$(ARCH).o \ -$(SRC)/psphotFixedPSF.$(ARCH).o \ -$(SRC)/psphotEnsemblePSF.$(ARCH).o \ -$(SRC)/psphotReapplyPSF.$(ARCH).o \ $(SRC)/psphotFitGalaxies.$(ARCH).o \ -$(SRC)/psphotOutput.$(ARCH).o \ -$(SRC)/psphotMarkPSF.$(ARCH).o \ +$(SRC)/psphotReplaceUnfit.$(ARCH).o \ +$(SRC)/psphotApResid.$(ARCH).o \ +$(SRC)/psphotOutput.$(ARCH).o + +SUPPORT = \ +$(SRC)/psphotModelGroupInit.$(ARCH).o \ +$(SRC)/psphotGrowthCurve.$(ARCH).o \ +$(SRC)/psphotEvalPSF.$(ARCH).o \ +$(SRC)/psphotEvalFLT.$(ARCH).o \ +$(SRC)/psphotSourceFits.$(ARCH).o \ $(SRC)/psphotSortBySN.$(ARCH).o \ $(SRC)/psphotDefinePixels.$(ARCH).o \ $(SRC)/psphotMagnitudes.$(ARCH).o \ -$(SRC)/psphotRadiusChecks.$(ARCH).o \ -$(SRC)/psphotReplaceUnfit.$(ARCH).o \ -$(SRC)/psphotEvalPSF.$(ARCH).o \ -$(SRC)/psphotEvalFLT.$(ARCH).o \ -$(SRC)/psphotFullFit.$(ARCH).o \ +$(SRC)/psphotRadiusChecks.$(ARCH).o + +PSMODULES = \ +$(SRC)/psLine.$(ARCH).o \ +$(SRC)/psSparse.$(ARCH).o \ +$(SRC)/psPolynomialUtils.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ $(SRC)/pmSourceContour.$(ARCH).o \ -$(SRC)/psLine.$(ARCH).o \ -$(SRC)/psModulesUtils.$(ARCH).o \ -$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \ -$(SRC)/pmSourceFitFixed.$(ARCH).o \ -$(SRC)/psSparse.$(ARCH).o \ -$(SRC)/psImageData.$(ARCH).o \ -$(SRC)/psphotModelTest.$(ARCH).o \ -$(SRC)/psphotImageBackground.$(ARCH).o \ -$(SRC)/psphotBasicDeblend.$(ARCH).o \ -$(SRC)/psphotModelGroupInit.$(ARCH).o \ -$(SRC)/pmModelFitSet.$(ARCH).o \ -$(SRC)/pmSourceFitSet.$(ARCH).o \ -$(SRC)/psphotSourceFits.$(ARCH).o \ -$(SRC)/psphotBlendFit.$(ARCH).o \ -$(SRC)/psphotApResid.$(ARCH).o \ -$(SRC)/psphotAssessPSF.$(ARCH).o \ -$(SRC)/psBicube.$(ARCH).o - -MODELS = \ -$(SRC)/models/pmModel_GAUSS.c \ -$(SRC)/models/pmModel_PGAUSS.c \ -$(SRC)/models/pmModel_QGAUSS.c \ -$(SRC)/models/pmModel_SGAUSS.c \ -$(SRC)/models/pmModel_TGAUSS.c - -$(SRC)/pmModelGroup.$(ARCH).o: $(MODELS) - -MODELTEST = \ -$(SRC)/modelTestFitSource.$(ARCH).o \ -$(SRC)/modelTestArguments.$(ARCH).o \ -$(SRC)/psphotModelGroupInit.$(ARCH).o \ -$(SRC)/psModulesUtils.$(ARCH).o \ -$(SRC)/psImageData.$(ARCH).o \ -$(SRC)/psphotSetup.$(ARCH).o \ -$(SRC)/psphotDefinePixels.$(ARCH).o \ -$(SRC)/psphotModelTest.$(ARCH).o +$(SRC)/pmModelFitSet.$(ARCH).o \ +$(SRC)/pmSourceFitSet.$(ARCH).o TEST = \ -$(SRC)/psphotTest.$(ARCH).o \ -$(SRC)/psphotTestArguments.$(ARCH).o \ -$(SRC)/psphotApResid.$(ARCH).o +$(SRC)/psphotTest.$(ARCH).o \ +$(SRC)/psBicube.$(ARCH).o \ +$(SRC)/psphotTestArguments.$(ARCH).o psphot: $(BIN)/psphot.$(ARCH) -$(BIN)/psphot.$(ARCH) : $(PSPHOT) +$(BIN)/psphot.$(ARCH) : $(PSPHOT) $(PSMODULES) $(SUPPORT) $(PSPHOT) : $(SRC)/psphot.h -modeltest.install: psphotModelTest.install -modeltest: $(BIN)/psphotModelTest.$(ARCH) -$(BIN)/psphotModelTest.$(ARCH) : $(MODELTEST) -$(MODELTEST): $(SRC)/psphot.h +test.install: psphotTest.install +test: $(BIN)/psphotTest.$(ARCH) +$(DESTBIN)/test : $(DESTBIN)/psphotTest -test-install: psphotTest.install -test: $(BIN)/psphotTest.$(ARCH) $(BIN)/psphotTest.$(ARCH) : $(TEST) $(TEST): $(SRC)/psphot.h Index: unk/psphot/src/pmSourceFitFixed.c =================================================================== --- /trunk/psphot/src/pmSourceFitFixed.c (revision 5992) +++ (revision ) @@ -1,90 +1,0 @@ -# include "psphot.h" -# define USECLIP 0 - -bool pmSourceFitFixed (pmSource *source, pmModel *model) { - - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - - psF32 Ro, ymodel; - - psVector *params = model->params; - psVector *dparams = model->dparams; - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - - psF32 Sky = params->data.F32[0]; - - // construct the coordinate and value entries - psVector *x = psVectorAlloc(100, PS_TYPE_F32); - psVector *y = psVectorAlloc(100, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(100, PS_TYPE_F32); - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - - x->n = y->n = yErr->n = 0; - psS32 count = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j]) { - continue; - } - - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - ymodel = modelFunc (NULL, params, coord); - - x->data.F32[count] = source->pixels->data.F32[i][j]; - y->data.F32[count] = ymodel; - - // this test modifies the weight based on deviation from the model flux - Ro = 1.0 + fabs (x->data.F32[count] - ymodel) / sqrt(PS_SQR(ymodel - Sky) + PS_SQR(Sky)); - yErr->data.F32[count] = sqrt (source->weight->data.F32[i][j] * Ro); - - psVectorExtend (x, 100, 1); - psVectorExtend (y, 100, 1); - psVectorExtend (yErr, 100, 1); - count++; - } - } - - psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD); - -# if (USECLIP) - psVector *mask = psVectorAlloc (count, PS_TYPE_U8); - psVectorInit (mask, 0); - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->clipIter = 3; - poly = psVectorClipFitPolynomial1D (poly, stats, mask, 0, y, yErr, x); -# else - poly = psVectorFitPolynomial1D (poly, NULL, 0, y, yErr, x); -# endif - - // fitted values are variations from current values - params->data.F32[0] += poly->coeff[0]; - params->data.F32[1] *= poly->coeff[1]; - - // XXX need to get errors back from psVectorClipFit... - dparams->data.F32[0] = 0; - dparams->data.F32[1] = 0; - - // XXX need to get chisq back from psVectorClipFit... - model->chisq = 0; - model->nIter = 0; - model->nDOF = 0; - - psFree(x); - psFree(y); - psFree(yErr); - -# if (USECLIP) - psFree(stats); - psFree(mask); -# endif - - return true; -} - Index: /trunk/psphot/src/pmSourceFitSet.c =================================================================== --- /trunk/psphot/src/pmSourceFitSet.c (revision 5992) +++ /trunk/psphot/src/pmSourceFitSet.c (revision 5993) @@ -91,5 +91,5 @@ } - // PSF model only fits first 4 parameters, FLT model fits all + // PSF model only fits first 4 parameters, EXT model fits all int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1; Index: unk/psphot/src/psBicube.c =================================================================== --- /trunk/psphot/src/psBicube.c (revision 5992) +++ (revision ) @@ -1,143 +1,0 @@ -# include "psphot.h" - -float poly3DTestValue (float x, float y, float z) { - float value = 5 + 3*x - 5*x*x - 4*x*y - 2*y*y + 2.5*z - 3*x*z - 6*z*z; - return (value); -} - - -float poly2DTestValue (float x, float y) { - float value = 5 + 3*x - 5*x*x - 4*x*y - 2*y*y; - return (value); -} - -float bicubeTestValue (int x, int y) { - float value = 5 - 5*x*x - 2*y*y + 3*x - 4*x*y; - return (value); -} - -bool bicubeTest () { - - psImage *test = psImageAlloc (3, 3, PS_TYPE_F32); - for (int ny = 0; ny < 3; ny++) { - for (int nx = 0; nx < 3; nx++) { - test->data.F32[ny][nx] = bicubeTestValue (nx - 1, ny - 1); - } - } - psPolynomial2D *poly = psImageBicubeFit (test, 1, 1); - psPlane peak = psImageBicubeMin (poly); - - fprintf (stderr, "done test: %f, %f\n", peak.x, peak.y); - return true; -} - -bool poly3DTest_F64 () { - - int Ncell = 50; - psVector *x = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64); - psVector *y = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64); - psVector *z = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64); - psVector *f = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F64); - - fprintf (stderr, "starting %s\n", __func__); - int cell = 0; - for (int nx = 0; nx < Ncell; nx++) { - for (int ny = 0; ny < Ncell; ny++) { - for (int nz = 0; nz < Ncell; nz++) { - x->data.F64[cell] = nx; - y->data.F64[cell] = ny; - z->data.F64[cell] = nz; - f->data.F64[cell] = poly3DTestValue (nx, ny, nz); - cell++; - } - } - } - fprintf (stderr, "built input data\n"); - - psPolynomial3D *poly = psPolynomial3DAlloc (2, 2, 2, PS_POLYNOMIAL_ORD); - - psTimerStart ("fit"); - poly = psVectorFitPolynomial3D (poly, NULL, 0, f, NULL, x, y, z); - fprintf (stderr, "fit 3d : %6.1f msec\n", 1000*psTimerMark ("fit")); - - 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 : %f\n", nx, ny, nz, poly->coeff[nx][ny][nz]); - } - } - } - return true; -} - -bool poly3DTest_F32 () { - - int Ncell = 50; - psVector *x = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32); - psVector *y = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32); - psVector *z = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32); - psVector *f = psVectorAlloc (Ncell*Ncell*Ncell, PS_TYPE_F32); - - fprintf (stderr, "starting %s\n", __func__); - int cell = 0; - for (int nx = 0; nx < Ncell; nx++) { - for (int ny = 0; ny < Ncell; ny++) { - for (int nz = 0; nz < Ncell; nz++) { - x->data.F32[cell] = nx; - y->data.F32[cell] = ny; - z->data.F32[cell] = nz; - f->data.F32[cell] = poly3DTestValue (nx, ny, nz); - cell++; - } - } - } - fprintf (stderr, "built input data\n"); - - psPolynomial3D *poly = psPolynomial3DAlloc (2, 2, 2, PS_POLYNOMIAL_ORD); - - psTimerStart ("fit"); - poly = psVectorFitPolynomial3D (poly, NULL, 0, f, NULL, x, y, z); - fprintf (stderr, "fit 3d : %6.1f msec\n", 1000*psTimerMark ("fit")); - - 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 : %f\n", nx, ny, nz, poly->coeff[nx][ny][nz]); - } - } - } - return true; -} - -bool poly2DTest () { - - int Ncell = 50; - psVector *x = psVectorAlloc (Ncell*Ncell, PS_TYPE_F64); - psVector *y = psVectorAlloc (Ncell*Ncell, PS_TYPE_F64); - psVector *f = psVectorAlloc (Ncell*Ncell, PS_TYPE_F64); - - fprintf (stderr, "starting __func__\n"); - int cell = 0; - for (int nx = 0; nx < Ncell; nx++) { - for (int ny = 0; ny < Ncell; ny++) { - x->data.F64[cell] = nx; - y->data.F64[cell] = ny; - f->data.F64[cell] = poly2DTestValue (nx, ny); - cell++; - } - } - fprintf (stderr, "built input data\n"); - - psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD); - - psTimerStart ("fit"); - poly = psVectorFitPolynomial2D (poly, NULL, 0, f, NULL, x, y); - fprintf (stderr, "fit 2d : %6.1f msec\n", 1000*psTimerMark ("fit")); - - for (int ny = 0; ny <= poly->nY; ny++) { - for (int nx = 0; nx <= poly->nX; nx++) { - fprintf (stderr, "%d %d : %f\n", nx, ny, poly->coeff[nx][ny]); - } - } - return true; -} Index: unk/psphot/src/psLibUtils.c =================================================================== --- /trunk/psphot/src/psLibUtils.c (revision 5992) +++ (revision ) @@ -1,65 +1,0 @@ -# include // for strncasecmp -# include // for write -# include -# include "psLibUtils.h" - -// XXX EAM : These two functions (psLogArguments and psTraceArguments) have been -// rolled together into a single function in psLib -// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) -// 2 = default, -v = 3, -vv = 4, -vvv = 5 -psS32 psLogArguments (int *argc, char **argv) { - - int N, level; - - // default log level is 2 - level = 2; - - // set in order, so that -vvv overrides -vv overrides -v - if ((N = psArgumentGet (*argc, argv, "-v"))) { - psArgumentRemove (N, argc, argv); - level = 3; - } - if ((N = psArgumentGet (*argc, argv, "-vv"))) { - psArgumentRemove (N, argc, argv); - level = 4; - } - if ((N = psArgumentGet (*argc, argv, "-vvv"))) { - psArgumentRemove (N, argc, argv); - level = 5; - } - - if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { - if (*argc < N + 2) { - psAbort ("psLogArguments", "USAGE: -logfmt (format)"); - } - psArgumentRemove (N, argc, argv); - psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid - psArgumentRemove (N, argc, argv); - } - - // set the level, return the level - psLogSetLevel (level); - return (level); -} - -// set trace levels by facility -psS32 psTraceArguments (int *argc, char **argv) { - - int N; - - // argument format is: -trace (facil) (level) - while ((N = psArgumentGet (*argc, argv, "-trace"))) { - if (*argc < N + 3) { - psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); - } - psArgumentRemove (N, argc, argv); - psTraceSetLevel (argv[N], atoi(argv[N+1])); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { - psTracePrintLevels (); - exit (2); - } - return (TRUE); -} Index: /trunk/psphot/src/psModulesUtils.c =================================================================== --- /trunk/psphot/src/psModulesUtils.c (revision 5992) +++ /trunk/psphot/src/psModulesUtils.c (revision 5993) @@ -123,2 +123,41 @@ return (new); } + +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) { + + int Xs, Xe, Ys, Ye; + int xi, xj, yi, yj; + int xIs, xJs, yIs, yJs; + int xIe, yIe; + float flux; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Ti = Mi->mask; + psImage *Tj = Mj->mask; + + Xs = PS_MAX (Pi->col0, Pj->col0); + Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); + + Ys = PS_MAX (Pi->row0, Pj->row0); + Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); + + xIs = Xs - Pi->col0; + xJs = Xs - Pj->col0; + yIs = Ys - Pi->row0; + yJs = Ys - Pj->row0; + + xIe = Xe - Pi->col0; + yIe = Ye - Pi->row0; + + flux = 0; + for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { + for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { + if (Ti->data.U8[yi][xi]) continue; + if (Tj->data.U8[yj][xj]) continue; + flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]; + } + } + return (flux); +} Index: /trunk/psphot/src/psPolynomialUtils.c =================================================================== --- /trunk/psphot/src/psPolynomialUtils.c (revision 5993) +++ /trunk/psphot/src/psPolynomialUtils.c (revision 5993) @@ -0,0 +1,168 @@ +# include +# include // for strcasecmp +# include // for unlink +# include + +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); +} + +psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) { + + int ix = x - image->col0; + int iy = y - image->row0; + + psF32 *Fm = &image->data.F32[iy - 1][ix]; + psF32 *Fo = &image->data.F32[iy + 0][ix]; + psF32 *Fp = &image->data.F32[iy + 1][ix]; + + double Fxm = Fm[-1] + Fo[-1] + Fp[-1]; + double Fxp = Fm[+1] + Fo[+1] + Fp[+1]; + double Fym = Fm[-1] + Fm[+0] + Fm[+1]; + double Fyp = Fp[-1] + Fp[+0] + Fp[+1]; + double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1]; + + psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD); + poly->mask[2][2] = 1; + poly->mask[1][2] = 1; + poly->mask[2][1] = 1; + + poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ; + + poly->coeff[1][0] = (Fxp - Fxm)/6.0; + poly->coeff[0][1] = (Fyp - Fym)/6.0; + + poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0; + poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0; + + poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0; + + return (poly); +} + +psPlane psImageBicubeMin (psPolynomial2D *poly) { + + psPlane min; + + min.xErr = min.yErr = 0; + + double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); + + min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det; + min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det; + return (min); +} + +// may need to extend this to change the keyname in the copy +bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key) { + + psMetadataItem *item = psMetadataLookup (in, key); + if (item == NULL) return false; + + psMetadataAddItem (out, item, PS_LIST_TAIL, PS_META_REPLACE); + return true; +} Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5992) +++ /trunk/psphot/src/psphot.c (revision 5993) @@ -67,12 +67,10 @@ case 3: - // fit extended objects with galaxy models psphotApplyPSF (imdata, config, sources, psf, sky); break; case 4: - // fit extended objects with galaxy models psphotApplyPSF (imdata, config, sources, psf, sky); - psphotFitGalaxies (imdata, config, sources, sky); + psphotFitExtended (imdata, config, sources, sky); break; } Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5992) +++ /trunk/psphot/src/psphot.h (revision 5993) @@ -24,29 +24,41 @@ psMetadata *psphotArguments (int *argc, char **argv); eamReadout *psphotSetup (psMetadata *config); +bool psphotModelTest (eamReadout *imdata, psMetadata *config); psStats *psphotImageStats (eamReadout *imdata, psMetadata *config); +psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky); psArray *psphotFindPeaks (eamReadout *imdata, psMetadata *config, psStats *sky); psArray *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks); bool psphotRoughClass (psArray *sources, psMetadata *config); +bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky); pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); -bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); void psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); -bool psphotSubtractPSF (pmSource *source); +// optional object analysis steps +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotReplaceUnfit (psArray *sources); +bool psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats); +bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf); + +// basic support functions +pmModel *pmModelCopy (pmModel *model); +pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); +psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); +void psphotModelGroupInit (void); int psphotSortBySN (const void **a, const void **b); int psphotSortByY (const void **a, const void **b); -int psphotSaveImage (psMetadata *header, psImage *image, char *filename); +bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf); +void psphotTestArguments (int *argc, char **argv); + +// functions to set the correct source pixels +bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type); +bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model); +bool psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type); +bool psphotCheckRadiusEXT (eamReadout *imdata, pmSource *source, pmModel *model); 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); - -bool pmSourceFitFixed (pmSource *source, pmModel *model); -pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); -pmModel *pmSourceSelectModel (pmSource *source); - -// eamReadout functions -eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); // output functions @@ -57,70 +69,45 @@ bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky); +bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); +bool pmModelWriteEXTs (psArray *sources, char *filename); +bool pmModelWriteNULLs (psArray *sources, char *filename); bool pmPeaksWriteText (psArray *sources, char *filename); bool pmMomentsWriteText (psArray *sources, char *filename); -bool pmModelWritePSFs (psArray *sources, psMetadata *config, char *filename, pmPSF *psf); -bool pmModelWriteFLTs (psArray *sources, char *filename); -bool pmModelWriteNULLs (psArray *sources, char *filename); +bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image); +bool psphotDumpMoments (psMetadata *config, psArray *sources); +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); +bool psphotUpdateHeader (psMetadata *header, psMetadata *config); int pmSourcesDophotType (pmSource *source); +bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key); -// psphotModelTest functions -bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -float psphotCrossProduct (pmSource *Mi, pmSource *Mj); -psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky); -bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +// PSF / DBL / EXT evaluation functions +bool psphotEvalPSF (pmSource *source, pmModel *model); +bool psphotEvalDBL (pmSource *source, pmModel *model); +bool psphotEvalEXT (pmSource *source, pmModel *model); -pmModel *pmModelCopy (pmModel *model); -psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); -bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky); +// functions to support simultaneous multi-source fitting +bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF); +bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF); +psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x); +bool pmModelFitSetInit (pmModelType type); -bool psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool psphotInitLimitsPSF (psMetadata *config); -bool psphotEvalPSF (pmSource *source, pmModel *model); -bool psphotEvalDBL (pmSource *source, pmModel *model); -bool psphotEvalFLT (pmSource *source, pmModel *model); -bool psphotInitRadiusPSF (psMetadata *config, psStats *sky, pmModelType type); -bool psphotCheckRadiusPSF (eamReadout *imdata, pmSource *source, pmModel *model); -bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type); -bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model); -bool psphotSamplePSFs (psMetadata *config, pmPSF *psf, psImage *image); -bool psphotReplaceUnfit (psArray *sources); -bool psphotDumpMoments (psMetadata *config, psArray *sources); -bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf); -bool psphotAssessPSF (eamReadout *imdata, psMetadata *config, pmPSF *psf); +// functions to support the source fitting process +bool psphotInitLimitsPSF (psMetadata *config); +bool psphotInitLimitsEXT (psMetadata *config, psStats *sky); +bool psphotFitBlend (eamReadout *imdata, pmSource *source); +bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources); +bool psphotFitPSF (eamReadout *imdata, pmSource *source); +pmModel *psphotFitEXT (eamReadout *imdata, pmSource *source); +psArray *psphotFitDBL (eamReadout *imdata, pmSource *source); -bool psphotWritePSF (pmPSF *psf, char *filename); -pmPSF *psphotReadPSF (char *filename); +// eamReadout functions +eamReadout *eamReadoutAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); -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); +// bicubic interpolation +psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); +psPlane psImageBicubeMin (psPolynomial2D *poly); -bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF); -bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF); -psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x); -bool pmModelFitSetInit (pmModelType type); - -bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); - -bool psphotInitLimitsFLT (psMetadata *config, psStats *sky); - -bool psphotFitPSF (eamReadout *imdata, pmSource *source); -bool psphotFitBlend (eamReadout *imdata, pmSource *source); -bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources); - -pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source); -psArray *psphotFitDBL (eamReadout *imdata, pmSource *source); - -psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); -psPlane psImageBicubeMin (psPolynomial2D *poly); -void psphotTestArguments (int *argc, char **argv); - -bool psphotGrowthCurve (eamReadout *imdata, pmPSF *psf); - +// optional mode for clip fit? psPolynomial4D *psVectorChiClipFitPolynomial4D( psPolynomial4D *poly, @@ -134,6 +121,2 @@ 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 5992) +++ /trunk/psphot/src/psphotApResid.c (revision 5993) @@ -1,111 +1,3 @@ # include "psphot.h" - -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 @@ -347,9 +239,11 @@ psf->nApResid = residStats->clippedNvalues; - /* - (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) - (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) - (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) - */ + // save results for later output + psMetadataAdd (config, PS_LIST_TAIL, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); + psMetadataAdd (config, PS_LIST_TAIL, "SKYSAT", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skySat); + psMetadataAdd (config, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (config, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (config, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->growth->apLoss); + psMetadataAdd (config, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid); psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot")); @@ -370,2 +264,8 @@ } +/* + (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) + (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) + (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) +*/ + Index: /trunk/psphot/src/psphotApplyPSF.c =================================================================== --- /trunk/psphot/src/psphotApplyPSF.c (revision 5992) +++ /trunk/psphot/src/psphotApplyPSF.c (revision 5993) @@ -25,9 +25,9 @@ // use the source moments, etc to guess basic model parameters - pmModel *modelFLT = pmSourceModelGuess (source, psf->type); + pmModel *modelEXT = pmSourceModelGuess (source, psf->type); // set PSF parameters for this model - pmModel *model = pmModelFromPSF (modelFLT, psf); - psFree (modelFLT); + pmModel *model = pmModelFromPSF (modelEXT, psf); + psFree (modelEXT); source->modelPSF = model; Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 5992) +++ /trunk/psphot/src/psphotArguments.c (revision 5993) @@ -123,5 +123,6 @@ psMetadataAdd (config, PS_LIST_TAIL, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); - config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE); + // config file values override defaults set here + config = psMetadataConfigParse (config, &Nfail, argv[3], TRUE); psMetadataItem *item = NULL; @@ -129,7 +130,7 @@ while ((item = psMetadataGetAndIncrement (iter)) != NULL) { psMetadataAddItem (config, item, PS_LIST_TAIL, PS_META_REPLACE); - // psFree (item); + psFree (item); } - // psFree (iter); + psFree (iter); // identify input image & optional weight & mask images Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 5992) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 5993) @@ -36,5 +36,5 @@ psArray *overlap = psArrayAlloc (100); - // XXX make sure this results in decreasing, not increasing, SN + // examine sources in decreasing SN order for (int i = sources->n - 1; i >= 0; i--) { N = index->data.U32[i]; @@ -73,10 +73,9 @@ // generate source contour (1/4 peak counts) if (overlap->n > 0) { - // XXX EAM : make the 1/4 user-defined. - // XXX keep threshold from dropping too low (N*sky.sigma) + // set the threshold based on user inputs threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; threshold = PS_MAX (threshold, minThreshold); - // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now + // generate a basic contour (set of x,y coordinates at-or-below flux level) psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold); if (contour == NULL) continue; @@ -97,23 +96,4 @@ if (xv->data.F32[j+1] < testSource->peak->x) break; - # if (0) - int xp0 = source->moments->x - source->pixels->col0; - int xp1 = source->peak->x - source->pixels->col0; - int xp2 = testSource->moments->x - source->pixels->col0; - int xp3 = testSource->peak->x - source->pixels->col0; - - int yp0 = source->moments->y - source->pixels->row0; - int yp1 = source->peak->y - source->pixels->row0; - int yp2 = testSource->moments->y - testSource->pixels->row0; - int yp3 = testSource->peak->y - testSource->pixels->row0; - - fprintf (f, "%d %d (%f, %f) : %d %d (%f, %f) vs %f\n", - source->peak->x, source->peak->y, - source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1], - testSource->peak->x, testSource->peak->y, - testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold - ); - # endif - testSource->mode |= PM_SOURCE_BLEND; Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 5992) +++ /trunk/psphot/src/psphotBlendFit.c (revision 5993) @@ -1,6 +1,6 @@ # include "psphot.h" -bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) -{ +// XXX I don't like this name +bool psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { psTimerStart ("psphot"); @@ -10,5 +10,5 @@ psphotInitLimitsPSF (config); - psphotInitLimitsFLT (config, sky); + psphotInitLimitsEXT (config, sky); psphotInitRadiusPSF (config, sky, psf->type); Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 5992) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 5993) @@ -2,7 +2,6 @@ // try PSF models and select best option +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) { -pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) -{ bool status; char *modelName; @@ -55,5 +54,5 @@ } psFree (iter); - // psFree (list); XXX EAM - is list freed with iter? + // psFree (list); XXX is list freed with iter? psFree (stars); @@ -78,5 +77,5 @@ for (int i = 0; i < try->sources->n; i++) { pmSource *source = try->sources->data[i]; - if (try->mask->data.U8[i] & PSFTRY_MASK_FLT_FAIL) { + if (try->mask->data.U8[i] & PSFTRY_MASK_EXT_FAIL) { source->mode &= ~PM_SOURCE_PSFSTAR; } @@ -86,4 +85,6 @@ pmPSF *psf = psMemCopy(try->psf); psFree (models); + + psMetadataAdd (config, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "number of stars used to make PSF", psf->nPSFstars); modelName = pmModelGetType (psf->type); Index: /trunk/psphot/src/psphotDefinePixels.c =================================================================== --- /trunk/psphot/src/psphotDefinePixels.c (revision 5992) +++ /trunk/psphot/src/psphotDefinePixels.c (revision 5993) @@ -49,9 +49,8 @@ if (extend) { // Grab a new subimage - // psFree (mySource->pixels); - // psFree (mySource->weight); - // psFree (mySource->mask); + 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); @@ -62,6 +61,2 @@ 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 5992) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 5993) @@ -1,49 +1,3 @@ # include "psphot.h" - -psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) { - - int ix = x - image->col0; - int iy = y - image->row0; - - psF32 *Fm = &image->data.F32[iy - 1][ix]; - psF32 *Fo = &image->data.F32[iy + 0][ix]; - psF32 *Fp = &image->data.F32[iy + 1][ix]; - - double Fxm = Fm[-1] + Fo[-1] + Fp[-1]; - double Fxp = Fm[+1] + Fo[+1] + Fp[+1]; - double Fym = Fm[-1] + Fm[+0] + Fm[+1]; - double Fyp = Fp[-1] + Fp[+0] + Fp[+1]; - double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1]; - - psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD); - poly->mask[2][2] = 1; - poly->mask[1][2] = 1; - poly->mask[2][1] = 1; - - poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ; - - poly->coeff[1][0] = (Fxp - Fxm)/6.0; - poly->coeff[0][1] = (Fyp - Fym)/6.0; - - poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0; - poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0; - - poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0; - - return (poly); -} - -psPlane psImageBicubeMin (psPolynomial2D *poly) { - - psPlane min; - - min.xErr = min.yErr = 0; - - double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]); - - min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det; - min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det; - return (min); -} bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { @@ -107,10 +61,10 @@ } - // XXX EAM : add option to use FLT or PSF form + // XXX EAM : add option to use EXT or PSF form // use the source moments, etc to guess basic model parameters - pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); + pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); if (inSource->mode & PM_SOURCE_SATSTAR) { - modelFLT->params->data.F32[2] = inSource->moments->x; - modelFLT->params->data.F32[3] = inSource->moments->y; + modelEXT->params->data.F32[2] = inSource->moments->x; + modelEXT->params->data.F32[3] = inSource->moments->y; } else { // peak-up on peak (for non-sat objects) @@ -123,10 +77,10 @@ psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", - modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy); + modelEXT->params->data.F32[2], modelEXT->params->data.F32[3], min.x + ix, min.y + iy); // if min point is too deviant, keep the old value if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { - modelFLT->params->data.F32[2] = min.x + ix; - modelFLT->params->data.F32[3] = min.y + iy; + modelEXT->params->data.F32[2] = min.x + ix; + modelEXT->params->data.F32[3] = min.y + iy; } psFree (bicube); @@ -134,6 +88,6 @@ // set PSF parameters for this model - pmModel *model = pmModelFromPSF (modelFLT, psf); - psFree (modelFLT); + pmModel *model = pmModelFromPSF (modelEXT, psf); + psFree (modelEXT); // save the original coords @@ -185,9 +139,9 @@ // diagonal element (auto-cross-product) - f = psphotCrossProduct (Mi, Mi); + f = pmSourceCrossProduct (Mi, Mi); psSparseMatrixElement (sparse, i, i, f); // find the image x model value - f = psphotCrossProduct (Fi, Mi); + f = pmSourceCrossProduct (Fi, Mi); psSparseVectorElement (sparse, i, f); @@ -203,5 +157,5 @@ // got an overlap; calculate cross-product and add to output array - f = psphotCrossProduct (Mi, Mj); + f = pmSourceCrossProduct (Mi, Mj); psSparseMatrixElement (sparse, j, i, f); } @@ -235,42 +189,2 @@ return true; } - - -float psphotCrossProduct (pmSource *Mi, pmSource *Mj) { - - int Xs, Xe, Ys, Ye; - int xi, xj, yi, yj; - int xIs, xJs, yIs, yJs; - int xIe, yIe; - float flux; - - psImage *Pi = Mi->pixels; - psImage *Pj = Mj->pixels; - - psImage *Ti = Mi->mask; - psImage *Tj = Mj->mask; - - Xs = PS_MAX (Pi->col0, Pj->col0); - Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); - - Ys = PS_MAX (Pi->row0, Pj->row0); - Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); - - xIs = Xs - Pi->col0; - xJs = Xs - Pj->col0; - yIs = Ys - Pi->row0; - yJs = Ys - Pj->row0; - - xIe = Xe - Pi->col0; - yIe = Ye - Pi->row0; - - flux = 0; - for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { - for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.U8[yi][xi]) continue; - if (Tj->data.U8[yj][xj]) continue; - flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]; - } - } - return (flux); -} Index: /trunk/psphot/src/psphotEvalFLT.c =================================================================== --- /trunk/psphot/src/psphotEvalFLT.c (revision 5992) +++ /trunk/psphot/src/psphotEvalFLT.c (revision 5993) @@ -1,9 +1,9 @@ # include "psphot.h" -bool psphotEvalFLT (pmSource *source, pmModel *model) +bool psphotEvalEXT (pmSource *source, pmModel *model) { int keep; - // do we actually have a valid FLT model? + // do we actually have a valid EXT model? if (model == NULL) { source->mode &= ~PM_SOURCE_FITTED; @@ -22,15 +22,14 @@ case PM_MODEL_OFFIMAGE: default: - psLogMsg ("psphot", 5, "GAL fail fit\n"); + psLogMsg ("psphot", 5, "EXT fail fit\n"); source->mode |= PM_SOURCE_FAIL; return false; } - // unless we prove otherwise, this object is a galaxy. - // XXX EAM : do we need to save the old type? - source->type = PM_SOURCE_GALAXY; + // unless we prove otherwise, this object is extended + source->type = PM_SOURCE_EXTENDED; - // the following source->mode information pertains to modelPSF: - source->mode |= PM_SOURCE_FLTMODEL; + // the following source->mode information pertains to modelEXT: + source->mode |= PM_SOURCE_EXTMODEL; // if the object has a fitted peak below 0, the fit did not converge cleanly @@ -44,5 +43,5 @@ // poor-quality fit; only keep if nothing else works... - psLogMsg ("psphot", 5, "GAL poor fit\n"); + psLogMsg ("psphot", 5, "EXT poor fit\n"); source->mode |= PM_SOURCE_POOR; return false; Index: /trunk/psphot/src/psphotEvalPSF.c =================================================================== --- /trunk/psphot/src/psphotEvalPSF.c (revision 5992) +++ /trunk/psphot/src/psphotEvalPSF.c (revision 5993) @@ -19,5 +19,5 @@ // saturated stars should fall outside (larger), but have peaks above SATURATION -// galaxies should be larger, cosmic rays smaller +// extended sources should be larger, cosmic rays smaller // we also reject objects with S/N too low or ChiSquare to high @@ -138,7 +138,7 @@ } - // object appears to be large, suspected galaxy + // object appears to be large, suspected extended source if ((nSx >= PSF_SHAPE_NSIGMA) || (nSy >= PSF_SHAPE_NSIGMA)) { - source->type = PM_SOURCE_GALAXY; + source->type = PM_SOURCE_EXTENDED; return false; } Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 5992) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 5993) @@ -1,5 +1,5 @@ # include "psphot.h" -bool psphotFitGalaxies (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) +bool psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats) { bool status; @@ -12,10 +12,10 @@ psTimerStart ("psphot"); - float GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); - float GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); - char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); + float EXT_MIN_SN = psMetadataLookupF32 (&status, config, "EXT_MIN_SN"); + float EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS"); + char *modelName = psMetadataLookupPtr (&status, config, "EXT_MODEL"); pmModelType modelType = pmModelSetType (modelName); - psphotInitRadiusFLT (config, skyStats, modelType); + psphotInitRadiusEXT (config, skyStats, modelType); for (int i = 0; i < sources->n; i++) { @@ -23,12 +23,12 @@ // only fit suspected extended sources - if (source->type != PM_SOURCE_GALAXY) continue; + if (source->type != PM_SOURCE_EXTENDED) continue; if (source->mode & PM_SOURCE_BLEND) continue; // skip possible extended sources below fit threshold - if (source->moments->SN < GAL_MIN_SN) continue; + if (source->moments->SN < EXT_MIN_SN) continue; - // recalculate the source moments using the larger galaxy moments radius - status = pmSourceMoments (source, GAL_MOMENTS_RAD); + // recalculate the source moments using the larger extended-source moments radius + status = pmSourceMoments (source, EXT_MOMENTS_RAD); if (!status) { source->mode |= PM_SOURCE_FAIL; // better choice? @@ -38,12 +38,12 @@ // use the source moments, etc to guess basic model parameters pmModel *model = pmSourceModelGuess (source, modelType); - source->modelFLT = model; + source->modelEXT = model; x = model->params->data.F32[2]; y = model->params->data.F32[3]; // sets the model radius (via source->model) and adjusts pixels as needed - psphotCheckRadiusFLT (imdata, source, model); + psphotCheckRadiusEXT (imdata, source, model); - // fit FLT (not PSF) model (set/unset the pixel mask) + // fit EXT (not PSF) model (set/unset the pixel mask) psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); pmSourceFitModel (source, model, false); @@ -61,6 +61,6 @@ } } - psLogMsg ("psphot", 3, "fit GAL models: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); - psLogMsg ("psphot", 4, "subtracted %d GAL models\n", Nsub); + psLogMsg ("psphot", 3, "fit EXT models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d EXT objects\n", Nsub); return (true); } Index: unk/psphot/src/psphotFixedPSF.c =================================================================== --- /trunk/psphot/src/psphotFixedPSF.c (revision 5992) +++ (revision ) @@ -1,78 +1,0 @@ -# include "psphot.h" - -// fit psf model to all objects with centroid fixed - -bool psphotFixedPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) -{ - bool status; - float x; - float y; - int Nfit = 0; - - psTimerStart ("psphot"); - - // we may set this differently here from the value used to mark likely saturated stars - float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); - float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); - - // set the object surface-brightness limit for fitted pixels - float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; - psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); - - // this function specifies the radius at this the model hits the given flux - pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); - - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // skip non-astronomical objects (very likely defects) - // XXX EAM : should we try these anyway? - if (source->type == PM_SOURCE_DEFECT) continue; - if (source->type == PM_SOURCE_SATURATED) continue; - - // use the source moments, etc to guess basic model parameters - pmModel *modelFLT = pmSourceModelGuess (source, psf->type); - - // set PSF parameters for this model - pmModel *model = pmModelFromPSF (modelFLT, psf); - // model->params->data.F32[2] = source->peak->x; - // model->params->data.F32[3] = source->peak->y; - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - psFree (modelFLT); - - // set the fit radius based on the object flux limit and the model - // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) - model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; - if (isnan(model->radius)) { - psAbort ("apply_psf_model", "error in radius"); - } - - psphotRedefinePixels (source, imdata, x, y, model->radius); - - // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitFixed (source, model); - psImageKeepCircle (source->mask, x, y, model->radius, "AND", ~PSPHOT_MASK_MARKED); - - if (!status || (model->params->data.F32[1] < 0)) { - psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (peak %f\n", x, y, model->params->data.F32[1]); - source->mode |= PM_SOURCE_FAIL; // better choice? - psFree (model); - continue; - } - - // model succeeded : tentatively keep it - source->modelPSF = model; - Nfit ++; - - // subtract object, leave local sky - pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); - source->mode |= PM_SOURCE_SUBTRACTED; - source->mode &= ~PM_SOURCE_TEMPSUB; - } - - psLogMsg ("psphot", 3, "fit fixed PSF models: %f sec for %d objects\n", psTimerMark ("psphot"), Nfit); - return (true); -} Index: /trunk/psphot/src/psphotFullFit.c =================================================================== --- /trunk/psphot/src/psphotFullFit.c (revision 5992) +++ /trunk/psphot/src/psphotFullFit.c (revision 5993) @@ -14,15 +14,15 @@ sources = psArraySort (sources, psphotSortBySN); - // galaxy model parameters - float GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); - float GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); + // extended source model parameters + float EXT_MIN_SN = psMetadataLookupF32 (&status, config, "EXT_MIN_SN"); + float EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS"); - // 'galaxy' model descriptions - char *modelNameFLT = psMetadataLookupPtr (&status, config, "GAL_MODEL"); - pmModelType modelTypeFLT = pmModelSetType (modelNameFLT); + // extended source model descriptions + char *modelNameEXT = psMetadataLookupPtr (&status, config, "EXT_MODEL"); + pmModelType modelTypeEXT = pmModelSetType (modelNameEXT); psphotInitLimitsPSF (config); psphotInitRadiusPSF (config, sky, psf->type); - psphotInitRadiusFLT (config, sky, modelTypeFLT); + psphotInitRadiusEXT (config, sky, modelTypeEXT); for (int i = 0; i < sources->n; i++) { @@ -82,34 +82,34 @@ // skip the source if we don't think it is extended - if (source->type != PM_SOURCE_GALAXY) goto subLINEAR; - if (source->moments->SN < GAL_MIN_SN) goto subLINEAR; + if (source->type != PM_SOURCE_EXTENDED) goto subLINEAR; + if (source->moments->SN < EXT_MIN_SN) goto subLINEAR; // add the double-PSF fit mode // how do we compare the results? chisq? - // recalculate the source moments using the larger galaxy moments radius - if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) goto subLINEAR; + // recalculate the source moments using the larger extended-source moments radius + if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) goto subLINEAR; // use the source moments, etc to guess basic model parameters - source->modelFLT = pmSourceModelGuess (source, modelTypeFLT); - pmModel *modelFLT = source->modelFLT; + source->modelEXT = pmSourceModelGuess (source, modelTypeEXT); + pmModel *modelEXT = source->modelEXT; - psphotCheckRadiusFLT (imdata, source, modelFLT); + psphotCheckRadiusEXT (imdata, source, modelEXT); - x = modelFLT->params->data.F32[2]; - y = modelFLT->params->data.F32[3]; + x = modelEXT->params->data.F32[2]; + y = modelEXT->params->data.F32[3]; - // fit FLT (not PSF) model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, modelFLT->radius, "OR", PSPHOT_MASK_MARKED); - pmSourceFitModel (source, modelFLT, false); - psImageKeepCircle (source->mask, x, y, modelFLT->radius, "AND", ~PSPHOT_MASK_MARKED); + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, modelEXT->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, modelEXT, false); + psImageKeepCircle (source->mask, x, y, modelEXT->radius, "AND", ~PSPHOT_MASK_MARKED); // track model evaluations - Niter += modelFLT[0].nIter; + Niter += modelEXT[0].nIter; Nfit++; - // does the FLT model succeed? - if (psphotEvalFLT (source, source->modelFLT)) { - pmSourceSubModel (source->pixels, source->mask, source->modelFLT, false, false); + // does the EXT model succeed? + if (psphotEvalEXT (source, source->modelEXT)) { + pmSourceSubModel (source->pixels, source->mask, source->modelEXT, false, false); source->mode |= PM_SOURCE_SUBTRACTED; source->mode &= ~PM_SOURCE_TEMPSUB; Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 5992) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 5993) @@ -21,6 +21,6 @@ break; - case PM_SOURCE_GALAXY: - model = source->modelFLT; + case PM_SOURCE_EXTENDED: + model = source->modelEXT; if (model == NULL) return NULL; radius = model->radius; Index: unk/psphot/src/psphotMarkPSF.c =================================================================== --- /trunk/psphot/src/psphotMarkPSF.c (revision 5992) +++ (revision ) @@ -1,132 +1,0 @@ -# include "psphot.h" - -// given a pmSource which has been fitted using modelPSF, evaluate the -// resulting fit: did the fit succeed? is this object PSF-like? is this object -// extended? return status is TRUE for a valid PSF, FALSE otherwise. - -// identify objects consistent with PSF shape/magnitude distribution -// we expect dparams[4],dparams[5] to have a scatter of: -// sigma_x / (S/N) * sqrt(2) -// 1 / (params[4],params[5])*(S/N) - -// sigma_x : 1 / SX -// dsx : 1 / (SX * SN) -// dsx_o = hypot (1/(SX*SN), MIN_DSX) - -// any objects which is consistent with the PSF should have -// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) -// this includes a minimum buffer (DS) for the brighter objects - -// saturated stars should fall outside (larger), but have peaks above SATURATE -// galaxies should be larger, cosmic rays smaller -// we also reject objects with S/N too low or ChiSquare to high - -// floor for DS value -// XXX EAM : add to configuration? -# define MIN_DS 0.01 - -bool psphotMarkPSF (pmSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE) -{ - int keep; - float dSX, dSY, SX, SY, SN; - float nSx, nSy, Chi; - - pmModel *model = source->modelPSF; - - // do we actually have a valid PSF model? - if (model == NULL) { - source->mode &= ~PM_SOURCE_FITTED; - return false; - } - - // did the model fit fail for one or another reason? - switch (model->status) { - case PM_MODEL_UNTRIED: - source->mode &= ~PM_SOURCE_FITTED; - return false; - case PM_MODEL_BADARGS: - source->mode |= PM_SOURCE_FAIL; - return false; - case PM_MODEL_SUCCESS: - break; - case PM_MODEL_NONCONVERGE: - case PM_MODEL_OFFIMAGE: - psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", model->params->data.F32[2], model->params->data.F32[3], model->nIter); - source->mode |= PM_SOURCE_FAIL; - return false; - default: - return false; - } - - // unless we prove otherwise, this object is a star. - source->type = PM_SOURCE_STAR; - - // if the object has fitted peak above saturation, label as SATSTAR - // this is a valid PSF object, but ignore the other quality tests - // remember: fit does not use saturated pixels (masked) - // XXX no extended object can saturate and stay extended... - if (model->params->data.F32[1] >= SATURATE) { - if (source->mode & PM_SOURCE_PSFSTAR) { - psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); - } - source->mode &= ~PM_SOURCE_PSFSTAR; - source->mode |= PM_SOURCE_SATSTAR; - return true; - } - - // if the object has a fitted peak below 0, the fit did not converge cleanly - if (model->params->data.F32[1] < 0) { - source->mode |= PM_SOURCE_FAIL; - return false; - } - - // if the source was predicted to be a SATSTAR, but it fitted below saturation, - // make a note to the user - if (source->mode & PM_SOURCE_SATSTAR) { - psLogMsg ("psphot", 5, "SATSTAR marked normal (fitted peak below saturation)\n"); - source->mode &= ~PM_SOURCE_SATSTAR; - } - - SN = model->params->data.F32[1]/model->dparams->data.F32[1]; - SX = model->params->data.F32[4]; - SY = model->params->data.F32[5]; - dSX = model->dparams->data.F32[4]; - dSY = model->dparams->data.F32[5]; - Chi = model->chisq / model->nDOF; - - // swing of sigma_x,y in sigmas - nSx = dSX / hypot (MIN_DS, 1 / (SX * SN)); - nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); - - // assign PM_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] - keep = TRUE; - keep &= (fabs(nSx) < shapeNsigma); - keep &= (fabs(nSy) < shapeNsigma); - keep &= (SN > minSN); - keep &= (Chi < maxChi); - if (keep) return true; - - // this source is not a star, unflag as PSFSTAR - // XXX : if this object was used to build the PSF, this flag should - // be set even if the object is not a star... - if (source->mode & PM_SOURCE_PSFSTAR) { - source->mode &= ~PM_SOURCE_PSFSTAR; - psLogMsg ("psphot", 5, "PSFSTAR demoted based on fit quality\n"); - } - - // object appears to be small, suspected defect - if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { - source->type = PM_SOURCE_DEFECT; - return false; - } - - // object appears to be large, suspected galaxy - if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { - source->type = PM_SOURCE_GALAXY; - return false; - } - - // poor-quality fit; only keep if nothing else works... - source->mode |= PM_SOURCE_POOR; - return false; -} Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 5992) +++ /trunk/psphot/src/psphotModelTest.c (revision 5993) @@ -6,4 +6,5 @@ bool status; int modelType; + unsigned int Nfail; float obsMag, fitMag, value; char name[64]; @@ -18,5 +19,5 @@ char *psfModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE"); if (!status) { - psfModeWord = psStringCopy ("FLT"); + psfModeWord = psStringCopy ("EXT"); } bool psfMode = !strcasecmp (psfModeWord, "PSF"); @@ -26,5 +27,6 @@ char *psfFile = psMetadataLookupPtr (&status, config, "PSF_INPUT_FILE"); if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied"); - psf = psphotReadPSF (psfFile); + psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, psfFile, FALSE); + psf = pmPSFfromMD (psfData); modelType = psf->type; } else { Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5992) +++ /trunk/psphot/src/psphotOutput.c (revision 5993) @@ -14,7 +14,11 @@ if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage); - if (psfFile != NULL) psphotWritePSF (psf, psfFile); psphotSamplePSFs (config, psf, imdata->image); + if (psfFile != NULL) { + psMetadata *psfData = pmPSFtoMD (NULL, psf); + psMetadataConfigWrite (psfData, psfFile); + psFree (psfData); + } if (outputFile == NULL) { psLogMsg ("output", 3, "no data output file selected"); @@ -166,21 +170,8 @@ float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); if (!status) ZERO_POINT = 25.0; - char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE"); // write necessary information to output header - psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); - 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, "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); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + psphotUpdateHeader (imdata->header, config); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); // create file, write-out header @@ -227,9 +218,9 @@ psLineAdd (line, "%2d ", type); psLineAdd (line, "%3.1f ", lsky); - psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal + psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); - psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' - psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' - psLineAdd (line, "%5.1f\n", 0); // should be theta + psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' + psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' + psLineAdd (line, "%5.1f\n", 0); // should be theta fwrite (line->line, 1, line->Nline, f); } @@ -256,23 +247,7 @@ float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS"); float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT"); - char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE"); // write necessary information to output header - psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); - 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, "SKYBIAS", PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias", psf->skyBias); - 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); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); - psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); - - // set NAXIS to 0 - mdi = psMetadataLookup (imdata->header, "NAXIS"); - mdi->data.S32 = 0; - mdi->type = PS_DATA_S32; + psphotUpdateHeader (imdata->header, config); table = psArrayAlloc (sources->n); @@ -299,10 +274,10 @@ psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", type); psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", lsky); - psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0); // should be 'Mgal' + psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0); psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", source->apMag + ZERO_POINT); - psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); // should be 'FHWM x' - psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); // should be 'FHWM y' - psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_TYPE_S32, "", 0); // should be theta - // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY", PS_DATA_STRING, "", NULL); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_STRING, "", "0"); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "123"); psArrayAdd (table, 100, row); @@ -310,28 +285,21 @@ } - // psAbort ("psphotOutput", "FITS table output not finished"); - - // psImage *image = psImageAlloc (1, 1, PS_TYPE_F32); - // XXX missing: theader = psFitsHeaderFromTable (NULL, table, "SMPFILE"); - // psFitsWriteImage (fits, imdata->header, image, 0); - + // write the correct number of sources in the header + psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n); + psMetadataAdd (imdata->header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true); + + // set NAXIS to 0 (we don't write out the data array) + mdi = psMetadataLookup (imdata->header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create the basic table header theader = psMetadataAlloc (); psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); - psImage *image; - - 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; - } - } - + // write out the header and table psFits *fits = psFitsOpen (filename, "w"); - // psFitsWriteHeader (imdata->header, fits); - // psFitsWriteTable (fits, theader, table); - - psFitsWriteImage (fits, NULL, image, 0); - psFitsWriteTable (fits, NULL, table); + psFitsWriteHeader (imdata->header, fits); + psFitsWriteTable (fits, theader, table); psFitsClose (fits); @@ -348,6 +316,6 @@ pmModelWritePSFs (sources, config, name, psf); - sprintf (name, "%s.flt.dat", filename); - pmModelWriteFLTs (sources, name); + sprintf (name, "%s.ext.dat", filename); + pmModelWriteEXTs (sources, name); sprintf (name, "%s.nul.dat", filename); @@ -417,5 +385,5 @@ // dump the sources to an output file -bool pmModelWriteFLTs (psArray *sources, char *filename) { +bool pmModelWriteEXTs (psArray *sources, char *filename) { double dPos, dMag; @@ -427,5 +395,5 @@ f = fopen (filename, "w"); if (f == NULL) { - psLogMsg ("pmModelWriteFLTs", 3, "can't open output file for moments%s\n", filename); + psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename); return false; } @@ -435,5 +403,5 @@ pmSource *source = (pmSource *) sources->data[i]; - if (source->type != PM_SOURCE_GALAXY) continue; + if (source->type != PM_SOURCE_EXTENDED) continue; model = pmSourceMagnitudes (source, NULL, 0.0); if (model == NULL) continue; @@ -491,7 +459,7 @@ source = sources->data[i]; - // skip these sources (in PSF or FLT) + // skip these sources (in PSF or EXT) if (source->type == PM_SOURCE_STAR) continue; - if (source->type == PM_SOURCE_GALAXY) continue; + if (source->type == PM_SOURCE_EXTENDED) continue; if (source->moments == NULL) { @@ -580,5 +548,5 @@ return (1); - case PM_SOURCE_GALAXY: + case PM_SOURCE_EXTENDED: return (2); @@ -597,5 +565,5 @@ } -psImage *pmModelPSFatXY (psImage *image, pmModel *modelFLT, pmPSF *psf, int x, int y, int dx, int dy) { +psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) { psRegion region = {x - dx, x + dx, y - dy, y + dy}; @@ -603,7 +571,7 @@ psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32); psImageInit (sample, 0); - modelFLT->params->data.F32[2] = x; - modelFLT->params->data.F32[3] = y; - pmModel *modelPSF = pmModelFromPSF (modelFLT, psf); + modelEXT->params->data.F32[2] = x; + modelEXT->params->data.F32[3] = y; + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); pmSourceAddModel (sample, NULL, modelPSF, false, false); psFree (modelPSF); @@ -625,300 +593,23 @@ if (output[0] == 0) return false; - pmModel *modelFLT = pmModelAlloc (psf->type); - modelFLT->params->data.F32[0] = 0; - modelFLT->params->data.F32[1] = 1; + pmModel *modelEXT = pmModelAlloc (psf->type); + modelEXT->params->data.F32[0] = 0; + modelEXT->params->data.F32[1] = 1; psFits *fits = psFitsOpen (output, "w"); - sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25); + sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25); psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, image->numRows - 25, 25, 25); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25); psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, 25, 25, 25); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25); psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelFLT, psf, 25, image->numRows - 25, 25, 25); + sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25); psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols / 2, image->numRows / 2, 25, 25); + sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25); psFitsWriteImage (fits, NULL, sample, 0); psFitsClose (fits); return (TRUE); -} - -psPolynomial2D *psPolynomial2DfromMD (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"); - - psPolynomial2D *poly = psPolynomial2DAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD); - - for (int nx = 0; nx < poly->nX + 1; nx++) { - for (int ny = 0; ny < poly->nY + 1; ny++) { - sprintf (keyword, "VAL_X%02d_Y%02d", nx, ny); - poly->coeff[nx][ny] = psMetadataLookupF32 (&status, folder, keyword); - if (!status) poly->mask[nx][ny] = 1; - - sprintf (keyword, "ERR_X%02d_Y%02d", nx, ny); - poly->coeffErr[nx][ny] = psMetadataLookupF32 (&status, folder, keyword); - } - } - return (poly); -} - -psPolynomial3D *psPolynomial3DfromMD (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"); - - psPolynomial3D *poly = psPolynomial3DAlloc (nXorder, nYorder, nZorder, 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++) { - sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz); - poly->coeff[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword); - if (!status) poly->mask[nx][ny][nz] = 1; - - sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz); - poly->coeffErr[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword); - } - } - } - return (poly); -} - -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, ...) { - - 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 2D 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); - - // place polynomial entries on folder - for (int nx = 0; nx < poly->nX + 1; nx++) { - for (int ny = 0; ny < poly->nY + 1; ny++) { - if (poly->mask[nx][ny]) continue; - psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny], nx, ny); - psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient error", poly->coeffErr[nx][ny], nx, ny); - } - } - return true; -} - -bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *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 3D 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); - - // 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++) { - if (poly->mask[nx][ny][nz]) continue; - psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz], nx, ny, nz); - psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz], nx, ny, nz); - } - } - } - return true; -} - -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) { - - psMetadata *psfdata = psMetadataAlloc (); - - char *modelName = pmModelGetType (psf->type); - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); - - int nPar = pmModelParameterCount (psf->type) ; - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); - - for (int i = 0; i < nPar - 4; i++) { - psPolynomial2D *poly = psf->params->data[i]; - psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i); - } - psPolynomial4DtoMD (psfdata, psf->ApTrend, "APTREND"); - - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid); - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); - - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq); - psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars); - - psMetadataConfigWrite (psfdata, filename); - psFree (psfdata); - - return true; -} - -pmPSF *psphotReadPSF (char *filename) { - - bool status; - unsigned int Nfail; - char keyword[80]; - - psMetadata *psfdata = psMetadataConfigParse (NULL, &Nfail, filename, FALSE); - - char *modelName = psMetadataLookupPtr (&status, psfdata, "PSF_MODEL_NAME"); - pmModelType type = pmModelSetType (modelName); - - pmPSF *psf = pmPSFAlloc (type); - - int nPar = psMetadataLookupS32 (&status, psfdata, "PSF_MODEL_NPAR"); - if (nPar != pmModelParameterCount (psf->type)) psAbort ("read PSF" , "mismatch model par count"); - - for (int i = 0; i < nPar - 4; i++) { - sprintf (keyword, "PSF_PAR%02d", i); - psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword); - psPolynomial2D *poly = psPolynomial2DfromMD (folder); - psFree (psf->params->data[i]); - psf->params->data[i] = poly; - } - sprintf (keyword, "APTREND"); - psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword); - psPolynomial4D *poly = psPolynomial4DfromMD (folder); - psFree (psf->ApTrend); - psf->ApTrend = poly; - - psf->ApResid = psMetadataLookupF32 (&status, psfdata, "PSF_AP_RESID"); - psf->dApResid = psMetadataLookupF32 (&status, psfdata, "PSF_dAP_RESID"); - psf->skyBias = psMetadataLookupF32 (&status, psfdata, "PSF_SKY_BIAS"); - - psf->chisq = psMetadataLookupF32 (&status, psfdata, "PSF_CHISQ"); - psf->nPSFstars = psMetadataLookupS32 (&status, psfdata, "PSF_NSTARS"); - - psFree (psfdata); - return (psf); } @@ -937,2 +628,31 @@ return true; } + +bool psphotUpdateHeader (psMetadata *header, psMetadata *config) { + + // write necessary information to output header + psMetadataItemTransfer (header, config, "NSTARS"); + psMetadataItemTransfer (header, config, "ZERO_PT"); + psMetadataItemTransfer (header, config, "APMIFIT"); + psMetadataItemTransfer (header, config, "dAPMIFIT"); + psMetadataItemTransfer (header, config, "SKYBIAS"); + psMetadataItemTransfer (header, config, "PHOTCODE"); + + // write necessary information to output header + psMetadataItemTransfer (header, config, "NPSFSTAR"); + psMetadataItemTransfer (header, config, "SKYBIAS"); + psMetadataItemTransfer (header, config, "SKYSAT"); + psMetadataItemTransfer (header, config, "APMIFIT"); + psMetadataItemTransfer (header, config, "DAPMIFIT"); + psMetadataItemTransfer (header, config, "NAPMIFIT"); + psMetadataItemTransfer (header, config, "APLOSS"); + + // XXX these need to be defined from elsewhere + psMetadataAdd (header, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "ANGLE", PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); + psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); + + return true; +} Index: /trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- /trunk/psphot/src/psphotRadiusChecks.c (revision 5992) +++ /trunk/psphot/src/psphotRadiusChecks.c (revision 5993) @@ -31,62 +31,35 @@ } -# 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 } +static float EXT_OUTER_RADIUS; +static float EXT_FIT_PADDING; +static float EXT_FLUX_LIMIT; +static pmModelRadius modelRadiusEXT; -static float GAL_OUTER_RADIUS; -static float GAL_FIT_PADDING; -static float GAL_FLUX_LIMIT; -static pmModelRadius modelRadiusGAL; - -bool psphotInitRadiusFLT (psMetadata *config, psStats *sky, pmModelType type) { +bool psphotInitRadiusEXT (psMetadata *config, psStats *sky, pmModelType type) { bool status; - float GAL_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA"); - GAL_OUTER_RADIUS = psMetadataLookupF32 (&status, config, "GAL_OUTER_RADIUS"); - GAL_FIT_PADDING = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING"); - GAL_FLUX_LIMIT = GAL_FIT_NSIGMA * sky->sampleStdev; + float EXT_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "EXT_FIT_NSIGMA"); + EXT_OUTER_RADIUS = psMetadataLookupF32 (&status, config, "EXT_OUTER_RADIUS"); + EXT_FIT_PADDING = psMetadataLookupF32 (&status, config, "EXT_FIT_PADDING"); + EXT_FLUX_LIMIT = EXT_FIT_NSIGMA * sky->sampleStdev; // this function specifies the radius at this the model hits the given flux - modelRadiusGAL = pmModelRadius_GetFunction (type); + modelRadiusEXT = pmModelRadius_GetFunction (type); return true; } -bool psphotCheckRadiusFLT (eamReadout *imdata, pmSource *source, pmModel *model) { +bool psphotCheckRadiusEXT (eamReadout *imdata, pmSource *source, pmModel *model) { // set the fit radius based on the object flux limit and the model - model->radius = modelRadiusGAL (model->params, GAL_FLUX_LIMIT) + GAL_FIT_PADDING; - if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); + model->radius = modelRadiusEXT (model->params, EXT_FLUX_LIMIT) + EXT_FIT_PADDING; + if (isnan(model->radius)) psAbort (__func__, "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) { - // (re)-allocate image, weight, mask arrays for each peak (square of radius OUTER) - // do I need to save the old pixel and mask data? - psphotDefinePixels (source, imdata, model->params->data.F32[2], model->params->data.F32[3], model->radius); - return true; - } - return false; -# endif } Index: unk/psphot/src/psphotReapplyPSF.c =================================================================== --- /trunk/psphot/src/psphotReapplyPSF.c (revision 5992) +++ (revision ) @@ -1,68 +1,0 @@ -# include "psphot.h" - -bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) -{ - bool status; - float x, y, r; - int Nfit = 0; - int Nsub = 0; - int Niter = 0; - - psTimerStart ("psphot"); - - // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, psphotSortBySN); - - // we may set this differently here from the value used to mark likely saturated stars - float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); - float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); - float PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); - float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); - - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // skip non-astronomical objects (very likely defects) - // XXX EAM : should we try these anyway? - if (source->mode |= PM_SOURCE_BLEND) continue; - if (source->type == PM_SOURCE_DEFECT) continue; - if (source->type == PM_SOURCE_SATURATED) continue; - - pmModel *model = source->modelPSF; - pmModel *save = pmModelCopy (model); - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - r = model->radius; - - // mask outside radius to limit the following actions to the source - psImageKeepCircle (source->mask, x, y, r, "OR", PSPHOT_MASK_MARKED); - - // replace object in image - pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false); - - // fit PSF model (set/unset the pixel mask) - pmSourceFitModel (source, model, true); - - // model succeeded : tentatively keep it - Niter += model[0].nIter; - Nfit ++; - - // is it a good PSF model? - if (psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION)) { - psFree (save); - Nsub ++; - } else { - psFree (model); - source->modelPSF = save; - } - pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); - - // clear the masked region - psImageKeepCircle (source->mask, x, y, r, "AND", ~PSPHOT_MASK_MARKED); - } - - psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); - psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub); - return (true); -} Index: /trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- /trunk/psphot/src/psphotReplaceUnfit.c (revision 5992) +++ /trunk/psphot/src/psphotReplaceUnfit.c (revision 5993) @@ -20,25 +20,2 @@ return true; } - -# if (0) -if (source->mode & PM_SOURCE_BLEND) continue; -if (source->type == PM_SOURCE_SATURATED) continue; - -if (source->type == PM_SOURCE_DEFECT) { - if (source->mode & PM_SOURCE_SUBTRACTED) goto addPSF; - continue; -} - -if (source->type == PM_SOURCE_GALAXY) { - if (source->mode & PM_SOURCE_FAIL) goto addPSF; - if (source->mode & PM_SOURCE_POOR) goto addPSF; - continue; -} - -// need to skip the successful fits -if (source->type == PM_SOURCE_STAR) { - if (source->mode & PM_SOURCE_FAIL) goto addPSF; - if (source->mode & PM_SOURCE_POOR) goto addPSF; - continue; -} -# endif Index: /trunk/psphot/src/psphotRoughClass.c =================================================================== --- /trunk/psphot/src/psphotRoughClass.c (revision 5992) +++ /trunk/psphot/src/psphotRoughClass.c (revision 5993) @@ -9,5 +9,5 @@ psfClump = pmSourcePSFClump (sources, config); - // group into STAR, COSMIC, GALAXY, SATURATED, etc. + // group into STAR, COSMIC, EXTENDED, SATURATED, etc. pmSourceRoughClass (sources, config, psfClump); Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 5992) +++ /trunk/psphot/src/psphotSourceFits.c (revision 5993) @@ -41,20 +41,20 @@ -static float GAL_MIN_SN; -static float GAL_MOMENTS_RAD; -static pmModelType modelTypeFLT; - -bool psphotInitLimitsFLT (psMetadata *config, psStats *sky) { +static float EXT_MIN_SN; +static float EXT_MOMENTS_RAD; +static pmModelType modelTypeEXT; + +bool psphotInitLimitsEXT (psMetadata *config, psStats *sky) { bool status; - // galaxy model parameters - GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); - GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); - - // 'galaxy' model descriptions - char *modelNameFLT = psMetadataLookupPtr (&status, config, "GAL_MODEL"); - modelTypeFLT = pmModelSetType (modelNameFLT); - psphotInitRadiusFLT (config, sky, modelTypeFLT); + // extended source model parameters + EXT_MIN_SN = psMetadataLookupF32 (&status, config, "EXT_MIN_SN"); + EXT_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "EXT_MOMENTS_RADIUS"); + + // extended source model descriptions + char *modelNameEXT = psMetadataLookupPtr (&status, config, "EXT_MODEL"); + modelTypeEXT = pmModelSetType (modelNameEXT); + psphotInitRadiusEXT (config, sky, modelTypeEXT); return true; @@ -63,6 +63,6 @@ bool psphotFitBlob (eamReadout *imdata, pmSource *source, psArray *sources) { - bool okFLT, okDBL; - float chiFLT, chiDBL; + bool okEXT, okDBL; + float chiEXT, chiDBL; // skip the source if we don't think it is extended @@ -70,8 +70,8 @@ if (source->type == PM_SOURCE_DEFECT) return false; if (source->type == PM_SOURCE_SATURATED) return false; - if (source->moments->SN < GAL_MIN_SN) return false; - - // recalculate the source moments using the larger galaxy moments radius - if (!pmSourceMoments (source, GAL_MOMENTS_RAD)) return false; + if (source->moments->SN < EXT_MIN_SN) return false; + + // recalculate the source moments using the larger extended-source moments radius + if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false; psTrace ("psphot.blend", 5, "trying blob...\n"); @@ -79,7 +79,7 @@ pmSource *tmpSrc = pmSourceAlloc (); - pmModel *FLT = psphotFitFLT (imdata, source); - okFLT = psphotEvalFLT (tmpSrc, FLT); - chiFLT = FLT->chisq / FLT->nDOF; + pmModel *EXT = psphotFitEXT (imdata, source); + okEXT = psphotEvalEXT (tmpSrc, EXT); + chiEXT = EXT->chisq / EXT->nDOF; psArray *DBL = psphotFitDBL (imdata, source); @@ -89,22 +89,22 @@ chiDBL = ONE->chisq / ONE->nDOF; - if (okFLT && okDBL) { - psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiFLT, chiDBL); + if (okEXT && okDBL) { + psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL); // XXX EAM : a bogus bias: need to examine this better - if (3*chiFLT > chiDBL) goto keepDBL; - goto keepFLT; - } - - if (okFLT && !okDBL) goto keepFLT; - if (!okFLT && okDBL) goto keepDBL; + if (3*chiEXT > chiDBL) goto keepDBL; + goto keepEXT; + } + + if (okEXT && !okDBL) goto keepEXT; + if (!okEXT && okDBL) goto keepDBL; return false; -keepFLT: - // sub FLT - pmSourceSubModel (source->pixels, source->mask, FLT, false, false); - psTrace ("psphot.blend", 5, "blob as FLT: %f %f\n", FLT->params->data.F32[2], FLT->params->data.F32[3]); +keepEXT: + // sub EXT + pmSourceSubModel (source->pixels, source->mask, EXT, false, false); + psTrace ("psphot.blend", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[2], EXT->params->data.F32[3]); // save new model - source->modelFLT = FLT; + source->modelEXT = EXT; source->mode |= PM_SOURCE_SUBTRACTED; source->mode &= ~PM_SOURCE_TEMPSUB; @@ -168,5 +168,5 @@ y = PSF->params->data.F32[3]; - // fit FLT (not PSF) model (set/unset the pixel mask) + // fit EXT (not PSF) model (set/unset the pixel mask) psImageKeepCircle (source->mask, x, y, PSF->radius, "OR", PSPHOT_MASK_MARKED); pmSourceFitSet (source, modelSet, true); @@ -176,22 +176,22 @@ } -pmModel *psphotFitFLT (eamReadout *imdata, pmSource *source) { +pmModel *psphotFitEXT (eamReadout *imdata, pmSource *source) { float x, y; // use the source moments, etc to guess basic model parameters - pmModel *FLT = pmSourceModelGuess (source, modelTypeFLT); - - psphotCheckRadiusFLT (imdata, source, FLT); - - x = FLT->params->data.F32[2]; - y = FLT->params->data.F32[3]; - - // fit FLT (not PSF) model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, FLT->radius, "OR", PSPHOT_MASK_MARKED); - pmSourceFitModel (source, FLT, false); - psImageKeepCircle (source->mask, x, y, FLT->radius, "AND", ~PSPHOT_MASK_MARKED); - - return (FLT); + pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); + + psphotCheckRadiusEXT (imdata, source, EXT); + + x = EXT->params->data.F32[2]; + y = EXT->params->data.F32[3]; + + // fit EXT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, EXT->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceFitModel (source, EXT, false); + psImageKeepCircle (source->mask, x, y, EXT->radius, "AND", ~PSPHOT_MASK_MARKED); + + return (EXT); } Index: unk/psphot/src/psphotSparseMatrix.c =================================================================== --- /trunk/psphot/src/psphotSparseMatrix.c (revision 5992) +++ (revision ) @@ -1,116 +1,0 @@ -# include "psphot.h" - -bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { - - bool status; - float x; - float y; - float Sky; - int Nfit = 0; - - psTimerStart ("psphot"); - - // source analysis is done in spatial order - sources = psArraySort (sources, psphotSortByY); - - 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"); - - // set the object surface-brightness limit for fitted pixels - float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; - psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); - - // vector for model input coordinates - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - - // this function specifies the radius at this the model hits the given flux - pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); - - // this is the function to get the model function - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - - // pre-calculate all model pixels - psArray *pixlist = psArrayAlloc (sources->n); - for (int k = 0; k < sources->n; k++) { - source = sources->data[k]; - - // use the source moments, etc to guess basic model parameters - pmModel *modelFLT = pmSourceModelGuess (source, psf->type); - - // set PSF parameters for this model - pmModel *model = pmModelFromPSF (modelFLT, psf); - // model->params->data.F32[2] = source->peak->x; - // model->params->data.F32[3] = source->peak->y; - model->params->data.F32[1] = 1.0; - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - psFree (modelFLT); - - // set the fit radius based on the object flux limit and the model - // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) - model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; - if (isnan(model->radius)) { - psAbort ("apply_psf_model", "error in radius"); - } - - // build the model image - nx = src->pixels->numCols; - ny = src->pixels->numRows; - im = psImageAlloc (nx, ny, PS_DATA_F32); - mk = psImageAlloc (nx, ny, PS_DATA_U8); - - // fill in the model pixel values - psImageKeepCircle (mk, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); - psVector *params = model->params; - for (int i = 0; i < nx; i++) { - for (int j = 0; j < ny; j++) { - if (mk->data.U8[j][i]) continue; - x->data.F32[0] = (float) i + pix->col0; - x->data.F32[1] = (float) j + pix->row0; - im->data.F32[j][i] = modelFunc (NULL, params, x); - } - } - pixlist->data[k] = im; - } - - // fill out the sparse matrix - psSparse *sparse = psSparseAlloc (sources->n, 100); - psVector *Aij = sparse->Aij; - psVector *Bfj = sparse->Bfj; - psVector *Si = sparse->Si; - psVector *Sj = sparse->Sj; - - int Ri, Rj; - - Nelem = 0; - for (i = 0; i < pixlist->n; i++) { - pf = ((psSource *)sources->data[k])->pixels; - pi = pixlist->data[i]; - - // diagonal element (auto-cross-product) - f = psphotCrossProduct (pi, pi); - psSparseMatrixElement (sparse, i, i, f); - - // find the image x model value - f = psphotCrossProduct (pi, pf); - psSparseVectorElement (sparse, i, f); - - // loop over all other stars following this one - for (j = i + 1; j < sources->n; j++) { - pj = ((pmSource *)sources->data[j])->pixels; - - if (pi->col0 + pi->numCols < pj->col0) break; - if (pj->col0 + pj->numCols < pi->col0) continue; - if (pi->row0 + pi->numRows < pj->row0) continue; - if (pj->row0 + pj->numRows < pi->row0) continue; - - // got an overlap; calculate cross-product and add to output array - f = psphotCrossProduct (pi, pj); - psSparseMatrixElement (sparse, j, i, f); - } - } - - psSparseResort (sparse); - return (match); -} Index: unk/psphot/src/psphotSubtractPSF.c =================================================================== --- /trunk/psphot/src/psphotSubtractPSF.c (revision 5992) +++ (revision ) @@ -1,29 +1,0 @@ -# include "psphot.h" - -// subtract model fits from image (only PSFSTAR && BRIGHTSTAR) -// need to control application of this with some thresholding - -bool psphotSubtractPSF (pmSource *source) -{ - pmModel *model; - - // only subtract successful fits - if (source->type == PM_SOURCE_SATSTAR) goto valid; - if (source->type == PM_SOURCE_PSFSTAR) goto valid; - if (source->type == PM_SOURCE_GOODSTAR) goto valid; - return (false); - -valid: - model = source->modelPSF; - if (model == NULL) { - psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y); - return (false); - } - - // subtract object, leave local sky - pmSourceSubModel (source->pixels, source->mask, model, false, false); - - // XXX EAM : amplify the weight matrix a la dophot? - - return (true); -} Index: /trunk/psphot/src/psphotTest.c =================================================================== --- /trunk/psphot/src/psphotTest.c (revision 5992) +++ /trunk/psphot/src/psphotTest.c (revision 5993) @@ -1,4 +1,3 @@ # include "psphot.h" -bool apResidSetPower (float value); void psExit (int status, char *process, char *format, ...) { @@ -16,51 +15,52 @@ int main (int argc, char **argv) { - // read in a file consisting of: - // x y radius fitmag apmag dmag + psMetadata *row; + psArray *table; - double fit; - double ap; + psMetadataItem *mdi; + psRegion region = {0,0,0,0}; // a region representing the entire array psphotTestArguments (&argc, argv); - FILE *f = fopen (argv[1], "r"); - if (f == NULL) psExit (2, "test", "no file\n"); + psFits *file = psFitsOpen (argv[1], "r"); + psMetadata *header = psFitsReadHeader (NULL, file); + psImage *image = psFitsReadImage (NULL, file, region, 0); + psFitsClose (file); - float power = atof (argv[2]); - apResidSetPower (power); + psMetadataConfigWrite (header, argv[2]); - 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); + // attempt to write image with NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create a test image + // psImage *tmpimage = psImageAlloc (10, 10, PS_DATA_F32); - 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) { + // create a test table + table = psArrayAlloc (10); + table->n = 0; - if (fabs(ap-fit) > 0.1) continue; - if (fit > -10) continue; + for (int i = 0; i < 10; i++) { + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "ROW", PS_DATA_S32, "", i); + psMetadataAdd (row, PS_LIST_TAIL, "FROW", PS_TYPE_F32, "", 0.1*i); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "test line"); + + table->data[i] = row; + } + table->n = 10; + + psMetadata *theader = psMetadataAlloc (); + psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); - apres->data.F64[Npt] = ap-fit; - rflux->data.F64[Npt] = pow(10.0, 0.4*fit); + psFits *fits = psFitsOpen (argv[3], "w"); + // psFitsWriteImage (fits, header, tmpimage, 0); + psFitsWriteHeader (header, fits); + psFitsWriteTable (fits, theader, table); + psFitsClose (fits); - 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; - - psphotApResidPolyFit (x, y, rad, rflux, apres, dmag); - + psFree (header); + psFree (image); exit (0); } Index: /trunk/psphot/src/psphotTestArguments.c =================================================================== --- /trunk/psphot/src/psphotTestArguments.c (revision 5992) +++ /trunk/psphot/src/psphotTestArguments.c (revision 5993) @@ -8,5 +8,5 @@ psArgumentVerbosity (argc, argv); - if (*argc != 3) usage (); + if (*argc != 4) usage (); return; @@ -15,5 +15,5 @@ static int usage () { - fprintf (stderr, "USAGE: psphotModelTest (apresid.dat) (power)\n"); + fprintf (stderr, "USAGE: psphotTest (input.fits) (output.hdr) (out.fits)\n"); exit (2); }