IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24529


Ignore:
Timestamp:
Jun 23, 2009, 7:53:02 AM (17 years ago)
Author:
eugene
Message:

updates to fix moments problem and to do multi-plane photometry

Location:
branches/eam_branches/20090522
Files:
5 added
16 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090522/ippconfig/recipes/psphot.config

    r23301 r24529  
    1111SAVE.PLOTS                          BOOL  FALSE
    1212
     13# format for output CMF files
     14OUTPUT.FORMAT                       STR   PS1_V1
     15
    1316# the zero point is used to set a basic scale for DVO
    1417# XXX it may not currently be read : double check this (EAM)
    1518ZERO_POINT                          F32   25.000          # zero point used by DVO
    1619ZERO_PT                             F32   25.000          # zero point used by DVO
    17 
    18 OUTPUT.FORMAT                       STR   PS1_V1
    1920
    2021# these parameter govern how the background is measured
    … …  
    5354MOMENTS_SY_MAX                      F32   50.0
    5455MOMENTS_AR_MAX                      F32   1.5             # maximum axial ratio: 1 / AR < (sx / sy) < AR
     56MOMENTS_MIN_PIXEL_SN                F32   1.0             # XXX caution on this: too high and we will excessively clip
     57MOMENTS_GAUSS_SIGMA                 F32   2.0             # XXX TEST THIS: sigma in pixels (should be >> seeing equiv)
    5558
    5659# basic object statistics

  • branches/eam_branches/20090522/ppSim/src/ppSimArguments.c

    r24099 r24529  
    169169    }
    170170    ppSimType type = ppSimTypeFromString (typeStr);
     171    // XXX handle error
    171172
    172173    ppSimArgToRecipeStr(&status, options, "FILTER", arguments, "-filter"); // Filter name

  • branches/eam_branches/20090522/ppSim/src/ppSimUtils.c

    r24099 r24529  
    117117    if (!strcasecmp (typeStr, "OBJECT")) return PPSIM_TYPE_OBJECT;
    118118    psAbort("Should never get here.");
     119    // XXX raise error instead of aborting
    119120}
    120121

  • branches/eam_branches/20090522/psModules/src/imcombine/pmPSFEnvelope.c

    r24212 r24529  
    298298        }
    299299
    300         if (!pmSourceMoments(source, maxRadius)) {
     300        // measure the source moments: tophat windowing, pixel S/N > 1.0
     301        if (!pmSourceMoments(source, maxRadius, 0.0, 1.0)) {
    301302            // Can't do anything about it; limp along as best we can
    302303            psErrorClear();

  • branches/eam_branches/20090522/psModules/src/objects/pmSource.c

    r23989 r24529  
    291291
    292292    bool status = true;                 // Status of MD lookup
    293     float PSF_CLUMP_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_SN_LIM");
     293    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    294294    if (!status) {
    295         PSF_CLUMP_SN_LIM = 0;
     295        PSF_SN_LIM = 0;
    296296    }
    297297    float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
    … …  
    342342            }
    343343
    344             if (src->moments->SN < PSF_CLUMP_SN_LIM) {
     344            if (src->moments->SN < PSF_SN_LIM) {
    345345                psTrace("psModules.objects", 10, "Rejecting source from clump because of low S/N (%f)\n",
    346346                        src->moments->SN);
    … …  
    450450            if (tmpSrc->moments == NULL)
    451451                continue;
    452             if (tmpSrc->moments->SN < PSF_CLUMP_SN_LIM)
     452            if (tmpSrc->moments->SN < PSF_SN_LIM)
    453453                continue;
    454454
    … …  
    534534    bool status;
    535535    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
    536     if (!status)
    537         PSF_SN_LIM = 20.0;
     536    if (!status) PSF_SN_LIM = 20.0;
    538537    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    539     if (!status)
    540         PSF_CLUMP_NSIGMA = 1.5;
    541     float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     538    if (!status) PSF_CLUMP_NSIGMA = 1.5;
     539
     540    // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    542541
    543542    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
    … …  
    582581            source->type = PM_SOURCE_TYPE_STAR;
    583582            source->mode |= PM_SOURCE_MODE_SATSTAR;
    584             // recalculate moments here with larger box?
    585             pmSourceMoments (source, INNER_RADIUS);
     583            // why do we recalculate moments here?
     584            // we already attempt to do this in psphotSourceStats
     585            // pmSourceMoments (source, INNER_RADIUS);
    586586            Nsatstar ++;
    587587            continue;

  • branches/eam_branches/20090522/psModules/src/objects/pmSource.h

    r23487 r24529  
    213213 */
    214214bool pmSourceMoments(
    215     pmSource *source,   ///< The input pmSource for which moments will be computed
    216     float radius   ///< Use a circle of pixels around the peak
     215    pmSource *source, ///< The input pmSource for which moments will be computed
     216    float radius,     ///< Use a circle of pixels around the peak
     217    float sigma,      ///< size of Gaussian window function (<= 0.0 -> skip window)
     218    float minSN       ///< minimum pixel significance
    217219);
    218220

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceIO_RAW.c

    r20937 r24529  
    254254    }
    255255
     256    fprintf (f, "# %5s %5s  %8s  %7s %7s  %6s %6s  %10s %7s %7s %7s  %4s %4s %5s\n",
     257             "x", "y", "peak", "Mx", "My", "Mxx", "Myy", "Sum", "Peak", "Sky", "SN", "nPix", "type", "mode");
     258
    256259    for (i = 0; i < sources->n; i++) {
    257260        source = sources->data[i];
    258261        if (source->moments == NULL)
    259262            continue;
    260         fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
     263        fprintf (f, "%5d %5d  %8.1f  %7.1f %7.1f  %6.3f %6.3f  %10.1f %7.1f %7.1f %7.1f  %4d %2d %#5x\n",
    261264                 source->peak->x, source->peak->y, source->peak->value,
    262265                 source->moments->Mx, source->moments->My,

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceMoments.c

    r21363 r24529  
    5757# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    5858
    59 bool pmSourceMoments(pmSource *source, psF32 radius)
     59bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN)
    6060{
    6161    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    8484    psF32 Y1 = 0.0;
    8585    psF32 R2 = PS_SQR(radius);
     86    psF32 minSN2 = PS_SQR(minSN);
     87    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
    8688
    8789    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    9799
    98100    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     101
     102        psF32 yDiff = row - yPeak;
     103        if (fabs(yDiff) > radius) continue;
    99104
    100105        psF32 *vPix = source->pixels->data.F32[row];
    … …  
    113118
    114119            psF32 xDiff = col - xPeak;
    115             psF32 yDiff = row - yPeak;
     120            if (fabs(xDiff) > radius) continue;
    116121
    117122            // radius is just a function of (xDiff, yDiff)
    … …  
    121126            psF32 wDiff = *vWgt;
    122127
    123             // XXX EAM : should this limit be user-defined?
    124             if (PS_SQR(pDiff) < wDiff) continue;
    125 
    126             Var += wDiff;
    127             Sum += pDiff;
    128 
    129             psF32 xWght = xDiff * pDiff;
    130             psF32 yWght = yDiff * pDiff;
    131 
    132             X1  += xWght;
    133             Y1  += yWght;
     128            // skip pixels below specified significance level
     129            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     130            if (pDiff < 0) continue;
     131
     132            if (sigma > 0.0) {
     133              // apply a pseudo-gaussian weight
     134
     135              // XXX a lot of extra flops; can we do pre-calculate?
     136              psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     137              assert (z >= 0.0);
     138              psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0));
     139              psF32 weight  = 1.0 / t;
     140
     141              // fprintf (stderr, "%6.1f %6.1f  %8.1f %8.1f  %5.3f  ", xDiff, yDiff, pDiff, wDiff, weight);
     142
     143              wDiff *= weight;
     144              pDiff *= weight;
     145            }
     146
     147            Var += wDiff;
     148            Sum += pDiff;
     149
     150            psF32 xWght = xDiff * pDiff;
     151            psF32 yWght = yDiff * pDiff;
     152
     153            X1  += xWght;
     154            Y1  += yWght;
     155
     156            // fprintf (stderr, " : %6.1f %6.1f  %8.1f %8.1f\n", X1, Y1, Sum, Var);
    134157
    135158            peakPixel = PS_MAX (*vPix, peakPixel);
    … …  
    188211    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    189212
     213        psF32 yDiff = row - yCM;
     214        if (fabs(yDiff) > radius) continue;
     215
    190216        psF32 *vPix = source->pixels->data.F32[row];
    191217        psF32 *vWgt = source->variance->data.F32[row];
    … …  
    203229
    204230            psF32 xDiff = col - xCM;
    205             psF32 yDiff = row - yCM;
     231            if (fabs(xDiff) > radius) continue;
    206232
    207233            // radius is just a function of (xDiff, yDiff)
    … …  
    215241            // XXX EAM : should this limit be user-defined?
    216242
    217             if (PS_SQR(pDiff) < wDiff) continue;
     243            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    218244            if (pDiff < 0) continue;
     245
     246            if (sigma > 0.0) {
     247              // apply a pseudo-gaussian weight
     248
     249              // XXX a lot of extra flops; can we do pre-calculate?
     250              psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     251              assert (z >= 0.0);
     252              psF32 t = 1.0 + z*(1.0 + z/2.0*(1.0 + z/3.0));
     253              psF32 weight  = 1.0 / t;
     254
     255              // fprintf (stderr, "%6.1f %6.1f  %8.1f %8.1f  %5.3f  ", xDiff, yDiff, pDiff, wDiff, weight);
     256
     257              wDiff *= weight;
     258              pDiff *= weight;
     259            }
    219260
    220261            Sum += pDiff;

  • branches/eam_branches/20090522/psModules/src/objects/pmSourceUtils.c

    r23187 r24529  
    9494    source->peak = pmPeakAlloc (xChip, yChip, Io, PM_PEAK_LONE);
    9595
    96     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    97     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    98     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    99     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
     96    float xReadout, yReadout;
    10097
    101     // XXX fix the binning : currently not selected from concepts
    102     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    103     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    104     int xBin = 1;
    105     int yBin = 1;
     98    // if we have information about the chip & cell, adjust the coordinates chip->cell->readout
     99    // otherwise, assume 0,0 offset and 1,1 parity
     100    if (cell) {
     101      int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
     102      int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
     103      int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
     104      int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    106105
    107     // Position on the cell
    108     float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
    109     float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
     106      // XXX fix the binning : currently not selected from concepts
     107      // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
     108      // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
     109      int xBin = 1;
     110      int yBin = 1;
    110111
    111     // Position on the readout
    112     // float xReadout = CELL_TO_READOUT(xCell, x0Readout);
    113     // float yReadout = CELL_TO_READOUT(yCell, y0Readout);
     112      // Position on the cell
     113      float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
     114      float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
     115
     116      // Position on the readout
     117      // float xReadout = CELL_TO_READOUT(xCell, x0Readout);
     118      // float yReadout = CELL_TO_READOUT(yCell, y0Readout);
     119      xReadout = xCell;
     120      yReadout = yCell;
     121    } else {
     122      xReadout = xChip;
     123      yReadout = yChip;
     124    }
    114125   
    115     pmSourceDefinePixels (source, readout, xCell, yCell, radius);
     126    pmSourceDefinePixels (source, readout, xReadout, yReadout, radius);
    116127
    117128    return (source);

  • branches/eam_branches/20090522/psModules/test/objects/tap_pmSource.c

    r21471 r24529  
    432432        psArray *sources = psArrayAlloc(NUM_SOURCES);
    433433        psMetadata *recipe = psMetadataAlloc();
    434         bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0);
     434        bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0);
    435435        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0);
    436436        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0);
    … …  
    452452        psArray *sources = psArrayAlloc(NUM_SOURCES);
    453453        psMetadata *recipe = psMetadataAlloc();
    454         bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0);
     454        bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0);
    455455        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0);
    456456        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0);
    … …  
    490490        }
    491491        psMetadata *recipe = psMetadataAlloc();
    492         bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_CLUMP_SN_LIM", 0, NULL, 0.0);
     492        bool rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "PSF_SN_LIM", 0, NULL, 0.0);
    493493        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SX_MAX", 0, NULL, 10.0);
    494494        rc = psMetadataAddF32(recipe, PS_LIST_HEAD, "MOMENTS_SY_MAX", 0, NULL, 10.0);

  • branches/eam_branches/20090522/psphot/src/Makefile.am

    r23808 r24529  
    1616
    1717# Force recompilation of psphotVersion.c, since it gets the version information
    18 psphotVersion.c: psphotVersionDefinitions.h
    19 psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
    20         -$(RM) psphotVersionDefinitions.h
    21         $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
    22 FORCE: ;
     18# psphotVersion.c: psphotVersionDefinitions.h
     19# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
     20#       -$(RM) psphotVersionDefinitions.h
     21#       $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
     22# FORCE: ;
    2323
    2424libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotTest
     27bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
    2828
    2929psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    3434psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    3535psphotTest_LDADD = libpsphot.la
     36
     37psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     38psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     39psphotMomentsStudy_LDADD = libpsphot.la
    3640
    3741psphot_SOURCES = \
    … …  
    5458psphotTest_SOURCES = \
    5559        psphotTest.c
     60
     61psphotMomentsStudy_SOURCES = \
     62        psphotMomentsStudy.c
    5663
    5764libpsphot_la_SOURCES = \

  • branches/eam_branches/20090522/psphot/src/psphotArguments.c

    r24144 r24529  
    105105    pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE);
    106106    if (config == NULL) {
    107         psError(PSPHOT_ERR_CONFIG, false, "Can't read site configuration");
     107      psErrorStackPrint(stderr, "Can't read site configuration");
    108108        exit(PS_EXIT_CONFIG_ERROR);
    109109    }
    … …  
    121121    if ((N = psArgumentGet (argc, argv, "-modeltest"))) {
    122122        if (argc<=N+2) {
    123           psError(PSPHOT_ERR_ARGUMENTS, true,
    124                   "Expected to see 2 more arguments; saw %d", argc - 1);
     123          psErrorStackPrint(stderr, "Expected to see 2 more arguments; saw %d", argc - 1);
    125124          exit(PS_EXIT_CONFIG_ERROR);
    126125        }
    … …  
    137136        if ((N = psArgumentGet (argc, argv, "-model"))) {
    138137            if (argc<=N+1) {
    139               psError(PSPHOT_ERR_ARGUMENTS, true,
    140                       "Expected to see 1 more argument; saw %d", argc - 1);
     138              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    141139              exit(PS_EXIT_CONFIG_ERROR);
    142140            }
    … …  
    149147        if ((N = psArgumentGet (argc, argv, "-fitmode"))) {
    150148            if (argc<=N+1) {
    151               psError(PSPHOT_ERR_ARGUMENTS, true,
    152                       "Expected to see 1 more argument; saw %d", argc - 1);
     149              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    153150              exit(PS_EXIT_CONFIG_ERROR);
    154151            }
    … …  
    159156        if ((N = psArgumentGet (argc, argv, "-fitset"))) {
    160157            if (argc<=N+1) {
    161               psError(PSPHOT_ERR_ARGUMENTS, true,
    162                       "Expected to see 1 more argument; saw %d", argc - 1);
     158              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    163159              exit(PS_EXIT_CONFIG_ERROR);
    164160            }
    … …  
    172168    if ((N = psArgumentGet (argc, argv, "-photcode"))) {
    173169        if (argc<=N+1) {
    174           psError(PSPHOT_ERR_ARGUMENTS, true,
    175                   "Expected to see 1 more argument; saw %d", argc - 1);
     170          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    176171          exit(PS_EXIT_CONFIG_ERROR);
    177172        }
    … …  
    190185    if ((N = psArgumentGet (argc, argv, "-break"))) {
    191186        if (argc<=N+1) {
    192           psError(PSPHOT_ERR_ARGUMENTS, true,
    193                   "Expected to see 1 more argument; saw %d", argc - 1);
     187          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    194188          exit(PS_EXIT_CONFIG_ERROR);
    195189        }
    … …  
    202196    if ((N = psArgumentGet (argc, argv, "-region"))) {
    203197        if (argc<=N+1) {
    204           psError(PSPHOT_ERR_ARGUMENTS, true,
    205                   "Expected to see 1 more argument; saw %d", argc - 1);
     198          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    206199          exit(PS_EXIT_CONFIG_ERROR);
    207200        }
    … …  
    214207    if ((N = psArgumentGet (argc, argv, "-chip"))) {
    215208        if (argc<=N+1) {
    216           psError(PSPHOT_ERR_ARGUMENTS, true,
    217                   "Expected to see 1 more argument; saw %d", argc - 1);
     209          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    218210          exit(PS_EXIT_CONFIG_ERROR);
    219211        }

  • branches/eam_branches/20090522/psphot/src/psphotModelTest.c

    r21183 r24529  
    139139
    140140    // get the source moments
    141     status = pmSourceMoments (source, mRADIUS);
     141    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0);
    142142    if (!status) psAbort("psSourceMoments error");
    143143    source->peak->value = source->moments->Peak;

  • branches/eam_branches/20090522/psphot/src/psphotSourceFits.c

    r24188 r24529  
    224224
    225225    // recalculate the source moments using the larger extended-source moments radius
    226     if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false;
     226    // at this stage, skip Gaussian windowing, and allow pixels with S/N > 0.5
     227    if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.5)) return false;
    227228
    228229    psTrace ("psphot", 5, "trying blob...\n");

  • branches/eam_branches/20090522/psphot/src/psphotSourceStats.c

    r23219 r24529  
    8989
    9090            if (!psThreadJobAddPending(job)) {
    91                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     91                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    9292                psFree (job);
    9393                return NULL;
    9494            }
    9595            psFree(job);
    96 
    97 # if (0)
    98                 int nfail = 0;
    99                 int nmoments = 0;
    100                 if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
    101                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    102                     return NULL;
    103                 }
    104                 Nfail += nfail;
    105                 Nmoments += nmoments;
    106 # endif
    10796        }
    10897
    10998        // wait for the threads to finish and manage results
    11099        if (!psThreadPoolWait (false)) {
    111             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     100            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
    112101            return NULL;
    113102        }
    … …  
    149138    psMetadata *recipe              = job->args->data[1];
    150139
    151     float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    152     if (!status) return false;
    153     float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    154     if (!status) return false;
    155     float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     140    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     141    if (!status) return false;
     142    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     143    if (!status) return false;
     144    float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     145    if (!status) return false;
     146    float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN");
     147    if (!status) return false;
     148    float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
    156149    if (!status) return false;
    157150
    … …  
    202195
    203196        // measure basic source moments
    204         status = pmSourceMoments (source, RADIUS);
     197        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
    205198        if (status) {
    206199            Nmoments ++;
    … …  
    212205        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    213206        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    214         status = pmSourceMoments (source, BIG_RADIUS);
     207        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
    215208        if (status) {
    216209            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    … …  
    292285
    293286        // measure basic source moments
    294         status = pmSourceMoments (source, RADIUS);
     287        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
    295288        if (status) {
    296289            Nmoments ++;
    … …  
    302295        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    303296        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    304         status = pmSourceMoments (source, BIG_RADIUS);
     297        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
    305298        if (status) {
    306299            Nmoments ++;

  • branches/eam_branches/20090522/psphot/src/psphotVisual.c

    r23259 r24529  
    471471    bool status;
    472472    Graphdata graphdata;
    473 
    474     if (!pmVisualIsVisual()) return true;
     473    KapaSection section;
     474
     475    // if (!pmVisualIsVisual()) return true;
    475476
    476477    if (kapa3 == -1) {
    … …  
    485486    KapaClearPlots (kapa3);
    486487    KapaInitGraph (&graphdata);
    487 
    488     // there are N regions: use the first (guaranteed to exist) to get the overall limits
    489     psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    490 
    491     float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    492     float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    493 
    494     // examine sources to set data range
    495     graphdata.xmin = -0.05;
    496     graphdata.ymin = -0.05;
    497     graphdata.xmax = 2.0*psfX;
    498     graphdata.ymax = 2.0*psfY;
    499     KapaSetLimits (kapa3, &graphdata);
    500 
    501     KapaSetFont (kapa3, "helvetica", 14);
    502     KapaBox (kapa3, &graphdata);
    503     KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    504     KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
    505 
     488    KapaSetFont (kapa3, "courier", 14);
     489
     490    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
     491
     492    // select the max psfX,Y values for the plot limits
     493    float Xmin = 0.0, Xmax = 0.0;
     494    float Ymin = 0.0, Ymax = 0.0;
     495    {
     496        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     497        for (int n = 0; n < nRegions; n++) {
     498
     499            char regionName[64];
     500            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     501            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     502
     503            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     504            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     505            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     506            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     507
     508            float X0 = psfX - 4.0*psfdX;
     509            float X1 = psfX + 4.0*psfdX;
     510            float Y0 = psfY - 4.0*psfdY;
     511            float Y1 = psfY + 4.0*psfdY;
     512
     513            if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
     514            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     515            if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
     516            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     517        }
     518    }
     519
     520    // storage vectors for data to be plotted
    506521    psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    507522    psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     523    psVector *mBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    508524    psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    509525    psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     526    psVector *mFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    510527
    511528    // construct the vectors
    … …  
    519536        xFaint->data.F32[nF] = source->moments->Mxx;
    520537        yFaint->data.F32[nF] = source->moments->Myy;
     538        mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum);
    521539        nF++;
    522540
    523541        // XXX make this a user-defined cutoff
    524         if (source->moments->SN < 50)
     542        if (source->moments->SN < SN_LIM)
    525543            continue;
    526544
    527545        xBright->data.F32[nB] = source->moments->Mxx;
    528546        yBright->data.F32[nB] = source->moments->Myy;
     547        mBright->data.F32[nB] = -2.5*log10(source->moments->Sum);
    529548        nB++;
    530549    }
    531550    xFaint->n = nF;
    532551    yFaint->n = nF;
     552    mFaint->n = nF;
    533553
    534554    xBright->n = nB;
    535555    yBright->n = nB;
     556    mBright->n = nB;
     557
     558    // three sections: MxxMyy, MagMxx, MagMyy
     559
     560    // first section: MxxMyy
     561    section.dx = 0.75;
     562    section.dy = 0.75;
     563    section.x  = 0.00;
     564    section.y  = 0.00;
     565    section.name = psStringCopy ("MxxMyy");
     566    KapaSetSection (kapa3, &section);
     567    psFree (section.name);
     568
     569    graphdata.color = KapaColorByName ("black");
     570    graphdata.xmin = Xmin;
     571    graphdata.ymin = Ymin;
     572    graphdata.xmax = Xmax;
     573    graphdata.ymax = Ymax;
     574    KapaSetLimits (kapa3, &graphdata);
     575
     576    KapaBox (kapa3, &graphdata);
     577    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
     578    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
    536579
    537580    graphdata.color = KapaColorByName ("black");
    … …  
    551594    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
    552595
    553     // XXX draw N circles to outline the clumps
     596    // second section: MagMyy
     597    section.dx = 0.75;
     598    section.dy = 0.25;
     599    section.x  = 0.00;
     600    section.y  = 0.80;
     601    section.name = psStringCopy ("MagMyy");
     602    KapaSetSection (kapa3, &section);
     603    psFree (section.name);
     604
     605    graphdata.color = KapaColorByName ("black");
     606    graphdata.xmin = -17.1;
     607    graphdata.xmax =  -7.9;
     608    graphdata.ymin = Ymin;
     609    graphdata.ymax = Ymax;
     610    KapaSetLimits (kapa3, &graphdata);
     611
     612    strcpy (graphdata.labels, "0210");
     613    KapaBox (kapa3, &graphdata);
     614    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP);
     615    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
     616
     617    graphdata.color = KapaColorByName ("black");
     618    graphdata.ptype = 0;
     619    graphdata.size = 0.3;
     620    graphdata.style = 2;
     621    KapaPrepPlot (kapa3, nF, &graphdata);
     622    KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
     623    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     624
     625    graphdata.color = KapaColorByName ("red");
     626    graphdata.ptype = 2;
     627    graphdata.size = 0.5;
     628    graphdata.style = 2;
     629    KapaPrepPlot (kapa3, nB, &graphdata);
     630    KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
     631    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     632
     633    // third section: MagMxx
     634    section.dx = 0.25;
     635    section.dy = 0.75;
     636    section.x  = 0.80;
     637    section.y  = 0.00;
     638    section.name = psStringCopy ("MagMxx");
     639    KapaSetSection (kapa3, &section);
     640    psFree (section.name);
     641
     642    graphdata.color = KapaColorByName ("black");
     643    graphdata.xmin = Xmin;
     644    graphdata.xmax = Xmax;
     645    graphdata.ymin =  -7.9;
     646    graphdata.ymax = -17.1;
     647    KapaSetLimits (kapa3, &graphdata);
     648
     649    strcpy (graphdata.labels, "2001");
     650    KapaBox (kapa3, &graphdata);
     651    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
     652    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP);
     653
     654    graphdata.color = KapaColorByName ("black");
     655    graphdata.ptype = 0;
     656    graphdata.size = 0.3;
     657    graphdata.style = 2;
     658    KapaPrepPlot (kapa3, nF, &graphdata);
     659    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
     660    KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
     661
     662    graphdata.color = KapaColorByName ("red");
     663    graphdata.ptype = 2;
     664    graphdata.size = 0.5;
     665    graphdata.style = 2;
     666    KapaPrepPlot (kapa3, nB, &graphdata);
     667    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
     668    KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
     669
     670    // draw N circles to outline the clumps
    554671    {
     672        KapaSelectSection (kapa3, "MxxMyy");
     673
    555674        // draw a circle centered on psfX,Y with size of the psf limit
    556675        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
    … …  
    562681        graphdata.color = KapaColorByName ("blue");
    563682        graphdata.style = 0;
     683
     684        graphdata.xmin = Xmin;
     685        graphdata.ymin = Ymin;
     686        graphdata.xmax = Xmax;
     687        graphdata.ymax = Ymax;
     688        KapaSetLimits (kapa3, &graphdata);
    564689
    565690        for (int n = 0; n < nRegions; n++) {
    … …  
    626751    psFree (xBright);
    627752    psFree (yBright);
     753    psFree (mBright);
    628754    psFree (xFaint);
    629755    psFree (yFaint);
     756    psFree (mFaint);
    630757
    631758    // pause and wait for user input:
Note: See TracChangeset for help on using the changeset viewer.