IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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/psModules/src/objects
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • 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);
Note: See TracChangeset for help on using the changeset viewer.