Index: /branches/eam_branches/ipp-20100621/psModules/src/config/pmConfig.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/config/pmConfig.c (revision 28434) +++ /branches/eam_branches/ipp-20100621/psModules/src/config/pmConfig.c (revision 28435) @@ -32,4 +32,5 @@ #include "pmConfig.h" +#include "pmVisualUtils.h" #ifdef HAVE_NEBCLIENT @@ -638,4 +639,23 @@ psArgumentVerbosity(argc, argv); // XXX: substitute the string for the default log level for "2". + } + + // Set the visualization levels + // argument format is: -visual (facil) (level) + while ((argNum = psArgumentGet(*argc, argv, "-visual"))) { + if ( (*argc < argNum + 3) ) { + psError(PS_ERR_IO, true, "-visual switch specified without facility and level."); + return NULL; + } + psArgumentRemove(argNum, argc, argv); + pmVisualSetLevel(argv[argNum], atoi(argv[argNum+1])); + psArgumentRemove(argNum, argc, argv); + psArgumentRemove(argNum, argc, argv); + } + if ((argNum = psArgumentGet(*argc, argv, "-visual-all"))) { + pmVisualSetLevel(".", 10); + } + if ((argNum = psArgumentGet(*argc, argv, "-visual-levels"))) { + pmVisualPrintLevels(stdout); } Index: /branches/eam_branches/ipp-20100621/psModules/src/extras/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/extras/Makefile.am (revision 28434) +++ /branches/eam_branches/ipp-20100621/psModules/src/extras/Makefile.am (revision 28435) @@ -9,4 +9,5 @@ pmKapaPlots.c \ pmVisual.c \ + pmVisualUtils.c \ ippStages.c @@ -17,4 +18,5 @@ pmKapaPlots.h \ pmVisual.h \ + pmVisualUtils.h \ ippDiffMode.h \ ippStages.h Index: /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisualUtils.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisualUtils.c (revision 28435) +++ /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisualUtils.c (revision 28435) @@ -0,0 +1,505 @@ +/** These utility functions manage the visual level information (equivalent to pmVisual levels) + * @author Eugene Magnier, IfA + * @date June 18, 2010 + */ + +/* Include Files */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +#include + +#define PM_UNKNOWN_VISUAL_LEVEL -9999 // we don't know this name's level +#define PM_DEFAULT_VISUAL_LEVEL -1 +#define PM_THE_OTHER_DEFAULT_VISUAL_LEVEL 0 // ??? +#define MAX_COMPONENT_LENGTH 1024 + +/** Basic structure for the component tree. A component is a string of the + form aaa.bbb.ccc, and may itself contain further subcomponents. The + Component structure doesn't in fact contain it's full name, but only the + last part. */ + +static pmVisualComponent* cRoot = NULL; // The root of the visual component tree + +/** static function prototypes **/ +static void componentFree(pmVisualComponent* comp); +static pmVisualComponent* componentAlloc(const char *name, int level); +static int getLevel(const char *facil); +static void initVisualLevels(void); +static void componentFree(pmVisualComponent* comp); +static pmVisualComponent* componentAlloc(const char *name, int level); +static void doGetVisualLevels(psMetadata *out, const pmVisualComponent* comp, psString parent, int defLevel); +static void doPrintVisualLevels(const pmVisualComponent* comp, FILE *output, psS32 depth, psS32 defLevel); +static psS32 doGetVisualLevel(const char *aname); +static bool componentAdd(const char *addNodeName, psS32 level); + +/***************************************************************************** +Set all visual levels at or below the specified node to zero. + *****************************************************************************/ +void pmVisualReset(pmVisualComponent* currentNode) +{ + psAssert(currentNode, "impossible"); + + psS32 i = 0; + + if (NULL == currentNode) { + return; + } + + currentNode->level = 0; + for (i = 0; i < currentNode->n; i++) { + if (!currentNode->subcomp[i]) { + psLogMsg("pmVisualReset", PS_LOG_WARN, _("Sub-component %d of node %s in visual tree is NULL."), i, currentNode->name); + } else { + pmVisualReset(currentNode->subcomp[i]); + } + } + return; +} + +/***************************************************************************** +Free all visual levels + *****************************************************************************/ +void pmVisualCleanup() +{ + psFree(cRoot); + cRoot = NULL; +} + +// psSetVisualLevel(): add the component named "comp" to the component tree, +int pmVisualSetLevel(const char *comp, // component of interest + int level) // desired visual level +{ + PS_ASSERT_STRING_NON_EMPTY(comp, 0); + + char *compName = NULL; + int prevLevel = -1; + + // If the root component tree does not exist, then initialize it. + if (cRoot == NULL) { + initVisualLevels(); + } + + // If the component name has no leading dot, then supply it. + if (comp[0] != '.') { + compName = (char *) psAlloc(10 + strlen(comp)); + strcpy(compName, "."); + compName = strcat(compName, comp); + } else { + compName = (char *) comp; + } + prevLevel = getLevel(compName); + // Add the new component to the component tree. + if ( !componentAdd(compName, level) ) { + psError(PS_ERR_UNKNOWN, false, + _("Failed to set visual level (%d) to '%s'."), + level, + compName); + + if (comp[0] != '.') { + psFree(compName); + } + // return false; + return -1; + } + + if (comp[0] != '.') { + psFree(compName); + } + + // return true; + return prevLevel; +} + +// Append the function name to the facility +// NB: declares TARGET! +#define FACILITY(TARGET, FUNC, FACIL) \ + size_t _facilLength = strlen(FACIL); /* Length of facility name */ \ + size_t _funcLength = strlen(FUNC); /* Length of function name */ \ + char TARGET[_facilLength + _funcLength + 2]; /* facility + the function name */ \ + strcpy(&TARGET[0], FACIL); \ + TARGET[_facilLength] = '.'; \ + strcpy(&TARGET[_facilLength + 1], FUNC); + +int p_pmVisualGetLevel(const char *func, const char *name) +{ + PS_ASSERT_STRING_NON_EMPTY(name, 0); + PS_ASSERT_STRING_NON_EMPTY(func, 0); + + FACILITY(facility, func, name); + return getLevel(facility); +} + +/** the main tool in the pmVisual system : check that the specified level matches the desired level **/ +bool p_pmVisualTestLevel(const char *func, const char *name, int level) +{ + PS_ASSERT_STRING_NON_EMPTY(name, false); + PS_ASSERT_STRING_NON_EMPTY(func, false); + + // return true if level is set to be shown + FACILITY(facility, func, name); + bool valid = level <= getLevel(facility); + + if (valid) { + psLogMsg ("psModules.visual", PS_LOG_DETAIL, "visualization %s (%d)\n", facility, level); + } + + return (valid); +} + +// psPrintVisualLevels(): Simply print all the visual levels in the visual level +void pmVisualPrintLevels(FILE *output) +{ + if (cRoot == NULL) { + return; + } + + doPrintVisualLevels(cRoot, output, 0, PM_THE_OTHER_DEFAULT_VISUAL_LEVEL); +} + +// generate a metadata of the visual levels (is this really needed?) +psMetadata *pmVisualLevels(void) +{ + if (cRoot == NULL) { + return psMetadataAlloc(); + } + + psMetadata *out = psMetadataAlloc();// Output metadata with the levels + doGetVisualLevels(out, cRoot, NULL, PM_THE_OTHER_DEFAULT_VISUAL_LEVEL); + + return out; +} + +/****************************** static functions ******************************/ + +/***************************************************************************** +componentAdd(): Adds the component named "addNodeName" to the root tree. + *****************************************************************************/ +static bool componentAdd(const char *addNodeName, psS32 level) +{ + psAssert(addNodeName, "impossible"); + + psS32 i = 0; // Loop index variable. + char name[MAX_COMPONENT_LENGTH]; // buffer for writeable copy. + char *firstComponent = NULL; // first component of name + pmVisualComponent* currentNode = cRoot; + psS32 nodeExists = 0; + + // XXX: Verify that this is the correct behavior. + if (strcmp("", addNodeName) == 0) { + psError(PS_ERR_BAD_PARAMETER_NULL,true, + _("Failed to add null component to visual tree.")); + return false; + } + + // Is this the root node? If so, simply set level and return. + if (strcmp(".", addNodeName) == 0) { + cRoot->level = level; + return true; + } + + if (addNodeName[0] != '.') { + psError(PS_ERR_BAD_PARAMETER_VALUE,true, + _("Failed to add '%s' to the root component tree; component must start with '.'."), + addNodeName); + return false; + } + + strncpy(name, addNodeName, MAX_COMPONENT_LENGTH); + char *pname = name+1; // Take off the period + // Iterate through the components of addNodeName. Strip off the first + // component of the name, find that in the root tree, or add it if it + // does not exist, then move to the next component in the name. + + while (pname != NULL) { + firstComponent = pname; + pname = strchr(firstComponent, '.'); + if (pname != NULL) { + *pname = '\0'; + pname++; + } + nodeExists = 0; + for (i = 0; i < currentNode->n; i++) { + if (strcmp(currentNode->subcomp[i]->name, firstComponent) == 0) { + currentNode = currentNode->subcomp[i]; + nodeExists = 1; + if (pname == NULL) { + currentNode->level = level; + } + } + } + + if (nodeExists == 0) { + currentNode->subcomp = psRealloc(currentNode->subcomp, + (currentNode->n + 1) * sizeof(pmVisualComponent* )); + psMemSetPersistent(currentNode->subcomp,true); + + currentNode->n = (currentNode->n) + 1; + + if (pname == NULL) { + // This is the final component to add. + currentNode->subcomp[(currentNode->n) - 1] = + componentAlloc(firstComponent, level); + } else { + // We are adding an intermediate component. The visual level is not defined. + // An undefined visual level inherits the visual level of it's parent. However, + // we do not set that specifically here since that would result in a static, + // one-time, type of behavior. + + currentNode->subcomp[(currentNode->n) - 1] = + componentAlloc(firstComponent, PM_DEFAULT_VISUAL_LEVEL); + } + currentNode = currentNode->subcomp[(currentNode->n) - 1]; + } + } + + return true; +} + +/***************************************************************************** + doGetVisualLevel() + This function recursively searches the root component tree for the + component named "name", which is supplied by a parameter. If it + finds that component, it returns the level of that component. + Otherwise, it returns ???. + + NOTE: We modified this so that the user may omit the leading "," in a + component name. Since the code was already implemented assuming the "." + was required, rather than change all that code, in this function, I + simply add a leading "." to the component name if there is none. + + Inputs: + name: + Outputs: + none + Returns: + The visual level of the "name" component. + *****************************************************************************/ +static psS32 doGetVisualLevel(const char *aname) +{ + psAssert(aname, "impossible"); + char name[strlen(aname) + 1]; // need a writeable copy: for strsep() + char *pname = name; + char *firstComponent = NULL; // first component of name + pmVisualComponent* currentNode = cRoot; + psS32 i = 0; + psS32 defaultLevel = 0; + + if (NULL == currentNode) { + return (PM_UNKNOWN_VISUAL_LEVEL); + } + + if (strcmp(".", aname) == 0) { + return (cRoot->level); + } + + if (aname[0] != '.') { + return (PM_UNKNOWN_VISUAL_LEVEL); + } + + defaultLevel = cRoot->level; + strcpy(name, aname); + pname = name+1; + while (pname != NULL) { + firstComponent = pname; + pname = strchr(firstComponent, '.'); + if (pname != NULL) { + *pname = '\0'; + pname++; + } + for (i = 0; i < currentNode->n; i++) { + if (NULL == currentNode->subcomp[i]) { + psLogMsg("p_pmVisualReset", PS_LOG_WARN, + _("Sub-component %d of node %s in visual tree is NULL."), + i, currentNode->name); + } + + if (strcmp(currentNode->subcomp[i]->name, firstComponent) == 0) { + currentNode = currentNode->subcomp[i]; + // For level inheritance purpose, we save the level of this + // component if it is not DEFAULT. + if (currentNode->level != PM_DEFAULT_VISUAL_LEVEL) { + defaultLevel = currentNode->level; + } + // Determine if this is the last component: + if (pname == NULL) { + if (currentNode->level != PM_DEFAULT_VISUAL_LEVEL) { + return (currentNode->level); + } else { + return(defaultLevel); + } + } + } + } + } + return(defaultLevel); +} + +/***************************************************************************** + getLevel() + Return a visual level of "name" in the root component tree. If the + exact string of components in "name" does not exist in the root + tree, we return the deepest level of the match. + *****************************************************************************/ +static int getLevel(const char *name) +{ + if (cRoot == NULL) { + return (PM_UNKNOWN_VISUAL_LEVEL); + } + + psS32 visualLevel; + + // If the component name has no leading dot, then supply it. + if (name[0] != '.') { + char compName[strlen(name) + 2]; + compName[0] = '.'; + strcpy(&compName[1], name); + + visualLevel = doGetVisualLevel(compName); + } else { + // Search the component root tree, determine the visual level. + visualLevel = doGetVisualLevel(name); + } + + // XXX: The default visual level is currently set at -1, which is not a + // valid visual level. This is convenient in determining whether or not + // a component should inherit the visual level from parent nodes. However, + // it's not clear that -1 should ever be returned by this function. + // The SDR is unclear on this point and we should probably request IfA + // comment. + if (visualLevel == PM_DEFAULT_VISUAL_LEVEL) { + visualLevel = PM_THE_OTHER_DEFAULT_VISUAL_LEVEL; + } + + return(visualLevel); +} + +/***************************************************************************** + doPrintVisualLevels() + This function recursively searches the component tree supplied by the + parameter "comp" and prints the name and level of each component. + Inputs: + comp: a node in the component tree. + level: the level of that node + Outputs: + none + Returns: + null + *****************************************************************************/ +static void doPrintVisualLevels(const pmVisualComponent* comp, + FILE *output, + psS32 depth, + psS32 defLevel) +{ + psAssert(comp, "impossible"); + + char line[1024]; + psS32 i = 0; + + if (comp->name[0] == '\0') { + return; + } else { + if (comp->level == PM_DEFAULT_VISUAL_LEVEL) { + sprintf(line,"%*s%-*s %d\n", depth, "", 20 - depth, comp->name, defLevel); + fwrite (line, 1, strlen(line), output); + } else { + sprintf(line, "%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); + fwrite (line, 1, strlen(line), output); + } + } + + for (i = 0; i < comp->n; i++) { + if (comp->level == PM_DEFAULT_VISUAL_LEVEL) { + doPrintVisualLevels(comp->subcomp[i], output, depth + 1, defLevel); + } else { + doPrintVisualLevels(comp->subcomp[i], output, depth + 1, comp->level); + } + } +} + +static void doGetVisualLevels(psMetadata *out, // Output metadata with the visual levels + const pmVisualComponent* comp, // Component to add + psString parent, // Name of parent level + int defLevel // Default level + ) +{ + if (comp->name[0] == '\0') { + return; + } + + psString name = psStringCopy(parent); // Name of this level + if (comp->name[0] == '.') { + psStringAppend(&name, "%s", comp->name + 1); + } else if (!parent) { + psStringAppend(&name, "%s", comp->name); + } else { + psStringAppend(&name, ".%s", comp->name); + } + + int level = (comp->level == PM_DEFAULT_VISUAL_LEVEL) ? defLevel : comp->level; // Level for component + if (name) { + psMetadataAddS32(out, PS_LIST_TAIL, name, 0, NULL, level); + } + for (int i = 0; i < comp->n; i++) { + doGetVisualLevels(out, comp->subcomp[i], name, level); + } + + psFree(name); + + return; +} + +/***************************************************************************** +componentAlloc(): allocate memory for a new node, and initialize members. + *****************************************************************************/ +static pmVisualComponent* componentAlloc(const char *name, int level) +{ + psAssert(name, "impossible"); + + pmVisualComponent* comp = psAlloc(sizeof(pmVisualComponent)); + psMemSetDeallocator(comp, (psFreeFunc) componentFree); + + comp->name = psStringCopy(name); + comp->level = level; + comp->n = 0; + comp->specified = false; + comp->subcomp = NULL; + return comp; +} + +/***************************************************************************** +componentFree(): free the current node in the root tree, and all children +nodes as well. + *****************************************************************************/ +static void componentFree(pmVisualComponent* comp) +{ + if (comp == NULL) { + return; + } + + if (comp->subcomp != NULL) { + for (psS32 i = 0; i < comp->n; i++) { + psFree(comp->subcomp[i]); + } + psFree(comp->subcomp); + } + psFree(comp->name); +} + +/***************************************************************************** +initVisualLevels(): simply initialize the component root tree. +*****************************************************************************/ +static void initVisualLevels(void) +{ + if (cRoot == NULL) { + cRoot = componentAlloc(".", PM_DEFAULT_VISUAL_LEVEL); + } +} Index: /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisualUtils.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisualUtils.h (revision 28435) +++ /branches/eam_branches/ipp-20100621/psModules/src/extras/pmVisualUtils.h (revision 28435) @@ -0,0 +1,33 @@ +/* @file pmVisualUtils.h + * @brief functions to create visual diagnostics with the help of 'kapa' + * @author Chris Beaumont, IfA + * + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_VISUAL_UTILS_H +#define PM_VISUAL_UTILS_H + +typedef struct pmVisualComponent +{ + const char *name; ///< last part of name of component + psS32 level; ///< visual level for this component + bool specified; ///< whether the component is specified + psS32 n; ///< number of subcomponents + struct pmVisualComponent* *subcomp; ///< next level of subcomponents +} +pmVisualComponent; + +psMetadata *pmVisualLevels(void); +void pmVisualPrintLevels(FILE *output); +int pmVisualSetLevel(const char *comp, int level); +void pmVisualCleanup(); +void pmVisualReset(pmVisualComponent* currentNode); + +int p_pmVisualGetLevel(const char *func, const char *name); +# define pmVisualGetLevel(facil) p_pmVisualGetLevel(__func__, facil) + +bool p_pmVisualTestLevel(const char *func, const char *name, int level); +# define pmVisualTestLevel(facil, level) p_pmVisualTestLevel(__func__, facil, level) + +#endif // PM_VISUAL_UTILS_H Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceVisual.c (revision 28434) +++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceVisual.c (revision 28435) @@ -13,4 +13,5 @@ #include #include "pmVisual.h" +#include "pmVisualUtils.h" // functions used to visualize the analysis as it goes @@ -34,5 +35,5 @@ Graphdata graphdata; - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.metric", 2)) return true; if (kapa1 == -1) { @@ -118,5 +119,5 @@ Graphdata graphdata; - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.subpix", 3)) return true; if (kapa1 == -1) { @@ -280,5 +281,5 @@ bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.fits", 2)) return true; if (kapa2 == -1) { @@ -360,5 +361,6 @@ bool pmSourceVisualShowModelFit (pmSource *source) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.fitresid", 2)) return true; + if (!source->pixels) return false; if (!source->modelFlux) return false; @@ -404,5 +406,6 @@ Graphdata graphdata; - if (!pmVisualIsVisual() || !plotPSF) return true; + if (!plotPSF) return true; + if (!pmVisualTestLevel("psphot.psf.resid", 2)) return true; if (kapa1 == -1) { Index: /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h =================================================================== --- /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h (revision 28434) +++ /branches/eam_branches/ipp-20100621/psModules/src/psmodules.h (revision 28435) @@ -10,4 +10,5 @@ #include #include +#include #include #include Index: /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am (revision 28434) +++ /branches/eam_branches/ipp-20100621/psphot/src/Makefile.am (revision 28435) @@ -20,5 +20,5 @@ -$(RM) psphotVersionDefinitions.h $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h -FORCE: ; +# FORCE: ; libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotArguments.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotArguments.c (revision 28434) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotArguments.c (revision 28435) @@ -176,10 +176,4 @@ } - // visual : interactive display mode - if ((N = psArgumentGet (argc, argv, "-visual"))) { - psArgumentRemove (N, &argc, argv); - pmVisualSetVisual(true); - } - // break : used from recipe throughout psphotReadout if ((N = psArgumentGet (argc, argv, "-break"))) { Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotPetrosianVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotPetrosianVisual.c (revision 28434) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotPetrosianVisual.c (revision 28435) @@ -1,4 +1,3 @@ # include "psphotInternal.h" -# define FORCE_VISUAL 0 // this function displays representative images as the psphot analysis progresses: @@ -54,6 +53,5 @@ Graphdata graphdata; - // return true; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.byangle", 2)) return true; if (kapa2 == -1) { @@ -101,6 +99,5 @@ Graphdata graphdata; - // return true; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.radii", 2)) return true; if (kapa == -1) { @@ -173,5 +170,5 @@ KapaSection section; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.stats", 2)) return true; if (kapa2 == -1) { @@ -311,5 +308,5 @@ Graphdata graphdata; - if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.petro.ellipse", 2)) return true; if (kapa == -1) { Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c (revision 28434) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotSourceStats.c (revision 28435) @@ -501,5 +501,7 @@ psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - psphotVisualPlotMoments (recipe, analysis, sources); + if (pmVisualTestLevel("psphot.moments.full", 2)) { + psphotVisualPlotMoments (recipe, analysis, sources); + } Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; Index: /branches/eam_branches/ipp-20100621/psphot/src/psphotVisual.c =================================================================== --- /branches/eam_branches/ipp-20100621/psphot/src/psphotVisual.c (revision 28434) +++ /branches/eam_branches/ipp-20100621/psphot/src/psphotVisual.c (revision 28435) @@ -66,9 +66,9 @@ int myKapa = psphotKapaChannel (channel); if (!(strcasecmp (overlay, "all"))) { - KiiEraseOverlay (myKapa, "red"); - KiiEraseOverlay (myKapa, "green"); - KiiEraseOverlay (myKapa, "blue"); - KiiEraseOverlay (myKapa, "yellow"); - return true; + KiiEraseOverlay (myKapa, "red"); + KiiEraseOverlay (myKapa, "green"); + KiiEraseOverlay (myKapa, "blue"); + KiiEraseOverlay (myKapa, "yellow"); + return true; } KiiEraseOverlay (myKapa, overlay); @@ -182,5 +182,5 @@ bool psphotVisualShowImage (pmReadout *readout) { - if (!pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image", 1)) return true; int kapa = psphotKapaChannel (1); @@ -199,5 +199,5 @@ pmReadout *backgnd; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image.backgnd", 2)) return true; int kapa = psphotKapaChannel (1); @@ -208,7 +208,7 @@ if (file->mode == PM_FPA_MODE_INTERNAL) { - backgnd = file->readout; + backgnd = file->readout; } else { - backgnd = pmFPAviewThisReadout (view, file->fpa); + backgnd = pmFPAviewThisReadout (view, file->fpa); } @@ -222,5 +222,5 @@ bool psphotVisualShowSignificance (psImage *image, float min, float max) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image.signif", 2)) return true; int kapa = psphotKapaChannel (1); @@ -233,4 +233,5 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowPeaks (pmDetections *detections) { @@ -238,5 +239,5 @@ KiiOverlay *overlay; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.peaks", 1)) return true; int kapa = psphotKapaChannel (1); @@ -252,15 +253,15 @@ for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - if (peak == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_BOX; - overlay[Noverlay].x = peak->xf; - overlay[Noverlay].y = peak->yf; - overlay[Noverlay].dx = 2.0; - overlay[Noverlay].dy = 2.0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmPeak *peak = peaks->data[i]; + if (peak == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 2.0; + overlay[Noverlay].dy = 2.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -272,4 +273,5 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowFootprints (pmDetections *detections) { @@ -277,5 +279,5 @@ KiiOverlay *overlay; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.footprints", 3)) return true; int kapa = psphotKapaChannel (1); @@ -292,82 +294,82 @@ for (int i = 0; i < footprints->n; i++) { - pmSpan *span = NULL; - - pmFootprint *footprint = footprints->data[i]; - if (footprint == NULL) continue; - if (footprint->spans == NULL) continue; - if (footprint->spans->n < 1) continue; - - // draw the top - // XXX need to allow top (and bottom) to have more than one span - span = footprint->spans->data[0]; - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = span->x0; - overlay[Noverlay].y = span->y; - overlay[Noverlay].dx = span->x1 - span->x0; - overlay[Noverlay].dy = 0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - int ys = span->y; - int x0s = span->x0; - int x1s = span->x1; - - // draw the outer span edges - for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int ye = span1->y; - int x0e = span1->x0; - int x1e = span1->x1; - - // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) - // find all of the spans in this row and generate x0e, x01: - for (int k = j + 1; k < footprint->spans->n; k++) { - pmSpan *span2 = footprint->spans->data[k]; - if (span2->y > span1->y) break; - x0e = PS_MIN (x0e, span2->x0); - x1e = PS_MAX (x1e, span2->x1); - j++; - } - - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = x0s; - overlay[Noverlay].y = ys; - overlay[Noverlay].dx = x0e - x0s; - overlay[Noverlay].dy = ye - ys; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = x1s; - overlay[Noverlay].y = ys; - overlay[Noverlay].dx = x1e - x1s; - overlay[Noverlay].dy = ye - ys; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - ys = ye; - x0s = x0e; - x1s = x1e; - } - - // draw the bottom - span = footprint->spans->data[footprint->spans->n - 1]; - overlay[Noverlay].type = KII_OVERLAY_LINE; - overlay[Noverlay].x = span->x0; - overlay[Noverlay].y = span->y; - overlay[Noverlay].dx = span->x1 - span->x0; - overlay[Noverlay].dy = 0; - overlay[Noverlay].angle = 0.0; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + pmSpan *span = NULL; + + pmFootprint *footprint = footprints->data[i]; + if (footprint == NULL) continue; + if (footprint->spans == NULL) continue; + if (footprint->spans->n < 1) continue; + + // draw the top + // XXX need to allow top (and bottom) to have more than one span + span = footprint->spans->data[0]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + int ys = span->y; + int x0s = span->x0; + int x1s = span->x1; + + // draw the outer span edges + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int ye = span1->y; + int x0e = span1->x0; + int x1e = span1->x1; + + // we cannot have two discontinuous spans on the top or bottom, right? (no, probably not right) + // find all of the spans in this row and generate x0e, x01: + for (int k = j + 1; k < footprint->spans->n; k++) { + pmSpan *span2 = footprint->spans->data[k]; + if (span2->y > span1->y) break; + x0e = PS_MIN (x0e, span2->x0); + x1e = PS_MAX (x1e, span2->x1); + j++; + } + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x0s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x0e - x0s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = x1s; + overlay[Noverlay].y = ys; + overlay[Noverlay].dx = x1e - x1s; + overlay[Noverlay].dy = ye - ys; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + ys = ye; + x0s = x0e; + x1s = x1e; + } + + // draw the bottom + span = footprint->spans->data[footprint->spans->n - 1]; + overlay[Noverlay].type = KII_OVERLAY_LINE; + overlay[Noverlay].x = span->x0; + overlay[Noverlay].y = span->y; + overlay[Noverlay].dx = span->x1 - span->x0; + overlay[Noverlay].dy = 0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); } @@ -379,4 +381,5 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowMoments (psArray *sources) { @@ -387,5 +390,5 @@ psEllipseAxes axes; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.moments", 2)) return true; int kapa = psphotKapaChannel (1); @@ -401,27 +404,27 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - pmMoments *moments = source->moments; - if (moments == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = moments->Mx; - overlay[Noverlay].y = moments->My; - - emoments.x2 = moments->Mxx; - emoments.xy = moments->Mxy; - emoments.y2 = moments->Myy; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - - overlay[Noverlay].dx = 2.0*axes.major; - overlay[Noverlay].dy = 2.0*axes.minor; - - overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; - - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->Mx; + overlay[Noverlay].y = moments->My; + + emoments.x2 = moments->Mxx; + emoments.xy = moments->Mxy; + emoments.y2 = moments->Myy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -439,5 +442,5 @@ KapaSection section; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.moments", 1)) return true; int myKapa = psphotKapaChannel (2); @@ -456,26 +459,26 @@ float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); - for (int n = 0; n < nRegions; n++) { - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); - float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); - float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - - float X0 = psfX - 4.0*psfdX; - float X1 = psfX + 4.0*psfdX; - float Y0 = psfY - 4.0*psfdY; - float Y1 = psfY + 4.0*psfdY; - - if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } - if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } - if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } - if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } - } + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + + float X0 = psfX - 4.0*psfdX; + float X1 = psfX + 4.0*psfdX; + float Y0 = psfY - 4.0*psfdY; + float Y1 = psfY + 4.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } } Xmin = PS_MAX(Xmin, -0.1); @@ -498,21 +501,21 @@ int nF = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->moments == NULL) - continue; - - xFaint->data.F32[nF] = source->moments->Mxx; - yFaint->data.F32[nF] = source->moments->Myy; - mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum); - nF++; - - // XXX make this a user-defined cutoff - if (source->moments->SN < SN_LIM) - continue; - - xBright->data.F32[nB] = source->moments->Mxx; - yBright->data.F32[nB] = source->moments->Myy; - mBright->data.F32[nB] = -2.5*log10(source->moments->Sum); - nB++; + pmSource *source = sources->data[i]; + if (source->moments == NULL) + continue; + + xFaint->data.F32[nF] = source->moments->Mxx; + yFaint->data.F32[nF] = source->moments->Myy; + mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum); + nF++; + + // XXX make this a user-defined cutoff + if (source->moments->SN < SN_LIM) + continue; + + xBright->data.F32[nB] = source->moments->Mxx; + yBright->data.F32[nB] = source->moments->Myy; + mBright->data.F32[nB] = -2.5*log10(source->moments->Sum); + nB++; } xFaint->n = nF; @@ -652,45 +655,45 @@ // draw N circles to outline the clumps { - KapaSelectSection (myKapa, "MxxMyy"); - - // draw a circle centered on psfX,Y with size of the psf limit - psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); - psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); - - int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); - float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); - - graphdata.color = KapaColorByName ("blue"); - graphdata.style = 0; - - graphdata.xmin = Xmin; - graphdata.ymin = Ymin; - graphdata.xmax = Xmax; - graphdata.ymax = Ymax; - KapaSetLimits (myKapa, &graphdata); - - for (int n = 0; n < nRegions; n++) { - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); - float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); - float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - float Rx = psfdX * PSF_CLUMP_NSIGMA; - float Ry = psfdY * PSF_CLUMP_NSIGMA; - - for (int i = 0; i < xLimit->n; i++) { - xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; - yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; - } - KapaPrepPlot (myKapa, xLimit->n, &graphdata); - KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); - KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); - } - psFree (xLimit); - psFree (yLimit); + KapaSelectSection (myKapa, "MxxMyy"); + + // draw a circle centered on psfX,Y with size of the psf limit + psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); + psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); + + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 0; + + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + float Rx = psfdX * PSF_CLUMP_NSIGMA; + float Ry = psfdY * PSF_CLUMP_NSIGMA; + + for (int i = 0; i < xLimit->n; i++) { + xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; + yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; + } + KapaPrepPlot (myKapa, xLimit->n, &graphdata); + KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); + } + psFree (xLimit); + psFree (yLimit); } @@ -721,28 +724,28 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - if (source->type != type) continue; - if (mode && !(source->mode & mode)) continue; - - pmMoments *moments = source->moments; - if (moments == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = moments->Mx; - overlay[Noverlay].y = moments->My; - - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - - overlay[Noverlay].dx = 2.0*axes.major; - overlay[Noverlay].dy = 2.0*axes.minor; - overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + if (source->type != type) continue; + if (mode && !(source->mode & mode)) continue; + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->Mx; + overlay[Noverlay].y = moments->My; + + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = 2.0*axes.major; + overlay[Noverlay].dy = 2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -753,7 +756,8 @@ } +// XXX : requires psphotVisualShowImage bool psphotVisualShowRoughClass (psArray *sources) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.size", 3)) return true; int myKapa = psphotKapaChannel (1); @@ -776,5 +780,5 @@ bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.model", 1)) return true; int myKapa = psphotKapaChannel (3); @@ -797,28 +801,28 @@ // generate a fake model at each of the 3x3 image grid positions for (int x = -2; x <= +2; x ++) { - for (int y = -2; y <= +2; y ++) { - // use the center of the center pixel of the image - float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5; - float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5; - - // assign the x and y coords to the image center - // create an object with center intensity of 1000 - modelRef->params->data.F32[PM_PAR_SKY] = 0; - modelRef->params->data.F32[PM_PAR_I0] = 1000; - modelRef->params->data.F32[PM_PAR_XPOS] = xc; - modelRef->params->data.F32[PM_PAR_YPOS] = yc; - - // create modelPSF from this model - pmModel *model = pmModelFromPSF (modelRef, psf); - if (!model) continue; - - // place the reference object in the image center - // no need to mask the source here - // XXX should we measure this for the analytical model only or the full model? - pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); - psFree (model); - } + for (int y = -2; y <= +2; y ++) { + // use the center of the center pixel of the image + float xc = (int)((0.5 + 0.225*x)*readout->image->numCols) + readout->image->col0 + 0.5; + float yc = (int)((0.5 + 0.225*y)*readout->image->numRows) + readout->image->row0 + 0.5; + + // assign the x and y coords to the image center + // create an object with center intensity of 1000 + modelRef->params->data.F32[PM_PAR_SKY] = 0; + modelRef->params->data.F32[PM_PAR_I0] = 1000; + modelRef->params->data.F32[PM_PAR_XPOS] = xc; + modelRef->params->data.F32[PM_PAR_YPOS] = yc; + + // create modelPSF from this model + pmModel *model = pmModelFromPSF (modelRef, psf); + if (!model) continue; + + // place the reference object in the image center + // no need to mask the source here + // XXX should we measure this for the analytical model only or the full model? + pmModelAddWithOffset (psfMosaic, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + pmModelAddWithOffset (funMosaic, NULL, model, PM_MODEL_OP_FUNC | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + pmModelAddWithOffset (resMosaic, NULL, model, PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1 | PM_MODEL_OP_CENTER, 0, x*DX, y*DY); + psFree (model); + } } @@ -842,5 +846,5 @@ bool status; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.stars", 2)) return true; int myKapa = psphotKapaChannel (3); @@ -871,32 +875,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - bool keep = false; - keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); - if (!keep) continue; - - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows - - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + bool keep = false; + keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); + if (!keep) continue; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } NY += DY; @@ -918,59 +922,59 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - bool keep = false; - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - nPSF ++; - keep = true; - } - if (!keep) continue; - - if (Xo + DX > NX) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - Xo = DX; - dY = DY; - } - } else { - // extend this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); - if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - Xo += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + bool keep = false; + if (source->mode & PM_SOURCE_MODE_PSFSTAR) { + nPSF ++; + keep = true; + } + if (!keep) continue; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + Xo = DX; + dY = DY; + } + } else { + // extend this row + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY, true); + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY, true); + if (!subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + Xo += DX; + dY = PS_MAX (dY, DY); + } } @@ -994,5 +998,5 @@ bool status; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.psf.sat", 3)) return true; int myKapa = psphotKapaChannel (3); @@ -1023,32 +1027,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // only show "real" saturated stars (not defects) - if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; - if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; - - // how does this subimage get placed into the output image? - // DX = source->pixels->numCols - // DY = source->pixels->numRows - - if (dX + DX > NX) { - // too wide for the rest of this row - if (dX == 0) { - // alone on this row - NY += DY; - dX = 0; - dY = 0; - } else { - // start the next row - NY += dY; - dX = DX; - dY = DY; - } - } else { - // extend this row - dX += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + // only show "real" saturated stars (not defects) + if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; + + // how does this subimage get placed into the output image? + // DX = source->pixels->numCols + // DY = source->pixels->numRows + + if (dX + DX > NX) { + // too wide for the rest of this row + if (dX == 0) { + // alone on this row + NY += DY; + dX = 0; + dY = 0; + } else { + // start the next row + NY += dY; + dX = DX; + dY = DY; + } + } else { + // extend this row + dX += DX; + dY = PS_MAX (dY, DY); + } } NY += DY; @@ -1068,46 +1072,46 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // only show "real" saturated stars (not defects) - if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; - if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; - nSAT ++; - - if (Xo + DX > NX) { - // too wide for the rest of this row - if (Xo == 0) { - // place source alone on this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - Yo += DY; - Xo = 0; - dY = 0; - } else { - // start the next row - Yo += dY; - Xo = 0; - - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - Xo = DX; - dY = DY; - } - } else { - // extend this row - bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); - if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - - Xo += DX; - dY = PS_MAX (dY, DY); - } + pmSource *source = sources->data[i]; + + // only show "real" saturated stars (not defects) + if (!(source->mode & PM_SOURCE_MODE_SATSTAR)) continue;; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue;; + nSAT ++; + + if (Xo + DX > NX) { + // too wide for the rest of this row + if (Xo == 0) { + // place source alone on this row + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + Yo += DY; + Xo = 0; + dY = 0; + } else { + // start the next row + Yo += dY; + Xo = 0; + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + Xo = DX; + dY = DY; + } + } else { + // extend this row + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + psphotMosaicSubimage (outsat, source, Xo, Yo, DX, DY, false); + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + Xo += DX; + dY = PS_MAX (dY, DY); + } } @@ -1151,19 +1155,19 @@ float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; for (int iy = 0; iy < source->pixels->numRows; iy++) { - for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { - rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; - // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; - Rb->data.F32[nb] = log10(rb->data.F32[nb]); - fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); - nb++; - } else { - rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; - // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; - Rg->data.F32[ng] = log10(rg->data.F32[ng]); - fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]); - ng++; - } - } + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; + Rb->data.F32[nb] = log10(rb->data.F32[nb]); + fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); + nb++; + } else { + rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; + Rg->data.F32[ng] = log10(rg->data.F32[ng]); + fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]); + ng++; + } + } } @@ -1358,5 +1362,5 @@ KapaSection section; // put the positive profile in one and the residuals in another? - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.profiles", 3)) return true; int myKapa = psphotKapaChannel (2); @@ -1394,22 +1398,22 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; - - psphotVisualPlotRadialProfile (myKapa, source, maskVal); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - if (key[0] == 'e') { - KapaClearPlots (myKapa); - } - if (key[0] == 's') { - break; - } + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + psphotVisualPlotRadialProfile (myKapa, source, maskVal); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + if (key[0] == 'e') { + KapaClearPlots (myKapa); + } + if (key[0] == 's') { + break; + } } @@ -1429,5 +1433,5 @@ return true; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.flags", 3)) return true; int myKapa = psphotKapaChannel (1); @@ -1445,62 +1449,62 @@ for (int i = 0; i < sources->n; i++) { - float Xo, Yo, Rmaj, Rmin, cs, sn; - - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - pmMoments *moments = source->moments; - if (0) { - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - Xo = moments->Mx; - Yo = moments->My; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - Rmaj = 2.0*axes.major; - Rmin = 2.0*axes.minor; - cs = cos(axes.theta); - sn = sin(axes.theta); - } else { - Rmaj = Rmin = 5.0; - cs = 1.0; - sn = 0.0; - Xo = source->peak->xf; - Yo = source->peak->yf; - } - - unsigned short int flagMask = 0x01; - for (int j = 0; j < 8; j++) { - if (source->mode & flagMask) { - overlayE[NoverlayE].type = KII_OVERLAY_LINE; - overlayE[NoverlayE].x = Xo; - overlayE[NoverlayE].y = Yo; - - float phi = j*M_PI/4.0; - overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; - overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; - overlayE[NoverlayE].angle = 0; - overlayE[NoverlayE].text = NULL; - NoverlayE ++; - CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100); - } - flagMask <<= 1; - - if (source->mode & flagMask) { - overlayO[NoverlayO].type = KII_OVERLAY_LINE; - overlayO[NoverlayO].x = Xo + 1; - overlayO[NoverlayO].y = Yo; - - float phi = j*M_PI/4.0; - overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; - overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; - overlayO[NoverlayO].angle = 0; - overlayO[NoverlayO].text = NULL; - NoverlayO ++; - CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100); - } - flagMask <<= 1; - } + float Xo, Yo, Rmaj, Rmin, cs, sn; + + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + pmMoments *moments = source->moments; + if (0) { + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + Xo = moments->Mx; + Yo = moments->My; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + Rmaj = 2.0*axes.major; + Rmin = 2.0*axes.minor; + cs = cos(axes.theta); + sn = sin(axes.theta); + } else { + Rmaj = Rmin = 5.0; + cs = 1.0; + sn = 0.0; + Xo = source->peak->xf; + Yo = source->peak->yf; + } + + unsigned short int flagMask = 0x01; + for (int j = 0; j < 8; j++) { + if (source->mode & flagMask) { + overlayE[NoverlayE].type = KII_OVERLAY_LINE; + overlayE[NoverlayE].x = Xo; + overlayE[NoverlayE].y = Yo; + + float phi = j*M_PI/4.0; + overlayE[NoverlayE].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; + overlayE[NoverlayE].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; + overlayE[NoverlayE].angle = 0; + overlayE[NoverlayE].text = NULL; + NoverlayE ++; + CHECK_REALLOCATE (overlayE, KiiOverlay, NOVERLAYE, NoverlayE, 100); + } + flagMask <<= 1; + + if (source->mode & flagMask) { + overlayO[NoverlayO].type = KII_OVERLAY_LINE; + overlayO[NoverlayO].x = Xo + 1; + overlayO[NoverlayO].y = Yo; + + float phi = j*M_PI/4.0; + overlayO[NoverlayO].dx = +Rmaj*cos(phi)*cs - Rmin*sin(phi)*sn; + overlayO[NoverlayO].dy = +Rmaj*cos(phi)*sn + Rmin*sin(phi)*cs; + overlayO[NoverlayO].angle = 0; + overlayO[NoverlayO].text = NULL; + NoverlayO ++; + CHECK_REALLOCATE (overlayO, KiiOverlay, NOVERLAYO, NoverlayO, 100); + } + flagMask <<= 1; + } } @@ -1531,33 +1535,33 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source == NULL) continue; - - if (mode) { - if (keep) { - if (!(source->mode & mode)) continue; - } else { - if (source->mode & mode) continue; - } - } - - pmMoments *moments = source->moments; - if (moments == NULL) continue; - - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = moments->Mx; - overlay[Noverlay].y = moments->My; - - emoments.x2 = moments->Mxx; - emoments.y2 = moments->Myy; - emoments.xy = moments->Mxy; - - axes = psEllipseMomentsToAxes (emoments, 20.0); - - overlay[Noverlay].dx = scale*2.0*axes.major; - overlay[Noverlay].dy = scale*2.0*axes.minor; - overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + if (mode) { + if (keep) { + if (!(source->mode & mode)) continue; + } else { + if (source->mode & mode) continue; + } + } + + pmMoments *moments = source->moments; + if (moments == NULL) continue; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = moments->Mx; + overlay[Noverlay].y = moments->My; + + emoments.x2 = moments->Mxx; + emoments.y2 = moments->Myy; + emoments.xy = moments->Mxy; + + axes = psEllipseMomentsToAxes (emoments, 20.0); + + overlay[Noverlay].dx = scale*2.0*axes.major; + overlay[Noverlay].dy = scale*2.0*axes.minor; + overlay[Noverlay].angle = axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; } @@ -1570,5 +1574,5 @@ bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.size", 2)) return true; int myKapa = psphotKapaChannel (1); @@ -1597,5 +1601,5 @@ KapaSection section; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.size", 2)) return true; int myKapa = psphotKapaChannel (2); @@ -1612,26 +1616,26 @@ float Ymin = 1000.0, Ymax = 0.0; { - int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); - for (int n = 0; n < nRegions; n++) { - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); - float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); - float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - - float X0 = psfX - 10.0*psfdX; - float X1 = psfX + 10.0*psfdX; - float Y0 = psfY - 10.0*psfdY; - float Y1 = psfY + 10.0*psfdY; - - if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } - if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } - if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } - if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } - } + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + + float X0 = psfX - 10.0*psfdX; + float X1 = psfX + 10.0*psfdX; + float Y0 = psfY - 10.0*psfdY; + float Y1 = psfY + 10.0*psfdY; + + if (isfinite(X0)) { Xmin = PS_MIN(Xmin, X0); } + if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); } + if (isfinite(Y0)) { Ymin = PS_MIN(Ymin, Y0); } + if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); } + } } Xmin = PS_MAX(Xmin, -0.1); @@ -1677,53 +1681,53 @@ int nCR = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->moments == NULL) continue; + pmSource *source = sources->data[i]; + if (source->moments == NULL) continue; // only plot the measured sources... - if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue; - - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { - xCR->data.F32[nCR] = source->moments->Mxx; - yCR->data.F32[nCR] = source->moments->Myy; - mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum); - sCR->data.F32[nCR] = source->extNsigma; - nCR++; - } - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - xSAT->data.F32[nSAT] = source->moments->Mxx; - ySAT->data.F32[nSAT] = source->moments->Myy; - mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum); - sSAT->data.F32[nSAT] = source->extNsigma; - nSAT++; - } - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - xEXT->data.F32[nEXT] = source->moments->Mxx; - yEXT->data.F32[nEXT] = source->moments->Myy; - mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum); - sEXT->data.F32[nEXT] = source->extNsigma; - nEXT++; - continue; - } - if (source->mode & PM_SOURCE_MODE_DEFECT) { - xDEF->data.F32[nDEF] = source->moments->Mxx; - yDEF->data.F32[nDEF] = source->moments->Myy; - mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum); - sDEF->data.F32[nDEF] = source->extNsigma; - nDEF++; - continue; - } - if (source->errMag > 0.1) { - xLOW->data.F32[nLOW] = source->moments->Mxx; - yLOW->data.F32[nLOW] = source->moments->Myy; - mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum); - sLOW->data.F32[nLOW] = source->extNsigma; - nLOW++; - continue; - } - xPSF->data.F32[nPSF] = source->moments->Mxx; - yPSF->data.F32[nPSF] = source->moments->Myy; - mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum); - sPSF->data.F32[nPSF] = source->extNsigma; - nPSF++; + if (!(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)) continue; + + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) { + xCR->data.F32[nCR] = source->moments->Mxx; + yCR->data.F32[nCR] = source->moments->Myy; + mCR->data.F32[nCR] = -2.5*log10(source->moments->Sum); + sCR->data.F32[nCR] = source->extNsigma; + nCR++; + } + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + xSAT->data.F32[nSAT] = source->moments->Mxx; + ySAT->data.F32[nSAT] = source->moments->Myy; + mSAT->data.F32[nSAT] = -2.5*log10(source->moments->Sum); + sSAT->data.F32[nSAT] = source->extNsigma; + nSAT++; + } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + xEXT->data.F32[nEXT] = source->moments->Mxx; + yEXT->data.F32[nEXT] = source->moments->Myy; + mEXT->data.F32[nEXT] = -2.5*log10(source->moments->Sum); + sEXT->data.F32[nEXT] = source->extNsigma; + nEXT++; + continue; + } + if (source->mode & PM_SOURCE_MODE_DEFECT) { + xDEF->data.F32[nDEF] = source->moments->Mxx; + yDEF->data.F32[nDEF] = source->moments->Myy; + mDEF->data.F32[nDEF] = -2.5*log10(source->moments->Sum); + sDEF->data.F32[nDEF] = source->extNsigma; + nDEF++; + continue; + } + if (source->errMag > 0.1) { + xLOW->data.F32[nLOW] = source->moments->Mxx; + yLOW->data.F32[nLOW] = source->moments->Myy; + mLOW->data.F32[nLOW] = -2.5*log10(source->moments->Sum); + sLOW->data.F32[nLOW] = source->extNsigma; + nLOW++; + continue; + } + xPSF->data.F32[nPSF] = source->moments->Mxx; + yPSF->data.F32[nPSF] = source->moments->Myy; + mPSF->data.F32[nPSF] = -2.5*log10(source->moments->Sum); + sPSF->data.F32[nPSF] = source->extNsigma; + nPSF++; } @@ -2029,45 +2033,45 @@ // draw N circles to outline the clumps { - KapaSelectSection (myKapa, "MxxMyy"); - - // draw a circle centered on psfX,Y with size of the psf limit - psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); - psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); - - int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); - float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); - - graphdata.color = KapaColorByName ("blue"); - graphdata.style = 0; - - graphdata.xmin = Xmin; - graphdata.ymin = Ymin; - graphdata.xmax = Xmax; - graphdata.ymax = Ymax; - KapaSetLimits (myKapa, &graphdata); - - for (int n = 0; n < nRegions; n++) { - - char regionName[64]; - snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); - psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); - - float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); - float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); - float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); - float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); - float Rx = psfdX * PSF_CLUMP_NSIGMA; - float Ry = psfdY * PSF_CLUMP_NSIGMA; - - for (int i = 0; i < xLimit->n; i++) { - xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; - yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; - } - KapaPrepPlot (myKapa, xLimit->n, &graphdata); - KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); - KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); - } - psFree (xLimit); - psFree (yLimit); + KapaSelectSection (myKapa, "MxxMyy"); + + // draw a circle centered on psfX,Y with size of the psf limit + psVector *xLimit = psVectorAlloc (120, PS_TYPE_F32); + psVector *yLimit = psVectorAlloc (120, PS_TYPE_F32); + + int nRegions = psMetadataLookupS32 (&status, analysis, "PSF.CLUMP.NREGIONS"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); + + graphdata.color = KapaColorByName ("blue"); + graphdata.style = 0; + + graphdata.xmin = Xmin; + graphdata.ymin = Ymin; + graphdata.xmax = Xmax; + graphdata.ymax = Ymax; + KapaSetLimits (myKapa, &graphdata); + + for (int n = 0; n < nRegions; n++) { + + char regionName[64]; + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n); + psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, regionName); + + float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); + float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); + float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); + float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); + float Rx = psfdX * PSF_CLUMP_NSIGMA; + float Ry = psfdY * PSF_CLUMP_NSIGMA; + + for (int i = 0; i < xLimit->n; i++) { + xLimit->data.F32[i] = Rx*cos(i*2.0*M_PI/120.0) + psfX; + yLimit->data.F32[i] = Ry*sin(i*2.0*M_PI/120.0) + psfY; + } + KapaPrepPlot (myKapa, xLimit->n, &graphdata); + KapaPlotVector (myKapa, xLimit->n, xLimit->data.F32, "x"); + KapaPlotVector (myKapa, yLimit->n, yLimit->data.F32, "y"); + } + psFree (xLimit); + psFree (yLimit); } @@ -2108,5 +2112,5 @@ bool psphotVisualShowResidualImage (pmReadout *readout) { - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.image.resid", 2)) return true; int myKapa = psphotKapaChannel (1); @@ -2124,5 +2128,5 @@ float lineX[2], lineY[2]; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.apresid", 1)) return true; int myKapa = psphotKapaChannel (2); @@ -2144,19 +2148,19 @@ int n = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source) continue; - if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!isfinite (source->apMag)) continue; - if (!isfinite (source->psfMag)) continue; - - x->data.F32[n] = source->psfMag; - y->data.F32[n] = source->apMag - source->psfMag; - dy->data.F32[n] = source->errMag; - graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); - graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); - graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); - graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); - - n++; + pmSource *source = sources->data[i]; + if (!source) continue; + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (!isfinite (source->apMag)) continue; + if (!isfinite (source->psfMag)) continue; + + x->data.F32[n] = source->psfMag; + y->data.F32[n] = source->apMag - source->psfMag; + dy->data.F32[n] = source->errMag; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + + n++; } x->n = y->n = dy->n = n; @@ -2242,5 +2246,5 @@ Graphdata graphdata; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.chisq", 1)) return true; int myKapa = psphotKapaChannel (2); @@ -2263,22 +2267,22 @@ int n = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source) continue; - if (source->type != PM_SOURCE_TYPE_STAR) continue; - if (!source->moments) continue; - if (!isfinite(source->moments->Sum)) continue; - if (!source->modelPSF) continue; - if (!isfinite(source->modelPSF->chisq)) continue; - - x->data.F32[n] = -2.5*log10(source->moments->Sum); - y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF; - graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); - graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); - graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); - graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); - - fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]); - - n++; + pmSource *source = sources->data[i]; + if (!source) continue; + if (source->type != PM_SOURCE_TYPE_STAR) continue; + if (!source->moments) continue; + if (!isfinite(source->moments->Sum)) continue; + if (!source->modelPSF) continue; + if (!isfinite(source->modelPSF->chisq)) continue; + + x->data.F32[n] = -2.5*log10(source->moments->Sum); + y->data.F32[n] = source->modelPSF->chisq / source->modelPSF->nDOF; + graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]); + graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]); + graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); + graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); + + fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]); + + n++; } x->n = y->n = n; @@ -2326,5 +2330,5 @@ KiiOverlay *overlay; - if (!DEBUG && !pmVisualIsVisual()) return true; + if (!pmVisualTestLevel("psphot.objects.petro", 2)) return true; int kapa = psphotKapaChannel (1); @@ -2336,32 +2340,32 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - if (!source) continue; - if (!source->extpars) continue; - if (!source->extpars->petProfile) continue; - - float petrosianRadius = source->extpars->petrosianRadius; + pmSource *source = sources->data[i]; + + if (!source) continue; + if (!source->extpars) continue; + if (!source->extpars->petProfile) continue; + + float petrosianRadius = source->extpars->petrosianRadius; psEllipseAxes *axes = &source->extpars->axes; - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - overlay[Noverlay].dx = 1.0*petrosianRadius; - overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major; - overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - - // overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - // overlay[Noverlay].x = source->peak->xf; - // overlay[Noverlay].y = source->peak->yf; - // overlay[Noverlay].dx = 2.0*petrosianRadius; - // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major; - // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; - // overlay[Noverlay].text = NULL; - // Noverlay ++; - // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 1.0*petrosianRadius; + overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major; + overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + // overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + // overlay[Noverlay].x = source->peak->xf; + // overlay[Noverlay].y = source->peak->yf; + // overlay[Noverlay].dx = 2.0*petrosianRadius; + // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major; + // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; + // overlay[Noverlay].text = NULL; + // Noverlay ++; + // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); }