Index: trunk/psModules/src/objects/Makefile.am =================================================================== --- trunk/psModules/src/objects/Makefile.am (revision 6511) +++ trunk/psModules/src/objects/Makefile.am (revision 6872) @@ -1,12 +1,26 @@ noinst_LTLIBRARIES = libpsmoduleobjects.la -libpsmoduleobjects_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) +libpsmoduleobjects_la_CPPFLAGS = $(SRCINC) $(PSMODULE_CFLAGS) -I../pslib/ libpsmoduleobjects_la_LDFLAGS = -release $(PACKAGE_VERSION) libpsmoduleobjects_la_SOURCES = \ - pmObjects.c \ - pmPSF.c \ - pmPSFtry.c \ - pmModelGroup.c \ - psEllipse.c + pmPeaks.c \ + pmMoments.c \ + pmModel.c \ + pmModelGroup.c \ + pmSource.c \ + pmSourceSky.c \ + pmSourceContour.c \ + pmSourceFitModel.c \ + pmSourcePhotometry.c \ + pmSourceIO.c \ + pmSourceIO_CMF.c \ + pmSourceIO_CMP.c \ + pmSourceIO_OBJ.c \ + pmSourceIO_SX.c \ + pmSourceIO_RAW.c \ + pmPSF.c \ + pmPSF_IO.c \ + pmPSFtry.c \ + pmGrowthCurve.c EXTRA_DIST = \ @@ -18,7 +32,16 @@ psmoduleincludedir = $(includedir) psmoduleinclude_HEADERS = \ - pmObjects.h \ - pmPSF.h \ - pmPSFtry.h \ - pmModelGroup.h \ - psEllipse.h + pmPeaks.h \ + pmMoments.h \ + pmModel.h \ + pmModelGroup.h \ + pmSource.h \ + pmSourceSky.h \ + pmSourceContour.h \ + pmSourceFitModel.h \ + pmSourcePhotometry.h \ + pmSourceIO.h \ + pmPSF.h \ + pmPSF_IO.h \ + pmPSFtry.h \ + pmGrowthCurve.h Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 6511) +++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c (revision 6872) @@ -26,13 +26,13 @@ if (deriv != NULL) { - // note difference from a pure gaussian: q = PAR[1]*r + psF32 *dPAR = deriv->data.F32; psF32 q = PAR[1]*r*r*t; - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; + dPAR[0] = +1.0; + dPAR[1] = +r; + dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + dPAR[4] = -2.0*q*px*X; + dPAR[5] = -2.0*q*py*Y; + dPAR[6] = -q*X*Y; } return(f); @@ -50,5 +50,5 @@ beta_lim[0][0].data.F32[0] = 1000; - beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[1] = 3e6; beta_lim[0][0].data.F32[2] = 5; beta_lim[0][0].data.F32[3] = 5; @@ -66,5 +66,5 @@ params_max[0][0].data.F32[0] = 1e5; - params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[1] = 1e8; params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! @@ -133,4 +133,5 @@ params[5] = 1.2 / moments->Sy; params[6] = 0.0; + return(true); } Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 6511) +++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 6872) @@ -26,21 +26,27 @@ psF32 py = PAR[5]*Y; psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; - - psF32 r = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); - psF32 f = PAR[1]*r + PAR[0]; + psF32 zp = pow(z,1.25); + + psF32 r = 1.0 / (1 + PAR[7]*z + z*zp); + // test: psF32 r = 1.0 / (1 + PAR[7]*z + PS_SQR(z)); + psF32 r1 = PAR[1]*r; + psF32 f = r1 + PAR[0]; if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7] + 2.25*pow(z, 1.25)); - - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); - deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X); - deriv->data.F32[4] = -2.0*q*px*X; - deriv->data.F32[5] = -2.0*q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -t*z; + psF32 t = r1*r; + psF32 q = t*(PAR[7] + 2.25*zp); + // test: psF32 q = t*(PAR[7] + 2*z); + + dPAR[0] = +1.0; + dPAR[1] = +r; + dPAR[2] = q*(2.0*px*PAR[4] + PAR[6]*Y); + dPAR[3] = q*(2.0*py*PAR[5] + PAR[6]*X); + dPAR[4] = -2.0*q*px*X; + dPAR[5] = -2.0*q*py*Y; + dPAR[6] = -q*X*Y; + dPAR[7] = -t*z; } return(f); @@ -58,5 +64,5 @@ beta_lim[0][0].data.F32[0] = 1000; - beta_lim[0][0].data.F32[1] = 10000; + beta_lim[0][0].data.F32[1] = 3e6; beta_lim[0][0].data.F32[2] = 5; beta_lim[0][0].data.F32[3] = 5; @@ -76,5 +82,5 @@ params_max[0][0].data.F32[0] = 1e5; - params_max[0][0].data.F32[1] = 1e6; + params_max[0][0].data.F32[1] = 1e8; params_max[0][0].data.F32[2] = 1e4; // this should be set by image dimensions! params_max[0][0].data.F32[3] = 1e4; // this should be set by image dimensions! @@ -102,4 +108,5 @@ params[6] = 0.0; params[7] = 1.0; + return(true); } @@ -119,9 +126,10 @@ // the area needs to be multiplied by the integral of f(z) norm = 0.0; - for (z = 0.005; z < 50; z += 0.01) { + for (z = 0.05; z < 50; z += 0.1) { f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + // test: f = 1.0 / (1 + PAR[7]*z + PS_SQR(z)); norm += f; } - norm *= 0.01; + norm *= 0.1; psF64 Flux = PAR[1] * Area * norm; @@ -150,6 +158,7 @@ // we can do this much better with intelligent choices here - for (z = 0.0; z < 20.0; z += dz) { + for (z = 0.0; z < 30.0; z += dz) { f = 1.0 / (1 + PAR[7]*z + pow(z, 2.25)); + // test: f = 1.0 / (1 + PAR[7]*z + PS_SQR(z)); if (f < limit) break; @@ -201,6 +210,6 @@ status &= ((dPAR[1]/PAR[1]) < 0.5); - if (status) - return true; - return false; -} + if (!status) + return false; + return true; +} Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 6511) +++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c (revision 6872) @@ -117,7 +117,7 @@ psVector *params = model->params; - EllipseAxes axes; - EllipseShape shape; - EllipseMoments moments; + psEllipseAxes axes; + psEllipseShape shape; + psEllipseMoments moments; moments.x2 = PS_SQR(source->moments->Sx); @@ -125,6 +125,6 @@ moments.xy = source->moments->Sxy; - axes = EllipseMomentsToAxes(moments); - shape = EllipseAxesToShape(axes); + axes = psEllipseMomentsToAxes(moments); + shape = psEllipseAxesToShape(axes); params->data.F32[0] = source->moments->Sky; Index: trunk/psModules/src/objects/models/pmModel_SGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_SGAUSS.c (revision 6511) +++ trunk/psModules/src/objects/models/pmModel_SGAUSS.c (revision 6872) @@ -17,5 +17,5 @@ # define SQ(A)((A)*(A)) -psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image); +psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image); psF64 p_psImageGetElementF64(psImage *a, int i, int j); @@ -102,5 +102,5 @@ // measure the flux for the elliptical contour -psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image) +psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image) { @@ -149,7 +149,7 @@ psF32 *params = model->params->data.F32; - EllipseAxes axes; - EllipseShape shape; - EllipseMoments moments; + psEllipseAxes axes; + psEllipseShape shape; + psEllipseMoments moments; moments.x2 = PS_SQR(sMoments->Sx); @@ -158,6 +158,6 @@ // solve the math to go from Moments To Shape - axes = EllipseMomentsToAxes(moments); - shape = EllipseAxesToShape(axes); + axes = psEllipseMomentsToAxes(moments); + shape = psEllipseAxesToShape(axes); params[0] = sMoments->Sky; @@ -199,7 +199,7 @@ float f1, f2; - EllipseAxes axes; - EllipseShape shape; - EllipseMoments moments; + psEllipseAxes axes; + psEllipseShape shape; + psEllipseMoments moments; moments.x2 = PS_SQR(sMoments->Sx); @@ -208,6 +208,6 @@ // solve the math to go from Moments To Shape - axes = EllipseMomentsToAxes(moments); - shape = EllipseAxesToShape(axes); + axes = psEllipseMomentsToAxes(moments); + shape = psEllipseAxesToShape(axes); params[0] = sMoments->Sky; @@ -265,9 +265,9 @@ psF64 pmModelRadius_SGAUSS (const psVector *params, psF64 flux) { - psF64 r, z, pr, f; + psF64 r, z = 0.0, pr, f; psF32 *PAR = params->data.F32; - EllipseAxes axes; - EllipseShape shape; + psEllipseAxes axes; + psEllipseShape shape; if (flux <= 0) @@ -283,5 +283,5 @@ shape.sxy = PAR[6]; - axes = EllipseShapeToAxes (shape); + axes = psEllipseShapeToAxes (shape); psF64 dr = 1.0 / axes.major; psF64 limit = flux / PAR[1]; @@ -327,6 +327,6 @@ psF32 dP; bool status; - EllipseAxes axes; - EllipseShape shape; + psEllipseAxes axes; + psEllipseShape shape; psF32 *PAR = model->params->data.F32; @@ -337,5 +337,5 @@ shape.sxy = PAR[6]; - axes = EllipseShapeToAxes (shape); + axes = psEllipseShapeToAxes (shape); dP = 0; Index: trunk/psModules/src/objects/models/pmModel_ZGAUSS.c =================================================================== --- trunk/psModules/src/objects/models/pmModel_ZGAUSS.c (revision 6511) +++ trunk/psModules/src/objects/models/pmModel_ZGAUSS.c (revision 6872) @@ -85,6 +85,6 @@ psF32 *PAR = params->data.F32; - EllipseAxes axes; - EllipseShape shape; + psEllipseAxes axes; + psEllipseShape shape; if (flux <= 0) @@ -100,5 +100,5 @@ shape.sxy = PAR[6]; - axes = EllipseShapeToAxes (shape); + axes = psEllipseShapeToAxes (shape); psF64 dr = 1.0 / axes.major; psF64 limit = flux / PAR[1]; @@ -124,7 +124,7 @@ psVector *params = model->params; - EllipseAxes axes; - EllipseShape shape; - EllipseMoments moments; + psEllipseAxes axes; + psEllipseShape shape; + psEllipseMoments moments; moments.x2 = PS_SQR(source->moments->Sx); @@ -132,6 +132,6 @@ moments.xy = source->moments->Sxy; - axes = EllipseMomentsToAxes(moments); - shape = EllipseAxesToShape(axes); + axes = psEllipseMomentsToAxes(moments); + shape = psEllipseAxesToShape(axes); params->data.F32[0] = source->moments->Sky; Index: trunk/psModules/src/objects/pmGrowthCurve.c =================================================================== --- trunk/psModules/src/objects/pmGrowthCurve.c (revision 6872) +++ trunk/psModules/src/objects/pmGrowthCurve.c (revision 6872) @@ -0,0 +1,130 @@ +/** @file pmGrowthCurve.c + * + * Measure the curve-of-growth for sources + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include "pmGrowthCurve.h" + +static void pmGrowthCurveFree (pmGrowthCurve *growth) +{ + + if (growth == NULL) + return; + + psFree (growth->radius); + psFree (growth->apMag); + return; +} + +pmGrowthCurve *pmGrowthCurveAlloc (psF32 minRadius, psF32 maxRadius, psF32 refRadius) +{ + + pmGrowthCurve *growth = psAlloc (sizeof(pmGrowthCurve)); + psMemSetDeallocator(growth, (psFreeFunc) pmGrowthCurveFree); + + // fractional pixel radii are not well defined; use integer steps + growth->radius = psVectorCreate (NULL, minRadius, maxRadius, 1.0, PS_TYPE_F32); + growth->apMag = psVectorAlloc (growth->radius->n, PS_TYPE_F32); + + // XXX may want to extend this to allow for a different refRadius; + growth->refRadius = refRadius; + growth->maxRadius = maxRadius; + growth->apLoss = 0.0; + growth->fitMag = 0.0; + return growth; +} + +psF32 pmGrowthCurveCorrect (pmGrowthCurve *growth, psF32 radius) +{ + + float apRad = psVectorInterpolate (growth->radius, growth->apMag, radius); + float apCor = growth->apRef - apRad; + return apCor; +} + +// return the last entry below or first entry above key value +int psVectorBracket (psVector *index, psF32 key, bool above) +{ + + int N; + int Nlo = 0; + int Nhi = index->n; + + if (above) { + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (index->data.F32[N] > key) { + Nhi = N; + } else { + Nlo = N - 1; + } + } + // at this point, index[Nhi] > key >= index[Nlo] + N = Nlo; + while ((index->data.F32[N] <= key) && (N < Nhi)) { + N++; + } + return (N); + } + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (index->data.F32[N] < key) { + Nlo = N; + } else { + Nhi = N + 1; + } + } + // at this point, index[Nhi] >= key > index[Nlo] + N = Nhi; + while ((index->data.F32[N] >= key) && (N > Nlo)) { + N--; + } + return (N); +} + +// search for the bins bounding key in index, interpolate the corresponding values +psF32 psVectorInterpolate (psVector *index, psVector *value, psF32 key) +{ + + int n0 = 0; + int n1 = 0; + + // extrapolate at ends + if (key < index->data.F32[0]) { + n0 = 0; + n1 = 1; + } + + // extrapolate at ends + if (key > index->data.F32[index->n-1]) { + n0 = index->n-2; + n1 = index->n-1; + } + + if (n1 == 0) { + n0 = psVectorBracket (index, key, FALSE); + n1 = n0 + 1; + } + + if (n0 == index->n-1) { + n1 = n0; + n0 = n1 - 1; + } + + float dy = value->data.F32[n1] - value->data.F32[n0]; + float dx = index->data.F32[n1] - index->data.F32[n0]; + float dX = key - index->data.F32[n0]; + float dY = dX * (dy/dx); + float result = value->data.F32[n0] + dY; + return result; +} + Index: trunk/psModules/src/objects/pmGrowthCurve.h =================================================================== --- trunk/psModules/src/objects/pmGrowthCurve.h (revision 6872) +++ trunk/psModules/src/objects/pmGrowthCurve.h (revision 6872) @@ -0,0 +1,24 @@ + +# ifndef PM_GROWTH_CURVE_H +# define PM_GROWTH_CURVE_H + +typedef struct +{ + psVector *radius; + psVector *apMag; + psF32 refRadius; + psF32 maxRadius; + psF32 fitMag; + psF32 apRef; // apMag[refRadius] + psF32 apLoss; // fitMag - apRef +} +pmGrowthCurve; + +pmGrowthCurve *pmGrowthCurveAlloc (psF32 minRadius, psF32 maxRadius, psF32 refRadius); +int psVectorBracket (psVector *index, psF32 key, bool above); +psF32 psVectorInterpolate (psVector *index, psVector *value, psF32 key); +psF32 pmGrowthCurveCorrect (pmGrowthCurve *growth, psF32 radius); + +// XXX psVectorBracket,Interpolate should be put in pslib + +# endif /* PM_GROWTH_CURVE_H */ Index: trunk/psModules/src/objects/pmModel.c =================================================================== --- trunk/psModules/src/objects/pmModel.c (revision 6872) +++ trunk/psModules/src/objects/pmModel.c (revision 6872) @@ -0,0 +1,225 @@ +/** @file pmModel.c + * + * Functions to define and manipulate object models + * + * @author GLG, MHPCC + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmModel.h" + +/* The following types and functions need to be known by pmModel.c and are defined in + pmModelGroup.h. The use of pmSource and other types in pmModelGroup.h prevent us from + directly including it here ***/ + +typedef psMinimizeLMChi2Func pmModelFunc; +typedef bool (*pmModelFitStatusFunc)(pmModel *model); +pmModelFunc pmModelFunc_GetFunction (pmModelType type); +pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); +int pmModelParameterCount(pmModelType type); + +static void modelFree(pmModel *tmp) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + psFree(tmp->params); + psFree(tmp->dparams); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); +} + +/****************************************************************************** +pmModelAlloc(): Allocate the pmModel structure, along with its parameters, +and initialize the type member. Initialize the params to 0.0. +XXX EAM: simplifying code with pmModelParameterCount +*****************************************************************************/ +pmModel *pmModelAlloc(pmModelType type) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); + + tmp->type = type; + tmp->chisq = 0.0; + tmp->nIter = 0; + tmp->radiusTMP = 0; + tmp->status = PM_MODEL_UNTRIED; + + psS32 Nparams = pmModelParameterCount(type); + if (Nparams == 0) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return(NULL); + } + + tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); + + for (psS32 i = 0; i < tmp->params->n; i++) { + tmp->params->data.F32[i] = 0.0; + tmp->dparams->data.F32[i] = 0.0; + } + + psMemSetDeallocator(tmp, (psFreeFunc) modelFree); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +// copy model to a new structure +pmModel *pmModelCopy (pmModel *model) +{ + + pmModel *new = pmModelAlloc (model->type); + + new->chisq = model->chisq; + new->nDOF = model->nDOF; + new->nIter = model->nIter; + new->status = model->status; + new->radiusTMP = model->radiusTMP; + + for (int i = 0; i < new->params->n; i++) { + new->params->data.F32[i] = model->params->data.F32[i]; + new->dparams->data.F32[i] = model->dparams->data.F32[i]; + } + + return (new); +} + +/****************************************************************************** + pmModelEval(source, level, row): evaluates the model function at the specified coords. + + NOTE: The coords are in subImage source->pixel coords, not image coords. + + XXX: Use static vectors for x (NO: needs to be thread safe) +*****************************************************************************/ +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_PTR_NON_NULL(model->params, false); + + // Allocate the x coordinate structure and convert row/col to image space. + // + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = (psF32) (col + image->col0); + x->data.F32[1] = (psF32) (row + image->row0); + psF32 tmpF; + pmModelFunc modelFunc; + + modelFunc = pmModelFunc_GetFunction (model->type); + tmpF = modelFunc (NULL, model->params, x); + psFree(x); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpF); +} + +bool AddOrSubModel(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky, + bool add + ) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + PS_ASSERT_PTR_NON_NULL(model, false); + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + psVector *params = model->params; + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + psS32 imageCol; + psS32 imageRow; + psF32 skyValue = params->data.F32[0]; + psF32 pixelValue; + + for (psS32 i = 0; i < image->numRows; i++) { + for (psS32 j = 0; j < image->numCols; j++) { + if ((mask != NULL) && mask->data.U8[i][j]) + continue; + + // Convert i/j to image coord space: + // 'center' option changes meaning of i,j: + if (center) { + imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; + imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; + } else { + imageCol = j + image->col0; + imageRow = i + image->row0; + } + + x->data.F32[0] = (float) imageCol; + x->data.F32[1] = (float) imageRow; + + // set the appropriate pixel value for this coordinate + if (sky) { + pixelValue = modelFunc (NULL, params, x); + } else { + pixelValue = modelFunc (NULL, params, x) - skyValue; + } + + + // add or subtract the value + if (add + ) { + image->data.F32[i][j] += pixelValue; + } + else { + image->data.F32[i][j] -= pixelValue; + } + } + } + psFree(x); + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmModelAdd(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = AddOrSubModel(image, mask, model, center, sky, true); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +/****************************************************************************** + *****************************************************************************/ +bool pmModelSub(psImage *image, + psImage *mask, + pmModel *model, + bool center, + bool sky) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psBool rc = AddOrSubModel(image, mask, model, center, sky, false); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +// check for success of model fit +bool pmModelFitStatus (pmModel *model) +{ + + bool status; + + pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); + status = statusFunc (model); + + return (status); +} + Index: trunk/psModules/src/objects/pmModel.h =================================================================== --- trunk/psModules/src/objects/pmModel.h (revision 6872) +++ trunk/psModules/src/objects/pmModel.h (revision 6872) @@ -0,0 +1,113 @@ +/** @file pmObjects.h + * + * Functions to define and manipulate object models + * + * @author GLG, MHPCC + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_MODEL_H +# define PM_MODEL_H + +// type of model carried by the pmModel structure +typedef int pmModelType; + +typedef enum { + PM_MODEL_UNTRIED, ///< model fit not yet attempted + PM_MODEL_SUCCESS, ///< model fit succeeded + PM_MODEL_NONCONVERGE, ///< model fit did not converge + PM_MODEL_OFFIMAGE, ///< model fit drove out of range + PM_MODEL_BADARGS ///< model fit called with invalid args +} pmModelStatus; + +/** pmModel data structure + * + * Every source may have two types of models: a PSF model and a EXT (extended-source) + * model. The PSF model represents the best fit of the image PSF to the specific + * object. In this case, the PSF-dependent parameters are specified for the + * object by the PSF, not by the fit. The EXT model represents the best fit of + * the given model to the object, with all shape parameters floating in the fit. + * + */ +typedef struct +{ + pmModelType type; ///< Model to be used. + psVector *params; ///< Paramater values. + psVector *dparams; ///< Parameter errors. + float chisq; ///< Fit chi-squared. + float chisqNorm; ///< re-normalized fit chi-squared. + int nDOF; ///< number of degrees of freedom + int nIter; ///< number of iterations to reach min + pmModelStatus status; ///< fit status + float radiusTMP; ///< fit radius actually used +} +pmModel; + +/** pmModelAlloc() + * + */ +pmModel *pmModelAlloc(pmModelType type); + +// copy model to a new structure +pmModel *pmModelCopy (pmModel *model); + +psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row); + +/** pmModelAdd() + * + * Add the given source model flux to/from the provided image. The boolean + * option center selects if the source is re-centered to the image center or if + * it is placed at its centroid location. The boolean option sky selects if the + * background sky is applied (TRUE) or not. The pixel range in the target image + * is at most the pixel range specified by the source.pixels image. The success + * status is returned. + * + */ +bool pmModelAdd( + psImage *image, ///< The output image (float) + psImage *mask, ///< The image pixel mask (valid == 0) + pmModel *model, ///< The input pmModel + bool center, ///< A boolean flag that determines whether pixels are centered + bool sky ///< A boolean flag that determines if the sky is subtracted +); + + +/** pmModelSub() + * + * Subtract the given source model flux to/from the provided image. The + * boolean option center selects if the source is re-centered to the image center + * or if it is placed at its centroid location. The boolean option sky selects if + * the background sky is applied (TRUE) or not. The pixel range in the target + * image is at most the pixel range specified by the source.pixels image. The + * success status is returned. + * + */ +bool pmModelSub( + psImage *image, ///< The output image (float) + psImage *mask, ///< The image pixel mask (valid == 0) + pmModel *model, ///< The input pmModel + bool center, ///< A boolean flag that determines whether pixels are centered + bool sky ///< A boolean flag that determines if the sky is subtracted +); + +/** pmModelFitStatus() + * + * This function wraps the call to the model-specific function returned by + * pmModelFitStatusFunc_GetFunction. The model-specific function examines the + * model parameters, parameter errors, Chisq, S/N, and other parameters available + * from model to decide if the particular fit was successful or not. + * + * XXX: Must code this. + * + */ +bool pmModelFitStatus( + pmModel *model ///< Model to be used +); + +# endif /* PM_MODEL_H */ Index: trunk/psModules/src/objects/pmModelGroup.c =================================================================== --- trunk/psModules/src/objects/pmModelGroup.c (revision 6511) +++ trunk/psModules/src/objects/pmModelGroup.c (revision 6872) @@ -1,8 +1,35 @@ -# include "pmModelGroup.h" - +/** @file pmModelGroup.c + * + * Functions to define and manipulate object model attributes + * + * @author GLG, MHPCC + * @author EAM, IfA + * + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" + +// XXX shouldn't these be defined for us in pslib.h ??? double hypot(double x, double y); double sqrt (double x); -#include "psEllipse.h" #include "models/pmModel_GAUSS.c" #include "models/pmModel_PGAUSS.c" @@ -25,5 +52,4 @@ if (modelGroup == NULL) return; - psFree (modelGroup); return; } @@ -60,4 +86,11 @@ } Nmodels = Nnew; + return; +} + +void pmModelGroupCleanup (void) +{ + + psFree (models); return; } @@ -153,2 +186,28 @@ return (models[type].name); } + +/****************************************************************************** + pmSourceModelGuess(source, model): This function allocates a new + pmModel structure based on the given modelType specified in the argument list. + The corresponding pmModelGuess function is returned, and used to + supply the values of the params array in the pmModel structure. + + XXX: Many parameters are based on the src->moments structure, which is in + image, not subImage coords. Therefore, the calls to the model evaluation + functions will be in image, not subImage coords. Remember this. +*****************************************************************************/ +pmModel *pmSourceModelGuess(pmSource *source, + pmModelType modelType) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + + pmModel *model = pmModelAlloc(modelType); + + pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); + modelGuessFunc(model, source); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(model); +} + Index: trunk/psModules/src/objects/pmModelGroup.h =================================================================== --- trunk/psModules/src/objects/pmModelGroup.h (revision 6511) +++ trunk/psModules/src/objects/pmModelGroup.h (revision 6872) @@ -9,16 +9,84 @@ * @author EAM, IfA * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-24 01:24:32 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * */ -#include "pmObjects.h" -#include "pmPSF.h" -/** - * - * This function returns the number of parameters used by the listed function. - * + +# ifndef PM_MODEL_GROUP_H +# define PM_MODEL_GROUP_H + +// This function is the model chi-square minimization function for this model. +typedef psMinimizeLMChi2Func pmModelFunc; + +// This function returns the integrated flux for the given model parameters. +typedef psF64 (*pmModelFlux)(const psVector *params); + + +// This function returns the radius at which the given model and parameters +// achieves the given flux. +typedef psF64 (*pmModelRadius)(const psVector *params, double flux); + +/* This function sets the model parameter limits vectors for the given model + */ +typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); + +/* This function provides the model guess parameters based on the details of + * the given source. + */ +typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); + + +/* This function constructs the PSF model for the given source based on the + * supplied psf and the EXT model for the object. + */ +typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf); + +/* This function returns the success / failure status of the given model fit + */ +typedef bool (*pmModelFitStatusFunc)(pmModel *model); + +/* Every model instance belongs to a class of models, defined by the value of + * the pmModelType type entry. Various functions need access to information about + * each of the models. Some of this information varies from model to model, and + * may depend on the current parameter values or other data quantities. In order + * to keep the code from requiring the information about each model to be coded + * into the low-level fitting routines, we define a collection of functions which + * allow us to abstract this type of model-dependent information. These generic + * functions take the model type and return the corresponding function pointer + * for the specified model. Each model is defined by creating this collection of + * specific functions, and placing them in a single file for each model. We + * define the following structure to carry the collection of information about + * the models. + */ +typedef struct +{ + char *name; + int nParams; + pmModelFunc modelFunc; + pmModelFlux modelFlux; + pmModelRadius modelRadius; + pmModelLimits modelLimits; + pmModelGuessFunc modelGuessFunc; + pmModelFromPSFFunc modelFromPSFFunc; + pmModelFitStatusFunc modelFitStatusFunc; +} +pmModelGroup; + +// allocate a pmModelGroup to hold nModels entries +pmModelGroup *pmModelGroupAlloc (int nModels); + +// initialize the internal (static) model group with the default models +void pmModelGroupInit (void); + +// free the internal (static) model group +void pmModelGroupCleanup (void); + +// add a new model to the internal (static) model group +void pmModelGroupAdd (pmModelGroup *model); + +/* This function returns the number of parameters used by the listed function. */ int pmModelParameterCount( @@ -27,8 +95,5 @@ -/** - * - * This function returns the user-space model names for the specified model type. - * +/* This function returns the user-space model names for the specified model type. */ char *pmModelGetType( @@ -46,62 +111,4 @@ ); - -#ifndef PM_MODEL_GROUP_H -#define PM_MODEL_GROUP_H - -/** - * - * This function is the model chi-square minimization function for this model. - * - */ -typedef psMinimizeLMChi2Func pmModelFunc; - - -/** - * - * This function returns the integrated flux for the given model parameters. - */ -typedef psF64 (*pmModelFlux)(const psVector *params); - - -/** - * - * This function returns the radius at which the given model and parameters - * achieves the given flux. - * - */ -typedef psF64 (*pmModelRadius)(const psVector *params, double flux); - -/** - * - * This function sets the model parameter limits vectors for the given model - * - */ -typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max); - -/** - * - * This function provides the model guess parameters based on the details of - * the given source. - * - */ -typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source); - - -/** - * - * This function constructs the PSF model for the given source based on the - * supplied psf and the FLT model for the object. - * - */ -typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf); - -/** - * - * This function returns the success / failure status of the given model fit - * - */ -typedef bool (*pmModelFitStatusFunc)(pmModel *model); - /** * @@ -110,5 +117,4 @@ * */ - /** @@ -177,44 +183,17 @@ +/** pmSourceModelGuess() + * + * Convert available data to an initial guess for the given model. This + * function allocates a pmModel entry for the pmSource structure based on the + * provided model selection. The method of defining the model parameter guesses + * are specified for each model below. The guess values are placed in the model + * parameters. The function returns TRUE on success or FALSE on failure. + * + */ +pmModel *pmSourceModelGuess( + pmSource *source, ///< The input pmSource + pmModelType model ///< The type of model to be created. +); - -/** - * - * Every model instance belongs to a class of models, defined by the value of - * the pmModelType type entry. Various functions need access to information about - * each of the models. Some of this information varies from model to model, and - * may depend on the current parameter values or other data quantities. In order - * to keep the code from requiring the information about each model to be coded - * into the low-level fitting routines, we define a collection of functions which - * allow us to abstract this type of model-dependent information. These generic - * functions take the model type and return the corresponding function pointer - * for the specified model. Each model is defined by creating this collection of - * specific functions, and placing them in a single file for each model. We - * define the following structure to carry the collection of information about - * the models. - * - */ -typedef struct -{ - char *name; - int nParams; - pmModelFunc modelFunc; - pmModelFlux modelFlux; - pmModelRadius modelRadius; - pmModelLimits modelLimits; - pmModelGuessFunc modelGuessFunc; - pmModelFromPSFFunc modelFromPSFFunc; - pmModelFitStatusFunc modelFitStatusFunc; -} -pmModelGroup; - -// allocate a pmModelGroup to hold nModels entries -pmModelGroup *pmModelGroupAlloc (int nModels); - -// initialize the internal (static) model group with the default models -void pmModelGroupInit (void); - -// add a new model to the internal (static) model group -void pmModelGroupAdd (pmModelGroup *model); - -# endif +# endif /* PM_MODEL_GROUP_H */ Index: trunk/psModules/src/objects/pmMoments.c =================================================================== --- trunk/psModules/src/objects/pmMoments.c (revision 6872) +++ trunk/psModules/src/objects/pmMoments.c (revision 6872) @@ -0,0 +1,39 @@ +/** @file pmMoments.c + * + * Functions defining the pmMoments structure + * + * @author GLG, MHPCC + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include "pslib.h" +#include "pmMoments.h" + +/****************************************************************************** +pmMomentsAlloc(): Allocate the pmMoments structure and initialize the members +to zero. +*****************************************************************************/ +pmMoments *pmMomentsAlloc() +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmMoments *tmp = (pmMoments *) psAlloc(sizeof(pmMoments)); + tmp->x = 0.0; + tmp->y = 0.0; + tmp->Sx = 0.0; + tmp->Sy = 0.0; + tmp->Sxy = 0.0; + tmp->Sum = 0.0; + tmp->Peak = 0.0; + tmp->Sky = 0.0; + tmp->nPixels = 0; + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + Index: trunk/psModules/src/objects/pmMoments.h =================================================================== --- trunk/psModules/src/objects/pmMoments.h (revision 6872) +++ trunk/psModules/src/objects/pmMoments.h (revision 6872) @@ -0,0 +1,45 @@ +/** @file pmMoments.h + * + * Definitions of the moments structure + * + * @author GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_MOMENTS_H +# define PM_MOMENTS_H + +/** pmMoments data structure + * + * One of the simplest measurements which can be made quickly for an object + * are the object moments. We specify a structure to carry the moment information + * for a specific source: + * + */ +typedef struct +{ + float x; ///< X-coord of centroid. + float y; ///< Y-coord of centroid. + float Sx; ///< x-second moment. + float Sy; ///< y-second moment. + float Sxy; ///< xy cross moment. + float Sum; ///< Pixel sum above sky (background). + float Peak; ///< Peak counts above sky. + float Sky; ///< Sky level (background). + float dSky; ///< local Sky variance + float SN; ///< approx signal-to-noise + int nPixels; ///< Number of pixels used. +} +pmMoments; + +/** pmMomentsAlloc() + * + */ +pmMoments *pmMomentsAlloc(); + +# endif Index: trunk/psModules/src/objects/pmObjects.h =================================================================== --- trunk/psModules/src/objects/pmObjects.h (revision 6511) +++ trunk/psModules/src/objects/pmObjects.h (revision 6872) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-24 01:24:32 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -27,535 +27,4 @@ #include #include "pslib.h" -#include "pmAstrometry.h" -/** - * In the object analysis process, we will use specific mask values to mark the - * image pixels. The following structure defines the relevant mask values. - * - * XXX: This is probably a bad solution: we will want to set mask values - * outside of the PSPHOT code. Perhaps we can set up a registered set of mask - * values with specific meanings that other functions can add to or define? - */ -typedef enum { - PSPHOT_MASK_CLEAR = 0x00, - PSPHOT_MASK_INVALID = 0x01, - PSPHOT_MASK_SATURATED = 0x02, - PSPHOT_MASK_MARKED = 0x08, -} psphotMaskValues; - - -/** pmPeakType - * - * A peak pixel may have several features which may be determined when the - * peak is found or measured. These are specified by the pmPeakType enum. - * PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate - * neighbors. The PM_PEAK_EDGE represents a peak pixel which touching the image - * edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific - * number of neighbors at the same value, within some tolarence: - * - */ -typedef enum { - PM_PEAK_LONE, ///< Isolated peak. - PM_PEAK_EDGE, ///< Peak on edge. - PM_PEAK_FLAT, ///< Peak has equal-value neighbors. - PM_PEAK_UNDEF ///< Undefined. -} pmPeakType; - - -/** pmPeak data structure - * - * A source has the capacity for several types of measurements. The - * simplest measurement of a source is the location and flux of the peak pixel - * associated with the source: - * - */ -typedef struct -{ - int x; ///< X-coordinate of peak pixel. - int y; ///< Y-coordinate of peak pixel. - float counts; ///< Value of peak pixel (above sky?). - pmPeakType class; ///< Description of peak. -} -pmPeak; - - -/** pmMoments data structure - * - * One of the simplest measurements which can be made quickly for an object - * are the object moments. We specify a structure to carry the moment information - * for a specific source: - * - */ -typedef struct -{ - float x; ///< X-coord of centroid. - float y; ///< Y-coord of centroid. - float Sx; ///< x-second moment. - float Sy; ///< y-second moment. - float Sxy; ///< xy cross moment. - float Sum; ///< Pixel sum above sky (background). - float Peak; ///< Peak counts above sky. - float Sky; ///< Sky level (background). - float SN; ///< approx signal-to-noise - int nPixels; ///< Number of pixels used. -} -pmMoments; - - -/** pmPSFClump data structure - * - * A collection of object moment measurements can be used to determine - * approximate object classes. The key to this analysis is the location and - * statistics (in the second-moment plane, - * - */ -typedef struct -{ - float X; - float dX; - float Y; - float dY; -} -pmPSFClump; - -// type of model carried by the pmModel structure -typedef int pmModelType; - -typedef enum { - PM_MODEL_UNTRIED, ///< model fit not yet attempted - PM_MODEL_SUCCESS, ///< model fit succeeded - PM_MODEL_NONCONVERGE, ///< model fit did not converge - PM_MODEL_OFFIMAGE, ///< model fit drove out of range - PM_MODEL_BADARGS ///< model fit called with invalid args -} pmModelStatus; - -/** pmModel data structure - * - * Every source may have two types of models: a PSF model and a FLT (floating) - * model. The PSF model represents the best fit of the image PSF to the specific - * object. In this case, the PSF-dependent parameters are specified for the - * object by the PSF, not by the fit. The FLT model represents the best fit of - * the given model to the object, with all parameters floating in the fit. - * - */ -typedef struct -{ - pmModelType type; ///< Model to be used. - psVector *params; ///< Paramater values. - psVector *dparams; ///< Parameter errors. - float chisq; ///< Fit chi-squared. - int nDOF; ///< number of degrees of freedom - int nIter; ///< number of iterations to reach min - int status; ///< fit status - float radius; ///< fit radius actually used -} -pmModel; - -/** pmSourceType enumeration - * - * A given source may be identified as most-likely to be one of several source - * types. The pmSource entry pmSourceType defines the current best-guess for this - * source. - * - * XXX: The values given below are currently illustrative and will require - * some modification as the source classification code is developed. (TBD) - * - */ -typedef enum { - PM_SOURCE_DEFECT, ///< a cosmic-ray - PM_SOURCE_SATURATED, ///< random saturated pixels - - PM_SOURCE_SATSTAR, ///< a saturated star - PM_SOURCE_PSFSTAR, ///< a PSF star - PM_SOURCE_GOODSTAR, ///< a good-quality star - - PM_SOURCE_POOR_FIT_PSF, ///< poor quality PSF fit - PM_SOURCE_FAIL_FIT_PSF, ///< failed to get a good PSF fit - PM_SOURCE_FAINTSTAR, ///< below S/N cutoff - - PM_SOURCE_GALAXY, ///< an extended object (galaxy) - PM_SOURCE_FAINT_GALAXY, ///< a galaxy below S/N cutoff - PM_SOURCE_DROP_GALAXY, ///< ? - PM_SOURCE_FAIL_FIT_GAL, ///< failed on the galaxy fit - PM_SOURCE_POOR_FIT_GAL, ///< poor quality galaxy fit - - PM_SOURCE_OTHER, ///< unidentified -} pmSourceType; - -/** pmSource data structure - * - * This source has the capacity for several types of measurements. The - * simplest measurement of a source is the location and flux of the peak pixel - * associated with the source: - * - */ -typedef struct -{ - pmPeak *peak; ///< Description of peak pixel. - psImage *pixels; ///< Rectangular region including object pixels. - psImage *weight; ///< Image variance. - psImage *mask; ///< Mask which marks pixels associated with objects. - pmMoments *moments; ///< Basic moments measure for the object. - pmModel *modelPSF; ///< PSF Model fit (parameters and type) - pmModel *modelFLT; ///< FLT (floating) Model fit (parameters and type). - pmSourceType type; ///< Best identification of object. - float apMag; - float fitMag; -} -pmSource; - - -/** pmPeakAlloc() - * - * @return pmPeak* newly allocated pmPeak with all internal pointers set to NULL - */ -pmPeak *pmPeakAlloc( - int x, ///< Row-coordinate in image space - int y, ///< Col-coordinate in image space - float counts, ///< The value of the peak pixel - pmPeakType class ///< The type of peak pixel -); - - -/** pmMomentsAlloc() - * - */ -pmMoments *pmMomentsAlloc(); - - -/** pmModelAlloc() - * - */ -pmModel *pmModelAlloc(pmModelType type); - - -/** pmSourceAlloc() - * - */ -pmSource *pmSourceAlloc(); - - -/** pmFindVectorPeaks() - * - * Find all local peaks in the given vector above the given threshold. A peak - * is defined as any element with a value greater than its two neighbors and with - * a value above the threshold. Two types of special cases must be addressed. - * Equal value elements: If an element has the same value as the following - * element, it is not considered a peak. If an element has the same value as the - * preceding element (but not the following), then it is considered a peak. Note - * that this rule (arbitrarily) identifies flat regions by their trailing edge. - * Edge cases: At start of the vector, the element must be higher than its - * neighbor. At the end of the vector, the element must be higher or equal to its - * neighbor. These two rules again places the peak associated with a flat region - * which touches the image edge at the image edge. The result of this function is - * a vector containing the coordinates (element number) of the detected peaks - * (type psU32). - * - */ -psVector *pmFindVectorPeaks( - const psVector *vector, ///< The input vector (float) - float threshold ///< Threshold above which to find a peak -); - - -/** pmFindImagePeaks() - * - * Find all local peaks in the given image above the given threshold. This - * function should find all row peaks using pmFindVectorPeaks, then test each row - * peak and exclude peaks which are not local peaks. A peak is a local peak if it - * has a higher value than all 8 neighbors. If the peak has the same value as its - * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors - * have an equal value, the peak is considered a valid peak. Note two points: - * first, the +x neighbor condition is already enforced by pmFindVectorPeaks. - * Second, these rules have the effect of making flat-topped regions have single - * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be - * set. The result of this function is an array of pmPeak entries. - * - */ -psArray *pmFindImagePeaks( - const psImage *image, ///< The input image where peaks will be found (float) - float threshold ///< Threshold above which to find a peak -); - - -/** pmCullPeaks() - * - * Eliminate peaks from the psList that have a peak value above the given - * maximum, or fall outside the valid region. - * - */ -psList *pmCullPeaks( - psList *peaks, ///< The psList of peaks to be culled - float maxValue, ///< Cull peaks above this value - const psRegion valid ///< Cull peaks otside this psRegion -); - - -/** pmPeaksSubset() - * - * Create a new peaks array, removing certain types of peaks from the input - * array of peaks based on the given criteria. Peaks should be eliminated if they - * have a peak value above the given maximum value limit or if the fall outside - * the valid region. The result of the function is a new array with a reduced - * number of peaks. - * - */ -psArray *pmPeaksSubset( - psArray *peaks, ///< Add comment. - float maxvalue, ///< Add comment. - const psRegion valid ///< Add comment. -); - - -/** pmSourceDefinePixels() - * - * Define psImage subarrays for the source located at coordinates x,y on the - * image set defined by readout. The pixels defined by this operation consist of - * a square window (of full width 2Radius+1) centered on the pixel which contains - * the given coordinate, in the frame of the readout. The window is defined to - * have limits which are valid within the boundary of the readout image, thus if - * the radius would fall outside the image pixels, the subimage is truncated to - * only consist of valid pixels. If readout->mask or readout->weight are not - * NULL, matching subimages are defined for those images as well. This function - * fails if no valid pixels can be defined (x or y less than Radius, for - * example). This function should be used to define a region of interest around a - * source, including both source and sky pixels. - * - * XXX: must code this. - * - */ -// XXX: Uncommenting the pmReadout causes compile errors. -bool pmSourceDefinePixels( - pmSource *mySource, ///< Add comment. - pmReadout *readout, ///< Add comment. - psF32 x, ///< Add comment. - psF32 y, ///< Add comment. - psF32 Radius ///< Add comment. -); - - -/** pmSourceLocalSky() - * - * Measure the local sky in the vicinity of the given source. The Radius - * defines the square aperture in which the moments will be measured. This - * function assumes the source pixels have been defined, and that the value of - * Radius here is smaller than the value of Radius used to define the pixels. The - * annular region not contained within the radius defined here is used to measure - * the local background in the vicinity of the source. The local background - * measurement uses the specified statistic passed in via the statsOptions entry. - * This function allocates the pmMoments structure. The resulting sky is used to - * set the value of the pmMoments.sky element of the provided pmSource structure. - * - */ -bool pmSourceLocalSky( - pmSource *source, ///< The input image (float) - psStatsOptions statsOptions, ///< The statistic used in calculating the background sky - float Radius ///< The inner radius of the square annulus to exclude -); - - -/** pmSourceMoments() - * - * Measure source moments for the given source, using the value of - * source.moments.sky provided as the local background value and the peak - * coordinates as the initial source location. The resulting moment values are - * applied to the source.moments entry, and the source is returned. The moments - * are measured within the given circular radius of the source.peak coordinates. - * The return value indicates the success (TRUE) of the operation. - * - */ -bool pmSourceMoments( - pmSource *source, ///< The input pmSource for which moments will be computed - float radius ///< Use a circle of pixels around the peak -); - - -/** pmSourcePSFClump() - * - * We use the source moments to make an initial, approximate source - * classification, and as part of the information needed to build a PSF model for - * the image. As long as the PSF shape does not vary excessively across the - * image, the sources which are represented by a PSF (the start) will have very - * similar second moments. The function pmSourcePSFClump searches a collection of - * sources with measured moments for a group with moments which are all very - * similar. The function returns a pmPSFClump structure, representing the - * centroid and size of the clump in the sigma_x, sigma_y second-moment plane. - * - * The goal is to identify and characterize the stellar clump within the - * sigma_x, sigma_y second-moment plane. To do this, an image is constructed to - * represent this plane. The units of sigma_x and sigma_y are in image pixels. A - * pixel in this analysis image represents 0.1 pixels in the input image. The - * dimensions of the image need only be 10 pixels. The peak pixel in this image - * (above a threshold of half of the image maximum) is found. The coordinates of - * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution. - * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then - * * used to calculate the median and standard deviation of the sigma_x, sigma_y - * values. These resulting values are returned via the pmPSFClump structure. - * - * The return value indicates the success (TRUE) of the operation. - * - * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD). - * XXX: Save the clump parameters on the Metadata (TBD) - * - */ -pmPSFClump pmSourcePSFClump( - psArray *source, ///< The input pmSource - psMetadata *metadata ///< Contains classification parameters -); - - -/** pmSourceRoughClass() - * - * Based on the specified data values, make a guess at the source - * classification. The sources are provides as a psArray of pmSource entries. - * Definable parameters needed to make the classification are provided to the - * routine with the psMetadata structure. The rules (in SDRS) refer to values which - * can be extracted from the metadata using the given keywords. Except as noted, - * the data type for these parameters are psF32. - * - */ -bool pmSourceRoughClass( - psArray *source, ///< The input pmSource - psMetadata *metadata, ///< Contains classification parameters - pmPSFClump clump ///< Statistics about the PSF clump -); - - -/** pmSourceModelGuess() - * - * Convert available data to an initial guess for the given model. This - * function allocates a pmModel entry for the pmSource structure based on the - * provided model selection. The method of defining the model parameter guesses - * are specified for each model below. The guess values are placed in the model - * parameters. The function returns TRUE on success or FALSE on failure. - * - */ -pmModel *pmSourceModelGuess( - pmSource *source, ///< The input pmSource - pmModelType model ///< The type of model to be created. -); - - -/** pmContourType - * - * Only one type is defined at present. - * - */ -typedef enum { - PS_CONTOUR_CRUDE, - PS_CONTOUR_UNKNOWN01, - PS_CONTOUR_UNKNOWN02 -} pmContourType; - - -/** pmSourceContour() - * - * Find points in a contour for the given source at the given level. If type - * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak, - * running along each pixel row until the level is crossed, then interpolating to - * the level coordinate for that row. This is done for each row, with the - * starting point determined by the midpoint of the previous row, until the - * starting point has a value below the contour level. The returned contour - * consists of two vectors giving the x and y coordinates of the contour levels. - * This function may be used as part of the model guess inputs. Other contour - * types may be specified in the future for more refined contours (TBD) - * - */ -psArray *pmSourceContour( - pmSource *source, ///< The input pmSource - const psImage *image, ///< The input image (float) (this arg should be removed) - float level, ///< The level of the contour - pmContourType mode ///< Currently this must be PS_CONTOUR_CRUDE -); - - -/** pmSourceFitModel() - * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. - * - */ -bool pmSourceFitModel( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or FLT? -); - - -/** pmModelFitStatus() - * - * This function wraps the call to the model-specific function returned by - * pmModelFitStatusFunc_GetFunction. The model-specific function examines the - * model parameters, parameter errors, Chisq, S/N, and other parameters available - * from model to decide if the particular fit was successful or not. - * - * XXX: Must code this. - * - */ -bool pmModelFitStatus( - pmModel *model ///< Add comment. -); - - -/** pmSourceAddModel() - * - * Add the given source model flux to/from the provided image. The boolean - * option center selects if the source is re-centered to the image center or if - * it is placed at its centroid location. The boolean option sky selects if the - * background sky is applied (TRUE) or not. The pixel range in the target image - * is at most the pixel range specified by the source.pixels image. The success - * status is returned. - * - */ -bool pmSourceAddModel( - psImage *image, ///< The output image (float) - psImage *mask, ///< The image pixel mask (valid == 0) - pmModel *model, ///< The input pmModel - bool center, ///< A boolean flag that determines whether pixels are centered - bool sky ///< A boolean flag that determines if the sky is subtracted -); - - -/** pmSourceSubModel() - * - * Subtract the given source model flux to/from the provided image. The - * boolean option center selects if the source is re-centered to the image center - * or if it is placed at its centroid location. The boolean option sky selects if - * the background sky is applied (TRUE) or not. The pixel range in the target - * image is at most the pixel range specified by the source.pixels image. The - * success status is returned. - * - */ -bool pmSourceSubModel( - psImage *image, ///< The output image (float) - psImage *mask, ///< The image pixel mask (valid == 0) - pmModel *model, ///< The input pmModel - bool center, ///< A boolean flag that determines whether pixels are centered - bool sky ///< A boolean flag that determines if the sky is subtracted -); - - -/** - * - * The function returns both the magnitude of the fit, defined as -2.5log(flux), - * where the flux is integrated under the model, theoretically from a radius of 0 - * to infinity. In practice, we integrate the model beyond 50sigma. The aperture magnitude is - * defined as -2.5log(flux) , where the flux is summed for all pixels which are - * not excluded by the aperture mask. The model flux is calculated by calling the - * model-specific function provided by pmModelFlux_GetFunction. - * - * XXX: must code this. - * - */ -bool pmSourcePhotometry( - float *fitMag, ///< integrated fit magnitude - float *obsMag, ///< aperture flux magnitude - pmModel *model, ///< model used for photometry - psImage *image, ///< image pixels to be used - psImage *mask ///< mask of pixels to ignore -); - /** @@ -563,4 +32,5 @@ * This function converts the source classification into the closest available * approximation to the Dophot classification scheme: + * XXX EAM : fix this to use current source classification scheme * * PM_SOURCE_DEFECT: 8 @@ -598,60 +68,4 @@ ); -/** pmSourceFitModel_v5() - * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. - * - */ -bool pmSourceFitModel_v5( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or FLT? -); - - -/** pmSourceFitModel_v7() - * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. - * - */ -bool pmSourceFitModel_v7( - pmSource *source, ///< The input pmSource - pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or FLT? -); - - -/** pmSourcePhotometry() - * - * XXX: Need descriptions - * - */ -bool pmSourcePhotometry( - float *fitMag, - float *obsMag, - pmModel *model, - psImage *image, - psImage *mask -); - -/** pmModelEval() - * - * XXX: Need descriptions - * - */ -psF32 pmModelEval( - pmModel *model, - psImage *image, - psS32 col, - psS32 row -); #endif Index: trunk/psModules/src/objects/pmPSF.h =================================================================== --- trunk/psModules/src/objects/pmPSF.h (revision 6511) +++ trunk/psModules/src/objects/pmPSF.h (revision 6872) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-10 19:53:40 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -15,5 +15,4 @@ # ifndef PM_PSF_H # define PM_PSF_H - /** pmPSF data structure @@ -33,12 +32,29 @@ pmModelType type; ///< PSF Model in use psArray *params; ///< Model parameters (psPolynomial2D) - float chisq; ///< PSF goodness statistic - float ApResid; ///< ??? - float dApResid; ///< ??? - float skyBias; ///< ??? + psPolynomial1D *ChiTrend; ///< Chisq vs flux fit (correction for systematic errors) + psPolynomial4D *ApTrend; ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst) + pmGrowthCurve *growth; ///< apMag vs Radius + float ApResid; ///< apMag - psfMag (for PSF stars) + float dApResid; ///< scatter of ApResid + float skyBias; ///< implied residual sky offset from ApResid fit + float skySat; ///< roll-over of ApResid fit + float chisq; ///< PSF goodness statistic (unused??) int nPSFstars; ///< number of stars used to measure PSF + int nApResid; ///< number of stars used to measure ApResid + bool poissonErrors; } pmPSF; +typedef enum { + PM_PSF_NONE, + PM_PSF_CONSTANT, + PM_PSF_SKYBIAS, + PM_PSF_SKYSAT, + PM_PSF_XY_LIN, + PM_PSF_XY_QUAD, + PM_PSF_SKY_XY_LIN, + PM_PSF_SKYSAT_XY_LIN, + PM_PSF_ALL +} pmPSF_ApTrendOptions; /** @@ -48,5 +64,6 @@ */ pmPSF *pmPSFAlloc( - pmModelType type ///< Add comment + pmModelType type, // type of model for PSF + bool poissonErrors ///< use poissonian errors or not? ); @@ -85,3 +102,5 @@ ); +bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option); + # endif Index: trunk/psModules/src/objects/pmPSF_IO.c =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.c (revision 6872) +++ trunk/psModules/src/objects/pmPSF_IO.c (revision 6872) @@ -0,0 +1,312 @@ +/** @file pmPSF_IO.c + * + * This file contains functions to read and write PSF models using the psMetadata Config file + * format. + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ + +#include +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAfile.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmPSF_IO.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf) +{ + + if (metadata == NULL) { + metadata = psMetadataAlloc (); + } + + char *modelName = pmModelGetType (psf->type); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); + + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors); + + int nPar = pmModelParameterCount (psf->type) ; + psMetadataAdd (metadata, 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 (metadata, poly, "PSF_PAR%02d", i); + } + psPolynomial4DtoMD (metadata, psf->ApTrend, "APTREND"); + + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); + + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars); + + return metadata; +} + +pmPSF *pmPSFfromMD (psMetadata *metadata) +{ + + bool status; + char keyword[80]; + + char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME"); + pmModelType type = pmModelSetType (modelName); + + bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS"); + if (!status) + poissonErrors = true; + + pmPSF *psf = pmPSFAlloc (type, poissonErrors); + + int nPar = psMetadataLookupS32 (&status, metadata, "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, metadata, keyword); + psPolynomial2D *poly = psPolynomial2DfromMD (folder); + psFree (psf->params->data[i]); + psf->params->data[i] = poly; + } + sprintf (keyword, "APTREND"); + psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); + psPolynomial4D *poly = psPolynomial4DfromMD (folder); + psFree (psf->ApTrend); + psf->ApTrend = poly; + + psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID"); + psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID"); + psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS"); + + psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ"); + psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS"); + + psFree (metadata); + return (psf); +} + +bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAWritePSFmodel (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipWritePSFmodel (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellWritePSFmodel (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutWritePSFmodel (readout, view, file); + return true; +} + +// read in all chip-level PSFmodel files for this FPA +bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + pmChipWritePSFmodel (chip, view, file); + } + return true; +} + +// read in all cell-level PSFmodel files for this chip +bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + pmCellWritePSFmodel (cell, view, file); + } + return true; +} + +// read in all readout-level PSFmodel files for this cell +bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + pmReadoutWritePSFmodel (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file) +{ + bool status; + char *filename; + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + + switch (file->type) { + case PM_FPA_FILE_PSF: + filename = pmFPAfileNameFromRule (file->filerule, file, view); + psMetadata *psfData = pmPSFtoMD (NULL, psf); + psMetadataConfigWrite (psfData, filename); + psFree (psfData); + psFree (filename); + return true; + + default: + fprintf (stderr, "warning: type mismatch\n"); + break; + } + return false; +} + + + +bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAReadPSFmodel (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipReadPSFmodel (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellReadPSFmodel (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutReadPSFmodel (readout, view, file); + return true; +} + +// read in all chip-level PSFmodel files for this FPA +bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + pmChipReadPSFmodel (chip, view, file); + } + return true; +} + +// read in all cell-level PSFmodel files for this chip +bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + pmCellReadPSFmodel (cell, view, file); + } + return true; +} + +// read in all readout-level PSFmodel files for this cell +bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + pmReadoutReadPSFmodel (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file) +{ + + int Nfail; + char *filename; + + switch (file->type) { + case PM_FPA_FILE_PSF: + filename = pmFPAfileNameFromRule (file->filerule, file, view); + + psMetadata *psfData = psMetadataConfigParse (NULL, &Nfail, filename, FALSE); + pmPSF *psf = pmPSFfromMD (psfData); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf); + + psFree (psf); + psFree (psfData); + psFree (filename); + + return true; + + default: + fprintf (stderr, "warning: type mismatch\n"); + break; + } + return false; +} + Index: trunk/psModules/src/objects/pmPSF_IO.h =================================================================== --- trunk/psModules/src/objects/pmPSF_IO.h (revision 6872) +++ trunk/psModules/src/objects/pmPSF_IO.h (revision 6872) @@ -0,0 +1,32 @@ +/** @file pmPSF.h + * + * This file contains typedefs for the Point-Spread Function and prototypes + * for functions that calculate the PSF. + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_PSF_IO_H +# define PM_PSF_IO_H + +psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf); +pmPSF *pmPSFfromMD (psMetadata *metadata); +bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file); +bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file); +bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file); +bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file); +bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file); + +bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file); +bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file); +bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file); +bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file); +bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file); + +# endif Index: trunk/psModules/src/objects/pmPSFtry.h =================================================================== --- trunk/psModules/src/objects/pmPSFtry.h (revision 6511) +++ trunk/psModules/src/objects/pmPSFtry.h (revision 6872) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-24 01:24:32 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -18,11 +18,11 @@ /** - * + * * This structure contains a pointer to the collection of sources which will * be used to test the PSF model form. It lists the pmModelType type of model * being tests, and contains an element to store the resulting psf * representation. In addition, this structure carries the complete collection of - * FLT (floating parameter) and PSF (fixed parameter) model fits to each of the - * sources modelFLT and modelPSF. It also contains a mask which is set by the + * EXT (floating parameter) and PSF (fixed parameter) model fits to each of the + * sources modelEXT and modelPSF. It also contains a mask which is set by the * model fitting and psf fitting steps. For each model, the value of the quality * metric is stored in the vector metric and the fitted instrumental magnitude is @@ -38,9 +38,9 @@ * ultimate metric to intercompare multiple types of PSF models is the value of * the aperture correction scatter. - * + * * XXX: There are many more members in the SDRS then in the prototype code. * I stuck with the prototype code. - * - * + * + * */ typedef struct @@ -48,5 +48,5 @@ pmPSF *psf; ///< Add comment. psArray *sources; ///< pointers to the original sources - psArray *modelFLT; ///< model fits, floating parameters + psArray *modelEXT; ///< model fits, floating parameters psArray *modelPSF; ///< model fits, PSF parameters psVector *mask; ///< Add comment. @@ -58,13 +58,13 @@ /** pmPSFtryMaskValues - * + * * The following datatype defines the masks used by the pmPSFtry analysis to * identify sources which should or should not be included in the analysis. - * + * */ -enum { +typedef enum { PSFTRY_MASK_CLEAR = 0x00, ///< Add comment. PSFTRY_MASK_OUTLIER = 0x01, ///< 1: outlier in psf polynomial fit (provided by psPolynomials) - PSFTRY_MASK_FLT_FAIL = 0x02, ///< 2: flt model failed to converge + PSFTRY_MASK_EXT_FAIL = 0x02, ///< 2: ext model failed to converge PSFTRY_MASK_PSF_FAIL = 0x04, ///< 3: psf model failed to converge PSFTRY_MASK_BAD_PHOT = 0x08, ///< 4: invalid source photometry @@ -74,33 +74,35 @@ /** pmPSFtryAlloc() - * + * * Allocate a pmPSFtry data structure. - * + * */ pmPSFtry *pmPSFtryAlloc( psArray *stars, ///< Add comment. - char *modelName ///< Add comment. + char *modelName, ///< Add comment. + bool poissonErrors // use poissonian or constant errors? ); /** pmPSFtryModel() - * + * * This function takes the input collection of sources and performs a complete * analysis to determine a PSF model of the given type (specified by model name). * The result is a pmPSFtry with the results of the analysis. - * + * */ pmPSFtry *pmPSFtryModel( psArray *sources, ///< Add comment. char *modelName, ///< Add comment. - float radius ///< Add comment. + float radius, ///< Add comment. + bool poissonErrors // use poissonian or constant errors? ); /** pmPSFtryMetric() - * + * * This function is used to measure the PSF model metric for the set of * results contained in the pmPSFtry structure. - * + * */ bool pmPSFtryMetric( Index: trunk/psModules/src/objects/pmPeaks.c =================================================================== --- trunk/psModules/src/objects/pmPeaks.c (revision 6872) +++ trunk/psModules/src/objects/pmPeaks.c (revision 6872) @@ -0,0 +1,542 @@ +/** @file pmPeaks.c + * + * This file defines functions to detect and manipulate peaks in images + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmPeaks.h" + +/****************************************************************************** +myListAddPeak(): A private function which allocates a psArray, if the list +argument is NULL, otherwise it adds the peak to that list. +XXX EAM : changed the output to psArray +XXX EAM : Switched row, col args +XXX EAM : NOTE: this was changed in the call, so the new code is consistent +*****************************************************************************/ +static psArray *myListAddPeak(psArray *list, + psS32 row, + psS32 col, + psF32 counts, + pmPeakType type) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); + + if (list == NULL) { + list = psArrayAlloc(100); + list->n = 0; + } + psArrayAdd(list, 100, tmpPeak); + psFree (tmpPeak); + // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) + + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(list); +} + +/****************************************************************************** +getRowVectorFromImage(): a private function which simply returns a +psVector containing the specified row of data from the psImage. + +XXX: Is there a better way to do this? +XXX EAM: does this really need to alloc a new vector??? +*****************************************************************************/ +static psVector *getRowVectorFromImage(psImage *image, + psU32 row) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + + psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); + for (psU32 col = 0; col < image->numCols ; col++) { + tmpVector->data.F32[col] = image->data.F32[row][col]; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(tmpVector); +} + +/****************************************************************************** +isItInThisRegion(): a private function which simply returns a +boolean denoting if specified coordinate is in the region. +XXX: Macro this. +*****************************************************************************/ +static bool isItInThisRegion(const psRegion valid, + psS32 x, + psS32 y) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + if ((x >= valid.x0) && + (x <= valid.x1) && + (y >= valid.y0) && + (y <= valid.y1)) { + psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); + return(true); + } + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return(false); +} + +/****************************************************************************** +pmPeakAlloc(): Allocate the pmPeak data structure and set appropriate members. +*****************************************************************************/ +pmPeak *pmPeakAlloc(psS32 x, + psS32 y, + psF32 counts, + pmPeakType type) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *tmp = (pmPeak *) psAlloc(sizeof(pmPeak)); + tmp->x = x; + tmp->y = y; + tmp->counts = counts; + tmp->type = type; + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +// psSort comparison function for peaks +int pmPeaksCompareAscend (const void **a, const void **b) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); + return (-1); + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); + return (+1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return (0); +} + +// psSort comparison function for peaks +int pmPeaksCompareDescend (const void **a, const void **b) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmPeak *A = *(pmPeak **)a; + pmPeak *B = *(pmPeak **)b; + + psF32 diff; + + diff = A->counts - B->counts; + if (diff < FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); + return (+1); + } else if (diff > FLT_EPSILON) { + psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); + return (-1); + } + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return (0); +} + +/****************************************************************************** +pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector +above the given threshold. Returns a vector of type PS_TYPE_U32 containing +the location (x value) of all peaks. + +XXX: What types should be supported? Only F32 is implemented. + +XXX: We currently step through the input vector twice; once to determine the +size of the output vector, then to set the values of the output vector. +Depending upon actual use, this may need to be optimized. +*****************************************************************************/ +psVector *pmFindVectorPeaks(const psVector *vector, + psF32 threshold) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_VECTOR_NON_NULL(vector, NULL); + PS_ASSERT_VECTOR_NON_EMPTY(vector, NULL); + PS_ASSERT_VECTOR_TYPE(vector, PS_TYPE_F32, NULL); + int count = 0; + int n = vector->n; + + // + // Special case: the input vector has a single element. + // + if (n == 1) { + psVector *tmpVector = NULL; + ; + if (vector->data.F32[0] > threshold) { + tmpVector = psVectorAlloc(1, PS_TYPE_U32); + tmpVector->data.U32[0] = 0; + } else { + tmpVector = psVectorAlloc(0, PS_TYPE_U32); + } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpVector); + } + + // + // Determine if first pixel is a peak + // + if ((vector->data.F32[0] > vector->data.F32[1]) && + (vector->data.F32[0] > threshold)) { + count++; + } + + // + // Determine if interior pixels are peaks + // + for (psU32 i = 1; i < n-1 ; i++) { + if ((vector->data.F32[i] > vector->data.F32[i-1]) && + (vector->data.F32[i] > vector->data.F32[i+1]) && + (vector->data.F32[i] > threshold)) { + count++; + } + } + + // + // Determine if last pixel is a peak + // + if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && + (vector->data.F32[n-1] > threshold)) { + count++; + } + + // + // We know how many peaks exist, so we now allocate a psVector to store + // those peaks. + // + psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32); + count = 0; + + // + // Determine if first pixel is a peak + // + if ((vector->data.F32[0] > vector->data.F32[1]) && + (vector->data.F32[0] > threshold)) { + tmpVector->data.U32[count++] = 0; + } + + // + // Determine if interior pixels are peaks + // + for (psU32 i = 1; i < (n-1) ; i++) { + if ((vector->data.F32[i] > vector->data.F32[i-1]) && + (vector->data.F32[i] > vector->data.F32[i+1]) && + (vector->data.F32[i] > threshold)) { + tmpVector->data.U32[count++] = i; + } + } + + // + // Determine if last pixel is a peak + // + if ((vector->data.F32[n-1] > vector->data.F32[n-2]) && + (vector->data.F32[n-1] > threshold)) { + tmpVector->data.U32[count++] = n-1; + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpVector); +} + + +/****************************************************************************** +pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage +above the given threshold. Returns a psArray containing location (x/y value) +of all peaks. + +XXX: I'm not convinced the peak type definition in the SDRS is mutually +exclusive. Some peaks can have multiple types. Edges for sure. Also, a +digonal line with the same value at each point will have a peak for every +point on that line. + +XXX: This does not work if image has either a single row, or a single column. + +XXX: In the output psArray elements, should we use the image row/column offsets? + Currently, we do not. +XXX EAM : this function needs to return peaks in *parent* coords + +XXX: Merge with CVS 1.20. This had the proper code for images with a single +row or column. + +*****************************************************************************/ +psArray *pmFindImagePeaks(const psImage *image, + psF32 threshold) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + if ((image->numRows == 1) || (image->numCols == 1)) { + psError(PS_ERR_UNKNOWN, true, "Currently, input image must have at least 2 rows and 2 columns."); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + } + psVector *tmpRow = NULL; + psU32 col = 0; + psU32 row = 0; + psArray *list = NULL; + + psU32 col0 = image->col0; + psU32 row0 = image->row0; + + // + // Find peaks in row 0 only. + // + row = 0; + tmpRow = getRowVectorFromImage((psImage *) image, row); + psVector *row1 = pmFindVectorPeaks(tmpRow, threshold); + // pmFindVectorPeaks returns coords in the vector, not corrected for col0 + + for (psU32 i = 0 ; i < row1->n ; i++ ) { + col = row1->data.U32[i]; + // + // Determine if pixel (0,0) is a peak. + // + if (col == 0) { + if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + } else if (col < (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + + } else if (col == (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + + } else { + psError(PS_ERR_UNKNOWN, true, "peak specified valid column range."); + } + } + psFree (tmpRow); + psFree (row1); + + // + // Exit if this image has a single row. + // + if (image->numRows == 1) { + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(list); + } + + // + // Find peaks in interior rows only. + // + for (row = 1 ; row < (image->numRows - 1) ; row++) { + tmpRow = getRowVectorFromImage((psImage *) image, row); + row1 = pmFindVectorPeaks(tmpRow, threshold); + + // Step through all local peaks in this row. + for (psU32 i = 0 ; i < row1->n ; i++ ) { + pmPeakType myType = PM_PEAK_UNDEF; + col = row1->data.U32[i]; + + if (col == 0) { + // If col==0, then we can not read col-1 pixels + if ((image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + myType = PM_PEAK_EDGE; + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); + } + } else if (col < (image->numCols - 1)) { + // This is an interior pixel + if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) { + if (image->data.F32[row][col] > threshold) { + if ((image->data.F32[row][col] > image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row+1][col]) && + (image->data.F32[row][col] > image->data.F32[row+1][col+1])) { + myType = PM_PEAK_LONE; + } + + if ((image->data.F32[row][col] == image->data.F32[row-1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row-1][col]) || + (image->data.F32[row][col] == image->data.F32[row-1][col+1]) || + (image->data.F32[row][col] == image->data.F32[row][col-1]) || + (image->data.F32[row][col] == image->data.F32[row][col+1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col-1]) || + (image->data.F32[row][col] == image->data.F32[row+1][col]) || + (image->data.F32[row][col] == image->data.F32[row+1][col+1])) { + myType = PM_PEAK_FLAT; + } + + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); + } + } + } else if (col == (image->numCols - 1)) { + // If col==numCols - 1, then we can not read col+1 pixels + if ((image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row+1][col])) { + myType = PM_PEAK_EDGE; + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], myType); + } + } else { + psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); + } + + } + psFree (tmpRow); + psFree (row1); + } + + // + // Find peaks in the last row only. + // + row = image->numRows - 1; + tmpRow = getRowVectorFromImage((psImage *) image, row); + row1 = pmFindVectorPeaks(tmpRow, threshold); + for (psU32 i = 0 ; i < row1->n ; i++ ) { + col = row1->data.U32[i]; + if (col == 0) { + if ( (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col+1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + } else if (col < (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] >= image->data.F32[row-1][col+1]) && + (image->data.F32[row][col] > image->data.F32[row][col-1]) && + (image->data.F32[row][col] >= image->data.F32[row][col+1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + + } else if (col == (image->numCols - 1)) { + if ( (image->data.F32[row][col] >= image->data.F32[row-1][col-1]) && + (image->data.F32[row][col] > image->data.F32[row-1][col]) && + (image->data.F32[row][col] > image->data.F32[row][col-1])) { + if (image->data.F32[row][col] > threshold) { + list = myListAddPeak(list, row + row0, col + col0, image->data.F32[row][col], PM_PEAK_EDGE); + } + } + } else { + psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range."); + } + } + psFree (tmpRow); + psFree (row1); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(list); +} + + +/****************************************************************************** +psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have +a peak value above the given maximum, or fall outside the valid region. + +XXX: Should the sky value be used when comparing the maximum? + +XXX: warning message if valid is NULL? + +XXX: changed API to create a NEW output psArray (should change name as well) + +XXX: Do we free the psList elements of those culled peaks? + +XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset? +*****************************************************************************/ +psList *pmCullPeaks(psList *peaks, + psF32 maxValue, + const psRegion valid) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(peaks, NULL); + + psListElem *tmpListElem = (psListElem *) peaks->head; + psS32 indexNum = 0; + + // printf("pmCullPeaks(): list size is %d\n", peaks->size); + while (tmpListElem != NULL) { + pmPeak *tmpPeak = (pmPeak *) tmpListElem->data; + if ((tmpPeak->counts > maxValue) || + (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) { + psListRemoveData(peaks, (psPtr) tmpPeak); + } + + indexNum++; + tmpListElem = tmpListElem->next; + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(peaks); +} + +// XXX EAM: I changed this to return a new, subset array +// rather than alter the existing one +// XXX: Fix the *valid pointer. +psArray *pmPeaksSubset( + psArray *peaks, + psF32 maxValue, + const psRegion valid) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(peaks, NULL); + + psArray *output = psArrayAlloc (200); + output->n = 0; + + psTrace (".pmObjects.pmCullPeaks", 3, "list size is %d\n", peaks->n); + + for (int i = 0; i < peaks->n; i++) { + pmPeak *tmpPeak = (pmPeak *) peaks->data[i]; + if (tmpPeak->counts > maxValue) + continue; + if (isItInThisRegion(valid, tmpPeak->x, tmpPeak->y)) + continue; + psArrayAdd (output, 200, tmpPeak); + } + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(output); +} Index: trunk/psModules/src/objects/pmPeaks.h =================================================================== --- trunk/psModules/src/objects/pmPeaks.h (revision 6872) +++ trunk/psModules/src/objects/pmPeaks.h (revision 6872) @@ -0,0 +1,143 @@ +/** @file pmPeaks.h + * + * The process of finding, measuring, and classifying astronomical sources on + * images is one of the critical tasks of the IPP or any astronomical software + * system. This file will define structures and functions related to the task + * of source detection and measurement. The elements defined in this section + * are generally low-level components which can be connected together to + * construct a complete object measurement suite. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_PEAKS_H +# define PM_PEAKS_H + +/** pmPeakType + * + * A peak pixel may have several features which may be determined when the + * peak is found or measured. These are specified by the pmPeakType enum. + * PM_PEAK_LONE represents a single pixel which is higher than its 8 immediate + * neighbors. The PM_PEAK_EDGE represents a peak pixel which touching the image + * edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific + * number of neighbors at the same value, within some tolarence: + * + */ +typedef enum { + PM_PEAK_LONE, ///< Isolated peak. + PM_PEAK_EDGE, ///< Peak on edge. + PM_PEAK_FLAT, ///< Peak has equal-value neighbors. + PM_PEAK_UNDEF ///< Undefined. +} pmPeakType; + + +/** pmPeak data structure + * + * A source has the capacity for several types of measurements. The + * simplest measurement of a source is the location and flux of the peak pixel + * associated with the source: + * + */ +typedef struct +{ + int x; ///< X-coordinate of peak pixel. + int y; ///< Y-coordinate of peak pixel. + float counts; ///< Value of peak pixel (above sky?). + pmPeakType type; ///< Description of peak. +} +pmPeak; + + +/** pmPeakAlloc() + * + * @return pmPeak* newly allocated pmPeak with all internal pointers set to NULL + */ +pmPeak *pmPeakAlloc( + int x, ///< Row-coordinate in image space + int y, ///< Col-coordinate in image space + float counts, ///< The value of the peak pixel + pmPeakType type ///< The type of peak pixel +); + + +/** pmFindVectorPeaks() + * + * Find all local peaks in the given vector above the given threshold. A peak + * is defined as any element with a value greater than its two neighbors and with + * a value above the threshold. Two types of special cases must be addressed. + * Equal value elements: If an element has the same value as the following + * element, it is not considered a peak. If an element has the same value as the + * preceding element (but not the following), then it is considered a peak. Note + * that this rule (arbitrarily) identifies flat regions by their trailing edge. + * Edge cases: At start of the vector, the element must be higher than its + * neighbor. At the end of the vector, the element must be higher or equal to its + * neighbor. These two rules again places the peak associated with a flat region + * which touches the image edge at the image edge. The result of this function is + * a vector containing the coordinates (element number) of the detected peaks + * (type psU32). + * + */ +psVector *pmFindVectorPeaks( + const psVector *vector, ///< The input vector (float) + float threshold ///< Threshold above which to find a peak +); + + +/** pmFindImagePeaks() + * + * Find all local peaks in the given image above the given threshold. This + * function should find all row peaks using pmFindVectorPeaks, then test each row + * peak and exclude peaks which are not local peaks. A peak is a local peak if it + * has a higher value than all 8 neighbors. If the peak has the same value as its + * +y neighbor or +x neighbor, it is NOT a local peak. If any other neighbors + * have an equal value, the peak is considered a valid peak. Note two points: + * first, the +x neighbor condition is already enforced by pmFindVectorPeaks. + * Second, these rules have the effect of making flat-topped regions have single + * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be + * set. The result of this function is an array of pmPeak entries. + * + */ +psArray *pmFindImagePeaks( + const psImage *image, ///< The input image where peaks will be found (float) + float threshold ///< Threshold above which to find a peak +); + + +/** pmCullPeaks() + * + * Eliminate peaks from the psList that have a peak value above the given + * maximum, or fall outside the valid region. + * + */ +psList *pmCullPeaks( + psList *peaks, ///< The psList of peaks to be culled + float maxValue, ///< Cull peaks above this value + const psRegion valid ///< Cull peaks otside this psRegion +); + + +/** pmPeaksSubset() + * + * Create a new peaks array, removing certain types of peaks from the input + * array of peaks based on the given criteria. Peaks should be eliminated if they + * have a peak value above the given maximum value limit or if the fall outside + * the valid region. The result of the function is a new array with a reduced + * number of peaks. + * + */ +psArray *pmPeaksSubset( + psArray *peaks, ///< Add comment. + float maxvalue, ///< Add comment. + const psRegion valid ///< Add comment. +); + +int pmPeaksCompareAscend (const void **a, const void **b); +int pmPeaksCompareDescend (const void **a, const void **b); + +# endif /* PM_PEAKS_H */ Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 6872) +++ trunk/psModules/src/objects/pmSource.c (revision 6872) @@ -0,0 +1,609 @@ +/** @file pmSource.c + * + * Functions to define and manipulate sources on images + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" + +static void sourceFree(pmSource *tmp) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + psFree(tmp->peak); + psFree(tmp->pixels); + psFree(tmp->weight); + psFree(tmp->mask); + psFree(tmp->moments); + psFree(tmp->modelPSF); + psFree(tmp->modelEXT); + psFree(tmp->blends); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); +} + +/****************************************************************************** +pmSourceAlloc(): Allocate the pmSource structure and initialize its members +to NULL. +*****************************************************************************/ +pmSource *pmSourceAlloc() +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmSource *tmp = (pmSource *) psAlloc(sizeof(pmSource)); + tmp->peak = NULL; + tmp->pixels = NULL; + tmp->weight = NULL; + tmp->mask = NULL; + tmp->moments = NULL; + tmp->blends = NULL; + tmp->modelPSF = NULL; + tmp->modelEXT = NULL; + tmp->type = PM_SOURCE_TYPE_UNKNOWN; + tmp->mode = PM_SOURCE_MODE_DEFAULT; + psMemSetDeallocator(tmp, (psFreeFunc) sourceFree); + + tmp->psfMag = 32.0; + tmp->extMag = 32.0; + tmp->errMag = 32.0; + tmp->apMag = 32.0; + tmp->pixWeight = 0.0; + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +// x,y are defined in the parent image coords of readout->image +bool pmSourceDefinePixels(pmSource *mySource, + const pmReadout *readout, + psF32 x, + psF32 y, + psF32 Radius) +{ + psRegion srcRegion; + + // Grab a subimage of the original image of size (2 * outerRadius). + srcRegion = psRegionForSquare (x, y, Radius); + srcRegion = psRegionForImage (readout->image, srcRegion); + + mySource->pixels = psImageSubset(readout->image, srcRegion); + mySource->weight = psImageSubset(readout->weight, srcRegion); + mySource->mask = psImageSubset(readout->mask, srcRegion); + mySource->region = srcRegion; + + return true; +} + +bool pmSourceRedefinePixels(pmSource *mySource, + const pmReadout *readout, + psF32 x, + psF32 y, + psF32 Radius) +{ + bool extend; + psRegion newRegion; + + if (Radius == 0) + return false; + + // check to see if new region is completely contained within old region + newRegion = psRegionForSquare (x, y, Radius); + newRegion = psRegionForImage (readout->image, newRegion); + + extend = false; + extend |= (int)(newRegion.x0) < (int)(mySource->region.x0); + extend |= (int)(newRegion.x1) > (int)(mySource->region.x1); + extend |= (int)(newRegion.y0) < (int)(mySource->region.y0); + extend |= (int)(newRegion.y1) > (int)(mySource->region.y1); + + extend |= (mySource->pixels == NULL); + extend |= (mySource->weight == NULL); + extend |= (mySource->mask == NULL); + + // extend = true; + if (extend) { + // Grab a new subimage + psFree (mySource->pixels); + psFree (mySource->weight); + psFree (mySource->mask); + + mySource->pixels = psImageSubset(readout->image, newRegion); + mySource->weight = psImageSubset(readout->weight, newRegion); + mySource->mask = psImageSubset(readout->mask, newRegion); + mySource->region = newRegion; + } + + return extend; +} + +/****************************************************************************** + pmSourcePSFClump(source, metadata): Find the likely PSF clump in the + sigma-x, sigma-y plane. return 0,0 clump in case of error. +*****************************************************************************/ + +// XXX EAM include a S/N cutoff in selecting the sources? +pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + # define NPIX 10 + # define SCALE 0.1 + + psArray *peaks = NULL; + pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; + pmPSFClump psfClump = emptyClump; + + PS_ASSERT_PTR_NON_NULL(sources, emptyClump); + PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); + + // find the sigmaX, sigmaY clump + { + psStats *stats = NULL; + psImage *splane = NULL; + int binX, binY; + bool status; + + psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); + if (!status) + SX_MAX = 10.0; + psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); + if (!status) + SY_MAX = 10.0; + + // construct a sigma-plane image + // psImageAlloc does zero the data + splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32); + for (int i = 0; i < splane->numRows; i++) + { + memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); + } + + // place the sources in the sigma-plane image (ignore 0,0 values?) + for (psS32 i = 0 ; i < sources->n ; i++) + { + pmSource *tmpSrc = (pmSource *) sources->data[i]; + if (tmpSrc == NULL) { + continue; + } + if (tmpSrc->moments == NULL) { + continue; + } + + // Sx,Sy are limited at 0. a peak at 0,0 is artificial + if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { + continue; + } + + // for the moment, force splane dimensions to be 10x10 image pix + binX = tmpSrc->moments->Sx/SCALE; + if (binX < 0) + continue; + if (binX >= splane->numCols) + continue; + + binY = tmpSrc->moments->Sy/SCALE; + if (binY < 0) + continue; + if (binY >= splane->numRows) + continue; + + splane->data.F32[binY][binX] += 1.0; + } + + // find the peak in this image + stats = psStatsAlloc (PS_STAT_MAX); + stats = psImageStats (stats, splane, NULL, 0); + peaks = pmFindImagePeaks (splane, stats[0].max / 2); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); + + psFree (splane); + psFree (stats); + + } + // XXX EAM : possible errors: + // 1) no peak in splane + // 2) no significant peak in splane + + // measure statistics on Sx, Sy if Sx, Sy within range of clump + { + pmPeak *clump; + psF32 minSx, maxSx; + psF32 minSy, maxSy; + psVector *tmpSx = NULL; + psVector *tmpSy = NULL; + psStats *stats = NULL; + + // XXX EAM : this lets us takes the single highest peak + psArraySort (peaks, pmPeaksCompareDescend); + clump = peaks->data[0]; + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); + + // define section window for clump + minSx = clump->x * SCALE - 0.2; + maxSx = clump->x * SCALE + 0.2; + minSy = clump->y * SCALE - 0.2; + maxSy = clump->y * SCALE + 0.2; + + tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); + tmpSx->n = 0; + tmpSy->n = 0; + + // XXX clip sources based on flux? + // create vectors with Sx, Sy values in window + for (psS32 i = 0 ; i < sources->n ; i++) + { + pmSource *tmpSrc = (pmSource *) sources->data[i]; + + if (tmpSrc->moments->Sx < minSx) + continue; + if (tmpSrc->moments->Sx > maxSx) + continue; + if (tmpSrc->moments->Sy < minSy) + continue; + if (tmpSrc->moments->Sy > maxSy) + continue; + tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; + tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; + tmpSx->n++; + tmpSy->n++; + if (tmpSx->n == tmpSx->nalloc) { + psVectorRealloc (tmpSx, tmpSx->nalloc + 100); + psVectorRealloc (tmpSy, tmpSy->nalloc + 100); + } + } + + // measures stats of Sx, Sy + stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); + + stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); + psfClump.X = stats->clippedMean; + psfClump.dX = stats->clippedStdev; + + stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); + psfClump.Y = stats->clippedMean; + psfClump.dY = stats->clippedStdev; + + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); + psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); + // these values should be pushed on the metadata somewhere + + psFree (stats); + psFree (peaks); + psFree (tmpSx); + psFree (tmpSy); + } + + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return (psfClump); +} + +/****************************************************************************** + pmSourceRoughClass(source, metadata): make a guess at the source + classification. + + XXX: push the clump info into the metadata? + + XXX: How can this function ever return FALSE? + + EAM: I moved S/N calculation to pmSourceMoments, using weight image +*****************************************************************************/ + +bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + + psBool rc = true; + + int Nsat = 0; + int Next = 0; + int Nstar = 0; + int Npsf = 0; + int Ncr = 0; + int Nsatstar = 0; + // psRegion allArray = psRegionSet (0, 0, 0, 0); + psRegion inner; + + // report stats on S/N values for star-like objects + psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); + starsn->n = 0; + + // check return status value (do these exist?) + bool status; + psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); + + // XXX allow clump size to be scaled relative to sigmas? + // make rough IDs based on clumpX,Y,DX,DY + for (psS32 i = 0 ; i < sources->n ; i++) { + + pmSource *tmpSrc = (pmSource *) sources->data[i]; + + tmpSrc->peak->type = 0; + + psF32 sigX = tmpSrc->moments->Sx; + psF32 sigY = tmpSrc->moments->Sy; + + // XXX EAM : can we use the value of SATURATE if mask is NULL? + // inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); + inner = psRegionForSquare (tmpSrc->peak->x, tmpSrc->peak->y, 2); + int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PM_SOURCE_MASK_SATURATED); + + // saturated star (size consistent with PSF or larger) + // Nsigma should be user-configured parameter + bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); + big = true; + if ((Nsatpix > 1) && big) { + tmpSrc->type = PM_SOURCE_TYPE_STAR; + tmpSrc->mode = PM_SOURCE_MODE_SATSTAR; + Nsatstar ++; + continue; + } + + // saturated object (not a star, eg bleed trails, hot pixels) + if (Nsatpix > 1) { + tmpSrc->type = PM_SOURCE_TYPE_SATURATED; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; + Nsat ++; + continue; + } + + // likely defect (too small to be stellar) (push out to 3 sigma) + // low S/N objects which are small are probably stellar + // only set candidate defects if + if ((sigX < 0.05) || (sigY < 0.05)) { + tmpSrc->type = PM_SOURCE_TYPE_DEFECT; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; + Ncr ++; + continue; + } + + // likely unsaturated extended source (too large to be stellar) + if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { + tmpSrc->type = PM_SOURCE_TYPE_EXTENDED; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; + Next ++; + continue; + } + + // the rest are probable stellar objects + starsn->data.F32[starsn->n] = tmpSrc->moments->SN; + starsn->n ++; + Nstar ++; + + // PSF star (within 1.5 sigma of clump center, S/N > limit) + psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); + if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { + tmpSrc->type = PM_SOURCE_TYPE_STAR; + tmpSrc->mode = PM_SOURCE_MODE_PSFSTAR; + Npsf ++; + continue; + } + + // random type of star + tmpSrc->type = PM_SOURCE_TYPE_STAR; + tmpSrc->mode = PM_SOURCE_MODE_DEFAULT; + } + + { + psStats *stats = NULL; + stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); + stats = psVectorStats (stats, starsn, NULL, NULL, 0); + psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); + psFree (stats); + psFree (starsn); + } + + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); + psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); + + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +/****************************************************************************** +pmSourceMoments(source, radius): this function takes a subImage defined in the +pmSource data structure, along with the peak location, and determines the +various moments associated with that peak. + +Requires the following to have been created: + pmSource + pmSource->peak + pmSource->pixels + pmSource->weight + pmSource->mask + +XXX: The peak calculations are done in image coords, not subImage coords. + +XXX EAM : this version clips input pixels on S/N +XXX EAM : this version returns false for several reasons +*****************************************************************************/ +# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) + +bool pmSourceMoments(pmSource *source, + psF32 radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false); + + // + // XXX: Verify the setting for sky if source->moments == NULL. + // + psF32 sky = 0.0; + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } else { + sky = source->moments->Sky; + } + + // + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // X2 = SUM (x - xc)^2 * (z - sky) + // XY = SUM (x - xc)*(y - yc)*(z - sky) + // + psF32 peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + psF32 Sum = 0.0; + psF32 Var = 0.0; + psF32 X1 = 0.0; + psF32 Y1 = 0.0; + psF32 X2 = 0.0; + psF32 Y2 = 0.0; + psF32 XY = 0.0; + psF32 x = 0; + psF32 y = 0; + psF32 R2 = PS_SQR(radius); + + psF32 xPeak = source->peak->x; + psF32 yPeak = source->peak->y; + psF32 xOff = source->pixels->col0 - source->peak->x; + psF32 yOff = source->pixels->row0 - source->peak->y; + + // XXX why do I get different results for these two methods of finding Sx? + // XXX Sx, Sy would be better measured if we clip pixels close to sky + // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? + // We loop through all pixels in this subimage (source->pixels), and for each + // pixel that is not masked, AND within the radius of the peak pixel, we + // proceed with the moments calculation. need to do two loops for a + // numerically stable result. first loop: get the sums. + // XXX EAM : mask == 0 is valid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + + psF32 *vPix = source->pixels->data.F32[row]; + psF32 *vWgt = source->weight->data.F32[row]; + psU8 *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row]; + + for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) { + if ((vMsk != NULL) && *vMsk) { + vMsk++; + continue; + } + + psF32 xDiff = col + xOff; + psF32 yDiff = row + yOff; + + // radius is just a function of (xDiff, yDiff) + if (!VALID_RADIUS(xDiff, yDiff, R2)) { + if (vMsk != NULL) + vMsk++; + continue; + } + + psF32 pDiff = *vPix - sky; + psF32 wDiff = *vWgt; + + // XXX EAM : check for valid S/N in pixel + // XXX EAM : should this limit be user-defined? + if (PS_SQR(pDiff) < wDiff) { + if (vMsk != NULL) + vMsk++; + continue; + } + + Var += wDiff; + Sum += pDiff; + + psF32 xWght = xDiff * pDiff; + psF32 yWght = yDiff * pDiff; + + X1 += xWght; + Y1 += yWght; + + XY += xDiff * yWght; + X2 += xDiff * xWght; + Y2 += yDiff * yWght; + + peakPixel = PS_MAX (*vPix, peakPixel); + numPixels++; + if (vMsk != NULL) + vMsk++; + } + } + + // if we have less than (1/4) of the possible pixels, force a retry + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 0.75*R2) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + psTrace (".psModules.pmSourceMoments", 5, + "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", + sky, Sum, X1, Y1, X2, Y2, XY, numPixels); + + // + // first moment X = X1/Sum + xc + // second moment X = sqrt (X2/Sum - (X1/Sum)^2) + // Sxy = XY / Sum + // + x = X1/Sum; + y = Y1/Sum; + if ((fabs(x) > radius) || (fabs(y) > radius)) { + psTrace (".psModules.pmSourceMoments", 3, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return (false); + } + + source->moments->x = x + xPeak; + source->moments->y = y + yPeak; + + // XXX EAM : Sxy needs to have x*y subtracted + source->moments->Sxy = XY/Sum - x*y; + source->moments->Sum = Sum; + source->moments->SN = Sum / sqrt(Var); + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + // XXX EAM : these values can be negative, so we need to limit the range + source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); + source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); + + psTrace (".psModules.pmSourceMoments", 4, + "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", + sky, Sum, source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, source->moments->Sxy); + + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return(true); +} + +pmModel *pmSourceSelectModel (pmSource *source) +{ + switch (source->type) { + case PM_SOURCE_TYPE_STAR: + return source->modelPSF; + + case PM_SOURCE_TYPE_EXTENDED: + return source->modelEXT; + + default: + return NULL; + } + return NULL; +} Index: trunk/psModules/src/objects/pmSource.h =================================================================== --- trunk/psModules/src/objects/pmSource.h (revision 6872) +++ trunk/psModules/src/objects/pmSource.h (revision 6872) @@ -0,0 +1,216 @@ +/** @file pmSource.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_H +# define PM_SOURCE_H + +/** + * In the object analysis process, we will use specific mask values to mark the + * image pixels. The following structure defines the relevant mask values. + * + * XXX: This is probably a bad solution: we will want to set mask values + * outside of the PSPHOT code. Perhaps we can set up a registered set of mask + * values with specific meanings that other functions can add to or define? + */ +typedef enum { + PM_SOURCE_MASK_CLEAR = 0x00, + PM_SOURCE_MASK_INVALID = 0x01, + PM_SOURCE_MASK_SATURATED = 0x02, + PM_SOURCE_MASK_MARKED = 0x08, +} psphotMaskValues; + +/** pmSourceType enumeration + * + * A given source may be identified as most-likely to be one of several source + * types. The pmSource entry pmSourceType defines the current best-guess for this + * source. + * + * XXX: The values given below are currently illustrative and will require + * some modification as the source classification code is developed. (TBD) + * + */ +typedef enum { + PM_SOURCE_TYPE_UNKNOWN, ///< a cosmic-ray + PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray + PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels + PM_SOURCE_TYPE_STAR, ///< a good-quality star + PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) +} pmSourceType; + +typedef enum { + PM_SOURCE_MODE_DEFAULT = 0x0000, ///< + PM_SOURCE_MODE_PSFMODEL = 0x0001, ///< + PM_SOURCE_MODE_EXTMODEL = 0x0002, ///< + PM_SOURCE_MODE_SUBTRACTED = 0x0004, ///< + PM_SOURCE_MODE_FITTED = 0x0008, ///< + PM_SOURCE_MODE_FAIL = 0x0010, ///< + PM_SOURCE_MODE_POOR = 0x0020, ///< + PM_SOURCE_MODE_PAIR = 0x0040, ///< + PM_SOURCE_MODE_PSFSTAR = 0x0080, ///< + PM_SOURCE_MODE_SATSTAR = 0x0100, ///< + PM_SOURCE_MODE_BLEND = 0x0200, ///< + PM_SOURCE_MODE_LINEAR = 0x0400, ///< + PM_SOURCE_MODE_TEMPSUB = 0x0800, ///< XXX get me a better name! +} pmSourceMode; + +/** pmSource data structure + * + * This source has the capacity for several types of measurements. The + * simplest measurement of a source is the location and flux of the peak pixel + * associated with the source: + * + */ +typedef struct +{ + pmPeak *peak; ///< Description of peak pixel. + psImage *pixels; ///< Rectangular region including object pixels. + psImage *weight; ///< Image variance. + psImage *mask; ///< Mask which marks pixels associated with objects. + pmMoments *moments; ///< Basic moments measure for the object. + pmModel *modelPSF; ///< PSF Model fit (parameters and type) + pmModel *modelEXT; ///< EXT (floating) Model fit (parameters and type). + pmSourceType type; ///< Best identification of object. + pmSourceMode mode; ///< Best identification of object. + psArray *blends; + float psfMag; + float extMag; + float errMag; + float apMag; + float pixWeight; // model-weighted coverage of valid pixels + psRegion region; // area on image covered by selected pixels +} +pmSource; + +/** pmPSFClump data structure + * + * A collection of object moment measurements can be used to determine + * approximate object classes. The key to this analysis is the location and + * statistics (in the second-moment plane, + * + */ +typedef struct +{ + float X; + float dX; + float Y; + float dY; +} +pmPSFClump; + + +/** pmSourceAlloc() + * + */ +pmSource *pmSourceAlloc(); + + +/** pmSourceDefinePixels() + * + * Define psImage subarrays for the source located at coordinates x,y on the + * image set defined by readout. The pixels defined by this operation consist of + * a square window (of full width 2Radius+1) centered on the pixel which contains + * the given coordinate, in the frame of the readout. The window is defined to + * have limits which are valid within the boundary of the readout image, thus if + * the radius would fall outside the image pixels, the subimage is truncated to + * only consist of valid pixels. If readout->mask or readout->weight are not + * NULL, matching subimages are defined for those images as well. This function + * fails if no valid pixels can be defined (x or y less than Radius, for + * example). This function should be used to define a region of interest around a + * source, including both source and sky pixels. + * + */ +bool pmSourceDefinePixels( + pmSource *mySource, ///< source to be re-defined + const pmReadout *readout, ///< base the source on this readout + psF32 x, ///< center coords of source + psF32 y, ///< center coords of source + psF32 Radius ///< size of box on source +); + +bool pmSourceRedefinePixels ( + pmSource *mySource, ///< source to be re-defined + const pmReadout *readout, ///< base the source on this readout + psF32 x, ///< center coords of source + psF32 y, ///< center coords of source + psF32 Radius ///< size of box on source +); + +/** pmSourcePSFClump() + * + * We use the source moments to make an initial, approximate source + * classification, and as part of the information needed to build a PSF model for + * the image. As long as the PSF shape does not vary excessively across the + * image, the sources which are represented by a PSF (the start) will have very + * similar second moments. The function pmSourcePSFClump searches a collection of + * sources with measured moments for a group with moments which are all very + * similar. The function returns a pmPSFClump structure, representing the + * centroid and size of the clump in the sigma_x, sigma_y second-moment plane. + * + * The goal is to identify and characterize the stellar clump within the + * sigma_x, sigma_y second-moment plane. To do this, an image is constructed to + * represent this plane. The units of sigma_x and sigma_y are in image pixels. A + * pixel in this analysis image represents 0.1 pixels in the input image. The + * dimensions of the image need only be 10 pixels. The peak pixel in this image + * (above a threshold of half of the image maximum) is found. The coordinates of + * this peak pixel represent the 2D mode of the sigma_x, sigma_y distribution. + * The sources with sigma_x, sigma_y within 0.2 pixels of this value are then + * * used to calculate the median and standard deviation of the sigma_x, sigma_y + * values. These resulting values are returned via the pmPSFClump structure. + * + * The return value indicates the success (TRUE) of the operation. + * + * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD). + * XXX: Save the clump parameters on the Metadata (TBD) + * + */ +pmPSFClump pmSourcePSFClump( + psArray *source, ///< The input pmSource + psMetadata *metadata ///< Contains classification parameters +); + + +/** pmSourceRoughClass() + * + * Based on the specified data values, make a guess at the source + * classification. The sources are provides as a psArray of pmSource entries. + * Definable parameters needed to make the classification are provided to the + * routine with the psMetadata structure. The rules (in SDRS) refer to values which + * can be extracted from the metadata using the given keywords. Except as noted, + * the data type for these parameters are psF32. + * + */ +bool pmSourceRoughClass( + psArray *source, ///< The input pmSource + psMetadata *metadata, ///< Contains classification parameters + pmPSFClump clump ///< Statistics about the PSF clump +); + + +/** pmSourceMoments() + * + * Measure source moments for the given source, using the value of + * source.moments.sky provided as the local background value and the peak + * coordinates as the initial source location. The resulting moment values are + * applied to the source.moments entry, and the source is returned. The moments + * are measured within the given circular radius of the source.peak coordinates. + * The return value indicates the success (TRUE) of the operation. + * + */ +bool pmSourceMoments( + pmSource *source, ///< The input pmSource for which moments will be computed + float radius ///< Use a circle of pixels around the peak +); + + +// select the model used for this source +pmModel *pmSourceSelectModel (pmSource *source); + +# endif /* PM_SOURCE_H */ Index: trunk/psModules/src/objects/pmSourceContour.c =================================================================== --- trunk/psModules/src/objects/pmSourceContour.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceContour.c (revision 6872) @@ -0,0 +1,446 @@ +/** @file pmSourceContour.c + * + * Functions to measure the local sky and sky variance for sources on images + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceContour.h" + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +# define LEFT false +# define RIGHT true + +// return the first coordinate at or below the threshold in the requested direction +static int findContourNeg( + psImage *image, + float threshold, + int x, + int y, + bool right) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = image->numCols - 1; + } else { + incr = -1; + lastColumn = 0; + } + + subCol = x; + + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value <= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +// return the last coordinate at or below the threshold in the requested direction +static int findContourPos( + psImage *image, + float threshold, + int x, + int y, + bool right, + int xEnd) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = PS_MIN (image->numCols - 1, xEnd); + } else { + incr = -1; + lastColumn = PS_MAX (0, xEnd); + } + + subCol = x; + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value >= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + if (subCol == x) { + return (subCol); + } + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +static psF32 findValue(pmSource *source, + psF32 level, + psU32 row, + psU32 col, + psU32 dir) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + // + // Convert coords to subImage space. + // + psU32 subRow = row - source->pixels->row0; + psU32 subCol = col - source->pixels->col0; + + // Ensure that the starting column is allowable. + if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { + psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + return(NAN); + } + if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); + return(NAN); + } + + // XXX EAM : i changed this to match pmModelEval above, but see + // XXX EAM the note below in pmSourceContour + psF32 oldValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); + if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(((psF32) (subCol + source->pixels->col0))); + } + + // + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + // + psS32 incr; + psS32 lastColumn; + if (dir == 0) { + incr = -1; + lastColumn = -1; + } else { + incr = 1; + lastColumn = source->pixels->numCols; + } + subCol+=incr; + + while (subCol != lastColumn) { + psF32 newValue = pmModelEval(source->modelEXT, source->pixels, subCol, subRow); + if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return((psF32) (subCol + source->pixels->col0)); + } + + if ((newValue <= level) && (level <= oldValue)) { + // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); + } + + if ((oldValue <= level) && (level <= newValue)) { + // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); + } + + subCol+=incr; + } + + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + return(NAN); +} + +/****************************************************************************** +new implementation of source contour function +*****************************************************************************/ +psArray *pmSourceContour (psImage *image, int xc, int yc, float threshold) +{ + + int xR, yR, x0, x1, x0s, x1s; + + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + + int x = xc - image->col0; + int y = yc - image->row0; + + // the requested point must be within the contour + if (image->data.F32[y][x] < threshold) + return NULL; + + // Ensure that the starting column is allowable. + if (x < 0) + return NULL; + if (y < 0) + return NULL; + if (x >= image->numCols) + return NULL; + if (y >= image->numRows) + return NULL; + + // Allocate data for x/y pairs. + psVector *xVec = psVectorAlloc(100, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(100, PS_TYPE_F32); + + // First row: find the left and right end-points + int Npt = 0; + + x0 = findContourNeg (image, threshold, x, y, LEFT); + x1 = findContourNeg (image, threshold, x, y, RIGHT); + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + y; + yVec->data.F32[Npt + 1] = image->row0 + y; + Npt += 2; + + x0s = x0; + x1s = x1; + + // look for contour outline above row + xR = x0s; + yR = y + 1; + while (yR < image->numRows) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt1; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR ++; + xR = x0; + } + +pt1: + // look for contour outline below row + xR = x0s; + x1 = x1s; + yR = y - 1; + while (yR >= 0) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt2; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + // fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR --; + } +pt2: + xVec->n = Npt; + yVec->n = Npt; + + // fprintf (stderr, "done\n"); + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) xVec; + tmpArray->data[1] = (psPtr *) yVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); +} + +/****************************************************************************** + pmSourceContour(src, img, level, mode): For an input subImage, and model, this + routine returns a psArray of coordinates that evaluate to the specified level. + + XXX: Probably should remove the "image" argument. + XXX: What type should the output coordinate vectors consist of? col,row? + XXX: Why a pmArray output? + XXX: doex x,y correspond with col,row or row/col? + XXX: What is mode? + XXX: The top, bottom of the contour is not correctly determined. + XXX EAM : this function is using the model for the contour, but it should + be using only the image counts +*****************************************************************************/ +psArray *pmSourceContour_Crude(pmSource *source, + const psImage *image, + psF32 level) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->modelEXT, false); + // XXX EAM : what is the purpose of modelPSF/modelEXT? + + // + // Allocate data for x/y pairs. + // + psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); + + // + // Start at the row with peak pixel, then decrement. + // + psS32 col = source->peak->x; + for (psS32 row = source->peak->y; row>= 0 ; row--) { + // XXX: yVec contain no real information. Do we really need it? + yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); + yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); + + // Starting at peak pixel, search leftwards for the column intercept. + psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; + + // Starting at peak pixel, search rightwards for the column intercept. + + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); + xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; + + // Set starting column for next row + col = (psS32) ((leftIntercept + rightIntercept) / 2.0); + } + // + // Start at the row (+1) with peak pixel, then increment. + // + col = source->peak->x; + for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { + // XXX: yVec contain no real information. Do we really need it? + yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); + yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); + + // Starting at peak pixel, search leftwards for the column intercept. + psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; + + // Starting at peak pixel, search rightwards for the column intercept. + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); + return(NULL); + //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); + } + xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; + + // Set starting column for next row + col = (psS32) ((leftIntercept + rightIntercept) / 2.0); + } + + // + // Allocate an array for result, store coord vectors there. + // + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) yVec; + tmpArray->data[1] = (psPtr *) xVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); +} Index: trunk/psModules/src/objects/pmSourceContour.h =================================================================== --- trunk/psModules/src/objects/pmSourceContour.h (revision 6872) +++ trunk/psModules/src/objects/pmSourceContour.h (revision 6872) @@ -0,0 +1,37 @@ +/** @file pmSourceContour.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_CONTOUR_H +# define PM_SOURCE_CONTOUR_H + +psArray *pmSourceContour (psImage *image, int xc, int yc, float threshold); + + +/** pmSourceContour() + * + * Find points in a contour for the given source at the given level. If type + * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak, + * running along each pixel row until the level is crossed, then interpolating to + * the level coordinate for that row. This is done for each row, with the + * starting point determined by the midpoint of the previous row, until the + * starting point has a value below the contour level. The returned contour + * consists of two vectors giving the x and y coordinates of the contour levels. + * This function may be used as part of the model guess inputs. Other contour + * types may be specified in the future for more refined contours (TBD) + * + */ +psArray *pmSourceContour_Crude( + pmSource *source, ///< The input pmSource + const psImage *image, ///< The input image (float) (this arg should be removed) + float level ///< The level of the contour +); + +# endif /* PM_SOURCE_PHOTOMETRY_H */ Index: trunk/psModules/src/objects/pmSourceFitModel.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceFitModel.c (revision 6872) @@ -0,0 +1,559 @@ +/** @file pmSourceFitModel.c + * + * fit single source models to image pixels + * + * @author EAM, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceFitModel.h" + +// save as static values so they may be set externally +static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; +static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; +static bool PM_PSF_POISSON_WEIGHTS = true; + +bool pmSourceFitModelInit (float nIter, float tol, bool poissonErrors) +{ + + PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; + PM_SOURCE_FIT_MODEL_TOLERANCE = tol; + PM_PSF_POISSON_WEIGHTS = poissonErrors; + + return true; +} + +bool pmSourceFitModel (pmSource *source, + pmModel *model, + pmSourceFitMode mode) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, 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); + + // XXX EAM : is it necessary for the mask & weight to exist? the + // tests below could be conditions (!NULL) + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + + psF32 dSky = source->moments->dSky; + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (PM_PSF_POISSON_WEIGHTS) { + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + } else { + yErr->data.F32[nPix] = 1.0 / dSky; + } + // XXX EAM : suggestion from RHL is to use the local sky as weight + // to avoid the bias from systematic errors + // here we would just use the source sky variance + nPix++; + } + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + // XXX EAM : the new minimization API supplies the constraints as a struct + // XXX the chisq if a fcn of source flux. the minimization criterion should + // probably be scaled somehow by flux. measure the trend? it depends on + // the about of systematic error in the models themselves... + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + psMinConstrain *constrain = psMinConstrainAlloc(); + + // set parameter mask based on fitting mode + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + psVectorInit (paramMask, 1); + + int nParams = 0; + switch (mode) { + case PM_SOURCE_FIT_NORM: + // NORM-only model fits only source normalization (Io) + nParams = 1; + paramMask->data.U8[1] = 0; + break; + case PM_SOURCE_FIT_PSF: + // PSF model only fits x,y,Io + nParams = 3; + paramMask->data.U8[1] = 0; + paramMask->data.U8[2] = 0; + paramMask->data.U8[3] = 0; + break; + case PM_SOURCE_FIT_EXT: + // EXT model fits all params (except sky) + nParams = params->n - 1; + psVectorInit (paramMask, 0); + paramMask->data.U8[0] = 1; + break; + case PM_SOURCE_FIT_PSF_AND_SKY: + nParams = 4; + psAbort ("pmSourceFitModel", "PSF + SKY not currently available"); + break; + case PM_SOURCE_FIT_EXT_AND_SKY: + nParams = params->n; + psAbort ("pmSourceFitModel", "EXT + SKY not currently available"); + break; + default: + psAbort ("pmSourceFitModel", "invalid fitting mode"); + } + constrain->paramMask = paramMask; + + if (nPix < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + psFree (myMin); + psFree (constrain); + psFree (paramMask); + return(false); + } + + // Set the parameter range checks + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + + fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + // set the model success or failure status + if (!fitStatus) { + model->status = PM_MODEL_NONCONVERGE; + } else { + model->status = PM_MODEL_SUCCESS; + } + + // models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + + source->mode |= PM_SOURCE_MODE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + + rc = (onPic && fitStatus); + psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); + return(rc); +} + +/********************* Source Model Set Functions ***************************/ + +// these static variables are used by one pass of pmSourceFitSet to store +// data for a model set. If we are going to make psphot thread-safe, these +// will have to go in a structure of their own. +// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, +// nPar = nSrc*(nOnePar - 1) + 1 +static pmModelFunc mFunc; +static int nPar; +static psVector *onePar; +static psVector *oneDeriv; + +bool pmModelFitSetInit (pmModelType type) +{ + + mFunc = pmModelFunc_GetFunction (type); + nPar = pmModelParameterCount (type); + + onePar = psVectorAlloc (nPar, PS_DATA_F32); + oneDeriv = psVectorAlloc (nPar, PS_DATA_F32); + + return true; +} + +void pmModelFitSetClear (void) +{ + + psFree (onePar); + psFree (oneDeriv); + return; +} + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x) +{ + + psF32 value; + psF32 model; + + psF32 *PAR = onePar->data.F32; + psF32 *dPAR = oneDeriv->data.F32; + + psF32 *pars = params->data.F32; + psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32; + + int nSrc = (params->n - 1) / (nPar - 1); + + PAR[0] = model = pars[0]; + for (int i = 0; i < nSrc; i++) { + int nOff = i*nPar - i; + for (int n = 1; n < nPar; n++) { + PAR[n] = pars[nOff + n]; + } + if (deriv == NULL) { + value = mFunc (NULL, onePar, x); + } else { + value = mFunc (oneDeriv, onePar, x); + for (int n = 1; n < nPar; n++) { + dpars[nOff + n] = dPAR[n]; + } + } + model += value; + } + if (deriv != NULL) { + dpars[0] = dPAR[0]*2.0; + } + return (model); +} + +/* +i: 0 1 2 +n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 +i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 +*/ + +bool pmSourceFitSet (pmSource *source, + psArray *modelSet, + pmSourceFitMode mode) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, 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); + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // base values on first model + pmModel *model = modelSet->data[0]; + + // set the static variables + pmModelFitSetInit (model->type); + + int nSrc = modelSet->n; + int nPar = model->params->n - 1; // number of object parameters (excluding sky) + + // define parameter vectors for source set + psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + + // define limits for a single-source model + psVector *oneDelta; + psVector *oneParMin; + psVector *oneParMax; + modelLimits (&oneDelta, &oneParMin, &oneParMax); + + psMinConstrain *constrain = psMinConstrainAlloc(); + constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8); + constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + + // set the parameter guesses and constraints for the multiple models + // first the values for the single sky parameter + params->data.F32[0] = model->params->data.F32[0]; + constrain->paramMin->data.F32[0] = oneParMin->data.F32[0]; + constrain->paramMax->data.F32[0] = oneParMax->data.F32[0]; + constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0]; + + // next, the values for the source parameters + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + for (int n = 1; n < nPar + 1; n++) { + params->data.F32[i*nPar + n] = model->params->data.F32[n]; + constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n]; + constrain->paramMin->data.F32[i*nPar + n] = oneParMin->data.F32[n]; + constrain->paramMax->data.F32[i*nPar + n] = oneParMax->data.F32[n]; + } + } + psFree (oneDelta); + psFree (oneParMin); + psFree (oneParMax); + + if (psTraceGetLevel(__func__) >= 5) { + for (int i = 0; i < params->n; i++) { + fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]); + } + } + + // set the parameter masks based on the fitting mode + int nParams = 0; + switch (mode) { + case PM_SOURCE_FIT_NORM: + // NORM-only model fits only source normalization (Io) + nParams = nSrc; + psVectorInit (constrain->paramMask, 1); + for (int i = 0; i < nSrc; i++) { + constrain->paramMask->data.U8[1 + i*nPar] = 0; + } + break; + case PM_SOURCE_FIT_PSF: + // PSF model only fits x,y,Io + nParams = 3*nSrc; + psVectorInit (constrain->paramMask, 1); + for (int i = 0; i < nSrc; i++) { + constrain->paramMask->data.U8[1 + i*nPar] = 0; + constrain->paramMask->data.U8[2 + i*nPar] = 0; + constrain->paramMask->data.U8[3 + i*nPar] = 0; + } + break; + case PM_SOURCE_FIT_EXT: + // EXT model fits all params (except sky) + nParams = nPar*nSrc; + psVectorInit (constrain->paramMask, 0); + constrain->paramMask->data.U8[0] = 1; + break; + case PM_SOURCE_FIT_PSF_AND_SKY: + nParams = 1 + 3*nSrc; + psAbort ("pmSourceFitModel", "PSF + SKY not currently available"); + break; + case PM_SOURCE_FIT_EXT_AND_SKY: + nParams = 1 + nPar*nSrc; + psAbort ("pmSourceFitModel", "EXT + SKY not currently available"); + break; + default: + psAbort ("pmSourceFitModel", "invalid fitting mode"); + } + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // local sky variance + psF32 dSky = source->moments->dSky; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (PM_PSF_POISSON_WEIGHTS) { + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + } else { + yErr->data.F32[nPix] = 1.0 / dSky; + } + nPix++; + } + } + if (nPix < nParams + 1) { + psTrace (__func__, 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree (dparams); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + return(false); + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (__func__, 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet); + + // parameter errors from the covariance matrix + for (int i = 0; i < dparams->n; i++) { + if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // get the Gauss-Newton distance for fixed model parameters + if (constrain->paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet); + for (int i = 0; i < dparams->n; i++) { + if (!constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + // assign back the parameters to the models + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + model->params->data.F32[0] = params->data.F32[0]; + for (int n = 1; n < nPar + 1; n++) { + model->params->data.F32[n] = params->data.F32[i*nPar + n]; + model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; + } + // save the resulting chisq, nDOF, nIter + // these are not unique for any one source + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // set the model success or failure status + model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE; + + // models can go insane: reject these + onPic &= (model->params->data.F32[2] >= source->pixels->col0); + onPic &= (model->params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (model->params->data.F32[3] >= source->pixels->row0); + onPic &= (model->params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + } + + source->mode |= PM_SOURCE_MODE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + psFree(params); + psFree(dparams); + + // free static memory used by pmModelFitSet + pmModelFitSetClear (); + + rc = (onPic && fitStatus); + psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF); + psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc); + return(rc); +} + Index: trunk/psModules/src/objects/pmSourceFitModel.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitModel.h (revision 6872) +++ trunk/psModules/src/objects/pmSourceFitModel.h (revision 6872) @@ -0,0 +1,69 @@ +/** @file pmSourceFitModel.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_FIT_MODEL_H +# define PM_SOURCE_FIT_MODEL_H + +typedef enum { + PM_SOURCE_FIT_NORM, + PM_SOURCE_FIT_PSF, + PM_SOURCE_FIT_EXT, + PM_SOURCE_FIT_PSF_AND_SKY, + PM_SOURCE_FIT_EXT_AND_SKY +} pmSourceFitMode; + +bool pmSourceFitModelInit( + float nIter, ///< max number of allowed iterations + float tol, ///< convergence criterion + bool poissonErrors // use poisson errors for fits? +); + +/** pmSourceFitModel() + * + * Fit the requested model to the specified source. The starting guess for the model is given + * by the input source.model parameter values. The pixels of interest are specified by the + * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image + * data. The function returns TRUE on success or FALSE on failure. + * + */ +bool pmSourceFitModel( + pmSource *source, ///< The input pmSource + pmModel *model, ///< model to be fitted + pmSourceFitMode mode ///< define parameters to be fitted +); + + +// initialize data for a group of object models +bool pmModelFitSetInit (pmModelType type); + +// clear data for a group of object models +void pmModelFitSetClear (void); + +// function used to fit a group of object models +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x); + +/** pmSourceFitSet() + * + * Fit the requested model to the specified source. The starting guess for the model is given + * by the input source.model parameter values. The pixels of interest are specified by the + * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image + * data. The function returns TRUE on success or FALSE on failure. + * + */ +bool pmSourceFitSet( + pmSource *source, ///< The input pmSource + psArray *modelSet, ///< model to be fitted + pmSourceFitMode mode ///< define parameters to be fitted +); + +# endif /* PM_SOURCE_FIT_MODEL_H */ Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 6872) @@ -0,0 +1,334 @@ +/** @file pmSourceFitSet.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceFitSet.h" + +// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, +// nPar = nSrc*(nOnePar - 1) + 1 + +static pmModelFunc mFunc; +static int nPar; +static psVector *onePar; +static psVector *oneDeriv; + +// save as static values so they may be set externally +static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; +static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; +static bool PM_PSF_POISSON_WEIGHTS = true; + +bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors) +{ + + PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; + PM_SOURCE_FIT_MODEL_TOLERANCE = tol; + + PM_PSF_POISSON_WEIGHTS = poissonErrors; + + return true; +} + +bool pmModelFitSetInit (pmModelType type) +{ + + mFunc = pmModelFunc_GetFunction (type); + nPar = pmModelParameterCount (type); + + onePar = psVectorAlloc (nPar, PS_DATA_F32); + oneDeriv = psVectorAlloc (nPar, PS_DATA_F32); + + return true; +} + +void pmModelFitSetClear (void) +{ + + psFree (onePar); + psFree (oneDeriv); + return; +} + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x) +{ + + psF32 value; + psF32 model; + + psF32 *PAR = onePar->data.F32; + psF32 *dPAR = oneDeriv->data.F32; + + psF32 *pars = params->data.F32; + psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32; + + int nSrc = (params->n - 1) / (nPar - 1); + + PAR[0] = model = pars[0]; + for (int i = 0; i < nSrc; i++) { + int nOff = i*nPar - i; + for (int n = 1; n < nPar; n++) { + PAR[n] = pars[nOff + n]; + } + if (deriv == NULL) { + value = mFunc (NULL, onePar, x); + } else { + value = mFunc (oneDeriv, onePar, x); + for (int n = 1; n < nPar; n++) { + dpars[nOff + n] = dPAR[n]; + } + } + model += value; + } + if (deriv != NULL) { + dpars[0] = dPAR[0]*2.0; + } + return (model); +} + +/* +i: 0 1 2 +n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 +i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 +*/ + +/* +# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 +# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 +*/ + +bool pmSourceFitSet (pmSource *source, + psArray *modelSet, + const bool PSF) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, 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); + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // base values on first model + pmModel *model = modelSet->data[0]; + + // set the static variables + pmModelFitSetInit (model->type); + + int nSrc = modelSet->n; + int nPar = model->params->n - 1; // number of object parameters (excluding sky) + + // define parameter vectors for source set + psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + + // define limits for single-source model + psVector *oneDelta; + psVector *oneParMin; + psVector *oneParMax; + modelLimits (&oneDelta, &oneParMin, &oneParMax); + + psMinConstrain *constrain = psMinConstrainAlloc(); + constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL; + + // all but the sky are allowed to vary independently (subject to PSF) + // set the parameters from the multiple models + // also, set limits based on single-source limits + params->data.F32[0] = model->params->data.F32[0]; + constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0]; + constrain->paramMin->data.F32[0] = oneParMin->data.F32[0]; + constrain->paramMax->data.F32[0] = oneParMax->data.F32[0]; + + // XXX all models skip the sky + constrain->paramMask->data.U8[0] = 1; + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + for (int n = 1; n < nPar + 1; n++) { + params->data.F32[i*nPar + n] = model->params->data.F32[n]; + dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n]; + constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n]; + constrain->paramMin->data.F32[i*nPar + n] = oneParMin->data.F32[n]; + constrain->paramMax->data.F32[i*nPar + n] = oneParMax->data.F32[n]; + if (PSF) { + constrain->paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1; + } + } + } + psFree (oneDelta); + psFree (oneParMin); + psFree (oneParMax); + + if (psTraceGetLevel(__func__) >= 5) { + for (int i = 0; i < params->n; i++) { + fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]); + } + } + + // PSF model only fits first 4 parameters, EXT model fits all + int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1; + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // local sky variance + psF32 dSky = source->moments->dSky; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (PM_PSF_POISSON_WEIGHTS) { + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + } else { + yErr->data.F32[nPix] = 1.0 / dSky; + } + nPix++; + } + } + if (nPix < nParams + 1) { + psTrace (__func__, 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree (dparams); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + return(false); + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (__func__, 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet); + + // parameter errors from the covariance matrix + for (int i = 0; i < dparams->n; i++) { + if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // get the Gauss-Newton distance for fixed model parameters + if (constrain->paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet); + for (int i = 0; i < dparams->n; i++) { + if (!constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + // assign back the parameters to the models + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + model->params->data.F32[0] = params->data.F32[0]; + for (int n = 1; n < nPar + 1; n++) { + model->params->data.F32[n] = params->data.F32[i*nPar + n]; + model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; + } + // save the resulting chisq, nDOF, nIter + // these are not unique for any one source + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // set the model success or failure status + model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE; + + // models can go insane: reject these + onPic &= (model->params->data.F32[2] >= source->pixels->col0); + onPic &= (model->params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (model->params->data.F32[3] >= source->pixels->row0); + onPic &= (model->params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + } + + source->mode |= PM_SOURCE_MODE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + psFree(params); + psFree(dparams); + + // free static memory used by pmModelFitSet + pmModelFitSetClear (); + + rc = (onPic && fitStatus); + psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF); + psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc); + return(rc); +} + Index: trunk/psModules/src/objects/pmSourceFitSet.h =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.h (revision 6872) +++ trunk/psModules/src/objects/pmSourceFitSet.h (revision 6872) @@ -0,0 +1,40 @@ +/** @file pmSourceFitset.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_FIT_SET_H +# define PM_SOURCE_FIT_SET_H + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x); + +bool pmModelFitSetInit (pmModelType type); + +void pmModelFitSetClear (void); + +bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors); + +/** pmSourceFitSet() + * + * Fit the requested model to the specified source. The starting guess for the + * model is given by the input source.model parameter values. The pixels of + * interest are specified by the source.pixelsand source.maskentries. This + * function calls psMinimizeLMChi2() on the image data. The function returns TRUE + * on success or FALSE on failure. + * + */ +bool pmSourceFitSet( + pmSource *source, ///< The input pmSource + psArray *modelSet, ///< model to be fitted + const bool PSF ///< Treat model as PSF or EXT? +); + +# endif /* PM_SOURCE_FIT_MODEL_H */ Index: trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO.c (revision 6872) @@ -0,0 +1,378 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psAdditionals.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// translations between psphot object types and dophot object types +int pmSourceDophotType (pmSource *source) +{ + + switch (source->type) { + + case PM_SOURCE_TYPE_DEFECT: + case PM_SOURCE_TYPE_SATURATED: + return (8); + + case PM_SOURCE_TYPE_STAR: + if (source->mode & PM_SOURCE_MODE_SATSTAR) + return (10); + if (source->mode & PM_SOURCE_MODE_POOR) + return (7); + if (source->mode & PM_SOURCE_MODE_FAIL) + return (4); + return (1); + + case PM_SOURCE_TYPE_EXTENDED: + return (2); + + default: + return (0); + } + return (0); +} + +// Given a FITS file pointer, read the table of object data +bool pmFPAviewWriteObjects (const pmFPAview *view, pmFPAfile *file) +{ + + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAWriteObjects (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipWriteObjects (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellWriteObjects (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutWriteObjects (readout, view, file); + return true; +} + +// read in all chip-level Objects files for this FPA +bool pmFPAWriteObjects (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + pmChipWriteObjects (chip, view, file); + } + return true; +} + +// read in all cell-level Objects files for this chip +bool pmChipWriteObjects (pmChip *chip, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + pmCellWriteObjects (cell, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmCellWriteObjects (pmCell *cell, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + pmReadoutWriteObjects (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmReadoutWriteObjects (pmReadout *readout, const pmFPAview *view, pmFPAfile *file) +{ + + bool status; + char *filename; + char *dataname; + char *headname; + pmHDU *hdu; + pmHDU *phu; + psMetadata *updates; + psMetadata *outhead; + + psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); + + switch (file->type) { + case PM_FPA_FILE_RAW: + filename = pmFPAfileNameFromRule (file->filerule, file, view); + pmSourcesWriteRAW (sources, filename); + psFree (filename); + break; + + case PM_FPA_FILE_OBJ: + filename = pmFPAfileNameFromRule (file->filerule, file, view); + pmSourcesWriteOBJ (sources, filename); + psFree (filename); + break; + + case PM_FPA_FILE_SX: + filename = pmFPAfileNameFromRule (file->filerule, file, view); + pmSourcesWriteSX (sources, filename); + psFree (filename); + break; + + case PM_FPA_FILE_CMP: + // a SPLIT format : only one header and object table per file + hdu = pmFPAviewThisHDU (view, file->fpa); + filename = pmFPAfileNameFromRule (file->filerule, file, view); + + // copy the header to an output header, add the output header data + outhead = psMetadataCopy (NULL, hdu->header); + + // check/fix first line (must be SIMPLE = F) + psMetadataItem *item = psMetadataGet (outhead, PS_LIST_HEAD); + if (strcmp (item->name, "SIMPLE") && strcmp (item->name, "XTENSION")) { + psErrorStackPrint(stderr, "invalid header: first line is neither SIMPLE nor XTENSION\n"); + exit(EXIT_FAILURE); + } + if (!strcmp (item->name, "XTENSION")) { + psMetadataRemoveIndex (outhead, PS_LIST_HEAD); + } + psMetadataAddBool (outhead, PS_LIST_HEAD, "SIMPLE", PS_META_REPLACE, "CMP file, not simple", false); + + // copy over the entries saved in the + updates = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.HEADER"); + psMetadataCopy (outhead, updates); + + pmSourcesWriteCMP (sources, filename, outhead); + psFree (outhead); + psFree (filename); + break; + + case PM_FPA_FILE_CMF: + // write a PHU? + // write a header? + // write the data + + // get the current header + hdu = pmFPAviewThisHDU (view, file->fpa); + + // if file does not yet have a PHU, write it to disk + if (file->phu == NULL) { + // get the corresponding phu + phu = pmFPAviewThisPHU (view, file->fpa); + + // CMF always has extensions + outhead = psMetadataCopy (NULL, phu->header); + psMetadataAdd (outhead, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true); + psFitsWriteHeaderNotImage (file->fits, outhead); + file->phu = phu->header; + psTrace ("pmFPAfile", 5, "wrote phu %s (type: %d)\n", file->filename, file->type); + psFree (outhead); + } + + // this this header block is new, write it to disk, + if (hdu->header != file->header) { + // determine name for header extension + headname = pmFPAfileNameFromRule (file->extxtra, file, view); + outhead = psMetadataCopy (NULL, hdu->header); + psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "extension name", headname); + psFitsWriteHeaderNotImage (file->fits, outhead); + psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type); + file->header = hdu->header; + psFree (outhead); + psFree (headname); + } + + // write this table to disk + dataname = pmFPAfileNameFromRule (file->extrule, file, view); + updates = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.HEADER"); + pmSourcesWriteCMF (file->fits, sources, updates, dataname); + psTrace ("pmFPAfile", 5, "wrote ext data %s (type: %d)\n", file->filename, file->type); + psFree (dataname); + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); + break; + } + return true; +} +// a MEF CMF file has: PHU, CELL-HEAD, TABLE, CELL-HEAD, TABLE, TABLE, TABLE... + +// Given a FITS file pointer, read the table of object data +bool pmFPAviewReadObjects (const pmFPAview *view, pmFPAfile *file) +{ + pmFPA *fpa = file->fpa; + + if (view->chip == -1) { + pmFPAReadObjects (fpa, view, file); + return true; + } + + if (view->chip >= fpa->chips->n) { + return false; + } + pmChip *chip = fpa->chips->data[view->chip]; + + if (view->cell == -1) { + pmChipReadObjects (chip, view, file); + return true; + } + + if (view->cell >= chip->cells->n) { + return false; + } + pmCell *cell = chip->cells->data[view->cell]; + + if (view->readout == -1) { + pmCellReadObjects (cell, view, file); + return true; + } + + if (view->readout >= cell->readouts->n) { + return false; + } + pmReadout *readout = cell->readouts->data[view->readout]; + + pmReadoutReadObjects (readout, view, file); + return true; +} + +// read in all chip-level Objects files for this FPA +bool pmFPAReadObjects (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < fpa->chips->n; i++) { + + pmChip *chip = fpa->chips->data[i]; + pmChipReadObjects (chip, view, file); + } + return true; +} + +// read in all cell-level Objects files for this chip +bool pmChipReadObjects (pmChip *chip, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < chip->cells->n; i++) { + + pmCell *cell = chip->cells->data[i]; + pmCellReadObjects (cell, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmCellReadObjects (pmCell *cell, const pmFPAview *view, pmFPAfile *file) +{ + + for (int i = 0; i < cell->readouts->n; i++) { + + pmReadout *readout = cell->readouts->data[i]; + pmReadoutReadObjects (readout, view, file); + } + return true; +} + +// read in all readout-level Objects files for this cell +bool pmReadoutReadObjects (pmReadout *readout, const pmFPAview *view, pmFPAfile *file) +{ + + bool status; + psArray *sources = NULL; + pmHDU *hdu; + + switch (file->type) { + case PM_FPA_FILE_OBJ: + fprintf (stderr, "warning: OBJ is not supported as an input object format\n"); + break; + + case PM_FPA_FILE_SX: + fprintf (stderr, "warning: SX is not supported as an input object format\n"); + break; + + case PM_FPA_FILE_CMP: + // a SPLIT format : only one header and object table per file + + // read in header, if not yet loaded + hdu = pmFPAviewThisHDU (view, file->fpa); + + char *filename = pmFPAfileNameFromRule (file->filerule, file, view); + file->fits = psFitsOpen (filename, "r"); + hdu->header = psFitsReadHeader (NULL, file->fits); + psFitsClose (file->fits); + sources = pmSourcesReadCMP (filename, hdu->header); + psFree (filename); + break; + + case PM_FPA_FILE_CMF: + // read in header, if not yet loaded + hdu = pmFPAviewThisHDU (view, file->fpa); + if (hdu->header == NULL) { + char *headname = pmFPAfileNameFromRule (file->extxtra, file, view); + psFitsMoveExtName (file->fits, headname); + hdu->header = psFitsReadHeader (NULL, file->fits); + psFree (headname); + } + + char *dataname = pmFPAfileNameFromRule (file->extrule, file, view); + psFitsMoveExtName (file->fits, dataname); + sources = pmSourcesReadCMF (file->fits, hdu->header); + psFree (dataname); + break; + + default: + fprintf (stderr, "warning: type mismatch\n"); + break; + } + status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); + return true; +} + + Index: trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- trunk/psModules/src/objects/pmSourceIO.h (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO.h (revision 6872) @@ -0,0 +1,44 @@ +/** @file pmSourceIO.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_IO_H +# define PM_SOURCE_IO_H + +int pmSourceDophotType (pmSource *source); + +bool pmSourcesWriteRAW (psArray *sources, char *filename); +bool pmSourcesWriteOBJ (psArray *sources, char *filename); +bool pmSourcesWriteSX (psArray *sources, char *filename); +bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header); +bool pmSourcesWriteCMF (psFits *fits, psArray *sources, psMetadata *header, char *extname); + +psArray *pmSourcesReadCMP (char *filename, psMetadata *header); +psArray *pmSourcesReadCMF (psFits *fits, psMetadata *header); + +bool pmSourcesWritePSFs (psArray *sources, char *filename); +bool pmSourcesWriteEXTs (psArray *sources, char *filename); +bool pmSourcesWriteNULLs (psArray *sources, char *filename); +bool pmMomentsWriteText (psArray *sources, char *filename); +bool pmPeaksWriteText (psArray *peaks, char *filename); + +bool pmFPAviewReadObjects (const pmFPAview *view, pmFPAfile *file); +bool pmFPAReadObjects (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file); +bool pmChipReadObjects (pmChip *chip, const pmFPAview *view, pmFPAfile *file); +bool pmCellReadObjects (pmCell *cell, const pmFPAview *view, pmFPAfile *file); +bool pmReadoutReadObjects (pmReadout *readout, const pmFPAview *view, pmFPAfile *file); + +bool pmFPAviewWriteObjects (const pmFPAview *view, pmFPAfile *file); +bool pmFPAWriteObjects (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file); +bool pmChipWriteObjects (pmChip *chip, const pmFPAview *view, pmFPAfile *file); +bool pmCellWriteObjects (pmCell *cell, const pmFPAview *view, pmFPAfile *file); +bool pmReadoutWriteObjects (pmReadout *readout, const pmFPAview *view, pmFPAfile *file); + +# endif /* PM_SOURCE_IO_H */ Index: trunk/psModules/src/objects/pmSourceIO_CMF.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMF.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO_CMF.c (revision 6872) @@ -0,0 +1,161 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// elixir-style FITS table output (header + table in 1st extension) +// this format consists of a header derived from the image header +// followed by a zero-size matrix, followed by the table data +bool pmSourcesWriteCMF (psFits *fits, psArray *sources, psMetadata *header, char *extname) +{ + + psArray *table; + psMetadata *row; + int i, type; + psF32 *PAR, *dPAR; + float lsky; + bool status; + psEllipseShape shape; + psEllipseAxes axes; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + + table = psArrayAlloc (sources->n); + table->n = 0; + + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = psEllipseShapeToAxes (shape); + + row = psMetadataAlloc (); + psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]); + psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", PS_MIN (99.0, source->psfMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", PS_MIN (99.0, source->extMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", PS_MIN (99.0, source->apMag + ZERO_POINT)); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", PS_MIN (999, 1000*source->errMag)); + psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", axes.major); + psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", axes.minor); + psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", axes.theta); + psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_U8, "", type); + psMetadataAdd (row, PS_LIST_TAIL, "WEIGHT", PS_DATA_U8, "", PS_MIN (255, PS_MAX(0, 255*source->pixWeight))); + psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "12"); + + psArrayAdd (table, 100, row); + psFree (row); + } + + // create the basic table header + psMetadataAdd (header, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING | PS_META_REPLACE, "extension name", extname); + + // XXX this is broken : does not write the existing header elements except EXTNAME + psFitsWriteTable (fits, header, table); + psTrace ("pmFPAfile", 5, "wrote ext data %s\n", extname); + + psFree (table); + return true; +} + +// read in a readout from the fits file +psArray *pmSourcesReadCMF (psFits *fits, psMetadata *header) +{ + + bool status; + psF32 *PAR, *dPAR; + psEllipseShape shape; + psEllipseAxes axes; + float lsky; + + // define PSF model type + int modelType = pmModelSetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelSetType (PSF_NAME); + } + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + psArray *table = psFitsReadTable (fits); + // validate a single row of the table (must match SMP) + + psArray *sources = psArrayAlloc (table->n); + + // convert the table to the pmSource entries + // XXX need to chooose PSF vs EXT, based on type? + for (int i = 0; i < table->n; i++) { + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + psMetadata *row = table->data[i]; + + lsky = psMetadataLookupF32 (&status, row, "LOG_SKY"); + PAR[0] = pow(lsky, 10.0); + PAR[2] = psMetadataLookupF32 (&status, row, "X_PIX"); + PAR[3] = psMetadataLookupF32 (&status, row, "Y_PIX"); + axes.major = psMetadataLookupF32 (&status, row, "FWHM_X"); + axes.minor = psMetadataLookupF32 (&status, row, "FWHM_Y"); + axes.theta = psMetadataLookupF32 (&status, row, "THETA"); + + shape = psEllipseAxesToShape (axes); + + PAR[4] = shape.sx; + PAR[5] = shape.sy; + PAR[6] = shape.sxy; + + source->psfMag = psMetadataLookupF32 (&status, row, "MAG_RAW") - ZERO_POINT; + source->extMag = psMetadataLookupF32 (&status, row, "MAG_GAL") - ZERO_POINT; + source->errMag = psMetadataLookupF32 (&status, row, "MAG_ERR"); + source->apMag = psMetadataLookupF32 (&status, row, "MAG_AP") - ZERO_POINT; + + sources->data[i] = source; + } + return (sources); +} Index: trunk/psModules/src/objects/pmSourceIO_CMP.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO_CMP.c (revision 6872) @@ -0,0 +1,250 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// elixir-style pseudo FITS table (header + ascii list) +bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header) +{ + + int i, type, nSrc; + psMetadataItem *mdi; + psF32 *PAR, *dPAR; + float lsky = 0; + bool status; + psEllipseShape shape; + psEllipseAxes axes; + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + // count the number of sources which will be written + for (i = nSrc = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + nSrc ++; + } + psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", nSrc); + + // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 + mdi = psMetadataLookup (header, "NAXIS"); + mdi->data.S32 = 0; + mdi->type = PS_DATA_S32; + + // create file, write-out header + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteHeader (header, fits); + psFitsClose (fits); + + // re-open, add data to end of file + FILE *f = fopen (filename, "a+"); + if (f == NULL) { + psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); + return false; + } + fseek (f, 0, SEEK_END); + + psLine *line = psLineAlloc (67); // 66 is imclean-defined line length + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + type = pmSourceDophotType (source); + lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); + + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = psEllipseShapeToAxes (shape); + + psLineInit (line); + psLineAdd (line, "%6.1f ", PAR[2]); + psLineAdd (line, "%6.1f ", PAR[3]); + psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->psfMag + ZERO_POINT)); + psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*source->errMag))); + psLineAdd (line, "%2d ", type); + psLineAdd (line, "%3.1f ", lsky); + psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->extMag + ZERO_POINT)); + psLineAdd (line, "%6.3f ", PS_MIN (99.0, source->apMag + ZERO_POINT)); + psLineAdd (line, "%6.2f ", axes.major); + psLineAdd (line, "%6.2f ", axes.minor); + psLineAdd (line, "%5.1f\n", axes.theta); + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + psFree (line); + return true; +} + +# define BYTES_STAR 66 +# define BLOCK 1000 + +// elixir-style pseudo FITS table (header + ascii list) +psArray *pmSourcesReadCMP (char *filename, psMetadata *header) +{ + + int Ninstar; + psF32 *PAR, *dPAR; + bool status; + psEllipseShape shape; + psEllipseAxes axes; + + // define PSF model type + int modelType = pmModelSetType ("PS_MODEL_GAUSS"); + + char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME"); + if (PSF_NAME != NULL) { + modelType = pmModelSetType (PSF_NAME); + } + + // find config information for output header + float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); + if (!status) + ZERO_POINT = 25.0; + + // how many lines in the header? + int nLines = header->list->n; + int nBytes = nLines * 80; + if (nBytes % 2880) { + int nBlock = 1 + (int)(nBytes / 2880); + nBytes = nBlock * 2880; + } + + // re-open, seek to end of header + FILE *f = fopen (filename, "r"); + if (f == NULL) { + psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename); + return NULL; + } + fseek (f, nBytes, SEEK_SET); + + // prepare array to store data + int nStars = psMetadataLookupS32 (&status, header, "NSTARS"); + psArray *sources = psArrayAlloc (nStars); + sources->n = 0; + + // we have fixed bytes / line : use that info + char *buffer = psAlloc (BYTES_STAR*nStars); + + int Nextra = 0; + while (true) { + /* load next data block */ + int Nbytes = BYTES_STAR * BLOCK - Nextra; + int nbytes = fread (&buffer[Nextra], 1, Nbytes, f); + if (nbytes == 0) { + goto done_load; + } + nbytes += Nextra; + + /* check line-by-line integrity */ + char *c = buffer; + char *c2 = NULL; + while (c < buffer + nbytes) { + for (c2 = c; *c2 == '\n'; c2++) + ; + if (c2 > c) { /* extra return chars */ + memmove (c, c2, (int)(buffer + nbytes - c2)); + int Nskip = c2 - c; + nbytes -= Nskip; + bzero (buffer + nbytes, Nskip); + psLogMsg (__func__, 4, "deleted %d extra return chars\n", Nskip); + } + c2 = strchr (c, '\n'); + if (c2 == (char *) NULL) { + goto done_check; + } + c2++; + if ((c2 - c) != BYTES_STAR) { /* bad line, delete it */ + memmove (c, c2, (int)(buffer + nbytes - c2)); + int Nskip = c2 - c; + nbytes -= Nskip; + bzero (buffer + nbytes, Nskip); + psLogMsg (__func__, 4, "deleted line, %d extra chars\n", Nskip); + } else { + c = c2; + } + } +done_check: + + /* extract data for stars */ + Ninstar = nbytes / BYTES_STAR; + Nextra = nbytes % BYTES_STAR; + for (int j = 0; j < Ninstar; j++) { + psString line = psStringNCopy (&buffer[j*BYTES_STAR], BYTES_STAR); + psList *list = psStringSplit (line, " "); + psArray *array = psListToArray (list); + + pmSource *source = pmSourceAlloc (); + pmModel *model = pmModelAlloc (modelType); + source->modelPSF = model; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + PAR[0] = pow (atof (array->data[5]), 10.0); + PAR[2] = atof (array->data[0]); + PAR[3] = atof (array->data[1]); + source->psfMag = atof (array->data[2]); + source->extMag = atof (array->data[6]); + source->errMag = atof (array->data[3]) / 1000.0; + source->apMag = atof (array->data[7]); + axes.major = atof (array->data[8]); + axes.minor = atof (array->data[9]); + axes.theta = atof (array->data[10]); + + shape = psEllipseAxesToShape (axes); + + PAR[4] = shape.sx; + PAR[5] = shape.sy; + PAR[6] = shape.sxy; + + // source->type = atof (array->data[4]); + + psArrayAdd (sources, 100, source); + } + } +done_load: + + // XXX free up temp memory in this block + // XXX if sources->n != nStars, give an error? + + fclose (f); + return (sources); +} Index: trunk/psModules/src/objects/pmSourceIO_OBJ.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO_OBJ.c (revision 6872) @@ -0,0 +1,94 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// dophot-style output list with fixed line width +bool pmSourcesWriteOBJ (psArray *sources, char *filename) +{ + + int type; + psF32 *PAR, *dPAR; + float dmag, apResid; + psEllipseShape shape; + psEllipseAxes axes; + + psTimerStart ("string"); + + psLine *line = psLineAlloc (104); // 104 is dophot-defined line length + + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg (__func__, 3, "can't open output file for output %s\n", filename); + return false; + } + + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + if ((source->apMag < 99.0) && (source->psfMag < 99.0)) { + apResid = source->apMag - source->psfMag; + } else { + apResid = 0.0; + } + + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = psEllipseShapeToAxes (shape); + + psLineInit (line); + psLineAdd (line, "%3d", type); + psLineAdd (line, "%8.2f", PAR[2]); + psLineAdd (line, "%8.2f", PAR[3]); + psLineAdd (line, "%8.3f", source->psfMag); + psLineAdd (line, "%6.3f", dmag); + psLineAdd (line, "%9.2f", PAR[0]); + psLineAdd (line, "%9.3f", axes.major); + psLineAdd (line, "%9.3f", axes.minor); + psLineAdd (line, "%7.2f", axes.theta); + psLineAdd (line, "%8.3f", source->extMag); + psLineAdd (line, "%8.3f", source->apMag); + psLineAdd (line, "%8.2f\n", apResid); + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + psFree (line); + fprintf (stderr, "%f seconds for %d objects\n", psTimerMark ("string"), (int)sources->n); + return true; +} Index: trunk/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO_RAW.c (revision 6872) @@ -0,0 +1,268 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourcePhotometry.h" +#include "pmSourceIO.h" + +/***** Text Output Methods *****/ +bool pmSourcesWriteRAW (psArray *sources, char *filename) +{ + + char *name = (char *) psAlloc (strlen(filename) + 10); + + sprintf (name, "%s.psf.dat", filename); + pmSourcesWritePSFs (sources, name); + + sprintf (name, "%s.ext.dat", filename); + pmSourcesWriteEXTs (sources, name); + + sprintf (name, "%s.nul.dat", filename); + pmSourcesWriteNULLs (sources, name); + + sprintf (name, "%s.mnt.dat", filename); + pmMomentsWriteText (sources, name); + + psFree (name); + return true; +} + +// write the PSF sources to an output file +bool pmSourcesWritePSFs (psArray *sources, char *filename) +{ + + double dPos; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg (__func__, 3, "can't open output file for PSFs: %s\n", filename); + return false; + } + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + if (source->type != PM_SOURCE_TYPE_STAR) + continue; + model = source->modelPSF; + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // dPos is positional error, dMag is mag error + dPos = hypot (dPAR[2], dPAR[3]); + + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->psfMag, source->errMag, dPos); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", PAR[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dPAR[j]); + } + fprintf (f, ": %8.4f %2d %#5x %7.3f %7.3f %7.1f %7.2f %4.2f %4d %2d\n", + source[0].apMag, source[0].type, source[0].mode, + log10(model[0].chisq/model[0].nDOF), + log10(model[0].chisqNorm/model[0].nDOF), + source[0].moments->SN, + model[0].radiusTMP, + source[0].pixWeight, + model[0].nDOF, + model[0].nIter); + } + fclose (f); + return true; +} + +// dump the sources to an output file +bool pmSourcesWriteEXTs (psArray *sources, char *filename) +{ + + double dPos; + int i, j; + FILE *f; + psF32 *PAR, *dPAR; + pmModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename); + return false; + } + + // write sources with models first + for (i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + if (source->type != PM_SOURCE_TYPE_EXTENDED) + continue; + + model = source->modelEXT; + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // dPos is shape error + // XXX these are hardwired for SGAUSS + dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5])); + + fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", + PAR[2], PAR[3], PAR[0], source->extMag, source->errMag, dPos); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", PAR[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dPAR[j]); + } + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.3f %7.1f %7.2f %4.2f %4d %2d\n", + source->apMag, + source[0].type, source[0].mode, + log10(model[0].chisq/model[0].nDOF), + log10(model[0].chisqNorm/model[0].nDOF), + source[0].moments->SN, + model[0].radiusTMP, + source[0].pixWeight, + model[0].nDOF, + model[0].nIter); + } + fclose (f); + return true; +} + +// dump the sources to an output file +bool pmSourcesWriteNULLs (psArray *sources, char *filename) +{ + + int i; + FILE *f; + pmMoments *moment = NULL; + pmSource *source = NULL; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); + return false; + } + + pmMoments *empty = pmMomentsAlloc (); + + // write sources with models first + for (i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // skip these sources (in PSF or EXT) + if (source->type == PM_SOURCE_TYPE_STAR) + continue; + if (source->type == PM_SOURCE_TYPE_EXTENDED) + continue; + + if (source->moments == NULL) { + moment = empty; + } else { + moment = source->moments; + } + + fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", + source->peak->x, source->peak->y, source->peak->counts, + source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, + source->moments->Sum, source->moments->Peak, + source->moments->Sky, source->moments->SN, + source->moments->nPixels, source->type); + } + fclose (f); + psFree (empty); + return true; +} + +// write the moments to an output file +bool pmMomentsWriteText (psArray *sources, char *filename) +{ + + int i; + FILE *f; + pmSource *source = NULL; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + source = sources->data[i]; + if (source->moments == NULL) + continue; + fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", + source->peak->x, source->peak->y, source->peak->counts, + source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, + source->moments->Sum, source->moments->Peak, + source->moments->Sky, source->moments->SN, + source->moments->nPixels, source->type, source->mode); + } + fclose (f); + return true; +} + +// write the peaks to an output file +bool pmPeaksWriteText (psArray *peaks, char *filename) +{ + + int i; + FILE *f; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); + return false; + } + + for (i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + if (peak == NULL) + continue; + fprintf (f, "%5d %5d %7.1f\n", + peak->x, peak->y, peak->counts); + } + fclose (f); + return true; +} + Index: trunk/psModules/src/objects/pmSourceIO_SX.c =================================================================== --- trunk/psModules/src/objects/pmSourceIO_SX.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceIO_SX.c (revision 6872) @@ -0,0 +1,86 @@ +/** @file pmSourceIO.c + * + * @author EAM, IfA + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "psLine.h" +#include "psEllipse.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" + +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourceIO.h" + +// elixir-mode / sextractor-style output list with fixed line width +bool pmSourcesWriteSX (psArray *sources, char *filename) +{ + + psF32 *PAR, *dPAR; + psEllipseShape shape; + psEllipseAxes axes; + + psLine *line = psLineAlloc (110); // 110 is sextractor line length + + FILE *f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg (__func__, 3, "can't open output file for output %s\n", filename); + return false; + } + + // write sources with models + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + pmModel *model = pmSourceSelectModel (source); + if (model == NULL) + continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // pmSourceSextractType (source, &type, &flags); + + shape.sx = PAR[4]; + shape.sy = PAR[5]; + shape.sxy = PAR[6]; + axes = psEllipseShapeToAxes (shape); + + psLineInit (line); + psLineAdd (line, "%5.2f", 0.0); // should be type + psLineAdd (line, "%11.3f", PAR[2]); + psLineAdd (line, "%11.3f", PAR[3]); + psLineAdd (line, "%9.4f", source->psfMag); + psLineAdd (line, "%9.4f", source->errMag); + psLineAdd (line, "%13.4f", PAR[0]); + psLineAdd (line, "%9.2f", axes.major); + psLineAdd (line, "%9.2f", axes.minor); + psLineAdd (line, "%6.1f", axes.theta); + psLineAdd (line, "%9.4f", source->extMag); + psLineAdd (line, "%9.4f", source->apMag); + psLineAdd (line, "%4d\n", 0); // should be flags + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + psFree (line); + return true; +} + +// XXX need to fix the FWHM / shape stuff, +// XXX make sure we are using the correct mags, etc Index: trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.c (revision 6872) +++ trunk/psModules/src/objects/pmSourcePhotometry.c (revision 6872) @@ -0,0 +1,387 @@ +/** @file pmSourcePhotometry.c + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmGrowthCurve.h" +#include "pmModel.h" +#include "pmPSF.h" +#include "pmSource.h" +#include "pmModelGroup.h" +#include "pmSourcePhotometry.h" + +# define DO_SKY 0 + +static float AP_MIN_SN = 0.0; + +bool pmSourceMagnitudesInit (psMetadata *config) +{ + + bool status; + + float limit = psMetadataLookupF32 (&status, config, "AP_MIN_SN"); + if (status) { + AP_MIN_SN = limit; + } + return true; +} + +/** + this function is used to calculate the three defined source magnitudes: + - apMag : only if S/N > AP_MIN_SN + : is optionally corrected for curve-of-growth if: + - the source is a STAR (PSF) + - the option is selected (how??) + - psfMag : all sources with non-NULL modelPSF + : is optionally corrected for aperture residual if: + - the source is a STAR (PSF) + - the option is selected (how??) + + - extMag : all sources with non-NULL modelEXT +**/ + +// XXX EAM : masked region should be (optionally) elliptical +bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode) +{ + + int status = false; + bool isPSF; + float x, y; + float rflux; + float SN; + pmModel *model; + + source->psfMag = 99.0; + source->extMag = 99.0; + source->errMag = 99.0; + source->apMag = 99.0; + + switch (source->type) { + case PM_SOURCE_TYPE_STAR: + model = source->modelPSF; + if (model == NULL) + return false; + isPSF = true; + break; + + case PM_SOURCE_TYPE_EXTENDED: + model = source->modelEXT; + if (model == NULL) + return false; + isPSF = false; + break; + + default: + return false; + } + + if (model->dparams->data.F32[1] > 0) { + SN = model->params->data.F32[1] / model->dparams->data.F32[1]; + source->errMag = 1.0 / SN; + } else { + SN = 0.0; + source->errMag = 0.0; + } + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // measure object model photometry + status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF); + status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS + if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf) { + // convert to the equivalent 2D model? + source->psfMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); + } + + if (SN < AP_MIN_SN) + return false; + + // replace source flux + // XXX test to see if source has been subtracted? + pmModelAdd (source->pixels, source->mask, model, false, false); + + // set aperture mask circle to model radius + psImageKeepCircle (source->mask, x, y, model->radiusTMP, "OR", PM_SOURCE_MASK_MARKED); + + // measure the weight of included pixels + if (mode & PM_SOURCE_PHOT_WEIGHT) { + pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->mask); + } + + // measure object aperture photometry + status = pmSourcePhotometryAper (&source->apMag, model, source->pixels, source->mask); + + // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS + if (status && isPSF && psf) { + if (psf->growth && (mode & PM_SOURCE_PHOT_GROWTH)) { + source->apMag += pmGrowthCurveCorrect (psf->growth, model->radiusTMP); + } + if (mode & PM_SOURCE_PHOT_APCORR) { + rflux = pow (10.0, 0.4*source->psfMag); + source->apMag -= PS_SQR(model->radiusTMP)*rflux * psf->skyBias + psf->skySat / rflux; + } + } + + // unmask aperture + psImageKeepCircle (source->mask, x, y, model->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED); + + // subtract object, leave local sky + pmModelSub (source->pixels, source->mask, model, false, false); + + return status; +} + +/* +aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) +(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) + +*/ + +// return source model magnitude +bool pmSourcePhotometryModel (float *fitMag, pmModel *model) +{ + + float fitSum = 0; + *fitMag = 99.0; + + if (model == NULL) { + return false; + } + + // measure fitMag + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + if (fitSum <= 0) + return false; + if (!isfinite(fitSum)) + return false; + *fitMag = -2.5*log10(fitSum); + + return (true); +} + +// return source aperture magnitude +bool pmSourcePhotometryAper (float *apMag, pmModel *model, psImage *image, psImage *mask) +{ + float apSum = 0; + float sky = 0; + *apMag = 99.0; + + if (model == NULL) { + return false; + } + + if (DO_SKY) { + sky = model->params->data.F32[0]; + } else { + sky = 0; + } + + // measure apMag + for (int ix = 0; ix < image->numCols; ix++) { + for (int iy = 0; iy < image->numRows; iy++) { + if (mask->data.U8[iy][ix]) + continue; + apSum += image->data.F32[iy][ix] - sky; + } + } + if (apSum <= 0) + return false; + + *apMag = -2.5*log10(apSum); + return (true); +} + +// return source aperture magnitude +bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask) +{ + float modelSum = 0; + float validSum = 0; + float sky = 0; + float value; + + int Xo, Yo, dP; + int dX, DX, NX; + int dY, DY, NY; + + *pixWeight = 0.0; + + if (model == NULL) + return false; + + // we only care about the value of the object model, not the local sky + if (DO_SKY) { + sky = model->params->data.F32[0]; + } else { + sky = 0; + } + + // the model function returns the source flux at a position + pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + psVector *params = model->params; + + Xo = params->data.F32[2]; + Yo = params->data.F32[3]; + + dX = Xo - image->col0; + dP = image->numCols - dX; + DX = PS_MAX(dX, dP); + NX = image->numCols; + + dY = Yo - image->row0; + dP = image->numRows - dY; + DY = PS_MAX(dY, dP); + NY = image->numRows; + + // measure modelSum and validSum + // XXX this does not work for sources near the edge: we need to measure for + // a full square region + for (int ix = -DX; ix < DX + 1; ix++) { + int mx = ix + dX; + for (int iy = -DY; iy < DY + 1; iy++) { + int my = iy + dY; + + coord->data.F32[0] = (psF32) (ix + Xo); + coord->data.F32[1] = (psF32) (iy + Yo); + + // for the full model, add all points + value = modelFunc (NULL, params, coord) - sky; + modelSum += value; + + // include count only the unmasked pixels within the image area + if (mx < 0) + continue; + if (my < 0) + continue; + if (mx >= NX) + continue; + if (my >= NY) + continue; + if (mask->data.U8[my][mx]) + continue; + + validSum += value; + } + } + if (validSum <= 0) + return false; + + *pixWeight = validSum / modelSum; + return (true); +} + +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, wt; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Wi = Mi->weight; + + 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; + + // note that this is addressing the same image pixels, + // though only if both are source not model images + 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; + wt = Wi->data.F32[yi][xi]; + if (wt > 0) { + flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt; + } + } + } + return (flux); +} + +float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) +{ + + int Xs, Xe, Ys, Ye; + int xi, xj, yi, yj; + int xIs, xJs, yIs, yJs; + int xIe, yIe; + float flux, wt; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Wi = Mi->weight; + + 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; + + // note that this is addressing the same image pixels, + // though only if both are source not model images + 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; + wt = Wi->data.F32[yi][xi]; + if (wt > 0) { + flux += 1.0 / wt; + } + } + } + return (flux); +} + Index: trunk/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- trunk/psModules/src/objects/pmSourcePhotometry.h (revision 6872) +++ trunk/psModules/src/objects/pmSourcePhotometry.h (revision 6872) @@ -0,0 +1,53 @@ +/** @file pmSourcePhotometry.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_PHOTOMETRY_H +# define PM_SOURCE_PHOTOMETRY_H + +/** + * + * The function returns both the magnitude of the fit, defined as -2.5log(flux), + * where the flux is integrated under the model, theoretically from a radius of 0 + * to infinity. In practice, we integrate the model beyond 50sigma. The aperture magnitude is + * defined as -2.5log(flux) , where the flux is summed for all pixels which are + * not excluded by the aperture mask. The model flux is calculated by calling the + * model-specific function provided by pmModelFlux_GetFunction. + * + * XXX: must code this. + * + */ + +typedef enum { + PM_SOURCE_PHOT_NONE = 0x0000, + PM_SOURCE_PHOT_GROWTH = 0x0001, + PM_SOURCE_PHOT_APCORR = 0x0002, + PM_SOURCE_PHOT_WEIGHT = 0x0004, +} pmSourcePhotometryMode; + +bool pmSourcePhotometryModel( + float *fitMag, ///< integrated fit magnitude + pmModel *model ///< model used for photometry +); + +bool pmSourcePhotometryAper( + float *apMag, ///< aperture flux magnitude + pmModel *model, ///< model used for photometry + psImage *image, ///< image pixels to be used + psImage *mask ///< mask of pixels to ignore +); + +bool pmSourceMagnitudesInit (psMetadata *config); +bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode); +float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); +float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj); +bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask); + +# endif /* PM_SOURCE_PHOTOMETRY_H */ Index: trunk/psModules/src/objects/pmSourceSky.c =================================================================== --- trunk/psModules/src/objects/pmSourceSky.c (revision 6872) +++ trunk/psModules/src/objects/pmSourceSky.c (revision 6872) @@ -0,0 +1,120 @@ +/** @file pmSourceSky.c + * + * Functions to measure the local sky and sky variance for sources on images + * + * @author GLG, MHPCC + * @author EAM, IfA: significant modifications. + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#include +#include +#include +#include "pslib.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModel.h" +#include "pmSource.h" +#include "pmSourceSky.h" + +/****************************************************************************** +pmSource *pmSourceLocalSky(source, statsOptions, Radius): this +routine creates a new pmSource.moments element if needed and sets pmSource.pmMoments.sky + +The sky value is set from the pixels in the square annulus surrounding the +peak pixel. + +The source.pixels and source.mask must already exist +*****************************************************************************/ + +bool pmSourceLocalSky( + pmSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_IMAGE_NON_NULL(source->pixels, false); + PS_ASSERT_IMAGE_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_INT_POSITIVE(Radius, false); + PS_ASSERT_INT_NONNEGATIVE(Radius, false); + + psImage *image = source->pixels; + psImage *mask = source->mask; + pmPeak *peak = source->peak; + psRegion srcRegion; + + srcRegion = psRegionForSquare(peak->x, peak->y, Radius); + srcRegion = psRegionForImage(mask, srcRegion); + + psImageMaskRegion(mask, srcRegion, "OR", PM_SOURCE_MASK_MARKED); + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion(mask, srcRegion, "AND", ~PM_SOURCE_MASK_MARKED); + + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + psFree(myStats); + + if (isnan(tmpF64)) { + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } + source->moments->Sky = (psF32) tmpF64; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return (true); +} + +// A complementary function to pmSourceLocalSky: calculate the local median variance +bool pmSourceLocalSkyVariance( + pmSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_IMAGE_NON_NULL(source->weight, false); + PS_ASSERT_IMAGE_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_INT_POSITIVE(Radius, false); + PS_ASSERT_INT_NONNEGATIVE(Radius, false); + + psImage *image = source->weight; + psImage *mask = source->mask; + pmPeak *peak = source->peak; + psRegion srcRegion; + + srcRegion = psRegionForSquare(peak->x, peak->y, Radius); + srcRegion = psRegionForImage(mask, srcRegion); + + psImageMaskRegion(mask, srcRegion, "OR", PM_SOURCE_MASK_MARKED); + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion(mask, srcRegion, "AND", ~PM_SOURCE_MASK_MARKED); + + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + psFree(myStats); + + if (isnan(tmpF64)) { + psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); + return(false); + } + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } + source->moments->dSky = (psF32) tmpF64; + psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); + return (true); +} Index: trunk/psModules/src/objects/pmSourceSky.h =================================================================== --- trunk/psModules/src/objects/pmSourceSky.h (revision 6872) +++ trunk/psModules/src/objects/pmSourceSky.h (revision 6872) @@ -0,0 +1,42 @@ +/** @file pmSourceSky.h + * + * @author EAM, IfA; GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:01:05 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +# ifndef PM_SOURCE_SKY_H +# define PM_SOURCE_SKY_H + +/** pmSourceLocalSky() + * + * Measure the local sky in the vicinity of the given source. The Radius + * defines the square aperture in which the moments will be measured. This + * function assumes the source pixels have been defined, and that the value of + * Radius here is smaller than the value of Radius used to define the pixels. The + * annular region not contained within the radius defined here is used to measure + * the local background in the vicinity of the source. The local background + * measurement uses the specified statistic passed in via the statsOptions entry. + * This function allocates the pmMoments structure. The resulting sky is used to + * set the value of the pmMoments.sky element of the provided pmSource structure. + * + */ +bool pmSourceLocalSky( + pmSource *source, ///< The input image (float) + psStatsOptions statsOptions, ///< The statistic used in calculating the background sky + float Radius ///< The inner radius of the square annulus to exclude +); + + +// A complementary function to pmSourceLocalSky: calculate the local sky variance +bool pmSourceLocalSkyVariance( + pmSource *source, ///< The input image (float) + psStatsOptions statsOptions, ///< The statistic used in calculating the background sky + float Radius ///< The inner radius of the square annulus to exclude +); + +# endif /* PM_SOURCE_PHOTOMETRY_H */