Index: /trunk/psModules/src/pmObjects.c =================================================================== --- /trunk/psModules/src/pmObjects.c (revision 3509) +++ /trunk/psModules/src/pmObjects.c (revision 3510) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-24 22:36:52 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-25 23:18:35 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -740,7 +740,7 @@ float SNlim, const psRegion *valid) -XXX: Waiting for sample code from IfA. - -XXX: Most code this. +XXX: The sigX and sigY stuff in the SDRS is unclear. + +XXX: How can this function ever return FALSE? *****************************************************************************/ #define SATURATE 0.0 @@ -760,7 +760,8 @@ psSource *tmpSrc = (psSource *) source->data[i]; PS_PTR_CHECK_NULL(tmpSrc->moments, false); + tmpSrc->peak->class = 0; if (tmpSrc->moments->Peak > SATURATE) { - tmpSrc->peak->class = PS_SOURCE_SATURATED; + tmpSrc->peak->class|= PS_SOURCE_SATURATED; } else { // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. @@ -773,1038 +774,1259 @@ PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain))); if (SN < FAINT_SN_LIM) { - tmpSrc->peak->class = PS_SOURCE_FAINTSTAR; + tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; } if (SN < PSF_SN_LIM) { - tmpSrc->peak->class = PS_SOURCE_FAINTSTAR; + tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; } + // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. + psF32 clumpX = 0.0; + psF32 clumpDX = 0.0; + psF32 clumpY = 0.0; + psF32 clumpYX = 0.0; + // XXX: The SDRS is not real clear on how to calculate sigX, sigY. + psF32 sigX = 0.0 + psF32 sigY = 0.0 + if ((fabs(sigX - clumpX) < clumpDX) && + (fabs(sigY - clumpY) < clumpDY)) { + tmpSrc->peak->class|= PS_SOURCE_PSFSTAR; + } + + if ((sigX < (clumpX - clumpDX)) && + (sigY < (clumpY - clumpDY))) + tmpSrc->peak->class|= PS_SOURCE_DEFECT; } - } - - return(rc); + + if ((sigX > (clumpX + clumpDX)) && + (sigY > (clumpY + clumpDY))) + tmpSrc->peak->class|= PS_SOURCE_GALAXY; + } + + if (tmpSrc->peak->class == 0) { + tmpSrc->peak->class|= PS_SOURCE_OTHER; + } } - - - -/****************************************************************************** -pmSourceSetPixelCircle(source, image, radius) - -XXX: Why boolean output? - -XXX: Why are we checking source->moments for NULL? Should the circle be - centered on the centroid or the peak? - -XXX: The circle will have a diameter of (1+radius). This is different from - the pmSourceSetLocal() function. - *****************************************************************************/ -bool pmSourceSetPixelCircle(psSource *source, - const psImage *image, - psF32 radius) -{ - PS_IMAGE_CHECK_NULL(image, false); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); - PS_PTR_CHECK_NULL(source, false); - // PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - PS_FLOAT_COMPARE(0.0, radius, false); - - // - // We define variables for code readability. - // - psS32 radiusS32 = (psS32) radius; - psS32 SubImageCenterRow = source->peak->y; - psS32 SubImageCenterCol = source->peak->x; - psS32 SubImageStartRow = SubImageCenterRow - radiusS32; - psS32 SubImageEndRow = SubImageCenterRow + radiusS32; - psS32 SubImageStartCol = SubImageCenterCol - radiusS32; - psS32 SubImageEndCol = SubImageCenterCol + radiusS32; - - if (SubImageStartRow < 0) { - psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", - SubImageStartRow); - return(false); - } - if (SubImageEndRow+1 >= image->numRows) { - psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", - SubImageEndRow); - return(false); - } - if (SubImageStartCol < 0) { - psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", - SubImageStartCol); - return(false); - } - if (SubImageEndCol+1 >= image->numCols) { - psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", - SubImageEndCol); - return(false); - } - - // XXX: Must recycle image. - if (source->pixels != NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); - psFree(source->pixels); - } - source->pixels = psImageSubset((psImage *) image, - SubImageStartCol, - SubImageStartRow, - SubImageEndCol+1, - SubImageEndRow+1); - - // XXX: Must recycle image. - if (source->mask != NULL) { - psFree(source->mask); - } - source->mask = psImageAlloc(1 + 2 * radiusS32, 1 + 2 * radiusS32, PS_TYPE_F32); - - // - // Loop through the subimage mask, initialize mask to 0 or 1. - // - for (psS32 row = 0 ; row < source->mask->numRows; row++) { - for (psS32 col = 0 ; col < source->mask->numCols; col++) { - - if (CheckRadius2((psF32) radiusS32, - (psF32) radiusS32, - radius, - (psF32) col, - (psF32) row)) { - source->mask->data.U8[row][col] = 1; - } else { - source->mask->data.U8[row][col] = 1; - } - } - } - - /* - for (psS32 row = SubImageCenterRow - radiusS32; row <= SubImageCenterRow + radiusS32; row++) { - for (psS32 col = SubImageCenterCol - radiusS32; col <= SubImageCenterCol + radiusS32; col++) { - if (CheckRadius(source->peak, radius, (psF32) col, (psF32) row)) { - source->mask->data.U8[row-SubImageCenterRow][col-SubImageCenterCol] = 1; - } - } - } - */ - - return(true); } - -/****************************************************************************** -pmSourceModelGuess(source, image, model): This function allocates a new -psModel structure and store it in the psSource data structure specified in the -argument list. The model type is specified in the argument list. The params -array in that psModel structure are allocated, and then set to the appropriate -values. This function returns true if everything was successful. - -XXX: Many of the initial parameters are set to 0.0 since I don't know what -the appropiate initial guesses are. - *****************************************************************************/ -bool pmSourceModelGuess(psSource *source, - const psImage *image, - psModelType model) -{ - PS_PTR_CHECK_NULL(source, false); - PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - PS_IMAGE_CHECK_NULL(image, false); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); - if (source->models != NULL) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n"); - psFree(source->models); - } - source->models = pmModelAlloc(PS_MODEL_UNDEFINED); - - switch (model) { - case PS_MODEL_GAUSS: - source->models->type = PS_MODEL_GAUSS; - source->models->Nparams = 7; - source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); - source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); - for (psS32 i = 0 ; i < 7 ; i++) { - source->models->params[i] = 0.0; - source->models->dparams[i] = 0.0; - } - source->models->params[0] = source->moments->Sky; - source->models->params[1] = source->peak->counts - source->moments->Sky; - source->models->params[2] = source->moments->x; - source->models->params[3] = source->moments->y; - source->models->params[4] = sqrt(2.0) / source->moments->Sx; - source->models->params[5] = sqrt(2.0) / source->moments->Sy; - source->models->params[6] = source->moments->Sxy; - source->models->chisq = 0.0; - return(true); - case PS_MODEL_PGAUSS: - source->models->type = PS_MODEL_PGAUSS; - source->models->Nparams = 7; - source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); - source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); - for (psS32 i = 0 ; i < 7 ; i++) { - source->models->params[i] = 0.0; - source->models->dparams[i] = 0.0; - } - source->models->params[0] = source->moments->Sky; - source->models->params[1] = source->peak->counts - source->moments->Sky; - source->models->params[2] = source->moments->x; - source->models->params[3] = source->moments->y; - source->models->params[4] = sqrt(2.0) / source->moments->Sx; - source->models->params[5] = sqrt(2.0) / source->moments->Sy; - source->models->params[6] = source->moments->Sxy; - source->models->chisq = 0.0; - return(true); - case PS_MODEL_TWIST_GAUSS: - source->models->type = PS_MODEL_TWIST_GAUSS; - source->models->Nparams = 11; - source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32)); - source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32)); - for (psS32 i = 0 ; i < 11 ; i++) { - source->models->params[i] = 0.0; - source->models->dparams[i] = 0.0; - } - - source->models->params[0] = source->moments->Sky; - source->models->params[1] = source->peak->counts - source->moments->Sky; - source->models->params[2] = source->moments->x; - source->models->params[3] = source->moments->y; - // XXX: What are these? - // source->models->params[4] = SxInner; - // source->models->params[5] = SyInner; - // source->models->params[6] = SxyInner; - // source->models->params[7] = SxOuter; - // source->models->params[8] = SyOuter; - // source->models->params[9] = SxyOuter; - // source->models->params[10] = N; - - source->models->chisq = 0.0; - return(true); - case PS_MODEL_WAUSS: - - source->models->params[0] = source->moments->Sky; - source->models->params[1] = source->peak->counts - source->moments->Sky; - source->models->params[2] = source->moments->x; - source->models->params[3] = source->moments->y; - source->models->params[4] = sqrt(2.0) / source->moments->Sx; - source->models->params[5] = sqrt(2.0) / source->moments->Sy; - source->models->params[6] = source->moments->Sxy; - // XXX: What are these? - // source->models->params[7] = B2; - // source->models->params[8] = B3; - - source->models->type = PS_MODEL_WAUSS; - source->models->Nparams = 9; - source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32)); - source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32)); - for (psS32 i = 0 ; i < 9 ; i++) { - source->models->params[i] = 0.0; - source->models->dparams[i] = 0.0; - } - source->models->chisq = 0.0; - return(true); - case PS_MODEL_SERSIC: - source->models->type = PS_MODEL_SERSIC; - source->models->Nparams = 8; - source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32)); - source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32)); - for (psS32 i = 0 ; i < 8 ; i++) { - source->models->params[i] = 0.0; - source->models->dparams[i] = 0.0; - } - - source->models->params[0] = source->moments->Sky; - source->models->params[1] = source->peak->counts - source->moments->Sky; - source->models->params[2] = source->moments->x; - source->models->params[3] = source->moments->y; - source->models->params[4] = sqrt(2.0) / source->moments->Sx; - source->models->params[5] = sqrt(2.0) / source->moments->Sy; - source->models->params[6] = source->moments->Sxy; - // XXX: What are these? - //source->models->params[7] = Nexp; - - source->models->chisq = 0.0; - return(true); - case PS_MODEL_SERSIC_CORE: - source->models->type = PS_MODEL_SERSIC_CORE; - source->models->Nparams = 12; - source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32)); - source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32)); - for (psS32 i = 0 ; i < 12 ; i++) { - source->models->params[i] = 0.0; - source->models->dparams[i] = 0.0; - } - - source->models->params[0] = source->moments->Sky; - source->models->params[1] = source->peak->counts - source->moments->Sky; - source->models->params[2] = source->moments->x; - source->models->params[3] = source->moments->y; - // XXX: What are these? - //source->models->params[4] SxInner; - //source->models->params[5] SyInner; - //source->models->params[6] SxyInner; - //source->models->params[7] Zd; - //source->models->params[8] SxOuter; - //source->models->params[9] SyOuter; - //source->models->params[10] = SxyOuter; - //source->models->params[11] = Nexp; - - source->models->chisq = 0.0; - return(true); - default: - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); - return(false); - } -} - - -/****************************************************************************** - - *****************************************************************************/ -psArray *pmSourceContour(psSource *source, - const psImage *image, - psF32 level, - psS32 mode) -{ - return(NULL); -} - -psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); -psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); -psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); - -//XXX: What should these values be? -#define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100 -#define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0 -/****************************************************************************** -pmSourceFitModel(source, image): must create the appropiate arguments to the -LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. - -XXX: should there be a mask value? - *****************************************************************************/ -bool pmSourceFitModel(psSource *source, - const psImage *image) -{ - PS_PTR_CHECK_NULL(source, false); - PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - PS_PTR_CHECK_NULL(source->pixels, false); - PS_PTR_CHECK_NULL(source->models, false); - PS_IMAGE_CHECK_NULL(image, false); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); - - psBool rc; - psS32 count = 0; - for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { - for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { - if (source->mask->data.U8[i][j] == 0) { - count++; - } - } - } - psArray *x = psArrayAlloc(count); - psVector *y = psVectorAlloc(count, PS_TYPE_F32); - for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { - for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { - if (source->mask->data.U8[i][j] == 0) { - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - // XXX: Should we use the subimage offsets here, or does it not matter? - coord->data.F32[0] = (psF32) (i); - coord->data.F32[1] = (psF32) (j); - x->data[count] = (psPtr *) coord; - y->data.F32[count] = source->pixels->data.F32[i][j]; - } - } - } - - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); - - psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); - - switch (source->models->type) { - case PS_MODEL_GAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); - break; - case PS_MODEL_PGAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); - break; - case PS_MODEL_TWIST_GAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); - break; - case PS_MODEL_WAUSS: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); - break; - case PS_MODEL_SERSIC: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); - break; - case PS_MODEL_SERSIC_CORE: - rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, - NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); - rc = false; - } - - psFree(x); - psFree(y); - psFree(myMin); - psFree(params); - return(rc); -} - -bool p_pmSourceAddOrSubModel(psImage *image, - psSource *source, - bool center, - psS32 flag) -{ - PS_PTR_CHECK_NULL(source, false); - PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - PS_PTR_CHECK_NULL(source->pixels, false); - PS_PTR_CHECK_NULL(source->models, false); - PS_IMAGE_CHECK_NULL(image, false); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); - - // XXX: make static, document. - psVector *deriv = psVectorAlloc(100, PS_TYPE_F32); - psVector *params = psVectorAlloc(100, PS_TYPE_F32); - psVector *x = psVectorAlloc(2, PS_TYPE_F32); - for (psS32 i = 0 ; i < source->models->Nparams ; i++) { - params->data.F32[i] = source->models->params[i]; - } - - for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { - for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { - psF32 pixelValue; - // XXX: Should you use offsets here? - // XXX: Make sure you have col/row order correct. - x->data.F32[0] = (float) j; - x->data.F32[1] = (float) i; - switch (source->models->type) { - case PS_MODEL_GAUSS: - pixelValue = pmMinLM_Gauss2D(deriv, params, x); - break; - case PS_MODEL_PGAUSS: - pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x); - break; - case PS_MODEL_TWIST_GAUSS: - pixelValue = pmMinLM_TwistGauss2D(deriv, params, x); - break; - case PS_MODEL_WAUSS: - pixelValue = pmMinLM_Wauss2D(deriv, params, x); - break; - case PS_MODEL_SERSIC: - pixelValue = pmMinLM_Sersic(deriv, params, x); - break; - case PS_MODEL_SERSIC_CORE: - pixelValue = pmMinLM_SersicCore(deriv, params, x); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); - psFree(x); - psFree(deriv); - psFree(params); - return(false); - } - if (flag == 1) { - pixelValue = -pixelValue; - } - - // XXX: Must figure out how to calculate the image coordinates and - // how to use the boolean "center" flag. - psS32 imageRow = 0; - psS32 imageCol = 0; - image->data.F32[imageRow][imageCol]+= pixelValue; - } - } - psFree(x); - psFree(deriv); - psFree(params); - return(true); -} - - - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceAddModel(psImage *image, - psSource *source, - bool center) -{ - return(p_pmSourceAddOrSubModel(image, source, center, 0)); -} - -/****************************************************************************** - *****************************************************************************/ -bool pmSourceSubModel(psImage *image, - psSource *source, - bool center) -{ - return(p_pmSourceAddOrSubModel(image, source, center, 1)); -} - - -// XXX: Put this is psConstants.h -#define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ -if (VEC1->n != N) { \ - psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ - "psVector %s has size %d, should be %d.", \ - #VEC1, VEC1->n, N); \ - return(RVAL); \ -} - - -/****************************************************************************** -pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, -y" coordinate pair. This function computes the gaussian, specified by the -parameters in "params" at that x,y point and returns the value. The -derivatives are also caculated and returned in the "deriv" argument. - - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = sqrt(2.0) / SigmaX; - params->data.F32[5] = sqrt(2.0) / SigmaY; - params->data.F32[6] = Sxy; - -XXX: Consider getting rid of the parameter checks since this might consume -a significant fraction of this function CPU time. - *****************************************************************************/ -psF32 pmMinLM_Gauss2D(psVector *deriv, - psVector *params, - psVector *x) -{ - PS_VECTOR_CHECK_NULL(deriv, NAN); - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); - PS_VECTOR_CHECK_NULL(params, NAN); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(params, 7, NAN); - PS_VECTOR_CHECK_NULL(x, NAN); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(x, 2, NAN); - - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[3]; - 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 r = exp(-z); - psF32 f = params->data.F32[1]*r + params->data.F32[0]; - - psF32 q = params->data.F32[1]*r; - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); - deriv->data.F32[3] = q*(2*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; - - return(f); -} - -/****************************************************************************** -p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above -function in a form that is usable in the LM minimization routines. - *****************************************************************************/ -psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, - psVector *params, - psArray *x) -{ - PS_IMAGE_CHECK_NULL(deriv, NULL); - PS_IMAGE_CHECK_EMPTY(deriv, NULL); - PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); - PS_PTR_CHECK_NULL(x, NULL); - if (deriv->numRows != x->n) { - psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); - } - psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); - // XXX: use static memory here. - psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); - - for (psS32 i = 0 ; i < x->n ; i++) { - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - tmpRow->data.F32[j] = deriv->data.F32[i][j]; - } - - tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, - params, - (psVector *) x->data[i]); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - - -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = sqrt(2) / SigmaX; - params->data.F32[5] = sqrt(2) / SigmaY; - params->data.F32[6] = Sxy; - *****************************************************************************/ -psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, +return(rc); + } + + + + /****************************************************************************** + pmSourceSetPixelCircle(source, image, radius) + + XXX: Why boolean output? + + XXX: Why are we checking source->moments for NULL? Should the circle be + centered on the centroid or the peak? + + XXX: The circle will have a diameter of (1+radius). This is different from + the pmSourceSetLocal() function. + *****************************************************************************/ + bool pmSourceSetPixelCircle(psSource *source, + const psImage *image, + psF32 radius) + { + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + PS_PTR_CHECK_NULL(source, false); + // PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_FLOAT_COMPARE(0.0, radius, false); + + // + // We define variables for code readability. + // + psS32 radiusS32 = (psS32) radius; + psS32 SubImageCenterRow = source->peak->y; + psS32 SubImageCenterCol = source->peak->x; + psS32 SubImageStartRow = SubImageCenterRow - radiusS32; + psS32 SubImageEndRow = SubImageCenterRow + radiusS32; + psS32 SubImageStartCol = SubImageCenterCol - radiusS32; + psS32 SubImageEndCol = SubImageCenterCol + radiusS32; + + if (SubImageStartRow < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", + SubImageStartRow); + return(false); + } + if (SubImageEndRow+1 >= image->numRows) { + psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", + SubImageEndRow); + return(false); + } + if (SubImageStartCol < 0) { + psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", + SubImageStartCol); + return(false); + } + if (SubImageEndCol+1 >= image->numCols) { + psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", + SubImageEndCol); + return(false); + } + + // XXX: Must recycle image. + if (source->pixels != NULL) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); + psFree(source->pixels); + } + source->pixels = psImageSubset((psImage *) image, + SubImageStartCol, + SubImageStartRow, + SubImageEndCol+1, + SubImageEndRow+1); + + // XXX: Must recycle image. + if (source->mask != NULL) { + psFree(source->mask); + } + source->mask = psImageAlloc(1 + 2 * radiusS32, + 1 + 2 * radiusS32, + PS_TYPE_F32); + + // + // Loop through the subimage mask, initialize mask to 0 or 1. + // + for (psS32 row = 0 ; row < source->mask->numRows; row++) { + for (psS32 col = 0 ; col < source->mask->numCols; col++) { + + if (CheckRadius2((psF32) radiusS32, + (psF32) radiusS32, + radius, + (psF32) col, + (psF32) row)) { + source->mask->data.U8[row][col] = 1; + } else { + source->mask->data.U8[row][col] = 0; + } + } + } + return(true); + } + + + /****************************************************************************** + pmSourceModelGuess(source, image, model): This function allocates a new + psModel structure and store it in the psSource data structure specified in the + argument list. The model type is specified in the argument list. The params + array in that psModel structure are allocated, and then set to the appropriate + values. This function returns true if everything was successful. + + XXX: Many of the initial parameters are set to 0.0 since I don't know what + the appropiate initial guesses are. + *****************************************************************************/ + bool pmSourceModelGuess(psSource *source, + const psImage *image, + psModelType model) + { + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + if (source->models != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n"); + psFree(source->models); + } + source->models = pmModelAlloc(PS_MODEL_UNDEFINED); + + switch (model) { + case PS_MODEL_GAUSS: + source->models->type = PS_MODEL_GAUSS; + source->models->Nparams = 7; + source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); + for (psS32 i = 0 ; i < 7 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + source->models->chisq = 0.0; + return(true); + case PS_MODEL_PGAUSS: + source->models->type = PS_MODEL_PGAUSS; + source->models->Nparams = 7; + source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); + for (psS32 i = 0 ; i < 7 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + source->models->chisq = 0.0; + return(true); + case PS_MODEL_TWIST_GAUSS: + source->models->type = PS_MODEL_TWIST_GAUSS; + source->models->Nparams = 11; + source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32)); + for (psS32 i = 0 ; i < 11 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + // XXX: What are these? + // source->models->params[4] = SxInner; + // source->models->params[5] = SyInner; + // source->models->params[6] = SxyInner; + // source->models->params[7] = SxOuter; + // source->models->params[8] = SyOuter; + // source->models->params[9] = SxyOuter; + // source->models->params[10] = N; + + source->models->chisq = 0.0; + return(true); + case PS_MODEL_WAUSS: + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + // XXX: What are these? + // source->models->params[7] = B2; + // source->models->params[8] = B3; + + source->models->type = PS_MODEL_WAUSS; + source->models->Nparams = 9; + source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32)); + for (psS32 i = 0 ; i < 9 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + source->models->chisq = 0.0; + return(true); + case PS_MODEL_SERSIC: + source->models->type = PS_MODEL_SERSIC; + source->models->Nparams = 8; + source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32)); + for (psS32 i = 0 ; i < 8 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + source->models->params[4] = sqrt(2.0) / source->moments->Sx; + source->models->params[5] = sqrt(2.0) / source->moments->Sy; + source->models->params[6] = source->moments->Sxy; + // XXX: What are these? + //source->models->params[7] = Nexp; + + source->models->chisq = 0.0; + return(true); + case PS_MODEL_SERSIC_CORE: + source->models->type = PS_MODEL_SERSIC_CORE; + source->models->Nparams = 12; + source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32)); + source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32)); + for (psS32 i = 0 ; i < 12 ; i++) { + source->models->params[i] = 0.0; + source->models->dparams[i] = 0.0; + } + + source->models->params[0] = source->moments->Sky; + source->models->params[1] = source->peak->counts - source->moments->Sky; + source->models->params[2] = source->moments->x; + source->models->params[3] = source->moments->y; + // XXX: What are these? + //source->models->params[4] SxInner; + //source->models->params[5] SyInner; + //source->models->params[6] SxyInner; + //source->models->params[7] Zd; + //source->models->params[8] SxOuter; + //source->models->params[9] SyOuter; + //source->models->params[10] = SxyOuter; + //source->models->params[11] = Nexp; + + source->models->chisq = 0.0; + return(true); + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + return(false); + } + } + + /****************************************************************************** + evalModel(source, level, row): a private function which evaluates the + source->model function at the specified coords. + + NOTE: The coords are in source->pixel coords, not image coords. + + XXX: reverse order of row,col args? + + XXX: This should probably be a public function. + + XXX: We don't care about the model derivatives here. So, for the individual + model evaluation functions, we should probably separate the derivative + calculation. + + XXX: Get rid of deriv. + + XXX: Use static vectors for x. + + XXX: Figure out if it's (row, col) or (col, row) for the model functions. + *****************************************************************************/ + psF32 evalModel(psSource *src, + psU32 row, + psU32 col) + { + // XXX: You won't need this once you separate the function evaluation. + // from the derivative calculation. + psVector *deriv = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); + // XXX: The following step will not be necessary if the models->params + // member is a psVector. + psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); + for (psS32 i = 0 ; i < src->models->Nparams ; i++) { + params->data.F32[i] = src->models->params[i]; + } + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + x->data.F32[0] = (psF32) (col + src->pixels->col0); + x->data.F32[1] = (psF32) (row + src->pixels->row0); + psF32 tmpF; + + + switch (src->models->type) { + case PS_MODEL_GAUSS: + tmpF = pmMinLM_Gauss2D(deriv, params, x); + break; + case PS_MODEL_PGAUSS: + tmpF = pmMinLM_PsuedoGauss2D(deriv, params, x); + break; + case PS_MODEL_TWIST_GAUSS: + tmpF = pmMinLM_TwistGauss2D(deriv, params, x); + break; + case PS_MODEL_WAUSS: + tmpF = pmMinLM_Wauss2D(deriv, params, x); + break; + case PS_MODEL_SERSIC: + tmpF = pmMinLM_Sersic(deriv, params, x); + break; + case PS_MODEL_SERSIC_CORE: + tmpF = pmMinLM_SersicCore(deriv, params, x); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + return(NAN); + } + + psFree(deriv); + psFree(params); + psFree(x); + return(tmpF); + } + + /****************************************************************************** + 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: The result is returned in subImage coords, not image coords. + *****************************************************************************/ + psF32 findValue(psSource *source, + psF32 level, + psU32 row, + psU32 col, + psU32 dir) + { + // Ensure that the starting column is allowable. + if (!((0 <= col) && (col < source->pixels->numCols))) { + psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); + return(NAN); + } + + psF32 oldValue = evalModel(source, row, col); + if (oldValue == level) { + return(((psF32) col)); + } + + // + // 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; + } + col+=incr; + + while (col != lastColumn) { + psF32 newValue = evalModel(source, row, col); + if (oldValue == level) { + return(col); + } + + if ((newValue <= level) && (level <= oldValue)) { + // This is simple linear interpolation. + return( ((psF32) col) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); + } + + if ((oldValue <= level) && (level <= newValue)) { + // This is simple linear interpolation. + return( ((psF32) col) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); + } + + col+=incr; + } + + return(NAN); + } + /****************************************************************************** + XXX: Probably should remove the "image" argument. + XXX: What type should the output coordinate vectors consist of? + XXX: Why a pmArray output? + XXX: doex x,y correspond with col,row or row/col? + *****************************************************************************/ + psArray *pmSourceContour(psSource *source, + const psImage *image, + psF32 level, + pmContourType mode) + { + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_PTR_CHECK_NULL(source->pixels, false); + PS_PTR_CHECK_NULL(source->models, false); + + // + // 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); + 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); + 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); + 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); + 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; + return(tmpArray); + } + + psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); + psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); + psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); + psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); + psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); + psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); + + //XXX: What should these values be? + #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100 + #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0 + /****************************************************************************** + pmSourceFitModel(source, image): must create the appropiate arguments to the + LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. + + XXX: should there be a mask value? + XXX: Probably should remove the "image" argument. + *****************************************************************************/ + bool pmSourceFitModel(psSource *source, + const psImage *image) + { + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_PTR_CHECK_NULL(source->pixels, false); + PS_PTR_CHECK_NULL(source->models, false); + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + + psBool rc; + psS32 count = 0; + for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { + for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { + if (source->mask->data.U8[i][j] == 0) { + count++; + } + } + } + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { + for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { + if (source->mask->data.U8[i][j] == 0) { + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + // XXX: Should we use the subimage offsets here, or does it not matter? + coord->data.F32[0] = (psF32) (i); + coord->data.F32[1] = (psF32) (j); + x->data[count] = (psPtr *) coord; + y->data.F32[count] = source->pixels->data.F32[i][j]; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); + + switch (source->models->type) { + case PS_MODEL_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); + break; + case PS_MODEL_PGAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); + break; + case PS_MODEL_TWIST_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); + break; + case PS_MODEL_WAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); + break; + case PS_MODEL_SERSIC: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); + break; + case PS_MODEL_SERSIC_CORE: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, + NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + rc = false; + } + + psFree(x); + psFree(y); + psFree(myMin); + psFree(params); + return(rc); + } + + bool p_pmSourceAddOrSubModel(psImage *image, + psSource *src, + bool center, + psS32 flag) + { + PS_PTR_CHECK_NULL(src, false); + PS_PTR_CHECK_NULL(src->moments, false); + PS_PTR_CHECK_NULL(src->peak, false); + PS_PTR_CHECK_NULL(src->pixels, false); + PS_PTR_CHECK_NULL(src->models, false); + PS_IMAGE_CHECK_NULL(image, false); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); + + // XXX: make static, document. + // We use a length of 100 here since it is larger enough for all models. + #define MAX_PARAMS + + psVector *deriv = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); + psVector *params = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); + psVector *x = psVectorAlloc(2, PS_TYPE_F32); + if (src->models->Nparams > MAX_PARAMS) { + psError(PS_ERR_UNKNOWN, true, "Internal error: increase MAX_PARAMS to %d", src->models->Nparams); + } + for (psS32 i = 0 ; i < src->models->Nparams ; i++) { + params->data.F32[i] = src->models->params[i]; + } + + for (psS32 i = 0 ; i < src->pixels->numRows ; i++) { + for (psS32 j = 0 ; j < src->pixels->numCols ; j++) { + psF32 pixelValue; + // XXX: Should you use offsets here? + // XXX: Make sure you have col/row order correct. + // XXX: Should you be adding the pixels for the entire subImage, + // or a radius of pixels around it? + + x->data.F32[0] = (float) j; + x->data.F32[1] = (float) i; + switch (src->models->type) { + case PS_MODEL_GAUSS: + pixelValue = pmMinLM_Gauss2D(deriv, params, x); + break; + case PS_MODEL_PGAUSS: + pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x); + break; + case PS_MODEL_TWIST_GAUSS: + pixelValue = pmMinLM_TwistGauss2D(deriv, params, x); + break; + case PS_MODEL_WAUSS: + pixelValue = pmMinLM_Wauss2D(deriv, params, x); + break; + case PS_MODEL_SERSIC: + pixelValue = pmMinLM_Sersic(deriv, params, x); + break; + case PS_MODEL_SERSIC_CORE: + pixelValue = pmMinLM_SersicCore(deriv, params, x); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + psFree(x); + psFree(deriv); + psFree(params); + return(false); + } + if (flag == 1) { + pixelValue = -pixelValue; + } + + // XXX: Must figure out how to calculate the image coordinates and + // how to use the boolean "center" flag. + psS32 imageRow = i + src->pixels->row0; + psS32 imageCol = j + src->pixles->col0; + + image->data.F32[imageRow][imageCol]+= pixelValue; + } + } + psFree(x); + psFree(deriv); + psFree(params); + return(true); + } + + + + /****************************************************************************** + *****************************************************************************/ + bool pmSourceAddModel(psImage *image, + psSource *src, + bool center) + { + return(p_pmSourceAddOrSubModel(image, src, center, 0)); + } + + /****************************************************************************** + *****************************************************************************/ + bool pmSourceSubModel(psImage *image, + psSource *src, + bool center) + { + return(p_pmSourceAddOrSubModel(image, src, center, 1)); + } + + + // XXX: Put this is psConstants.h + #define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ + if (VEC1->n != N) { \ + psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ + "psVector %s has size %d, should be %d.", \ + #VEC1, VEC1->n, N); \ + return(RVAL); \ + } + + + /****************************************************************************** + pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, + y" coordinate pair. This function computes the gaussian, specified by the + parameters in "params" at that x,y point and returns the value. The + derivatives are also caculated and returned in the "deriv" argument. + + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2.0) / SigmaX; + params->data.F32[5] = sqrt(2.0) / SigmaY; + params->data.F32[6] = Sxy; + + XXX: Consider getting rid of the parameter checks since this might consume + a significant fraction of this function CPU time. + *****************************************************************************/ + psF32 pmMinLM_Gauss2D(psVector *deriv, psVector *params, psVector *x) -{ - PS_VECTOR_CHECK_NULL(deriv, NAN); - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); - PS_VECTOR_CHECK_NULL(params, NAN); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(params, 7, NAN); - PS_VECTOR_CHECK_NULL(x, NAN); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(x, 2, NAN); - - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[3]; - 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 = 1 + z + 0.5*z*z; - psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ - psF32 f = params->data.F32[1]*r + params->data.F32[0]; - - // - // note difference from a pure gaussian: q = params->data.F32[1]*r - // - - psF32 q = params->data.F32[1]*r*r*t; - 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; - - return(f); -} - -/****************************************************************************** -p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the -above function in a form that is usable in the LM minimization routines. - *****************************************************************************/ -psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, + { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 7, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + 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 r = exp(-z); + psF32 f = params->data.F32[1]*r + params->data.F32[0]; + + psF32 q = params->data.F32[1]*r; + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); + deriv->data.F32[3] = q*(2*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; + + return(f); + } + + /****************************************************************************** + p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above + function in a form that is usable in the LM minimization routines. + *****************************************************************************/ + psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x) -{ - PS_IMAGE_CHECK_NULL(deriv, NULL); - PS_IMAGE_CHECK_EMPTY(deriv, NULL); - PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); - PS_PTR_CHECK_NULL(x, NULL); - if (deriv->numRows != x->n) { - psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); - } - psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); - // XXX: use static memory here. - psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); - - for (psS32 i = 0 ; i < x->n ; i++) { - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - tmpRow->data.F32[j] = deriv->data.F32[i][j]; - } - - tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow, - params, - (psVector *) x->data[i]); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - - - - -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = Sx; - params->data.F32[5] = Sy; - params->data.F32[6] = Sxy; - params->data.F32[7] = B2; - params->data.F32[8] = B3; - *****************************************************************************/ -psF32 pmMinLM_Wauss2D(psVector *deriv, - psVector *params, - psVector *x) -{ - PS_VECTOR_CHECK_NULL(deriv, NAN); - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); - PS_VECTOR_CHECK_NULL(params, NAN); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(params, 9, NAN); - PS_VECTOR_CHECK_NULL(x, NAN); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(x, 2, NAN); - - 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]; - - // - // 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); -} - -/****************************************************************************** -p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above -function in a form that is usable in the LM minimization routines. - *****************************************************************************/ -psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, - psVector *params, - psArray *x) -{ - PS_IMAGE_CHECK_NULL(deriv, NULL); - PS_IMAGE_CHECK_EMPTY(deriv, NULL); - PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); - PS_PTR_CHECK_NULL(x, NULL); - if (deriv->numRows != x->n) { - psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); - } - psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); - // XXX: use static memory here. - psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); - - for (psS32 i = 0 ; i < x->n ; i++) { - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - tmpRow->data.F32[j] = deriv->data.F32[i][j]; - } - - tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow, - params, - (psVector *) x->data[i]); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - - - - - - -// XXX: What should these be? -#define FFACTOR 1.0 -#define FSCALE 1.0 -/****************************************************************************** - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = SxInner; - params->data.F32[5] = SyInner; - params->data.F32[6] = SxyInner; - params->data.F32[7] = SxOuter; - params->data.F32[8] = SyOuter; - params->data.F32[9] = SxyOuter; - params->data.F32[10] = N; - *****************************************************************************/ -psF32 pmMinLM_TwistGauss2D(psVector *deriv, + { + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); + } + + + /****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = sqrt(2) / SigmaX; + params->data.F32[5] = sqrt(2) / SigmaY; + params->data.F32[6] = Sxy; + *****************************************************************************/ + psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, + psVector *params, + psVector *x) + { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 7, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + 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 = 1 + z + 0.5*z*z; + psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ + psF32 f = params->data.F32[1]*r + params->data.F32[0]; + + // + // note difference from a pure gaussian: q = params->data.F32[1]*r + // + + psF32 q = params->data.F32[1]*r*r*t; + 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; + + return(f); + } + + /****************************************************************************** + p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the + above function in a form that is usable in the LM minimization routines. + *****************************************************************************/ + psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) + { + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); + } + + + + + /****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = Sx; + params->data.F32[5] = Sy; + params->data.F32[6] = Sxy; + params->data.F32[7] = B2; + params->data.F32[8] = B3; + *****************************************************************************/ + psF32 pmMinLM_Wauss2D(psVector *deriv, + psVector *params, + psVector *x) + { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 9, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + 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]; + + // + // 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); + } + + /****************************************************************************** + p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above + function in a form that is usable in the LM minimization routines. + *****************************************************************************/ + psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) + { + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); + } + + + + + + + // XXX: What should these be? + #define FFACTOR 1.0 + #define FSCALE 1.0 + /****************************************************************************** + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = SxInner; + params->data.F32[5] = SyInner; + params->data.F32[6] = SxyInner; + params->data.F32[7] = SxOuter; + params->data.F32[8] = SyOuter; + params->data.F32[9] = SxyOuter; + params->data.F32[10] = N; + *****************************************************************************/ + psF32 pmMinLM_TwistGauss2D(psVector *deriv, + psVector *params, + psVector *x) + { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 11, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psF32 X = x->data.F32[0] - params->data.F32[2]; + psF32 Y = x->data.F32[1] - params->data.F32[3]; + psF32 px1 = params->data.F32[4]*X; + psF32 py1 = params->data.F32[5]*Y; + psF32 px2 = params->data.F32[7]*X; + psF32 py2 = params->data.F32[8]*Y; + psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; + psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; + psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); + + + psF32 f = params->data.F32[5]*r + params->data.F32[6]; + psF32 q1 = params->data.F32[5]*PS_SQR(r); + psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); + deriv->data.F32[0] = +1.0; + deriv->data.F32[1] = +r; + deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); + deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); + + // + // These fudge factors impede the growth of params->data.F32[4] beyond + // params->data.F32[7]. + // + psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); + psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; + deriv->data.F32[4] = -2.0*q1*px1*X*f2; + + // + // These fudge factors impede the growth of params->data.F32[5] beyond + // params->data.F32[8]. + // + f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); + f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; + deriv->data.F32[5] = -2.0*q1*py1*Y*f2; + deriv->data.F32[6] = -q1*X*Y; + deriv->data.F32[7] = -2.0*q2*px2*X; + deriv->data.F32[8] = -2.0*q2*py2*Y; + deriv->data.F32[9] = -q2*X*Y; + deriv->data.F32[10] = -q1*log(z2); + + return(f); + } + + /****************************************************************************** + p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above + function in a form that is usable in the LM minimization routines. + *****************************************************************************/ + psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, + psVector *params, + psArray *x) + { + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); + } + + + + /****************************************************************************** + float Sersic() + params->data.F32[0] = So; + params->data.F32[1] = Zo; + params->data.F32[2] = Xo; + params->data.F32[3] = Yo; + params->data.F32[4] = Sx; + params->data.F32[5] = Sy; + params->data.F32[6] = Sxy; + params->data.F32[7] = Nexp; + *****************************************************************************/ + psF32 pmMinLM_Sersic(psVector *deriv, psVector *params, psVector *x) -{ - PS_VECTOR_CHECK_NULL(deriv, NAN); - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); - PS_VECTOR_CHECK_NULL(params, NAN); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(params, 11, NAN); - PS_VECTOR_CHECK_NULL(x, NAN); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(x, 2, NAN); - - psF32 X = x->data.F32[0] - params->data.F32[2]; - psF32 Y = x->data.F32[1] - params->data.F32[3]; - psF32 px1 = params->data.F32[4]*X; - psF32 py1 = params->data.F32[5]*Y; - psF32 px2 = params->data.F32[7]*X; - psF32 py2 = params->data.F32[8]*Y; - psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; - psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; - psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); - - - psF32 f = params->data.F32[5]*r + params->data.F32[6]; - psF32 q1 = params->data.F32[5]*PS_SQR(r); - psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); - deriv->data.F32[0] = +1.0; - deriv->data.F32[1] = +r; - deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); - deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); - - // - // These fudge factors impede the growth of params->data.F32[4] beyond - // params->data.F32[7]. - // - psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); - psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; - deriv->data.F32[4] = -2.0*q1*px1*X*f2; - - // - // These fudge factors impede the growth of params->data.F32[5] beyond - // params->data.F32[8]. - // - f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); - f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; - deriv->data.F32[5] = -2.0*q1*py1*Y*f2; - deriv->data.F32[6] = -q1*X*Y; - deriv->data.F32[7] = -2.0*q2*px2*X; - deriv->data.F32[8] = -2.0*q2*py2*Y; - deriv->data.F32[9] = -q2*X*Y; - deriv->data.F32[10] = -q1*log(z2); - - return(f); -} - -/****************************************************************************** -p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above -function in a form that is usable in the LM minimization routines. - *****************************************************************************/ -psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, + { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 8, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); + return(0.0); + } + /****************************************************************************** + p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above + function in a form that is usable in the LM minimization routines. + *****************************************************************************/ + psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x) -{ - PS_IMAGE_CHECK_NULL(deriv, NULL); - PS_IMAGE_CHECK_EMPTY(deriv, NULL); - PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); - PS_PTR_CHECK_NULL(x, NULL); - if (deriv->numRows != x->n) { - psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); - } - psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); - // XXX: use static memory here. - psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); - - for (psS32 i = 0 ; i < x->n ; i++) { - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - tmpRow->data.F32[j] = deriv->data.F32[i][j]; - } - - tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow, - params, - (psVector *) x->data[i]); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - - - -/****************************************************************************** - float Sersic() - params->data.F32[0] = So; - params->data.F32[1] = Zo; - params->data.F32[2] = Xo; - params->data.F32[3] = Yo; - params->data.F32[4] = Sx; - params->data.F32[5] = Sy; - params->data.F32[6] = Sxy; - params->data.F32[7] = Nexp; - *****************************************************************************/ -psF32 pmMinLM_Sersic(psVector *deriv, - psVector *params, - psVector *x) -{ - PS_VECTOR_CHECK_NULL(deriv, NAN); - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); - PS_VECTOR_CHECK_NULL(params, NAN); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(params, 8, NAN); - PS_VECTOR_CHECK_NULL(x, NAN); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(x, 2, NAN); - - psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); - return(0.0); -} -/****************************************************************************** -p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above -function in a form that is usable in the LM minimization routines. - *****************************************************************************/ -psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, + { + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); + } + + /****************************************************************************** + float SersicBulge() + params->data.F32[0] So; + params->data.F32[1] Zo; + params->data.F32[2] Xo; + params->data.F32[3] Yo; + params->data.F32[4] SxInner; + params->data.F32[5] SyInner; + params->data.F32[6] SxyInner; + params->data.F32[7] Zd; + params->data.F32[8] SxOuter; + params->data.F32[9] SyOuter; + params->data.F32[10] = SxyOuter; + params->data.F32[11] = Nexp; + *****************************************************************************/ + psF32 pmMinLM_SersicCore(psVector *deriv, psVector *params, - psArray *x) -{ - PS_IMAGE_CHECK_NULL(deriv, NULL); - PS_IMAGE_CHECK_EMPTY(deriv, NULL); - PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); - PS_PTR_CHECK_NULL(x, NULL); - if (deriv->numRows != x->n) { - psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); - } - psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); - // XXX: use static memory here. - psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); - - for (psS32 i = 0 ; i < x->n ; i++) { - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - tmpRow->data.F32[j] = deriv->data.F32[i][j]; - } - - tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow, - params, - (psVector *) x->data[i]); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - -/****************************************************************************** - float SersicBulge() - params->data.F32[0] So; - params->data.F32[1] Zo; - params->data.F32[2] Xo; - params->data.F32[3] Yo; - params->data.F32[4] SxInner; - params->data.F32[5] SyInner; - params->data.F32[6] SxyInner; - params->data.F32[7] Zd; - params->data.F32[8] SxOuter; - params->data.F32[9] SyOuter; - params->data.F32[10] = SxyOuter; - params->data.F32[11] = Nexp; - *****************************************************************************/ -psF32 pmMinLM_SersicCore(psVector *deriv, - psVector *params, - psVector *x) -{ - PS_VECTOR_CHECK_NULL(deriv, NAN); - PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); - PS_VECTOR_CHECK_NULL(params, NAN); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(params, 12, NAN); - PS_VECTOR_CHECK_NULL(x, NAN); - PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); - PS_VECTOR_CHECK_SIZE(x, 2, NAN); - - psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); - return(0.0); -} -/****************************************************************************** -p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above -function in a form that is usable in the LM minimization routines. - *****************************************************************************/ -psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, - psVector *params, - psArray *x) -{ - PS_IMAGE_CHECK_NULL(deriv, NULL); - PS_IMAGE_CHECK_EMPTY(deriv, NULL); - PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); - PS_PTR_CHECK_NULL(x, NULL); - if (deriv->numRows != x->n) { - psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); - } - psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); - // XXX: use static memory here. - psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); - - for (psS32 i = 0 ; i < x->n ; i++) { - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - tmpRow->data.F32[j] = deriv->data.F32[i][j]; - } - - tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow, - params, - (psVector *) x->data[i]); - - for (psS32 j = 0 ; j < tmpRow->n ; j++) { - deriv->data.F32[i][j] = tmpRow->data.F32[j]; - } - } - - psFree(tmpRow); - return(tmpVec); -} - - - -/****************************************************************************** - *****************************************************************************/ -psF32 pmMinLM_PsuedoSersic(psVector *deriv, - psVector *params, - psVector *x) -{ - return(0.0); -} + psVector *x) + { + PS_VECTOR_CHECK_NULL(deriv, NAN); + PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(params, 12, NAN); + PS_VECTOR_CHECK_NULL(x, NAN); + PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_SIZE(x, 2, NAN); + + psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); + return(0.0); + } + /****************************************************************************** + p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above + function in a form that is usable in the LM minimization routines. + *****************************************************************************/ + psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, + psVector *params, + psArray *x) + { + PS_IMAGE_CHECK_NULL(deriv, NULL); + PS_IMAGE_CHECK_EMPTY(deriv, NULL); + PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); + PS_PTR_CHECK_NULL(x, NULL); + if (deriv->numRows != x->n) { + psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); + } + psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); + // XXX: use static memory here. + psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); + + for (psS32 i = 0 ; i < x->n ; i++) { + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + tmpRow->data.F32[j] = deriv->data.F32[i][j]; + } + + tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow, + params, + (psVector *) x->data[i]); + + for (psS32 j = 0 ; j < tmpRow->n ; j++) { + deriv->data.F32[i][j] = tmpRow->data.F32[j]; + } + } + + psFree(tmpRow); + return(tmpVec); + } + + + + /****************************************************************************** + *****************************************************************************/ + psF32 pmMinLM_PsuedoSersic(psVector *deriv, + psVector *params, + psVector *x) + { + return(0.0); + } Index: /trunk/psModules/src/pmObjects.h =================================================================== --- /trunk/psModules/src/pmObjects.h (revision 3509) +++ /trunk/psModules/src/pmObjects.h (revision 3510) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-24 22:36:52 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-25 23:18:35 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -65,4 +65,5 @@ } psModelType; +// XXX: It will be better if params, and dparams, were psVectors. typedef struct { @@ -173,8 +174,12 @@ /****************************************************************************** *****************************************************************************/ +typedef enum { + PS_CONTOUR_CRUDE, +} pmContourType; + psArray *pmSourceContour(psSource *source, const psImage *image, psF32 level, - psS32 mode); + pmContourType mode); /******************************************************************************