Index: trunk/psModules/src/objects/pmObjects.c =================================================================== --- trunk/psModules/src/objects/pmObjects.c (revision 5735) +++ trunk/psModules/src/objects/pmObjects.c (revision 5765) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-07 21:35:35 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:14:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -68,39 +68,4 @@ } -/****************************************************************************** -pmModelAlloc(): Allocate the pmModel structure, along with its parameters, -and initialize the type member. Initialize the params to 0.0. -XXX EAM: simplifying code with pmModelParameterCount -*****************************************************************************/ -pmModel *pmModelAlloc(pmModelType type) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); - - tmp->type = type; - tmp->chisq = 0.0; - tmp->nIter = 0; - psS32 Nparams = pmModelParameterCount(type); - if (Nparams == 0) { - psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); - return(NULL); - } - - tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); - tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); - - for (psS32 i = 0; i < tmp->params->n; i++) { - tmp->params->data.F32[i] = 0.0; - tmp->dparams->data.F32[i] = 0.0; - } - - psMemSetDeallocator(tmp, (psFreeFunc) modelFree); - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmp); -} - -/****************************************************************************** -XXX EAM : we can now free these pixels - memory ref is incremented now -*****************************************************************************/ static void sourceFree(pmSource *tmp) { @@ -115,4 +80,217 @@ psTrace(__func__, 4, "---- %s() end ----\n", __func__); } + +/****************************************************************************** +getRowVectorFromImage(): a private function which simply returns a +psVector containing the specified row of data from the psImage. + +XXX: Is there a better way to do this? +XXX EAM: does this really need to alloc a new vector??? +*****************************************************************************/ +static psVector *getRowVectorFromImage(psImage *image, + psU32 row) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + + psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); + for (psU32 col = 0; col < image->numCols ; col++) { + tmpVector->data.F32[col] = image->data.F32[row][col]; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(tmpVector); +} + +/****************************************************************************** +myListAddPeak(): A private function which allocates a psArray, if the list +argument is NULL, otherwise it adds the peak to that list. +XXX EAM : changed the output to psArray +XXX EAM : Switched row, col args +XXX EAM : NOTE: this was changed in the call, so the new code is consistent +*****************************************************************************/ +static psArray *myListAddPeak(psArray *list, + psS32 row, + psS32 col, + psF32 counts, + pmPeakType type) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); + + if (list == NULL) { + list = psArrayAlloc(100); + list->n = 0; + } + psArrayAdd(list, 100, tmpPeak); + psFree (tmpPeak); + // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) + + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(list); +} + + +/****************************************************************************** +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. +XXX: this is rather inefficient - at least compute and compare against radius^2 +*****************************************************************************/ +static bool checkRadius2(psF32 xCenter, + psF32 yCenter, + psF32 radius, + psF32 x, + psF32 y) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + /// XXX EAM should compare with hypot (x,y) for speed + if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { + return(true); + } + + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return(false); +} + +// XXX: Macro this. +static bool isItInThisRegion(const psRegion valid, + psS32 x, + psS32 y) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + if ((x >= valid.x0) && + (x <= valid.x1) && + (y >= valid.y0) && + (y <= valid.y1)) { + psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); + return(true); + } + psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); + return(false); +} + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +static psF32 findValue(pmSource *source, + psF32 level, + psU32 row, + psU32 col, + psU32 dir) +{ + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + // + // Convert coords to subImage space. + // + psU32 subRow = row - source->pixels->row0; + psU32 subCol = col - source->pixels->col0; + + // Ensure that the starting column is allowable. + if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { + psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + return(NAN); + } + if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); + return(NAN); + } + + // XXX EAM : i changed this to match pmModelEval above, but see + // XXX EAM the note below in pmSourceContour + psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); + if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return(((psF32) (subCol + source->pixels->col0))); + } + + // + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + // + psS32 incr; + psS32 lastColumn; + if (dir == 0) { + incr = -1; + lastColumn = -1; + } else { + incr = 1; + lastColumn = source->pixels->numCols; + } + subCol+=incr; + + while (subCol != lastColumn) { + psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); + if (oldValue == level) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return((psF32) (subCol + source->pixels->col0)); + } + + if ((newValue <= level) && (level <= oldValue)) { + // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); + } + + if ((oldValue <= level) && (level <= newValue)) { + // This is simple linear interpolation. + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); + } + + subCol+=incr; + } + + psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); + return(NAN); +} + +/****************************************************************************** +pmModelAlloc(): Allocate the pmModel structure, along with its parameters, +and initialize the type member. Initialize the params to 0.0. +XXX EAM: simplifying code with pmModelParameterCount +*****************************************************************************/ +pmModel *pmModelAlloc(pmModelType type) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel)); + + tmp->type = type; + tmp->chisq = 0.0; + tmp->nIter = 0; + psS32 Nparams = pmModelParameterCount(type); + if (Nparams == 0) { + psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); + return(NULL); + } + + tmp->params = psVectorAlloc(Nparams, PS_TYPE_F32); + tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32); + + for (psS32 i = 0; i < tmp->params->n; i++) { + tmp->params->data.F32[i] = 0.0; + tmp->dparams->data.F32[i] = 0.0; + } + + psMemSetDeallocator(tmp, (psFreeFunc) modelFree); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmp); +} + +/****************************************************************************** +XXX EAM : we can now free these pixels - memory ref is incremented now +*****************************************************************************/ /****************************************************************************** @@ -240,53 +418,4 @@ } -/****************************************************************************** -getRowVectorFromImage(): a private function which simply returns a -psVector containing the specified row of data from the psImage. - -XXX: Is there a better way to do this? -XXX EAM: does this really need to alloc a new vector??? -*****************************************************************************/ -static psVector *getRowVectorFromImage(psImage *image, - psU32 row) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - for (psU32 col = 0; col < image->numCols ; col++) { - tmpVector->data.F32[col] = image->data.F32[row][col]; - } - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(tmpVector); -} - -/****************************************************************************** -myListAddPeak(): A private function which allocates a psArray, if the list -argument is NULL, otherwise it adds the peak to that list. -XXX EAM : changed the output to psArray -XXX EAM : Switched row, col args -XXX EAM : NOTE: this was changed in the call, so the new code is consistent -*****************************************************************************/ -static psArray *myListAddPeak(psArray *list, - psS32 row, - psS32 col, - psF32 counts, - pmPeakType type) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - pmPeak *tmpPeak = pmPeakAlloc(col, row, counts, type); - - if (list == NULL) { - list = psArrayAlloc(100); - list->n = 0; - } - psArrayAdd(list, 100, tmpPeak); - psFree (tmpPeak); - // XXX EAM : is this free appropriate? (does psArrayAdd increment memory counter?) - - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(list); -} /****************************************************************************** @@ -503,21 +632,4 @@ } -// XXX: Macro this. -static bool isItInThisRegion(const psRegion valid, - psS32 x, - psS32 y) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - if ((x >= valid.x0) && - (x <= valid.x1) && - (y >= valid.y0) && - (y <= valid.y1)) { - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return(true); - } - psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); - return(false); -} - /****************************************************************************** @@ -627,12 +739,4 @@ *****************************************************************************/ - - - - - - - - bool pmSourceLocalSky( pmSource *source, @@ -673,27 +777,4 @@ psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); return (true); -} - -/****************************************************************************** -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. -XXX: this is rather inefficient - at least compute and compare against radius^2 -*****************************************************************************/ -static bool checkRadius2(psF32 xCenter, - psF32 yCenter, - psF32 radius, - psF32 x, - psF32 y) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - /// XXX EAM should compare with hypot (x,y) for speed - if ((PS_SQR(x - xCenter) + PS_SQR(y - yCenter)) < PS_SQR(radius)) { - return(true); - } - - psTrace(__func__, 4, "---- %s(false) end ----\n", __func__); - return(false); } @@ -1350,90 +1431,4 @@ /****************************************************************************** -findValue(source, level, row, col, dir): a private function which determines -the column coordinate of the model function which has the value "level". If -dir equals 0, then you loop leftwards from the peak pixel, otherwise, -rightwards. - -XXX: reverse order of row,col args? - -XXX: Input row/col are in image coords. - -XXX: The result is returned in image coords. -*****************************************************************************/ -static psF32 findValue(pmSource *source, - psF32 level, - psU32 row, - psU32 col, - psU32 dir) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - // - // Convert coords to subImage space. - // - psU32 subRow = row - source->pixels->row0; - psU32 subCol = col - source->pixels->col0; - - // Ensure that the starting column is allowable. - if (!((0 <= subCol) && (subCol < source->pixels->numCols))) { - psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - if (!((0 <= subRow) && (subRow < source->pixels->numRows))) { - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - psError(PS_ERR_UNKNOWN, true, "Starting row outside subImage range"); - return(NAN); - } - - // XXX EAM : i changed this to match pmModelEval above, but see - // XXX EAM the note below in pmSourceContour - psF32 oldValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); - if (oldValue == level) { - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return(((psF32) (subCol + source->pixels->col0))); - } - - // - // We define variables incr and lastColumn so that we can use the same loop - // whether we are stepping leftwards, or rightwards. - // - psS32 incr; - psS32 lastColumn; - if (dir == 0) { - incr = -1; - lastColumn = -1; - } else { - incr = 1; - lastColumn = source->pixels->numCols; - } - subCol+=incr; - - while (subCol != lastColumn) { - psF32 newValue = pmModelEval(source->modelFLT, source->pixels, subCol, subRow); - if (oldValue == level) { - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return((psF32) (subCol + source->pixels->col0)); - } - - if ((newValue <= level) && (level <= oldValue)) { - // This is simple linear interpolation. - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); - } - - if ((oldValue <= level) && (level <= newValue)) { - // This is simple linear interpolation. - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return( ((psF32) (subCol + source->pixels->col0)) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); - } - - subCol+=incr; - } - - psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__); - return(NAN); -} - -/****************************************************************************** pmSourceContour(src, img, level, mode): For an input subImage, and model, this routine returns a psArray of coordinates that evaluate to the specified level. @@ -1933,3 +1928,29 @@ } - +bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) +{ + + float obsSum = 0; + float fitSum = 0; + float sky = model->params->data.F32[0]; + + pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + + for (int ix = 0; ix < image->numCols; ix++) { + for (int iy = 0; iy < image->numRows; iy++) { + if (mask->data.U8[iy][ix]) + continue; + obsSum += image->data.F32[iy][ix] - sky; + } + } + if (obsSum <= 0) + return false; + if (fitSum <= 0) + return false; + + *fitMag = -2.5*log10(fitSum); + *obsMag = -2.5*log10(obsSum); + return (true); +} +