IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29443


Ignore:
Timestamp:
Oct 16, 2010, 10:49:43 AM (16 years ago)
Author:
eugene
Message:

split centroid from rest of moments, do not remeasure centroid for undetected or supplied sources

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100823/psModules/src/objects/pmSourceMoments.c

    r29124 r29443  
    6464void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
    6565
     66bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal);
     67
     68// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
     69
    6670bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
    6771{
    … …  
    7175    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    7276
    73     // use sky from moments if defined, 0.0 otherwise
    74 
    75     // XXX this value comes from the sky model at the source center, and tends to over-estimate
    76     // the sky in the vicinity of bright sources.  we are better off assuming the model worked
    77     // well:
     77    // this function assumes the sky has been well-subtracted for the image
    7878    psF32 sky = 0.0;
     79
    7980    if (source->moments == NULL) {
    8081      source->moments = pmMomentsAlloc();
    8182    }
    82     // XXX if (source->moments == NULL) {
    83     // XXX     source->moments = pmMomentsAlloc();
    84     // XXX } else {
    85     // XXX     sky = source->moments->Sky;
    86     // XXX }
    87 
    88     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    89     // Sum = SUM (z - sky)
    90     // X1  = SUM (x - xc)*(z - sky)
    91     // .. etc
    92 
    93     psF32 peakPixel = -PS_MAX_F32;
    94     psS32 numPixels = 0;
     83
    9584    psF32 Sum = 0.0;
    9685    psF32 Var = 0.0;
    97     psF32 X1 = 0.0;
    98     psF32 Y1 = 0.0;
    9986    psF32 R2 = PS_SQR(radius);
    10087    psF32 minSN2 = PS_SQR(minSN);
    … …  
    10996    // (int) so they can be used in the image index below.
    11097
    111     int xOff  = source->peak->x;
    112     int yOff  = source->peak->y;
    113     int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
    114     int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
    115 
    116     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
    117     // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
    118     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    119     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    120     // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    121 
    122     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    123 
    124         psF32 yDiff = row - yPeak;
    125         if (fabs(yDiff) > radius) continue;
    126 
    127         psF32 *vPix = source->pixels->data.F32[row];
    128         psF32 *vWgt = source->variance->data.F32[row];
    129         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    130 
    131         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    132             if (vMsk) {
    133                 if (*vMsk & maskVal) {
    134                     vMsk++;
    135                     continue;
    136                 }
    137                 vMsk++;
    138             }
    139             if (isnan(*vPix)) continue;
    140 
    141             psF32 xDiff = col - xPeak;
    142             if (fabs(xDiff) > radius) continue;
    143 
    144             // radius is just a function of (xDiff, yDiff)
    145             if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
    146 
    147             psF32 pDiff = *vPix - sky;
    148             psF32 wDiff = *vWgt;
    149 
    150             // skip pixels below specified significance level.  for a PSFs, this
    151             // over-weights the wings of bright stars compared to those of faint stars.
    152             // for the estimator used for extended source analysis (where the window
    153             // function is allowed to be arbitrarily large), we need to clip to avoid
    154             // negative second moments.
    155             if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
    156             if ((minSN > 0.0) && (pDiff < 0)) continue; //
    157 
    158             // Apply a Gaussian window function.  Be careful with the window function.  S/N
    159             // weighting over weights the sky for faint sources
    160             if (sigma > 0.0) {
    161                 // XXX a lot of extra flops; can we pre-calculate?
    162                 psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
    163                 assert (z >= 0.0);
    164                 psF32 weight  = exp(-z);
    165 
    166                 wDiff *= weight;
    167                 pDiff *= weight;
    168             }
    169 
    170             Var += wDiff;
    171             Sum += pDiff;
    172 
    173             psF32 xWght = xDiff * pDiff;
    174             psF32 yWght = yDiff * pDiff;
    175 
    176             X1  += xWght;
    177             Y1  += yWght;
    178 
    179             peakPixel = PS_MAX (*vPix, peakPixel);
    180             numPixels++;
    181         }
    182     }
    183 
    184     // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
    185     int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
    186 
    187     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    188     if ((numPixels < minPixels) || (Sum <= 0)) {
    189         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
    190         return (false);
    191     }
    192 
    193     // calculate the first moment.
    194     float Mx = X1/Sum;
    195     float My = Y1/Sum;
    196     if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
    197         psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
    198         return (false);
    199     }
    200 
    201     psTrace ("psModules.objects", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
    202 
    203     // add back offset of peak in primary image
    204     // also offset from pixel index to pixel coordinate
    205     // (the calculation above uses pixel index instead of coordinate)
    206     // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
    207     source->moments->Mx = Mx + xOff + 0.5;
    208     source->moments->My = My + yOff + 0.5;
    209 
    210     source->moments->Sum = Sum;
    211     source->moments->SN  = Sum / sqrt(Var);
    212     source->moments->Peak = peakPixel;
    213     source->moments->nPixels = numPixels;
     98    if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal)) {
     99        return false;
     100    }
    214101
    215102    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
    … …  
    234121    Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
    235122
    236     // center of mass in subimage.  Note: the calculation below uses pixel index, so we do not
    237     // correct xCM, yCM to pixel coordinate here.
    238     psF32 xCM = Mx + xPeak; // coord of peak in subimage
    239     psF32 yCM = My + yPeak; // coord of peak in subimage
     123    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     124    // xCM, yCM from pixel coords to pixel index here.
     125    psF32 xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
     126    psF32 yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
    240127
    241128    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    263150            // radius is just a function of (xDiff, yDiff)
    264151            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    265             psF32 r  = sqrt(r2);
    266             if (r > radius) continue;
     152            if (r2 > R2) continue;
    267153
    268154            psF32 fDiff = *vPix - sky;
    … …  
    274160            // stars.
    275161            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    276             // if (pDiff < 0) continue;
     162            if ((minSN > 0.0) && (pDiff < 0)) continue; //
    277163
    278164            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    279165            // weighting over weights the sky for faint sources
    280166            if (sigma > 0.0) {
    281                 // XXX a lot of extra flops; can we do pre-calculate?
    282                 // XXX we were re-calculating r2 (maybe the compiler caught this?)
    283                 // psF32 z  = (PS_SQR(xDiff) + PS_SQR(yDiff))*rsigma2;
    284                 psF32 z  = r2 * rsigma2;
     167                psF32 z = r2 * rsigma2;
    285168                assert (z >= 0.0);
    286169                psF32 weight  = exp(-z);
    … …  
    292175            Sum += pDiff;
    293176
    294 # if (1)
    295 # if (0)
    296             if (fDiff < 0) continue;
    297 # endif
     177            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
     178            psF32 r = sqrt(r2);
    298179            psF32 rf = r * fDiff;
    299180            psF32 rh = sqrt(r) * fDiff;
    300181            psF32 rs = fDiff;
    301 # else
    302             psF32 rf = r * pDiff;
    303             psF32 rh = sqrt(r) * pDiff;
    304             psF32 rs = pDiff;
    305 # endif
    306182
    307183            psF32 x = xDiff * pDiff;
    … …  
    363239
    364240    // Calculate the Kron magnitude (make this block optional?)
    365     // float radKron = 2.5*source->moments->Mrf;
    366241    float radKinner = 1.0*source->moments->Mrf;
    367242    float radKron   = 2.5*source->moments->Mrf;
    … …  
    397272            // radKron is just a function of (xDiff, yDiff)
    398273            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    399             psF32 r  = sqrt(r2);
    400274
    401275            psF32 pDiff = *vPix - sky;
    … …  
    407281            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    408282
     283            psF32 r  = sqrt(r2);
    409284            if (r < radKron) {
    410285                Sum += pDiff;
    … …  
    434309
    435310    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",
    436              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);
     311             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, source->moments->nPixels);
    437312
    438313    return(true);
    439314}
     315
     316bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal) {
     317
     318    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
     319    // Sum = SUM (z - sky)
     320    // X1  = SUM (x - xc)*(z - sky)
     321    // .. etc
     322
     323    psF32 sky = 0.0;
     324
     325    psF32 peakPixel = -PS_MAX_F32;
     326    psS32 numPixels = 0;
     327    psF32 Sum = 0.0;
     328    psF32 Var = 0.0;
     329    psF32 X1 = 0.0;
     330    psF32 Y1 = 0.0;
     331    psF32 R2 = PS_SQR(radius);
     332    psF32 minSN2 = PS_SQR(minSN);
     333    psF32 rsigma2 = 0.5 / PS_SQR(sigma);
     334
     335    int xOff  = source->peak->x;
     336    int yOff  = source->peak->y;
     337    int xPeak = source->peak->x - source->pixels->col0; // coord of peak in subimage
     338    int yPeak = source->peak->y - source->pixels->row0; // coord of peak in subimage
     339
     340    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     341    // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     342    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     343    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     344    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
     345
     346    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     347
     348        psF32 yDiff = row - yPeak;
     349        if (fabs(yDiff) > radius) continue;
     350
     351        psF32 *vPix = source->pixels->data.F32[row];
     352        psF32 *vWgt = source->variance->data.F32[row];
     353        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     354
     355        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     356            if (vMsk) {
     357                if (*vMsk & maskVal) {
     358                    vMsk++;
     359                    continue;
     360                }
     361                vMsk++;
     362            }
     363            if (isnan(*vPix)) continue;
     364
     365            psF32 xDiff = col - xPeak;
     366            if (fabs(xDiff) > radius) continue;
     367
     368            // radius is just a function of (xDiff, yDiff)
     369            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     370            if (r2 > R2) continue;
     371
     372            psF32 pDiff = *vPix - sky;
     373            psF32 wDiff = *vWgt;
     374
     375            // skip pixels below specified significance level.  for a PSFs, this
     376            // over-weights the wings of bright stars compared to those of faint stars.
     377            // for the estimator used for extended source analysis (where the window
     378            // function is allowed to be arbitrarily large), we need to clip to avoid
     379            // negative second moments.
     380            if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
     381            if ((minSN > 0.0) && (pDiff < 0)) continue; //
     382
     383            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     384            // weighting over weights the sky for faint sources
     385            if (sigma > 0.0) {
     386                psF32 z  = r2*rsigma2;
     387                assert (z >= 0.0);
     388                psF32 weight  = exp(-z);
     389
     390                wDiff *= weight;
     391                pDiff *= weight;
     392            }
     393
     394            Var += wDiff;
     395            Sum += pDiff;
     396
     397            psF32 xWght = xDiff * pDiff;
     398            psF32 yWght = yDiff * pDiff;
     399
     400            X1  += xWght;
     401            Y1  += yWght;
     402
     403            peakPixel = PS_MAX (*vPix, peakPixel);
     404            numPixels++;
     405        }
     406    }
     407
     408    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
     409    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
     410
     411    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     412    if ((numPixels < minPixels) || (Sum <= 0)) {
     413        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
     414        return (false);
     415    }
     416
     417    // calculate the first moment.
     418    float Mx = X1/Sum;
     419    float My = Y1/Sum;
     420    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
     421        psTrace ("psModules.objects", 3, "large centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
     422        return (false);
     423    }
     424
     425    psTrace ("psModules.objects", 5, "sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", sky, Mx, My, Sum, X1, Y1, numPixels);
     426
     427    // add back offset of peak in primary image
     428    // also offset from pixel index to pixel coordinate
     429    // (the calculation above uses pixel index instead of coordinate)
     430    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
     431
     432    // we only update the centroid if the position is not supplied from elsewhere
     433    bool skipCentroid = false;
     434    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
     435    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
     436
     437    if (skipCentroid) {
     438        source->moments->Mx = source->peak->xf;
     439        source->moments->My = source->peak->yf;
     440    } else {
     441        source->moments->Mx = Mx + xOff + 0.5;
     442        source->moments->My = My + yOff + 0.5;
     443    }
     444
     445    source->moments->Sum = Sum;
     446    source->moments->SN  = Sum / sqrt(Var);
     447    source->moments->Peak = peakPixel;
     448    source->moments->nPixels = numPixels;
     449
     450    return true;
     451}
Note: See TracChangeset for help on using the changeset viewer.