IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 2, 2009, 3:11:32 PM (17 years ago)
Author:
eugene
Message:

check in changes from genes development branch : extensive changes to moments calculation, psf model generation, aperture residuals

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSourceMoments.c

    r24577 r25754  
    3636#include "pmSource.h"
    3737
    38 bool pmSourceMoments_Old(pmSource *source, psF32 radius);
    39 
    4038/******************************************************************************
    4139pmSourceMoments(source, radius): this function takes a subImage defined in the
    … …  
    5048    pmSource->mask (optional)
    5149
    52 XXX: The peak calculations are done in image coords, not subImage coords.
    53 
    54 this version clips input pixels on S/N
    55 XXX EAM : this version returns false for several reasons
     50The peak calculations are done in image coords, not subImage coords.
     51
     52this version optionally clips input pixels on S/N
    5653*****************************************************************************/
    5754# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    … …  
    6461    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    6562
    66     // XXX supply sky in a different way?
     63    // use sky from moments if defined, 0.0 otherwise
    6764    psF32 sky = 0.0;
    6865    if (source->moments == NULL) {
    … …  
    9188    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
    9289    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
    93     // minimal round-off error in the sums
    94 
    95     psF32 xOff  = source->peak->x;
    96     psF32 yOff  = source->peak->y;
    97     psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
    98     psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
     90    // minimal round-off error in the sums.  since these values are subtracted just to minimize
     91    // the dynamic range and are added back below, the exact value does not matter. these are
     92    // (int) so they can be used in the image index below.
     93
     94    int xOff  = source->peak->x;
     95    int yOff  = source->peak->y;
     96    int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
     97    int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
     98
     99    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     100    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     101    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     102    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     103    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    99104
    100105    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    126131            psF32 wDiff = *vWgt;
    127132
    128             // skip pixels below specified significance level
     133            // skip pixels below specified significance level.  this is allowed, but should be
     134            // avoided -- the over-weights the wings of bright stars compared to those of faint
     135            // stars.
    129136            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    130137            if (pDiff < 0) continue;
    131138
     139            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     140            // weighting over weights the sky for faint sources
    132141            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;
     142                // XXX a lot of extra flops; can we pre-calculate?
     143                psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     144                assert (z >= 0.0);
     145                psF32 weight  = exp(-z);
     146
     147                wDiff *= weight;
     148                pDiff *= weight;
    145149            }
    146150
    … …  
    154158            Y1  += yWght;
    155159
    156             // fprintf (stderr, " : %6.1f %6.1f  %8.1f %8.1f\n", X1, Y1, Sum, Var);
    157 
    158             peakPixel = PS_MAX (*vPix, peakPixel);
    159             numPixels++;
    160         }
     160            peakPixel = PS_MAX (*vPix, peakPixel);
     161            numPixels++;
     162        }
    161163    }
    162164
    … …  
    164166    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    165167    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    166         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
    167         return (false);
     168        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
     169        return (false);
    168170    }
    169171
    … …  
    172174    float My = Y1/Sum;
    173175    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
    174         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
    175         return (false);
     176        psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
     177        return (false);
    176178    }
    177179
    … …  
    179181
    180182    // add back offset of peak in primary image
    181     source->moments->Mx = Mx + xOff;
    182     source->moments->My = My + yOff;
     183    // also offset from pixel index to pixel coordinate
     184    // (the calculation above uses pixel index instead of coordinate)
     185    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
     186    source->moments->Mx = Mx + xOff + 0.5;
     187    source->moments->My = My + yOff + 0.5;
    183188
    184189    source->moments->Sum = Sum;
    … …  
    205210    Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
    206211
    207     // center of mass in subimage
     212    // center of mass in subimage.  Note: the calculation below uses pixel index, so we do not
     213    // correct xCM, yCM to pixel coordinate here.
    208214    psF32 xCM = Mx + xPeak; // coord of peak in subimage
    209215    psF32 yCM = My + yPeak; // coord of peak in subimage
    … …  
    211217    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    212218
    213         psF32 yDiff = row - yCM;
     219        psF32 yDiff = row - yCM;
    214220        if (fabs(yDiff) > radius) continue;
    215221
    216         psF32 *vPix = source->pixels->data.F32[row];
    217         psF32 *vWgt = source->variance->data.F32[row];
    218         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    219 
    220         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    221             if (vMsk) {
    222                 if (*vMsk) {
    223                     vMsk++;
    224                     continue;
    225                 }
    226                 vMsk++;
    227             }
    228             if (isnan(*vPix)) continue;
    229 
    230             psF32 xDiff = col - xCM;
     222        psF32 *vPix = source->pixels->data.F32[row];
     223        psF32 *vWgt = source->variance->data.F32[row];
     224        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     225
     226        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     227            if (vMsk) {
     228                if (*vMsk) {
     229                    vMsk++;
     230                    continue;
     231                }
     232                vMsk++;
     233            }
     234            if (isnan(*vPix)) continue;
     235
     236            psF32 xDiff = col - xCM;
    231237            if (fabs(xDiff) > radius) continue;
    232238
    233             // radius is just a function of (xDiff, yDiff)
    234             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    235             psF32 r  = sqrt(r2);
    236             if (r > radius) continue;
    237 
    238             psF32 pDiff = *vPix - sky;
    239             psF32 wDiff = *vWgt;
    240 
    241             // XXX EAM : should this limit be user-defined?
    242 
    243             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    244             if (pDiff < 0) continue;
    245 
     239            // radius is just a function of (xDiff, yDiff)
     240            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     241            psF32 r  = sqrt(r2);
     242            if (r > radius) continue;
     243
     244            psF32 pDiff = *vPix - sky;
     245            psF32 wDiff = *vWgt;
     246
     247            // skip pixels below specified significance level.  this is allowed, but should be
     248            // avoided -- the over-weights the wings of bright stars compared to those of faint
     249            // stars.
     250            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     251            if (pDiff < 0) continue;
     252
     253            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     254            // weighting over weights the sky for faint sources
    246255            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;
     256                // XXX a lot of extra flops; can we do pre-calculate?
     257                psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
     258                assert (z >= 0.0);
     259                psF32 weight  = exp(-z);
     260
     261                wDiff *= weight;
     262                pDiff *= weight;
    259263            }
    260264
    261             Sum += pDiff;
    262 
    263             psF32 x = xDiff * pDiff;
    264             psF32 y = yDiff * pDiff;
    265 
    266             psF32 xx = xDiff * x;
    267             psF32 xy = xDiff * y;
    268             psF32 yy = yDiff * y;
    269 
    270             psF32 xxx = xDiff * xx / r;
    271             psF32 xxy = xDiff * xy / r;
    272             psF32 xyy = xDiff * yy / r;
    273             psF32 yyy = yDiff * yy / r;
    274 
    275             psF32 xxxx = xDiff * xxx / r2;
    276             psF32 xxxy = xDiff * xxy / r2;
    277             psF32 xxyy = xDiff * xyy / r2;
    278             psF32 xyyy = xDiff * yyy / r2;
    279             psF32 yyyy = yDiff * yyy / r2;
    280 
    281             XX  += xx;
    282             XY  += xy;
    283             YY  += yy;
    284 
    285             XXX  += xxx;
    286             XXY  += xxy;
    287             XYY  += xyy;
    288             YYY  += yyy;
    289 
    290             XXXX  += xxxx;
    291             XXXY  += xxxy;
    292             XXYY  += xxyy;
    293             XYYY  += xyyy;
    294             YYYY  += yyyy;
    295         }
     265            Sum += pDiff;
     266
     267            psF32 x = xDiff * pDiff;
     268            psF32 y = yDiff * pDiff;
     269
     270            psF32 xx = xDiff * x;
     271            psF32 xy = xDiff * y;
     272            psF32 yy = yDiff * y;
     273
     274            psF32 xxx = xDiff * xx / r;
     275            psF32 xxy = xDiff * xy / r;
     276            psF32 xyy = xDiff * yy / r;
     277            psF32 yyy = yDiff * yy / r;
     278
     279            psF32 xxxx = xDiff * xxx / r2;
     280            psF32 xxxy = xDiff * xxy / r2;
     281            psF32 xxyy = xDiff * xyy / r2;
     282            psF32 xyyy = xDiff * yyy / r2;
     283            psF32 yyyy = yDiff * yyy / r2;
     284
     285            XX  += xx;
     286            XY  += xy;
     287            YY  += yy;
     288
     289            XXX  += xxx;
     290            XXY  += xxy;
     291            XYY  += xyy;
     292            YYY  += yyy;
     293
     294            XXXX  += xxxx;
     295            XXXY  += xxxy;
     296            XXYY  += xxyy;
     297            XYYY  += xyyy;
     298            YYYY  += yyyy;
     299        }
    296300    }
    297301
    … …  
    312316
    313317    if (source->moments->Mxx < 0) {
    314       fprintf (stderr, "error: neg second moment??\n");
     318        fprintf (stderr, "error: neg second moment??\n");
    315319    }
    316320    if (source->moments->Myy < 0) {
    317       fprintf (stderr, "error: neg second moment??\n");
     321        fprintf (stderr, "error: neg second moment??\n");
    318322    }
    319323
    320324    psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
    321              source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
    322              source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
    323              source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     325             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
     326             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
     327             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
    324328
    325329    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  sky: %f  Npix: %d\n",
    326              source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, numPixels);
    327 
    328     if (source->moments->Mxx < 0) {
    329         fprintf (stderr, "error: neg second moment??\n");
    330     }
    331     if (source->moments->Myy < 0) {
    332         fprintf (stderr, "error: neg second moment??\n");
    333     }
    334 
    335     // XXX TEST:
    336     // pmSourceMoments_Old (source, radius);
     330             source->peak->xf, source->peak->yf, source->peak->flux, source->peak->SN, source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, numPixels);
     331
    337332    return(true);
    338333}
    339 
    340 
    341 bool pmSourceMoments_Old(pmSource *source, psF32 radius)
    342 {
    343     PS_ASSERT_PTR_NON_NULL(source, false);
    344     PS_ASSERT_PTR_NON_NULL(source->peak, false);
    345     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    346     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    347 
    348     psF32 sky = 0.0;
    349     if (source->moments == NULL) {
    350         source->moments = pmMomentsAlloc();
    351     } else {
    352         sky = source->moments->Sky;
    353     }
    354 
    355     // Sum = SUM (z - sky)
    356     // X1  = SUM (x - xc)*(z - sky)
    357     // X2  = SUM (x - xc)^2 * (z - sky)
    358     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
    359     psF32 peakPixel = -PS_MAX_F32;
    360     psS32 numPixels = 0;
    361     psF32 Sum = 0.0;
    362     psF32 Var = 0.0;
    363     psF32 X1 = 0.0;
    364     psF32 Y1 = 0.0;
    365     psF32 X2 = 0.0;
    366     psF32 Y2 = 0.0;
    367     psF32 XY = 0.0;
    368     psF32 x  = 0;
    369     psF32 y  = 0;
    370     psF32 R2 = PS_SQR(radius);
    371 
    372     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    373     // image coordinates.  the source->pixels image has an offset relative to its parent of
    374     // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
    375     // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
    376     // minimal round-off error in the sums
    377 
    378     psF32 xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
    379     psF32 yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
    380 
    381     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    382 
    383         psF32 *vPix = source->pixels->data.F32[row];
    384         psF32 *vWgt = source->variance->data.F32[row];
    385         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    386 
    387         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    388             if (vMsk) {
    389                 if (*vMsk) {
    390                     vMsk++;
    391                     continue;
    392                 }
    393                 vMsk++;
    394             }
    395             if (isnan(*vPix)) continue;
    396 
    397             psF32 xDiff = col - xPeak;
    398             psF32 yDiff = row - yPeak;
    399 
    400             // radius is just a function of (xDiff, yDiff)
    401             if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
    402 
    403             psF32 pDiff = *vPix - sky;
    404             psF32 wDiff = *vWgt;
    405 
    406             // XXX EAM : should this limit be user-defined?
    407             if (PS_SQR(pDiff) < wDiff) continue;
    408 
    409             Var += wDiff;
    410             Sum += pDiff;
    411 
    412             psF32 xWght = xDiff * pDiff;
    413             psF32 yWght = yDiff * pDiff;
    414 
    415             X1  += xWght;
    416             Y1  += yWght;
    417 
    418             X2  += xDiff * xWght;
    419             XY  += xDiff * yWght;
    420             Y2  += yDiff * yWght;
    421 
    422             peakPixel = PS_MAX (*vPix, peakPixel);
    423             numPixels++;
    424         }
    425     }
    426 
    427     // if we have less than (1/4) of the possible pixels, force a retry
    428     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    429     if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    430         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
    431         return (false);
    432     }
    433 
    434     psTrace ("psModules.objects", 4, "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
    435              sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    436 
    437     //
    438     // first moment X  = X1/Sum + xc
    439     // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    440     // Sxy             = XY / Sum
    441     //
    442     x = X1/Sum;
    443     y = Y1/Sum;
    444     if ((fabs(x) > radius) || (fabs(y) > radius)) {
    445         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n",
    446                  source->peak->x, source->peak->y);
    447         return (false);
    448     }
    449 
    450 # if (PS_TRACE_ON)
    451     float Sxx = PS_MAX(X2/Sum - PS_SQR(x), 0);
    452     float Sxy = XY / Sum;
    453     float Syy = PS_MAX(Y2/Sum - PS_SQR(y), 0);
    454 
    455     psTrace ("psModules.objects", 3,
    456              "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
    457              sky, Sum, x, y, Sxx, Sxy, Syy);
    458 # endif
    459 
    460     return(true);
    461 }
    462 
Note: See TracChangeset for help on using the changeset viewer.