Index: anches/eam_branches/20090715/psModules/src/objects/models/pmModel_SGAUSS.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/models/pmModel_SGAUSS.c (revision 25305) +++ (revision ) @@ -1,389 +1,0 @@ -/****************************************************************************** - * this file defines the SGAUSS source shape model (XXX need a better name!). Note that these - * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' - * statements of their own. The models use a psVector to represent the set of parameters, with - * the sequence used to specify the meaning of the parameter. The meaning of the parameters - * may thus vary depending on the specifics of the model. All models which are used a PSF - * representations share a few parameters, for which # define names are listed in pmModel.h: - - power-law with fitted slope and outer tidal radius - 1 / (1 + z^N + kz^4) - - * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future - * PM_PAR_I0 1 - central intensity - * PM_PAR_XPOS 2 - X center of object - * PM_PAR_YPOS 3 - Y center of object - * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) - * PM_PAR_SXY 6 - X*Y term of elliptical contour - * PM_PAR_7 7 - slope of power-law component (N) - * PM_PAR_8 8 - amplitude of the tidal cutoff (k) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - the code for guessing the power-law slope based on the radial profile - is either too slow or does not work well. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_SGAUSS -# define PM_MODEL_FLUX pmModelFlux_SGAUSS -# define PM_MODEL_GUESS pmModelGuess_SGAUSS -# define PM_MODEL_LIMITS pmModelLimits_SGAUSS -# define PM_MODEL_RADIUS pmModelRadius_SGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_SGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_SGAUSS - -psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image); - -psF32 PM_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 *PAR = params->data.F32; - - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y), 1e-8); - // note that if z -> 0, dPAR[7] -> -inf - // also z^(PAR[7]-1) -> Inf - - psF32 pr = z*PAR[8]; - psF32 pr3 = pr*pr*pr; - psF32 p = pow(z, PAR[7] - 1.0); - psF32 r = 1.0 / (1 + z*p + pr*pr3); - psF32 f = PAR[1]*r + PAR[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7]*p + 4*PAR[8]*pr3); - - 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] = -2*t*log(z)*z*p; - deriv->data.F32[8] = -2*t*4*z*pr3; - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (9, PS_TYPE_F32); - *params_min = psVectorAlloc (9, PS_TYPE_F32); - *params_max = psVectorAlloc (9, PS_TYPE_F32); - - beta_lim[0][0].data.F32[0] = 1000; - beta_lim[0][0].data.F32[1] = 10000; - beta_lim[0][0].data.F32[2] = 5; - beta_lim[0][0].data.F32[3] = 5; - beta_lim[0][0].data.F32[4] = 0.5; - beta_lim[0][0].data.F32[5] = 0.5; - beta_lim[0][0].data.F32[6] = 0.5; - beta_lim[0][0].data.F32[7] = 0.5; - beta_lim[0][0].data.F32[8] = 0.05; - - params_min[0][0].data.F32[0] = -1000; - params_min[0][0].data.F32[1] = 0; - params_min[0][0].data.F32[2] = -100; - params_min[0][0].data.F32[3] = -100; - params_min[0][0].data.F32[4] = 0.01; - params_min[0][0].data.F32[5] = 0.01; - params_min[0][0].data.F32[6] = -5.0; - params_min[0][0].data.F32[7] = 0.5; - params_min[0][0].data.F32[8] = 0.001; - - params_max[0][0].data.F32[0] = 1e5; - params_max[0][0].data.F32[1] = 1e6; - 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! - params_max[0][0].data.F32[4] = 2.0; - params_max[0][0].data.F32[5] = 2.0; - params_max[0][0].data.F32[6] = +3.0; - params_max[0][0].data.F32[7] = 5.0; - params_max[0][0].data.F32[8] = 0.5; - - return (TRUE); -} - -bool PM_MODEL_GUESS (pmModel *model, pmSource *source) -{ - - pmMoments *sMoments = source->moments; - pmPeak *peak = source->peak; - psF32 *params = model->params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(sMoments->Sx); - moments.y2 = PS_SQR(sMoments->Sy); - moments.xy = sMoments->Sxy; - - // solve the math to go from Moments To Shape - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params[0] = sMoments->Sky; - params[1] = sMoments->Peak - sMoments->Sky; - params[2] = peak->x; - params[3] = peak->y; - params[4] = 1.0 / shape.sx; - params[5] = 1.0 / shape.sy; - params[6] = shape.sxy; - - # if (0) - - f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - axes.major *= 2.0; - axes.minor *= 2.0; - f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - - if (f1 > f2) { - params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); - } else { - params[7] = 0.5; - } - # endif - - params[7] = 1.8; - params[8] = 0.1; - - - return(true); -} - -psF64 PM_MODEL_FLUX (const psVector *params) -{ - float f, norm, z; - - psF32 *PAR = params->data.F32; - - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); - // Area is equivalent to 2 pi sigma^2 - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - for (z = 0.005; z < 50; z += 0.01) { - psF32 pr = PAR[8]*z; - f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr))); - norm += f; - } - norm *= 0.01; - - psF64 Flux = PAR[1] * Area * norm; - - return(Flux); -} - -// XXX need to define the radius along the major axis -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - psF64 r, z = 0.0, pr, f; - psF32 *PAR = params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - - if (flux <= 0) - return (1.0); - if (PAR[1] <= 0) - return (1.0); - if (flux >= PAR[1]) - return (1.0); - - // convert Sx,Sy,Sxy to major/minor axes - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (shape); - psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; - - // XXX : we can do this faster with an intelligent starting choice - for (r = 0.0; r < 20.0; r += dr) { - z = PS_SQR(r); - pr = PAR[8]*z; - f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(pr))); - if (f < limit) - break; - } - psF64 radius = 2.0 * axes.major * sqrt (z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 9; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - // XXX: Verify this (from EAM change) - //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]); - out[i] = psPolynomial2DEval(poly, out[2], out[3]); - } - return(true); -} - -bool PM_MODEL_FIT_STATUS (pmModel *model) -{ - - psF32 dP; - bool status; - psEllipseAxes axes; - psEllipseShape shape; - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (shape); - - dP = 0; - dP += PS_SQR(dPAR[4] / PAR[4]); - dP += PS_SQR(dPAR[5] / PAR[5]); - dP += PS_SQR(dPAR[7] / PAR[7]); - dP = sqrt (dP); - - status = true; - status &= (dP < 0.5); - status &= (PAR[1] > 0); - status &= ((dPAR[1]/PAR[1]) < 0.5); - status &= (fabs(PAR[8]) < 0.5); - status &= (dPAR[8] < 0.1); - status &= (axes.major > 1.41); - status &= (axes.minor > 1.41); - status &= ((axes.major / axes.minor) < 5.0); - status &= (PAR[7] > 0.5); - - if (status) - return true; - return false; -} - -// measure the flux for the elliptical contour -psF64 psImageEllipseContour (psEllipseAxes axes, double xc, double yc, psImage *image) -{ - - double t, dt, ct, st, xo, yo, value; - int N, Nt, x, y; - - // choose dt to uniformly divide contour, with ~1 pix spacing at most - dt = asin (1 / axes.minor); - Nt = (int)(2*M_PI / dt) + 1; - dt = 2*M_PI / Nt; - - ct = cos(axes.theta); - st = sin(axes.theta); - xo = xc - image->col0; - yo = yc - image->row0; - - psVector *contour = psVectorAlloc (Nt, PS_TYPE_F32); - contour->n = contour->nalloc; - for (t = 0, N = 0; (t < 2*M_PI) && (N < Nt); t += dt) { - x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo; - y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo; - value = p_psImageGetElementF64(image, x, y); - if (isfinite(value)) { - contour->data.F32[N] = value; - N++; - } - } - contour->n = N; - // accept every pixel: double counting is not so problematic here... - - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - if (!psVectorStats (stats, contour, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - value = stats->sampleMedian; - - psFree (stats); - psFree (contour); - - return (value); -} - -// XXX EAM : test version using flux contours to guess slope -bool PM_MODEL_GUESS_HARD (pmModel *model, pmSource *source) -{ - - pmMoments *sMoments = source->moments; - pmPeak *peak = source->peak; - psF32 *params = model->params->data.F32; - float f1, f2; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(sMoments->Sx); - moments.y2 = PS_SQR(sMoments->Sy); - moments.xy = sMoments->Sxy; - - // solve the math to go from Moments To Shape - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params[0] = sMoments->Sky; - params[1] = sMoments->Peak - sMoments->Sky; - params[2] = peak->x; - params[3] = peak->y; - params[4] = 1.0 / shape.sx; - params[5] = 1.0 / shape.sy; - params[6] = shape.sxy; - - f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - axes.major *= 2.0; - axes.minor *= 2.0; - f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels); - - if (f1 > f2) { - params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0))); - } else { - params[7] = 0.5; - } - params[8] = 0.1; - - return(true); -} - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_FIT_STATUS Index: anches/eam_branches/20090715/psModules/src/objects/models/pmModel_TGAUSS.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/models/pmModel_TGAUSS.c (revision 25305) +++ (revision ) @@ -1,202 +1,0 @@ -/****************************************************************************** - * this file defines the TGAUSS source shape model (XXX need a better name!). Note that these - * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' - * statements of their own. The models use a psVector to represent the set of parameters, with - * the sequence used to specify the meaning of the parameter. The meaning of the parameters - * may thus vary depending on the specifics of the model. All models which are used a PSF - * representations share a few parameters, for which # define names are listed in pmModel.h: - - power-law with fixed slope and fitted amplitude - 1 / (1 + z + kz^2.2) - - * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future - * PM_PAR_I0 1 - central intensity - * PM_PAR_XPOS 2 - X center of object - * PM_PAR_YPOS 3 - Y center of object - * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) - * PM_PAR_SXY 6 - X*Y term of elliptical contour - * PM_PAR_7 7 - amplitude of high-order component (k) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_TGAUSS -# define PM_MODEL_FLUX pmModelFlux_TGAUSS -# define PM_MODEL_GUESS pmModelGuess_TGAUSS -# define PM_MODEL_LIMITS pmModelLimits_TGAUSS -# define PM_MODEL_RADIUS pmModelRadius_TGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_TGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_TGAUSS - -psF64 PS_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 *PAR = params->data.F32; - - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; - - psF32 p = pow(z, 1.2); - psF32 r = 1.0 / (1 + z + PAR[7]*z*p); - psF32 f = PAR[1]*r + PAR[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(1 + PAR[7]*2.2*p); - - 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*p; - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (8, PS_TYPE_F32); - *params_min = psVectorAlloc (8, PS_TYPE_F32); - *params_max = psVectorAlloc (8, PS_TYPE_F32); - - beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; - beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; - beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_7] = 0.5; - - params_min[0][0].data.F32[PM_PAR_SKY] = -1000; - params_min[0][0].data.F32[PM_PAR_I0] = 0; - params_min[0][0].data.F32[PM_PAR_XPOS] = -100; - params_min[0][0].data.F32[PM_PAR_YPOS] = -100; - params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; - params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; - params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; - params_min[0][0].data.F32[PM_PAR_7] = 0.1; - - params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; - params_max[0][0].data.F32[PM_PAR_I0] = 1e8; - params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; - params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; - params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; - params_max[0][0].data.F32[PM_PAR_7] = 10.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = 1.0/source->moments->Sx; - params->data.F32[5] = 1.0/source->moments->Sy; - // params->data.F32[6] = source->moments->Sxy; - params->data.F32[6] = 0.0; - params->data.F32[7] = 5.0; - return(true); -} - -psF64 PS_MODEL_FLUX (const psVector *params) -{ - psF64 A1 = 1 / PS_SQR(params->data.F32[4]); - psF64 A2 = 1 / PS_SQR(params->data.F32[5]); - psF64 A3 = params->data.F32[6]; - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); - // Area is equivalent to 2 pi sigma^2 - - psF64 Flux = params->data.F32[1] * Area; - - return(Flux); -} - -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - if (flux <= 0) - return (1.0); - if (params->data.F32[1] <= 0) - return (1.0); - if (flux >= params->data.F32[1]) - return (1.0); - - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 8; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } - return(true); -} - -bool PM_MODEL_FIT_STATUS (pmModel *model) -{ - - psF32 dP; - bool status; - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - dP = 0; - dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); - dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); - dP = sqrt (dP); - - status = true; - status &= (dP < 0.5); - status &= (PAR[PM_PAR_I0] > 0); - status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); - - return status; -} - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_FIT_STATUS Index: anches/eam_branches/20090715/psModules/src/objects/models/pmModel_WAUSS.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/models/pmModel_WAUSS.c (revision 25305) +++ (revision ) @@ -1,198 +1,0 @@ -/****************************************************************************** - * this file defines the WAUSS source shape model (XXX need a better name!). Note that these - * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' - * statements of their own. The models use a psVector to represent the set of parameters, with - * the sequence used to specify the meaning of the parameter. The meaning of the parameters - * may thus vary depending on the specifics of the model. All models which are used a PSF - * representations share a few parameters, for which # define names are listed in pmModel.h: - - power-law with fitted linear term - 1 / (1 + Az + Bz^2 + z^3/6) - - * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future - * PM_PAR_I0 1 - central intensity - * PM_PAR_XPOS 2 - X center of object - * PM_PAR_YPOS 3 - Y center of object - * PM_PAR_SXX 4 - X^2 term of elliptical contour (SigmaX / sqrt(2)) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) - * PM_PAR_SXY 6 - X*Y term of elliptical contour - * PM_PAR_7 7 - amplitude of the linear component (A) - * PM_PAR_8 8 - amplitude of the quadratic component (B) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_WAUSS -# define PM_MODEL_FLUX pmModelFlux_WAUSS -# define PM_MODEL_GUESS pmModelGuess_WAUSS -# define PM_MODEL_LIMITS pmModelLimits_WAUSS -# define PM_MODEL_RADIUS pmModelRadius_WAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_WAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_WAUSS - -psF64 PS_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[2]; - psF32 px = params->data.F32[4]*X; - psF32 py = params->data.F32[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; - psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); - psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ - psF32 f = params->data.F32[1]*r + params->data.F32[0]; - - if (deriv != NULL) { - // note difference from gaussian: q = params->data.F32[1]*r - psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); - deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[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] = -100.0*params->data.F32[1]*r*r*t; - deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; - // The values of 100 dampen the swing of params->data.F32[7,8] */ - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (8, PS_TYPE_F32); - *params_min = psVectorAlloc (8, PS_TYPE_F32); - *params_max = psVectorAlloc (8, PS_TYPE_F32); - - beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; - beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; - beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_7] = 0.5; - - params_min[0][0].data.F32[PM_PAR_SKY] = -1000; - params_min[0][0].data.F32[PM_PAR_I0] = 0; - params_min[0][0].data.F32[PM_PAR_XPOS] = -100; - params_min[0][0].data.F32[PM_PAR_YPOS] = -100; - params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; - params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; - params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; - params_min[0][0].data.F32[PM_PAR_7] = 0.1; - - params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; - params_max[0][0].data.F32[PM_PAR_I0] = 1e8; - params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; - params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; - params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; - params_max[0][0].data.F32[PM_PAR_7] = 10.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = sqrt(2.0) / source->moments->Sx; - params->data.F32[5] = sqrt(2.0) / source->moments->Sy; - params->data.F32[6] = source->moments->Sxy; - // XXX: What are these? - // params->data.F32[7] = B2; - // params->data.F32[8] = B3; - return(true); -} - -// this is probably wrong since it uses the gauss integral 2 pi sigma^2 -psF64 PS_MODEL_FLUX (const psVector *params) -{ - psF64 A1 = 1 / PS_SQR(params->data.F32[4]); - psF64 A2 = 1 / PS_SQR(params->data.F32[5]); - psF64 A3 = params->data.F32[6]; - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3)); - // Area is equivalent to 2 pi sigma^2 - - psF64 Flux = params->data.F32[1] * Area; - - return(Flux); -} - -// return the radius which yields the requested flux -psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - if (flux <= 0) - return (1.0); - if (params->data.F32[1] <= 0) - return (1.0); - if (flux >= params->data.F32[1] <= 0) - return (1.0); - - psF64 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); - psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); - return (radius); -} - -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 9; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } - return(true); -} - -bool PM_MODEL_FIT_STATUS (pmModel *model) -{ - - psF32 dP; - bool status; - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - dP = 0; - dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); - dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); - dP = sqrt (dP); - - status = true; - status &= (dP < 0.5); - status &= (PAR[PM_PAR_I0] > 0); - status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); - - return status; -} - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_FIT_STATUS Index: anches/eam_branches/20090715/psModules/src/objects/models/pmModel_ZGAUSS.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/models/pmModel_ZGAUSS.c (revision 25305) +++ (revision ) @@ -1,251 +1,0 @@ -/****************************************************************************** - * this file defines the ZGAUSS source shape model (XXX need a better name!). Note that these - * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include' - * statements of their own. The models use a psVector to represent the set of parameters, with - * the sequence used to specify the meaning of the parameter. The meaning of the parameters - * may thus vary depending on the specifics of the model. All models which are used a PSF - * representations share a few parameters, for which # define names are listed in pmModel.h: - - power-law with fitted slope and tidal cutoff - 1 / (1 + z^N + (Az)^4) - - * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future - * PM_PAR_I0 1 - central intensity - * PM_PAR_XPOS 2 - X center of object - * PM_PAR_YPOS 3 - Y center of object - * PM_PAR_SXX 4 - X^2 term of elliptical contour (sqrt(2) / SigmaX) - * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) - * PM_PAR_SXY 6 - X*Y term of elliptical contour - * PM_PAR_7 7 - slope of power-law component (N) - *****************************************************************************/ - -/*** - XXXX the model in this file needs to be tested more carefully. - fix up the code to follow conventions in the other model function files. -***/ - -XXX broken code - -# define PM_MODEL_FUNC pmModelFunc_ZGAUSS -# define PM_MODEL_FLUX pmModelFlux_ZGAUSS -# define PM_MODEL_GUESS pmModelGuess_ZGAUSS -# define PM_MODEL_LIMITS pmModelLimits_ZGAUSS -# define PM_MODEL_RADIUS pmModelRadius_ZGAUSS -# define PM_MODEL_FROM_PSF pmModelFromPSF_ZGAUSS -# define PM_MODEL_FIT_STATUS pmModelFitStatus_ZGAUSS - -# define PAR8 0.1 - -psF64 PS_MODEL_FUNC (psVector *deriv, - const psVector *params, - const psVector *x) -{ - psF32 *PAR = params->data.F32; - - psF32 X = x->data.F32[0] - PAR[2]; - psF32 Y = x->data.F32[1] - PAR[3]; - psF32 px = PAR[4]*X; - psF32 py = PAR[5]*Y; - psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y; - - psF32 pr = PAR8*z; - psF32 p = pow(z, PAR[7] - 1.0); - psF32 r = 1.0 / (1 + z*p + SQ(SQ(pr))); - psF32 f = PAR[1]*r + PAR[0]; - - if (deriv != NULL) { - // note difference from a pure gaussian: q = params->data.F32[1]*r - psF32 t = PAR[1]*r*r; - psF32 q = t*(PAR[7]*p + 4*PAR8*pr*pr*pr); - - 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] = -q*px*X; - deriv->data.F32[5] = -q*py*Y; - deriv->data.F32[6] = -q*X*Y; - deriv->data.F32[7] = -t*log(z)*z*p; - } - return(f); -} - -bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) -{ - - *beta_lim = psVectorAlloc (8, PS_TYPE_F32); - *params_min = psVectorAlloc (8, PS_TYPE_F32); - *params_max = psVectorAlloc (8, PS_TYPE_F32); - - beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; - beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; - beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; - beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; - beta_lim[0][0].data.F32[PM_PAR_7] = 0.5; - - params_min[0][0].data.F32[PM_PAR_SKY] = -1000; - params_min[0][0].data.F32[PM_PAR_I0] = 0; - params_min[0][0].data.F32[PM_PAR_XPOS] = -100; - params_min[0][0].data.F32[PM_PAR_YPOS] = -100; - params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; - params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; - params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; - params_min[0][0].data.F32[PM_PAR_7] = 0.1; - - params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; - params_max[0][0].data.F32[PM_PAR_I0] = 1e8; - params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! - params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; - params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; - params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; - params_max[0][0].data.F32[PM_PAR_7] = 10.0; - - return (TRUE); -} - -bool PS_MODEL_GUESS (psModel *model, psSource *source) -{ - - psVector *params = model->params; - - psEllipseAxes axes; - psEllipseShape shape; - psEllipseMoments moments; - - moments.x2 = PS_SQR(source->moments->Sx); - moments.y2 = PS_SQR(source->moments->Sy); - moments.xy = source->moments->Sxy; - - axes = psEllipseMomentsToAxes(moments); - shape = psEllipseAxesToShape(axes); - - params->data.F32[0] = source->moments->Sky; - params->data.F32[1] = source->peak->counts - source->moments->Sky; - params->data.F32[2] = source->moments->x; - params->data.F32[3] = source->moments->y; - params->data.F32[4] = 1.0 / shape.sx; - params->data.F32[5] = 1.0 / shape.sy; - params->data.F32[6] = shape.sxy; - params->data.F32[7] = 1.9; - return(true); -} - -psF64 PS_MODEL_FLUX (const psVector *params) -{ - float f, norm, z; - - psF32 *PAR = params->data.F32; - - psF64 A1 = PS_SQR(PAR[4]); - psF64 A2 = PS_SQR(PAR[5]); - psF64 A3 = PS_SQR(PAR[6]); - psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3); - // Area is equivalent to 2 pi sigma^2 - - // the area needs to be multiplied by the integral of f(z) - norm = 0.0; - psF32 pr = PAR8*z; - for (z = 0.005; z < 50; z += 0.01) { - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); - norm += f; - } - norm *= 0.01; - - psF64 Flux = PAR[1] * Area * norm; - - return(Flux); -} - -// XXX need to define the radius along the major axis -// define this function so it never returns Inf or NaN -// return the radius which yields the requested flux -psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) -{ - psF64 r, z, pr, f; - psF32 *PAR = params->data.F32; - - psEllipseAxes axes; - psEllipseShape shape; - - if (flux <= 0) - return (1.0); - if (PAR[1] <= 0) - return (1.0); - if (flux >= PAR[1]) - return (1.0); - - // convert Sx,Sy,Sxy to major/minor axes - shape.sx = 1.0 / PAR[4]; - shape.sy = 1.0 / PAR[5]; - shape.sxy = PAR[6]; - - axes = psEllipseShapeToAxes (shape); - psF64 dr = 1.0 / axes.major; - psF64 limit = flux / PAR[1]; - - // XXX : we can do this faster with an intelligent starting choice - for (r = 0.0; r < 20.0; r += dr) { - z = SQ(r); - pr = PAR8*z; - f = 1.0 / (1 + pow(z, PAR[7]) + SQ(SQ(pr))); - if (f < limit) - break; - } - psF64 radius = 2.0 * axes.major * sqrt (z); - if (isnan(radius)) { - fprintf (stderr, "error in code\n"); - } - return (radius); -} - -bool PS_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, const pmPSF *psf) -{ - - psF32 *out = modelPSF->params->data.F32; - psF32 *in = modelFLT->params->data.F32; - - out[0] = in[0]; - out[1] = in[1]; - out[2] = in[2]; - out[3] = in[3]; - - for (int i = 4; i < 8; i++) { - psPolynomial2D *poly = psf->params->data[i-4]; - out[i] = Polynomial2DEval (poly, out[2], out[3]); - } - return(true); -} - -bool PM_MODEL_FIT_STATUS (pmModel *model) -{ - - psF32 dP; - bool status; - - psF32 *PAR = model->params->data.F32; - psF32 *dPAR = model->dparams->data.F32; - - dP = 0; - dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]); - dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]); - dP = sqrt (dP); - - status = true; - status &= (dP < 0.5); - status &= (PAR[PM_PAR_I0] > 0); - status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); - - return status; -} - -# undef PM_MODEL_FUNC -# undef PM_MODEL_FLUX -# undef PM_MODEL_GUESS -# undef PM_MODEL_LIMITS -# undef PM_MODEL_RADIUS -# undef PM_MODEL_FROM_PSF -# undef PM_MODEL_FIT_STATUS Index: anches/eam_branches/20090715/psModules/src/objects/pmModelGroup.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmModelGroup.c (revision 25305) +++ (revision ) @@ -1,204 +1,0 @@ -/** @file pmModelGroup.c - * - * Functions to define and manipulate object model attributes - * - * @author GLG, MHPCC - * @author EAM, IfA - * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-12-08 02:51:14 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmSpan.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmGrowthCurve.h" -#include "pmResiduals.h" -#include "pmModel.h" -#include "pmModelGroup.h" -#include "pmErrorCodes.h" - -// XXX shouldn't these be defined for us in pslib.h ??? -double hypot(double x, double y); -double sqrt (double x); - -# include "models/pmModel_GAUSS.c" -# include "models/pmModel_PGAUSS.c" -# include "models/pmModel_QGAUSS.c" -# include "models/pmModel_RGAUSS.c" -# include "models/pmModel_SERSIC.c" - -static pmModelGroup defaultModels[] = { - {"PS_MODEL_GAUSS", 7, pmModelFunc_GAUSS, pmModelFlux_GAUSS, pmModelRadius_GAUSS, pmModelLimits_GAUSS, pmModelGuess_GAUSS, pmModelFromPSF_GAUSS, pmModelFitStatus_GAUSS}, - {"PS_MODEL_PGAUSS", 7, pmModelFunc_PGAUSS, pmModelFlux_PGAUSS, pmModelRadius_PGAUSS, pmModelLimits_PGAUSS, pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelFitStatus_PGAUSS}, - {"PS_MODEL_QGAUSS", 8, pmModelFunc_QGAUSS, pmModelFlux_QGAUSS, pmModelRadius_QGAUSS, pmModelLimits_QGAUSS, pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelFitStatus_QGAUSS}, - {"PS_MODEL_RGAUSS", 8, pmModelFunc_RGAUSS, pmModelFlux_RGAUSS, pmModelRadius_RGAUSS, pmModelLimits_RGAUSS, pmModelGuess_RGAUSS, pmModelFromPSF_RGAUSS, pmModelFitStatus_RGAUSS}, - {"PS_MODEL_SERSIC", 8, pmModelFunc_SERSIC, pmModelFlux_SERSIC, pmModelRadius_SERSIC, pmModelLimits_SERSIC, pmModelGuess_SERSIC, pmModelFromPSF_SERSIC, pmModelFitStatus_SERSIC} - }; - -static pmModelGroup *models = NULL; -static int Nmodels = 0; - -static void ModelGroupFree (pmModelGroup *modelGroup) -{ - - if (modelGroup == NULL) - return; - return; -} - -pmModelGroup *pmModelGroupAlloc (int nModels) -{ - - pmModelGroup *modelGroup = (pmModelGroup *) psAlloc (nModels * sizeof(pmModelGroup)); - psMemSetDeallocator(modelGroup, (psFreeFunc) ModelGroupFree); - return (modelGroup); -} - -bool psMemCheckModelGroup(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) ModelGroupFree); -} - -void pmModelGroupAdd (pmModelGroup *model) -{ - - if (models == NULL) { - pmModelGroupInit (); - } - - Nmodels ++; - models = (pmModelGroup *) psRealloc (models, Nmodels*sizeof(pmModelGroup)); - models[Nmodels-1] = model[0]; - return; -} - -bool pmModelGroupInit (void) -{ - - // if we do not need to init, return false; - if (models != NULL) - return false; - - int Nnew = sizeof (defaultModels) / sizeof (pmModelGroup); - - models = pmModelGroupAlloc (Nnew); - for (int i = 0; i < Nnew; i++) { - models[i] = defaultModels[i]; - } - Nmodels = Nnew; - return true; -} - -void pmModelGroupCleanup (void) -{ - psFree (models); - models = NULL; - return; -} - -pmModelFunc pmModelFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFunc); -} - -pmModelFlux pmModelFlux_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFlux); -} - -pmModelRadius pmModelRadius_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelRadius); -} - -pmModelLimits pmModelLimits_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelLimits); -} - -pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelGuessFunc); -} - -pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFitStatusFunc); -} - -pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].modelFromPSFFunc); -} - -psS32 pmModelParameterCount (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (0); - } - return (models[type].nParams); -} - -psS32 pmModelSetType (char *name) -{ - for (int i = 0; i < Nmodels; i++) { - if (!strcmp(models[i].name, name)) { - return (i); - } - } - return (-1); -} - -char *pmModelGetType (pmModelType type) -{ - if ((type < 0) || (type >= Nmodels)) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return (NULL); - } - return (models[type].name); -} Index: anches/eam_branches/20090715/psModules/src/objects/pmModelGroup.h =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmModelGroup.h (revision 25305) +++ (revision ) @@ -1,201 +1,0 @@ -/* @file pmModelGroup.h - * - * The object model function types are desined to allow for the flexible addition of new object - * models. Every object model, with parameters represented by pmModel, has an associated set of - * functions which provide necessary support operations. A set of abstract functions allow the - * programmer to select the approriate function or property for a specific named object 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. - * - * @author EAM, IfA - * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-10 01:09:20 $ - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ - -# ifndef PM_MODEL_GROUP_H -# define PM_MODEL_GROUP_H - -/// @addtogroup Objects Object Detection / Analysis Functions -/// @{ - -// This function is the model chi-square minimization function for this model. -typedef psMinimizeLMChi2Func pmModelFunc; - -// This function sets the parameter limits for this model. -typedef psMinimizeLMLimitFunc pmModelLimits; - -// 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 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 sets the model parameters based on the PSF for a given coordinate and central -// intensity -typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io); - -// This function returns the success / failure status of the given model fit -typedef bool (*pmModelFitStatusFunc)(pmModel *model); - -typedef struct -{ - char *name; - int nParams; - pmModelFunc modelFunc; - pmModelFlux modelFlux; - pmModelRadius modelRadius; - pmModelLimits modelLimits; - pmModelGuessFunc modelGuessFunc; - pmModelFromPSFFunc modelFromPSFFunc; - pmModelParamsFromPSF modelParamsFromPSF; - pmModelFitStatusFunc modelFitStatusFunc; -} -pmModelGroup; - -// allocate a pmModelGroup to hold nModels entries -pmModelGroup *pmModelGroupAlloc (int nModels); - -bool psMemCheckModelGroup(psPtr ptr); - -// initialize the internal (static) model group with the default models -bool 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( - pmModelType type ///< Add comment. -); - - -/* This function returns the user-space model names for the specified model type. - */ -char *pmModelGetType( - pmModelType type ///< Add comment. -); - - -/** - * - * This function returns the internal model type code for the user-space model names. - * - */ -pmModelType pmModelSetType( - char *name ///< Add comment. -); - -/** - * - * Each of the function types above has a corresponding function which returns - * the function given the model type: - * - */ - -/** - * - * pmModelFunc is the function used to determine the value of the model at a - * specific coordinate, and is the one used by psMinimizeLMChi2. - * - */ -pmModelFunc pmModelFunc_GetFunction (pmModelType type); - - -/** - * - * pmModelFlux returns the total integrated flux for the given input parameters. - * - */ -pmModelFlux pmModelFlux_GetFunction (pmModelType type); - - -/** - * - * pmModelRadius returns the scaling radius at which the flux of the model - * matches the specified flux. This presumes that the model is a function of an - * elliptical contour. - * - */ -pmModelRadius pmModelRadius_GetFunction (pmModelType type); - - -/** - * - * pmModelLimits sets the parameter limit vectors for the function. - * - */ -pmModelLimits pmModelLimits_GetFunction (pmModelType type); - - -/** - * - * pmModelGuessFunc generates an initial guess for the model based on the - * provided source statistics (moments and pixel values as needed). - * - */ -pmModelGuessFunc pmModelGuessFunc_GetFunction (pmModelType type); - - -/** - * - * pmModelFromPSFFunc takes as input a representation of the psf and a value - * for the model and fills in the PSF parameters of the model. The input - * primarily relies upon the centroid coordinates of the input model, though the - * normalization may potentially be used. - * - */ -pmModelFromPSFFunc pmModelFromPSFFunc_GetFunction (pmModelType type); - - -/** - * - * pmModelFitStatusFunc returns a true or false values based on the success or - * failure of a model fit. The success is determined by quantities such as the - * chisq or the signal-to-noise. - * - */ -pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type); - - -/** 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. -); - -/// @} -# endif /* PM_MODEL_GROUP_H */