Index: trunk/psModules/src/pmObjects.c =================================================================== --- trunk/psModules/src/pmObjects.c (revision 3510) +++ trunk/psModules/src/pmObjects.c (revision 3511) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-25 23:18:35 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-25 23:22:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -762,4 +762,12 @@ tmpSrc->peak->class = 0; + psF32 sigX = 0.0; + psF32 sigY = 0.0; + // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. + psF32 clumpX = 0.0; + psF32 clumpDX = 0.0; + psF32 clumpY = 0.0; + psF32 clumpDY = 0.0; + if (tmpSrc->moments->Peak > SATURATE) { tmpSrc->peak->class|= PS_SOURCE_SATURATED; @@ -779,1254 +787,1246 @@ 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; + 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; } if ((sigX > (clumpX + clumpDX)) && - (sigY > (clumpY + clumpDY))) + (sigY > (clumpY + clumpDY))) { tmpSrc->peak->class|= PS_SOURCE_GALAXY; - } - - if (tmpSrc->peak->class == 0) { - tmpSrc->peak->class|= PS_SOURCE_OTHER; - } -} -} - -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, + } + + if (tmpSrc->peak->class == 0) { + tmpSrc->peak->class|= PS_SOURCE_OTHER; + } + } + + 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 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, 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, +{ + 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_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, +{ + 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, 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_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_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, +{ + 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, 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); - } + 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); +}