Index: trunk/psModules/src/pmObjects.c =================================================================== --- trunk/psModules/src/pmObjects.c (revision 3514) +++ trunk/psModules/src/pmObjects.c (revision 3576) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-03-26 03:44:31 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-03-31 00:22:49 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -17,7 +17,8 @@ #include "psConstants.h" #include "pmObjects.h" -/****************************************************************************** - *****************************************************************************/ - + +/****************************************************************************** +pmPeakAlloc(): Allocate the psPeak data structure and set appropriate members. + *****************************************************************************/ psPeak *pmPeakAlloc(psS32 x, psS32 y, @@ -34,4 +35,8 @@ } +/****************************************************************************** +pmMomentsAlloc(): Allocate the psMoments structure and initialize the members +to zero. + *****************************************************************************/ psMoments *pmMomentsAlloc() { @@ -56,4 +61,8 @@ } +/****************************************************************************** +pmModelAlloc(): Allocate the psModel structure, along with its parameters, +and initialize the type member. Initialize the params to 0.0. + *****************************************************************************/ psModel *pmModelAlloc(psModelType type) { @@ -61,8 +70,45 @@ tmp->type = type; - tmp->Nparams = 0; - tmp->params = NULL; - tmp->dparams = NULL; tmp->chisq = 0.0; + switch (type) { + case PS_MODEL_GAUSS: + tmp->Nparams = 7; + tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); + tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); + break; + case PS_MODEL_PGAUSS: + tmp->Nparams = 7; + tmp->params = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); + tmp->dparams = (psF32 *) psAlloc(7 * sizeof(PS_TYPE_F32)); + break; + case PS_MODEL_TWIST_GAUSS: + tmp->Nparams = 11; + tmp->params = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32)); + tmp->dparams = (psF32 *) psAlloc(11 * sizeof(PS_TYPE_F32)); + break; + case PS_MODEL_WAUSS: + tmp->Nparams = 9; + tmp->params = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32)); + tmp->dparams = (psF32 *) psAlloc(9 * sizeof(PS_TYPE_F32)); + break; + case PS_MODEL_SERSIC: + tmp->Nparams = 8; + tmp->params = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32)); + tmp->dparams = (psF32 *) psAlloc(8 * sizeof(PS_TYPE_F32)); + break; + case PS_MODEL_SERSIC_CORE: + tmp->Nparams = 12; + tmp->params = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32)); + tmp->dparams = (psF32 *) psAlloc(12 * sizeof(PS_TYPE_F32)); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); + return(NULL); + } + + for (psS32 i = 0 ; i < tmp->Nparams ; i++) { + tmp->params[i] = 0.0; + tmp->dparams[i] = 0.0; + } p_psMemSetDeallocator(tmp, (psFreeFcn) p_psModelFree); @@ -70,7 +116,8 @@ } -// XXX: We don't free pixels and mask since that caused a memory error. -// We might need to increase the reference counter and decrease it here. -// +/****************************************************************************** +XXX: We don't free pixels and mask since that caused a memory error. +We might need to increase the reference counter and decrease it here. + *****************************************************************************/ static void p_psSourceFree(psSource *tmp) { @@ -82,4 +129,8 @@ } +/****************************************************************************** +pmSourceAlloc(): Allocate the psSource structure and initialize its members +to NULL. + *****************************************************************************/ psSource *pmSourceAlloc() { @@ -195,4 +246,7 @@ /****************************************************************************** +p_psGetRowVectorFromImage(): a private function which simply returns a +psVector containing the specified row of data from the psImage. + XXX: Is there a better way to do this? *****************************************************************************/ @@ -210,4 +264,8 @@ } +/****************************************************************************** +MyListAddPeak(): A private function which allocates a psList, if the list +argument is NULL, otherwise it adds the peak to that list. + *****************************************************************************/ // XXX: Switch row, col args? psList *MyListAddPeak(psList *list, @@ -504,4 +562,7 @@ // We define variables for code readability. // + // XXX: Since the peak->xy coords are in image, not subImage coords, + // these variables should be renamed for clarity (imageCenterRow, etc). + // psS32 SubImageCenterRow = peak->y; psS32 SubImageCenterCol = peak->x; @@ -612,7 +673,5 @@ XXX: macro this for performance. - -XXX: should arguments be (y, x) order? - *****************************************************************************/ + *****************************************************************************/ bool CheckRadius(psPeak *peak, psF32 radius, @@ -627,4 +686,10 @@ } +/****************************************************************************** +bool CheckRadius2(): private function which simply determines if the (x, y) +point is within the radius of the specified peak. + +XXX: macro this for performance. + *****************************************************************************/ bool CheckRadius2(psF32 xCenter, psF32 yCenter, @@ -641,5 +706,7 @@ /****************************************************************************** -pmSourceMoments(source, radius) +pmSourceMoments(source, radius): this function takes a subImage defined in the +psSource data structure, along with the peak location, and determines the +various moments associated with that peak. Requires the following to have been created: @@ -647,4 +714,6 @@ psSource->peak psSource->pixels + +XXX: The peak calculations are done in image coords, not subImage coords. XXX: mask values? @@ -733,11 +802,10 @@ /****************************************************************************** -pmSourceRoughClass(source, saturate, SNlim, valid): make a guess at the source +pmSourceRoughClass(source, metadata): make a guess at the source classification. -psSource *pmSourceRoughClass(psSource *source, - psF32 saturate, - float SNlim, - const psRegion *valid) +XXX: This is not useable code, as of the release date. There remains a fair +bit of coding to be completed. + XXX: The sigX and sigY stuff in the SDRS is unclear. @@ -838,4 +906,7 @@ // We define variables for code readability. // + // XXX: Since the peak->xy coords are in image, not subImage coords, + // these variables should be renamed for clarity (imageCenterRow, etc). + // psS32 radiusS32 = (psS32) radius; psS32 SubImageCenterRow = source->peak->y; @@ -910,11 +981,17 @@ /****************************************************************************** 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. +psModel structure and stores 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. + +XXX: The image argument is redundant. + +XXX: Many parameters are based on the src->moments structure, which is in +image, not subImage coords. Therefore, the calls to the model evaluation +functions will be in image, not subImage coords. Remember this. *****************************************************************************/ bool pmSourceModelGuess(psSource *source, @@ -931,16 +1008,8 @@ psFree(source->models); } - source->models = pmModelAlloc(PS_MODEL_UNDEFINED); + source->models = pmModelAlloc(model); 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; @@ -950,15 +1019,7 @@ 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; @@ -968,16 +1029,7 @@ 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; @@ -992,9 +1044,7 @@ // 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; @@ -1007,25 +1057,7 @@ // 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; @@ -1037,17 +1069,7 @@ // 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; @@ -1063,7 +1085,6 @@ //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"); @@ -1074,15 +1095,22 @@ /****************************************************************************** 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. +source->model function at the specified coords. The coords are subImage, not +image coords. + +NOTE: The coords are in subImage source->pixel coords, not image coords. XXX: reverse order of row,col args? -XXX: This should probably be a public function. +XXX: rename all coords in this file such that their name defines whether +the coords is in subImage or image space. + +XXX: This should probably be a public psModules function. XXX: Use static vectors for x. XXX: Figure out if it's (row, col) or (col, row) for the model functions. + +XXX: For a while, the first psVectorAlloc() was generating a seg fault during +testing. Try to reproduce that and debug. *****************************************************************************/ psF32 evalModel(psSource *src, @@ -1090,15 +1118,22 @@ psU32 col) { + PS_PTR_CHECK_NULL(src, false); + PS_PTR_CHECK_NULL(src->models, false); + PS_PTR_CHECK_NULL(src->models->params, false); + // 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); + // member is a psVector. Suggest to IfA. + psVector *params = psVectorAlloc(7, PS_TYPE_F32); for (psS32 i = 0 ; i < src->models->Nparams ; i++) { params->data.F32[i] = src->models->params[i]; } + + // + // Allocate the x coordinate structure and convert row/col to image space. + // psVector *x = psVectorAlloc(2, PS_TYPE_F32); x->data.F32[0] = (psF32) (col + src->pixels->col0); x->data.F32[1] = (psF32) (row + src->pixels->row0); psF32 tmpF; - switch (src->models->type) { @@ -1139,5 +1174,7 @@ XXX: reverse order of row,col args? -XXX: The result is returned in subImage coords, not image coords. +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. *****************************************************************************/ psF32 findValue(psSource *source, @@ -1147,13 +1184,23 @@ psU32 dir) { + // + // Convert coords to subImage space. + // + psU32 subRow = row - source->pixels->row0; + psU32 subCol = col - source->pixels->col0; + // Ensure that the starting column is allowable. - if (!((0 <= col) && (col < source->pixels->numCols))) { + if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); return(NAN); } - - psF32 oldValue = evalModel(source, row, col); + if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { + psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); + return(NAN); + } + + psF32 oldValue = evalModel(source, subRow, subCol); if (oldValue == level) { - return(((psF32) col)); + return(((psF32) (subCol + source->pixels->col0))); } @@ -1171,32 +1218,38 @@ lastColumn = source->pixels->numCols; } - col+=incr; - - while (col != lastColumn) { - psF32 newValue = evalModel(source, row, col); + subCol+=incr; + + while (subCol != lastColumn) { + psF32 newValue = evalModel(source, subRow, subCol); if (oldValue == level) { - return(col); + return((psF32) (subCol + source->pixels->col0)); } if ((newValue <= level) && (level <= oldValue)) { // This is simple linear interpolation. - return( ((psF32) col) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); + return( ((psF32) (subCol + source->pixels->col0)) + ((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( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); + } + + subCol+=incr; } return(NAN); } -/****************************************************************************** + +/****************************************************************************** +pmSourceContour(src, img, level, mode): For an input subImage, and model, this +routine returns a psArray of coordinates that evaluate to the specified level. + XXX: Probably should remove the "image" argument. -XXX: What type should the output coordinate vectors consist of? +XXX: What type should the output coordinate vectors consist of? col,row? XXX: Why a pmArray output? XXX: doex x,y correspond with col,row or row/col? +XXX: What is momde? +XXX: The top, bottom of the contour is not correctly determined. *****************************************************************************/ psArray *pmSourceContour(psSource *source, @@ -1206,4 +1259,5 @@ { PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(image, false); PS_PTR_CHECK_NULL(source->moments, false); PS_PTR_CHECK_NULL(source->peak, false); @@ -1228,8 +1282,22 @@ // Starting at peak pixel, search leftwards for the column intercept. psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + return(NULL); + } xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; // Starting at peak pixel, search rightwards for the column intercept. + psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + return(NULL); + } + //printf("HERE the intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; @@ -1248,8 +1316,20 @@ // Starting at peak pixel, search leftwards for the column intercept. psF32 leftIntercept = findValue(source, level, row, col, 0); + if (isnan(leftIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + return(NULL); + } xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; // Starting at peak pixel, search rightwards for the column intercept. psF32 rightIntercept = findValue(source, level, row, col, 1); + if (isnan(rightIntercept)) { + psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); + psFree(xVec); + psFree(yVec); + return(NULL); + } xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; @@ -1294,5 +1374,4 @@ PS_IMAGE_CHECK_NULL(image, false); PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); - psBool rc; psS32 count = 0; @@ -1306,13 +1385,15 @@ psArray *x = psArrayAlloc(count); psVector *y = psVectorAlloc(count, PS_TYPE_F32); + psS32 tmpCnt = 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) { 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]; + // XXX: Convert i/j to image space: + coord->data.F32[0] = (psF32) (i + source->pixels->row0); + coord->data.F32[1] = (psF32) (j + source->pixels->col0); + x->data[tmpCnt] = (psPtr *) coord; + y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + tmpCnt++; } } @@ -1326,4 +1407,5 @@ switch (source->models->type) { case PS_MODEL_GAUSS: + rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); @@ -1374,13 +1456,6 @@ 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 100 - - psVector *params = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); + psVector *params = psVectorAlloc(src->models->Nparams, 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]; @@ -1390,11 +1465,14 @@ 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; + // Convert i/j to imace coord space: + // XXX: Make sure you have col/row order correct. + psS32 imageRow = i + src->pixels->row0; + psS32 imageCol = j + src->pixels->col0; + + x->data.F32[0] = (float) imageCol; + x->data.F32[1] = (float) imageRow; switch (src->models->type) { case PS_MODEL_GAUSS: @@ -1428,6 +1506,4 @@ // 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->pixels->col0; image->data.F32[imageRow][imageCol]+= pixelValue; @@ -1489,4 +1565,8 @@ XXX: I added the following. Must conform with IfA. If deriv==NULL, then we simply don't perform the derivative calculations. + +XXX: It is not clear whether the subImage coords, or the image coords should +be used when calling these functions. I don't think it really matters, as +long as we are consistent. *****************************************************************************/ psF32 pmMinLM_Gauss2D(psVector *deriv, @@ -1553,7 +1633,8 @@ } + psVector *tmpVec2 = (psVector *) x->data[i]; tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, params, - (psVector *) x->data[i]); + tmpVec2); for (psS32 j = 0 ; j < tmpRow->n ; j++) {