Index: /branches/rel9/psModules/Makefile.am =================================================================== --- /branches/rel9/psModules/Makefile.am (revision 5767) +++ /branches/rel9/psModules/Makefile.am (revision 5768) @@ -6,5 +6,9 @@ pkgconfig_DATA= psmodule.pc -EXTRA_DIST = Doxyfile.in psmodule-config.in psmodule.pc.in +EXTRA_DIST = \ + Doxyfile.in \ + psmodule-config.in \ + psmodule.pc.in \ + autogen.sh if DOXYGEN Index: /branches/rel9/psModules/configure.ac =================================================================== --- /branches/rel9/psModules/configure.ac (revision 5767) +++ /branches/rel9/psModules/configure.ac (revision 5768) @@ -1,5 +1,5 @@ AC_PREREQ(2.59) -AC_INIT([psmodule],[0.9.0],[http://pan-starrs.ifa.hawaii.edu/bugzilla]) +AC_INIT([psmodule],[0.9.1],[http://pan-starrs.ifa.hawaii.edu/bugzilla]) AC_CONFIG_SRCDIR([psmodule.pc.in]) @@ -37,11 +37,24 @@ AC_SUBST([AM_CFLAGS]) +dnl ------------------- PERL options --------------------- + AC_ARG_WITH(perl, + [AS_HELP_STRING(--with-perl=FILE,Specify location of PERL executable.)], + [AC_CHECK_PROG(PERL, $withval, $withval)], + [AC_CHECK_PROG(PERL, perl, `which perl`)]) + if test "$PERL" == "" + then + AC_MSG_ERROR([PERL is required. Use --with-perl to specify its install location.]) + fi + AC_SUBST(PERL,$PERL) + SRCPATH='${top_srcdir}/src' -SRCDIRS="astrom camera config detrend imcombine imsubtract objects photom" -SRCINC=`echo "${SRCDIRS=}" | ${SED} "s|\(\\w\+\)|-I\${SRCPATH=}/\1|g"` -SRCSUBLIBS=`echo "${SRCDIRS=}" | ${SED} "s|\(\\w\+\)|\1/libpsmodule\1.la|g"` -AC_SUBST([SRCSUBLIBS],${SRCSUBLIBS=}) +SRCDIRS="astrom config detrend imcombine imsubtract objects" +# escape two escapes at this level so \\ gets passed to the shell and \ to perl +SRCINC=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|-I\\\\${SRCPATH=}/\1|g"` +SRCINC="-I${SRCPATH=} ${SRCINC=}" +SRCSUBLIBS=`echo "${SRCDIRS=}" | ${PERL} -pe "s|(\w+)|\1/libpsmodule\1.la|g"` +AC_SUBST(SRCSUBLIBS,${SRCSUBLIBS=}) +AC_SUBST(SRCINC,${SRCINC=}) AC_SUBST([SRCDIRS],${SRCDIRS=}) -AC_SUBST([SRCINC],${SRCINC=}) dnl doxygen ------------------------------------------------------------------- @@ -56,5 +69,5 @@ if test -z ${PSLIB_CONFIG} ; then - PKG_CHECK_MODULES([PSLIB], [pslib >= 0.8.0]) + PKG_CHECK_MODULES([PSLIB], [pslib >= 0.9.0]) else AC_CHECK_FILE($PSLIB_CONFIG,[], @@ -79,5 +92,4 @@ src/Makefile src/astrom/Makefile - src/camera/Makefile src/config/Makefile src/detrend/Makefile @@ -85,8 +97,6 @@ src/imsubtract/Makefile src/objects/Makefile - src/photom/Makefile test/Makefile test/astrom/Makefile - test/camera/Makefile test/config/Makefile test/detrend/Makefile @@ -94,5 +104,4 @@ test/imsubtract/Makefile test/objects/Makefile - test/photom/Makefile Doxyfile psmodule-config Index: /branches/rel9/psModules/src/astrom/pmAstrometry.c =================================================================== --- /branches/rel9/psModules/src/astrom/pmAstrometry.c (revision 5767) +++ /branches/rel9/psModules/src/astrom/pmAstrometry.c (revision 5768) @@ -13,6 +13,6 @@ * XXX: Should we implement non-linear cell->chip transforms? * -* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-05 21:28:55 $ +* @version $Revision: 1.10.2.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-12 21:52:17 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -799,2 +799,76 @@ return(cellCoord); } + +/***************************************************************************** + *****************************************************************************/ +bool pmFPASelectChip( + pmFPA *fpa, + int chipNum) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + if ((fpa->chips == NULL) || (chipNum >= fpa->chips->n)) { + return(false); + } + psBool rc = true; + + for (psS32 i = 0 ; i < fpa->chips->n ; i++) { + pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; + if (tmpChip == NULL) { + rc = false; + } else { + if (i == chipNum) { + tmpChip->valid = true; + } else { + tmpChip->valid = false; + } + } + } + + return(rc); +} + + +/***************************************************************************** +XXX: The SDRS is ambiguous on a few things: + Whether or not the other chips should be set valid=true. + Should we return the number of chip valid=true before or after they're set, + *****************************************************************************/ +/** + * + * pmFPAExcludeChip shall set valid to false only for the specified chip + * number (chipNum). In the event that the specified chip number does not exist + * within the fpa, the function shall generate a warning, and perform no action. + * The function shall return the number of chips within the fpa that have valid + * set to true. + * + */ +int pmFPAExcludeChip( + pmFPA *fpa, + int chipNum) +{ + PS_ASSERT_PTR_NON_NULL(fpa, false); + + if (fpa->chips == NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n"); + return(0); + } + if ((chipNum >= fpa->chips->n) || (NULL == (pmChip *) fpa->chips->data[chipNum])) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum); + return(0); + } + + psS32 numChips = 0; + for (psS32 i = 0 ; i < fpa->chips->n ; i++) { + pmChip *tmpChip = (pmChip *) fpa->chips->data[i]; + if (tmpChip != NULL) { + if (i == chipNum) { + tmpChip->valid = false; + } else { + tmpChip->valid = true; + numChips++; + } + } + } + + return(numChips); +} Index: /branches/rel9/psModules/src/astrom/pmAstrometry.h =================================================================== --- /branches/rel9/psModules/src/astrom/pmAstrometry.h (revision 5767) +++ /branches/rel9/psModules/src/astrom/pmAstrometry.h (revision 5768) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-05 21:28:55 $ +* @version $Revision: 1.4.2.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-12 21:52:17 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -423,3 +423,32 @@ ); +/** + * + * pmFPASelectChip shall set valid to true for the specified chip number + * (chipNum), and all other chips shall have valid set to false. In the event + * that the specified chip number does not exist within the fpa, the function + * shall return false. + * + */ +bool pmFPASelectChip( + pmFPA *fpa, + int chipNum +); + +/** + * + * pmFPAExcludeChip shall set valid to false only for the specified chip + * number (chipNum). In the event that the specified chip number does not exist + * within the fpa, the function shall generate a warning, and perform no action. + * The function shall return the number of chips within the fpa that have valid + * set to true. + * + */ +int pmFPAExcludeChip( + pmFPA *fpa, + int chipNum +); + + #endif // #ifndef PS_ASTROMETRY_H + Index: /branches/rel9/psModules/src/config/pmConfig.c =================================================================== --- /branches/rel9/psModules/src/config/pmConfig.c (revision 5767) +++ /branches/rel9/psModules/src/config/pmConfig.c (revision 5768) @@ -3,6 +3,6 @@ * @author PAP, IfA * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-15 20:09:03 $ + * @version $Revision: 1.4.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:18 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -28,5 +28,5 @@ const char *description) // Description of file { - int numBadLines = 0; + unsigned int numBadLines = 0; psLogMsg(__func__, PS_LOG_INFO, "Loading %s configuration from file %s\n", @@ -359,5 +359,5 @@ } else if (cameraItem->type == PS_DATA_STRING) { psTrace(__func__, 5, "Reading camera configuration for %s...\n", cameraItem->name); - int badLines = 0; // Number of bad lines in reading camera configuration + unsigned int badLines = 0; // Number of bad lines in reading camera configuration camera = psMetadataConfigParse(NULL, &badLines, cameraItem->data.V, true); if (badLines > 0) { Index: /branches/rel9/psModules/src/imcombine/pmImageCombine.c =================================================================== --- /branches/rel9/psModules/src/imcombine/pmImageCombine.c (revision 5767) +++ /branches/rel9/psModules/src/imcombine/pmImageCombine.c (revision 5768) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-28 20:43:52 $ + * @version $Revision: 1.1.10.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:20 $ * * XXX: pmRejectPixels() has a known bug with the pmImageTransform() call. @@ -35,15 +35,15 @@ XXX: Allocate a dummy psStats structure so that we don't destroy away its data. *****************************************************************************/ -psImage *pmCombineImages(psImage *combine, ///< Combined image (output) - psArray **questionablePixels, ///< Array of rejection masks - const psArray *images, ///< Array of input images - const psArray *errors, ///< Array of input error images - const psArray *masks, ///< Array of input masks - psU32 maskVal, ///< Mask value - const psPixels *pixels, ///< Pixels to combine - psS32 numIter, ///< Number of rejection iterations - psF32 sigmaClip, ///< Number of standard deviations at which to reject - const psStats *stats ///< Statistics to use in the combination - ) +psImage *pmCombineImages( + psImage *combine, ///< Combined image (output) + psArray **questionablePixels, ///< Array of rejection masks + const psArray *images, ///< Array of input images + const psArray *errors, ///< Array of input error images + const psArray *masks, ///< Array of input masks + psU32 maskVal, ///< Mask value + const psPixels *pixels, ///< Pixels to combine + psS32 numIter, ///< Number of rejection iterations + psF32 sigmaClip, ///< Number of standard deviations at which to reject + const psStats *stats) ///< Statistics to use in the combination { Index: /branches/rel9/psModules/src/objects/Makefile.am =================================================================== --- /branches/rel9/psModules/src/objects/Makefile.am (revision 5767) +++ /branches/rel9/psModules/src/objects/Makefile.am (revision 5768) @@ -8,6 +8,4 @@ pmPSFtry.c \ pmModelGroup.c \ - psModulesUtils.c \ - psLibUtils.c \ psEllipse.c @@ -24,5 +22,3 @@ pmPSFtry.h \ pmModelGroup.h \ - psModulesUtils.h \ - psLibUtils.h \ psEllipse.h Index: /branches/rel9/psModules/src/objects/pmObjects.c =================================================================== --- /branches/rel9/psModules/src/objects/pmObjects.c (revision 5767) +++ /branches/rel9/psModules/src/objects/pmObjects.c (revision 5768) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-20 23:06:24 $ + * @version $Revision: 1.3.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:22 $ * * 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); } @@ -1107,5 +1188,5 @@ bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); if ((Nsatpix > 1) && big) { - tmpSrc->type |= PM_SOURCE_SATSTAR; + tmpSrc->type = PM_SOURCE_SATSTAR; Nsatstar ++; continue; @@ -1114,5 +1195,5 @@ // saturated object (not a star, eg bleed trails, hot pixels) if (Nsatpix > 1) { - tmpSrc->type |= PM_SOURCE_SATURATED; + tmpSrc->type = PM_SOURCE_SATURATED; Nsat ++; continue; @@ -1123,5 +1204,5 @@ // only set candidate defects if if ((sigX < 0.05) || (sigY < 0.05)) { - tmpSrc->type |= PM_SOURCE_DEFECT; + tmpSrc->type = PM_SOURCE_DEFECT; Ncr ++; continue; @@ -1130,5 +1211,5 @@ // likely unsaturated galaxy (too large to be stellar) if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { - tmpSrc->type |= PM_SOURCE_GALAXY; + tmpSrc->type = PM_SOURCE_GALAXY; Ngal ++; continue; @@ -1143,5 +1224,5 @@ psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); if ((SN > PSF_SN_LIM) && (radius < 1.5)) { - tmpSrc->type |= PM_SOURCE_PSFSTAR; + tmpSrc->type = PM_SOURCE_PSFSTAR; Npsf ++; continue; @@ -1149,5 +1230,5 @@ // random type of star - tmpSrc->type |= PM_SOURCE_OTHER; + tmpSrc->type = PM_SOURCE_OTHER; } @@ -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); +} + Index: /branches/rel9/psModules/src/objects/pmObjects.h =================================================================== --- /branches/rel9/psModules/src/objects/pmObjects.h (revision 5767) +++ /branches/rel9/psModules/src/objects/pmObjects.h (revision 5768) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-10 19:53:40 $ + * @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:22 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -36,5 +36,5 @@ * values with specific meanings that other functions can add to or define? */ -enum { +typedef enum { PSPHOT_MASK_CLEAR = 0x00, PSPHOT_MASK_INVALID = 0x01, @@ -620,3 +620,29 @@ ); + +/** pmSourcePhotometry() + * + * XXX: Need descriptions + * + */ +bool pmSourcePhotometry( + float *fitMag, + float *obsMag, + pmModel *model, + psImage *image, + psImage *mask +); + +/** pmModelEval() + * + * XXX: Need descriptions + * + */ +psF32 pmModelEval( + pmModel *model, + psImage *image, + psS32 col, + psS32 row +); + #endif Index: /branches/rel9/psModules/src/objects/pmPSF.c =================================================================== --- /branches/rel9/psModules/src/objects/pmPSF.c (revision 5767) +++ /branches/rel9/psModules/src/objects/pmPSF.c (revision 5768) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-10 19:53:40 $ + * @version $Revision: 1.1.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:22 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -18,5 +18,4 @@ #include -#include "psLibUtils.h" #include "pmObjects.h" #include "pmPSF.h" @@ -94,6 +93,5 @@ for (int i = 0; i < psf->params->n; i++) { // XXX EAM : make this a user-defined value? - // XXX EAM : future version (0.7.0?) psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1); - psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1); + psf->params->data[i] = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); } Index: /branches/rel9/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/rel9/psModules/src/objects/pmPSFtry.c (revision 5767) +++ /branches/rel9/psModules/src/objects/pmPSFtry.c (revision 5768) @@ -1,6 +1,17 @@ +/** @file pmPSFtry.c + * + * XXX: need description of file purpose + * + * @author EAM, IfA + * + * @version $Revision: 1.1.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:22 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + # include -# include "psLibUtils.h" # include "pmObjects.h" -# include "psModulesUtils.h" # include "pmPSF.h" # include "pmPSFtry.h" @@ -41,5 +52,5 @@ type = pmModelSetType (modelName); test->psf = pmPSFAlloc (type); - test->sources = psMemCopy(sources); + test->sources = psMemIncrRefCounter(sources); test->modelFLT = psArrayAlloc (sources->n); test->modelPSF = psArrayAlloc (sources->n); @@ -79,123 +90,102 @@ int Npsf = 0; - pmPSFtry *try - = pmPSFtryAlloc (sources, modelName); + pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName); // stage 1: fit an independent model (freeModel) to all sources psTimerStart ("fit"); - for (int i = 0; i < try - ->sources->n; i++) { - - pmSource *source = try - ->sources->data[i]; - pmModel *model = pmSourceModelGuess (source, try - ->psf->type); - x = source->peak->x; - y = source->peak->y; - - // set temporary object mask and fit object - // fit model as FLT, not PSF - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, model, false); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - // exclude the poor fits - if (!status) { - try - ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; - psFree (model); - continue; - } - try - ->modelFLT->data[i] = model; - Nflt ++; - } + for (int i = 0; i < psfTry->sources->n; i++) { + + pmSource *source = psfTry->sources->data[i]; + pmModel *model = pmSourceModelGuess (source, psfTry->psf->type); + x = source->peak->x; + y = source->peak->y; + + // set temporary object mask and fit object + // fit model as FLT, not PSF + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, model, false); + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + // exclude the poor fits + if (!status) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL; + psFree (model); + continue; + } + psfTry->modelFLT->data[i] = model; + Nflt ++; + } psLogMsg ("psphot.psftry", 4, "fit flt: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n); // make this optional? - // DumpModelFits (try->modelFLT, "modelsFLT.dat"); + // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat"); // stage 2: construct a psf (pmPSF) from this collection of model fits - pmPSFFromModels (try - ->psf, try - ->modelFLT, try - ->mask); + pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask); // stage 3: refit with fixed shape parameters psTimerStart ("fit"); - for (int i = 0; i < try - ->sources->n; i++) { - // masked for: bad model fit, outlier in parameters - if (try - ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - - pmSource *source = try - ->sources->data[i]; - pmModel *modelFLT = try - ->modelFLT->data[i]; - - // set shape for this model based on PSF - pmModel *modelPSF = pmModelFromPSF (modelFLT, try - ->psf); - x = source->peak->x; - y = source->peak->y; - - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); - status = pmSourceFitModel (source, modelPSF, true); - - // skip poor fits - if (!status) { - try - ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; - psFree (modelPSF); - goto next_source; - } - - // otherwise, save the resulting model - try - ->modelPSF->data[i] = modelPSF; - - // XXX : use a different aperture radius from the fit radius? - // XXX : use a different estimator for the local sky? - // XXX : pass 'source' as input? - if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { - try - ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; - goto next_source; - } - - try - ->metric->data.F64[i] = obsMag - fitMag; - try - ->fitMag->data.F64[i] = fitMag; - Npsf ++; + for (int i = 0; i < psfTry->sources->n; i++) { + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + + pmSource *source = psfTry->sources->data[i]; + pmModel *modelFLT = psfTry->modelFLT->data[i]; + + // set shape for this model based on PSF + pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf); + x = source->peak->x; + y = source->peak->y; + + psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED); + status = pmSourceFitModel (source, modelPSF, true); + + // skip poor fits + if (!status) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL; + psFree (modelPSF); + goto next_source; + } + + // otherwise, save the resulting model + psfTry->modelPSF->data[i] = modelPSF; + + // XXX : use a different aperture radius from the fit radius? + // XXX : use a different estimator for the local sky? + // XXX : pass 'source' as input? + if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) { + psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT; + goto next_source; + } + + psfTry->metric->data.F64[i] = obsMag - fitMag; + psfTry->fitMag->data.F64[i] = fitMag; + Npsf ++; next_source: - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); - - } + psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED); + + } psLogMsg ("psphot.psftry", 4, "fit psf: %f sec for %d sources\n", psTimerMark ("fit"), sources->n); psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n); // make this optional - // DumpModelFits (try->modelPSF, "modelsPSF.dat"); + // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat"); // XXX this function wants aperture radius for pmSourcePhotometry - pmPSFtryMetric (try - , RADIUS); + pmPSFtryMetric (psfTry, RADIUS); psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", - modelName, try - ->psf->ApResid, try - ->psf->dApResid, try - ->psf->skyBias); - - return (try - ); + modelName, + psfTry->psf->ApResid, + psfTry->psf->dApResid, + psfTry->psf->skyBias); + + return (psfTry); } -bool pmPSFtryMetric (pmPSFtry *try - , float RADIUS) +bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS) { @@ -203,5 +193,5 @@ int nKeep, nSkip; - // the measured (aperture - fit) magnitudes (dA == try->metric) + // the measured (aperture - fit) magnitudes (dA == psfTry->metric) // depend on both the true ap-fit (dAo) and the bias in the sky measurement: // dA = dAo + dsky/flux @@ -213,18 +203,14 @@ // rflux = ten(0.4*fitMag); - psVector *rflux = psVectorAlloc (try - ->sources->n, PS_TYPE_F64); - for (int i = 0; i < try - ->sources->n; i++) { - if (try - ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue; - rflux->data.F64[i] = pow(10.0, 0.4*try - ->fitMag->data.F64[i]); - } + psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); + for (int i = 0; i < psfTry->sources->n; i++) { + if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL) + continue; + rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]); + } // find min and max of (1/flux): psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - psVectorStats (stats, rflux, NULL, try - ->mask, PSFTRY_MASK_ALL); + psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL); // build binned versions of rflux, metric @@ -240,29 +226,26 @@ for (int i = 0; i < daBin->n; i++) { - psVector *tmp = psVectorAlloc (try - ->sources->n, PS_TYPE_F64); + psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64); tmp->n = 0; // accumulate data within bin range - for (int j = 0; j < try - ->sources->n; j++) { - // masked for: bad model fit, outlier in parameters - if (try - ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue; - - // skip points with extreme dA values - if (fabs(try - ->metric->data.F64[j]) > 0.5) continue; - - // skip points outside of this bin - if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) - continue; - if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) - continue; - - tmp->data.F64[tmp->n] = try - ->metric->data.F64[j]; - tmp->n ++; - } + for (int j = 0; j < psfTry->sources->n; j++) { + // masked for: bad model fit, outlier in parameters + if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL) + continue; + + // skip points with extreme dA values + if (fabs(psfTry->metric->data.F64[j]) > 0.5) + continue; + + // skip points outside of this bin + if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) + continue; + if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) + continue; + + tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j]; + tmp->n ++; + } // is this a valid point? @@ -319,10 +302,7 @@ stats = psVectorStats (stats, daResid, NULL, maskB, 1); - try - ->psf->ApResid = poly->coeff[0]; - try - ->psf->dApResid = stats->clippedStdev; - try - ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); + psfTry->psf->ApResid = poly->coeff[0]; + psfTry->psf->dApResid = stats->clippedStdev; + psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS)); psFree (rflux); Index: /branches/rel9/psModules/src/objects/pmPSFtry.h =================================================================== --- /branches/rel9/psModules/src/objects/pmPSFtry.h (revision 5767) +++ /branches/rel9/psModules/src/objects/pmPSFtry.h (revision 5768) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-10-10 19:53:40 $ + * @version $Revision: 1.1.8.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:22 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -105,7 +105,6 @@ */ bool pmPSFtryMetric( - pmPSFtry *try - , ///< Add comment. - float RADIUS ///< Add comment. + pmPSFtry *psfTry, ///< Add comment. + float RADIUS ///< Add comment. ); Index: anches/rel9/psModules/src/objects/psLibUtils.c =================================================================== --- /branches/rel9/psModules/src/objects/psLibUtils.c (revision 5767) +++ (revision ) @@ -1,509 +1,0 @@ -# include // for strncasecmp -# include -# include "psLibUtils.h" - -// XXX EAM : this is NOT included in the C99 headers ?? -FILE *fdopen(int fildes, const char *mode); - -// XXX EAM : these utility functions should be added back into PSLib - -static psHash *timers = NULL; - -/* XXX: remove -bool psTimerClear (char *name) { - - bool status; - - if (name == NULL) return false; - - status = psHashRemove (timers, name); - return (status); -} -*/ - -void psTimerFree () -{ - - psFree (timers); - p_psTimeFinalize(); - return; -} - -// start/restart a named timer -bool psTimerStart (char *name) -{ - - psTime *start; - - if (timers == NULL) - timers = psHashAlloc (16); - - start = psTimeGetNow (PS_TIME_UTC); - psHashAdd (timers, name, start); - psFree (start); - return (TRUE); -} - -// get current elapsed time on named timer -psF64 psTimerMark (char *name) -{ - - psTime *start; - psTime *mark; - psF64 delta; - - if (timers == NULL) - return (0); - - start = psHashLookup (timers, name); - if (start == NULL) - return (0); - - mark = psTimeGetNow (PS_TIME_UTC); - delta = psTimeDelta (mark, start); - psFree (mark); - // psFree (start); -- XXX is psHashLookup not psMemCopying? - - return (delta); -} - -/* XXX: remove -// find the location of the specified argument -int psArgumentGet (int argc, char **argv, char *arg) { - - int i; - - for (i = 0; i < argc; i++) { - if (!strcmp(argv[i], arg)) - return (i); - } - - return ((int)NULL); -} -*/ - -/* XXX: remove -// remove the specified argument (by location) -int psArgumentRemove (int N, int *argc, char **argv) { - - int i; - - if ((N != (int)NULL) && (N != 0)) { - (*argc)--; - for (i = N; i < *argc; i++) { - argv[i] = argv[i+1]; - } - } - - return (N); -} -*/ - -// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) -// 2 = default, -v = 3, -vv = 4, -vvv = 5 -psS32 psLogArguments (int *argc, char **argv) -{ - - int N, level; - - // default log level is 2 - level = 2; - - // set in order, so that -vvv overrides -vv overrides -v - if ((N = psArgumentGet (*argc, argv, "-v"))) { - psArgumentRemove (N, argc, argv); - level = 3; - } - if ((N = psArgumentGet (*argc, argv, "-vv"))) { - psArgumentRemove (N, argc, argv); - level = 4; - } - if ((N = psArgumentGet (*argc, argv, "-vvv"))) { - psArgumentRemove (N, argc, argv); - level = 5; - } - - if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { - if (*argc < N + 2) { - psAbort ("psLogArguments", "USAGE: -logfmt (format)"); - } - psArgumentRemove (N, argc, argv); - psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid - psArgumentRemove (N, argc, argv); - } - - // set the level, return the level - psLogSetLevel (level); - return (level); -} - -// set trace levels by facility -psS32 psTraceArguments (int *argc, char **argv) -{ - - int N; - - // argument format is: -trace (facil) (level) - while ((N = psArgumentGet (*argc, argv, "-trace"))) { - if (*argc < N + 3) { - psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); - } - psArgumentRemove (N, argc, argv); - psTraceSetLevel (argv[N], atoi(argv[N+1])); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { - psTracePrintLevels (); - exit (2); - } - return (TRUE); -} - -# if (0) - // alternate implementation of this function from pmObjects.c - // now defined in psLib SDRS as psImageRow - psVector *psGetRowVectorFromImage(psImage *image, - psU32 row) -{ - PS_ASSERT_IMAGE_NON_NULL(image, NULL); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32)); - return(tmpVector); -} -# endif - -// XXX EAM : this is now in psLib -void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma) -{ - - int Nx, Ny, Npixel, Nrange; - float factor, g, s; - psVector *temp; - - // relevant terms - Nrange = sigma*Nsigma + 0.5; - Npixel = 2*Nrange + 1; - factor = -0.5/(sigma*sigma); - - Nx = image->numCols; - Ny = image->numRows; - - // generate gaussian - psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32); - for (int i = -Nrange; i < Nrange + 1; i++) { - gaussnorm->data.F32[i+Nrange] = exp (factor*i*i); - } - psF32 *gauss = &gaussnorm->data.F32[Nrange]; - - // smooth in X direction - temp = psVectorAlloc (Nx, PS_TYPE_F32); - for (int j = 0; j < Ny; j++) { - psF32 *vi = image->data.F32[j]; - psF32 *vo = temp->data.F32; - for (int i = 0; i < Nx; i++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (i+n < 0) - continue; - if (i+n >= Nx) - continue; - s += gauss[n]*vi[i+n]; - g += gauss[n]; - } - vo[i] = s / g; - } - memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32)); - } - psFree (temp); - - // smooth in Y direction - temp = psVectorAlloc (image->numRows, PS_TYPE_F32); - for (int i = 0; i < Nx; i++) { - psF32 *vo = temp->data.F32; - psF32 **vi = image->data.F32; - for (int j = 0; j < Ny; j++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (j+n < 0) - continue; - if (j+n >= Ny) - continue; - s += gauss[n]*vi[j+n][i]; - g += gauss[n]; - } - vo[j] = s / g; - } - // replace temp in image - for (int j = 0; j < Ny; j++) { - vi[j][i] = vo[j]; - } - } - psFree (temp); - psFree (gaussnorm); -} - -bool psImageInit (psImage *image,...) -{ - - va_list argp; - psU8 vU8; - psF32 vF32; - psF64 vF64; - - if (image == NULL) - return (false); - - va_start (argp, image); - - switch (image->type.type) { - case PS_TYPE_U8: - vU8 = va_arg (argp, psU32); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.U8[iy][ix] = vU8; - } - } - break; - - case PS_TYPE_F32: - vF32 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F32[iy][ix] = vF32; - } - } - return (true); - - case PS_TYPE_F64: - vF64 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F64[iy][ix] = vF64; - } - } - return (true); - - default: - psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type); - return (false); - } - return (false); -} - -/* XXX: remove -// count number of pixels with given mask value -int psImageCountPixelMask (psImage *mask, psU8 value) -{ - int Npixels = 0; - - for (int i = 0; i < mask->numRows; i++) { - for (int j = 0; j < mask->numCols; j++) { - if (mask->data.U8[i][j] & value) { - Npixels ++; - } - } - } - return (Npixels); -} -*/ - -// define a square region centered on the given coordinate -// XXX EAM : this is now in psLib -# if (0) - psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) -{ - psRegion *region; - region = psRegionAlloc (x - radius, x + radius + 1, - y - radius, y + radius + 1); - return (region); -} -# endif - -// set actual region based on image parameters: -// compensate for negative upper limits -// XXX this is inconsistent: the coordindates should always be in the parent -// frame, which means the negative values should subtract from Nx,Ny of -// the parent, not the child. but, we don't carry the dimensions of the -// parent in the psImage container. for now, use the child Nx,Ny -// force range to be on this subimage -// XXX EAM : this needs to be changed to use psRegion rather than psRegion* -// XXX EAM : this is now in psLib -# if (0) - psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) -{ - - // x0,y0, x1,y1 are in *parent* units - - if (out == NULL) { - out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1); - } else { - *out = *in; - } - // XXX these are probably wrong (see above) - if (out->x1 <= 0) { - out->x1 = image->col0 + image->numCols + out->x1; - } - if (out->y1 <= 0) { - out->y1 = image->row0 + image->numRows + out->y1; - } - - // force the lower-limits to be on the child - out->x0 = PS_MAX(image->col0, out->x0); - out->y0 = PS_MAX(image->row0, out->y0); - - // force the upper-limits to be on the child - out->x1 = PS_MIN(image->col0 + image->numCols, out->x1); - out->y1 = PS_MIN(image->row0 + image->numRows, out->y1); - return (out); -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) -{ - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) - continue; - if (ix + image->col0 >= region->x1) - continue; - if (iy + image->row0 < region->y0) - continue; - if (iy + image->row0 >= region->y1) - continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area not contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) -{ - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) - goto maskit; - if (ix + image->col0 >= region->x1) - goto maskit; - if (iy + image->row0 < region->y0) - goto maskit; - if (iy + image->row0 >= region->y1) - goto maskit; - continue; -maskit: - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) -{ - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 > R2) - continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -// mask the area contained by the region -// the region is defined wrt the parent image -// XXX EAM : this is now in psLib -# if (0) - void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) -{ - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 < R2) - continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} -# endif - -/* XXX: remove -psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) { - - int nBin = (upper - lower) / delta; - - psVector *out = psVectorAlloc (nBin, type); - - for (int i = 0; i < nBin; i++) { - out->data.F64[i] = lower + i * delta; - } - - return (out); -} -*/ - -// XXX EAM a utility function -bool p_psVectorPrintRow (int fd, psVector *a, char *name) -{ - - FILE *f; - f = fdopen(fd, "a+"); - fprintf (f, "vector: %s\n", name); - - for (int i = 0; i < a[0].n; i++) { - fprintf (f, "%f ", p_psVectorGetElementF64(a, i)); - } - fprintf (f, "\n"); - fclose(f); - return (true); -} -// XXX EAM is the use of fdopen here a good way to do this? Index: anches/rel9/psModules/src/objects/psLibUtils.h =================================================================== --- /branches/rel9/psModules/src/objects/psLibUtils.h (revision 5767) +++ (revision ) @@ -1,59 +1,0 @@ - -# ifndef PS_LIB_UTILS -# define PS_LIB_UTILS - -// XXX EAM : psEllipse needs to be be included in psLib -# include "psEllipse.h" - -// structure to carry a dynamic string -typedef struct -{ - int NLINE; - int Nline; - char *line; -} -psLine; - -# define psMemCopy(A)(psMemIncrRefCounter((A))) - -// XXX EAM : my version using varience instead of stdev -bool psMinimizeGaussNewtonDelta_EAM (psVector *delta, - const psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yErr, - psMinimizeLMChi2Func func); - -// minimize : using varience vs sigma and parameter limits -psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, psVector *Beta, psVector *Params, const psImage *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda); -psBool psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func); -psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda); - -// psLib extra utilities -bool psTimerStart (char *name); // added to SDRS -//XXX: I removed this: bool psTimerClear (char *name); // added to SDRS -psF64 psTimerMark (char *name); // added to SDRS -void psTimerFree (); // added to SDRS (as psTimerStop) -psS32 psLogArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) -psS32 psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity) -//XXX: remove: int psArgumentGet (int argc, char **argv, char *arg); // added to SDRS -//XXX: remove: int psArgumentRemove (int N, int *argc, char **argv); // added to SDRS -//XXX: remove: psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); // added to SDRS -// psVector *psGetRowVectorFromImage(psImage *image, psU32 row); // added to SDRS (as psImageRow) -//XXX: remove: int psImageCountPixelMask (psImage *mask, psU8 value); // added to SDRS - -// basic image functions -bool psImageInit (psImage *image,...); // added to SDRS (v.15) -void psImageSmooth_EAM (psImage *image, float sigma, float Nsigma); // added to SDRS (v.15) - -// psLine functions -- keep out for now? -psLine *psLineAlloc (int Nline); -bool psLineInit (psLine *line); -bool psLineAdd (psLine *line, char *format, ...); - -// not included in the .h file -- these are fairly lame, keep out? -bool p_psVectorPrint (int fd, psVector *a, char *name); -bool p_psVectorPrintRow (int fd, psVector *a, char *name); - -# endif Index: anches/rel9/psModules/src/objects/psModulesUtils.c =================================================================== --- /branches/rel9/psModules/src/objects/psModulesUtils.c (revision 5767) +++ (revision ) @@ -1,95 +1,0 @@ -# include "psModulesUtils.h" - -// XXXX: Must figure out the right names for these defines, then replace! -#define PS_META_STR 0 -#define PS_META_F32 0 - - - -// extract config informatin from config data or from image header, if specified -// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme -psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) -{ - - char *source; - char *keyword; - psF32 value; - psMetadataItem *item; - - // look for the entry in the config collection - item = psMetadataLookup (config, name); - if (item == NULL) { - psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name); - value = psMetadataLookupF32 (status, header, name); - if (*status == false) { - psAbort ("pmConfigLookupF32", "no key %s in header", name); - } - *status = true; - return (value); - } - - // I'm either expecting a string, with the name "HD:keyword"... - if (item->type == PS_META_STR) { - source = item->data.V; - if (!strncasecmp (source, "HD:", 3)) { - keyword = &source[3]; - value = psMetadataLookupF32 (status, header, keyword); - if (*status == false) { - psAbort ("pmConfigLookupF32", "no key %s in config", name); - } - *status = true; - // psFree (item); - return (value); - } - } - - // ... or a value (F32?) - if (item->type == PS_META_F32) { - value = item->data.F32; - // psFree (item); - return (value); - } - - *status = false; - psError(PS_ERR_UNKNOWN, true, "invalid item"); - return (0); -} - -bool pmModelFitStatus (pmModel *model) -{ - - bool status; - - pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); - status = statusFunc (model); - - return (status); -} - -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); -} - Index: anches/rel9/psModules/src/objects/psModulesUtils.h =================================================================== --- /branches/rel9/psModules/src/objects/psModulesUtils.h (revision 5767) +++ (revision ) @@ -1,20 +1,0 @@ - -#ifndef PS_MODULE_UTILS -#define PS_MODULE_UTILS - -#include // for strcasecmp -#include -#include "psLibUtils.h" -#include "pmObjects.h" -#include "pmModelGroup.h" - -// psModule extra utilities -// bool pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF); -bool pmModelFitStatus (pmModel *model); -int pmSourceDophotType (pmSource *source); -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask); - -// XXX: unify with paul's image/header/metadata functions -psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); - -#endif Index: /branches/rel9/psModules/test/astrom/tst_pmAstrometry01.c =================================================================== --- /branches/rel9/psModules/test/astrom/tst_pmAstrometry01.c (revision 5767) +++ /branches/rel9/psModules/test/astrom/tst_pmAstrometry01.c (revision 5768) @@ -1,18 +1,21 @@ /** @file tst_psAstrometry01.c -* -* @brief This code will test the pmFPA hierarchy transform code in psAstrometry.[ch] -* -* @author GLG, MHPCC -* -* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-05 21:28:55 $ -* -* XXX: Add tests were the coordinate does not transform to any legitimate cell -* or chip, or FPA, or whatever. -* -* XXX: For each function, add tests for bad input parameters, as well as failed transforms. -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ + * + * @brief This code will test the pmFPA hierarchy transform code in psAstrometry.[ch] + * + * @author GLG, MHPCC + * + * @version $Revision: 1.2.2.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-12 21:52:24 $ + * + * XXX: Add tests were the coordinate does not transform to any legitimate cell + * or chip, or FPA, or whatever. + * + * XXX: For each function, add tests for bad input parameters, as well as failed transforms. + * + * XXX: Must test pmFPASelectChip() and pmFPAExcludeChip(). + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + #include "config.h" #include @@ -23,8 +26,10 @@ static psS32 test3( void ); static psS32 test4( void ); +static psS32 test5( void ); testDescription tests[] = { {test3, -3, "pmAstrometry focal plane transformations", 0, false}, {test4, -3, "pmCheckParents()", 0, false}, + {test5, -3, "pmFPASelectChip() and pmFPAExcludeChip()", 0, false}, {NULL} }; @@ -584,2 +589,83 @@ return(testStatus); } + +/****************************************************************************** +test5(): This routine wil test the pmFPASelectChip() and pmFPAExcludeChip() +functions. We generate an pmFPA hierarchy, then set the ->valid members with +those routines, then verify. + *****************************************************************************/ +psS32 test5( void ) +{ + psS32 testStatus = 0; + pmChip *tmpChip = NULL; + + // + // Generate a pmFPA hierarchy. + // + pmFPA *tmpFPA = genSystem(); + + // + // We test the ->valid member for each chip. + // + for (psS32 i = 0 ; i < tmpFPA->chips->n ; i++) { + tmpChip = (pmChip *) tmpFPA->chips->data[i]; + if ((tmpChip == NULL) || (tmpChip->valid != false)) { + printf("TEST ERROR: Could not properly generate an FPA hierarchy.\n"); + testStatus = 1; + } + } + + // + // Exclude chip number 0, include all others, then test return value + // + psS32 numChips = pmFPAExcludeChip(tmpFPA, 0); + if (numChips != (NUM_CHIPS-1)) { + printf("TEST ERROR: pmFPAExcludeChip() did not return the correct number of chips.\n"); + testStatus = 2; + } + + // + // We test the ->valid member for each chip. + // + tmpChip = (pmChip *) tmpFPA->chips->data[0]; + if (tmpChip->valid != false) { + printf("TEST ERROR: pmFPAExcludeChip() did not set the proper chip->valid to FALSE.\n"); + testStatus = 3; + } + for (psS32 i = 1 ; i < tmpFPA->chips->n ; i++) { + pmChip *tmpChip = (pmChip *) tmpFPA->chips->data[i]; + if (tmpChip->valid != true) { + printf("TEST ERROR: pmFPAExcludeChip() did not set the proper chip->valids to FALSE.\n"); + testStatus = 4; + } + } + + + // + // Include chip number 0, exclude all others, then test return value + // + psBool tmpBool = pmFPASelectChip(tmpFPA, 0); + if (tmpBool != true) { + printf("TEST ERROR: pmFPASelectChip() returned FALSE.\n"); + testStatus = 5; + } + + // + // We test the ->valid member for each chip. + // + tmpChip = (pmChip *) tmpFPA->chips->data[0]; + if (tmpChip->valid != true) { + printf("TEST ERROR: pmFPASelectChip() did not set the proper chip->valid to FALSE.\n"); + testStatus = 6; + } + for (psS32 i = 1 ; i < tmpFPA->chips->n ; i++) { + pmChip *tmpChip = (pmChip *) tmpFPA->chips->data[i]; + if (tmpChip->valid != false) { + printf("TEST ERROR: pmFPASelectChip() did not set the proper chip->valids to FALSE.\n"); + testStatus = 7; + } + } + + psFree(tmpFPA); + return(testStatus); +}