IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 29, 2011, 1:26:12 PM (15 years ago)
Author:
eugene
Message:

improvements (?) to the kron flux analysis : attempt to window neighbors based on the symmetry of the object

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r32202 r32216  
    8080    }
    8181
    82     if (source->moments->nPixels != 0) {
    83         fprintf (stderr, "remeasure moments: %f,%f\n", source->peak->xf, source->peak->yf);
    84     }
    85 
    8682    float Sum = 0.0;
    8783    float Var = 0.0;
    … …  
    120116    // Xn  = SUM (x - xc)^n * (z - sky)
    121117
    122     float RFW = 0.0;
    123     float RHW = 0.0;
     118    float RFa = 0.0;
     119    float RSa = 0.0;
    124120
    125121    float RF = 0.0;
    … …  
    154150    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    155151
     152    // calculate the higher-order moments using Xo,Yo
    156153    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    157154
    … …  
    205202            Sum += pDiff;
    206203
    207             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
    208204            float r = sqrt(r2);
    209             float rf = r * fDiff;
    210             float rh = sqrt(r) * fDiff;
    211             float rs = fDiff;
    212 
    213             float rfw = r * pDiff;
    214             float rhw = sqrt(r) * pDiff;
    215205
    216206            float x = xDiff * pDiff;
    … …  
    232222            float yyyy = yDiff * yyy / r2;
    233223
    234             RF  += rf;
    235             RH  += rh;
    236             RS  += rs;
    237 
    238             RFW  += rfw;
    239             RHW  += rhw;
    240 
    241224            XX  += xx;
    242225            XY  += xy;
    … …  
    253236            XYYY  += xyyy;
    254237            YYYY  += yyyy;
     238
     239            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
     240            // XXX float r = sqrt(r2);
     241            // XXX float rf = r * fDiff;
     242            // XXX float rh = sqrt(r) * fDiff;
     243            // XXX float rs = fDiff;
     244            // XXX
     245            // XXX float rfw = r * pDiff;
     246            // XXX float rhw = sqrt(r) * pDiff;
     247            // XXX
     248            // XXX RF  += rf;
     249            // XXX RH  += rh;
     250            // XXX RS  += rs;
     251            // XXX
     252            // XXX RFW  += rfw;
     253            // XXX RHW  += rhw;
    255254        }
    256255    }
    257 
    258     source->moments->Mrf = RF/RS;
    259     source->moments->Mrh = RH/RS;
    260 
    261256    source->moments->Mxx = XX/Sum;
    262257    source->moments->Mxy = XY/Sum;
    … …  
    273268    source->moments->Mxyyy = XYYY/Sum;
    274269    source->moments->Myyyy = YYYY/Sum;
     270
     271# define TEST_X1 167
     272# define TEST_Y1 299
     273# define TEST_X2 180
     274# define TEST_Y2 300
     275    if ((fabs(Xo - TEST_X1) < 3) && (fabs(Yo - TEST_Y1) < 3)) {
     276        fprintf (stderr, "test obj 1\n");
     277    }
     278    if ((fabs(Xo - TEST_X2) < 3) && (fabs(Yo - TEST_Y2) < 3)) {
     279        fprintf (stderr, "test obj 2\n");
     280    }
     281
     282    float **vPix = source->pixels->data.F32;
     283    float **vWgt = source->variance->data.F32;
     284    psImageMaskType  **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     285
     286    // calculate the 1st radial moment (for kron flux) -- symmetrical averaging
     287    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     288
     289        float yDiff = row - yCM;
     290        if (fabs(yDiff) > radius) continue;
     291
     292        // coordinate of mirror pixel
     293        int yFlip = yCM - yDiff;
     294        if (yFlip < 0) continue;
     295        if (yFlip >= source->pixels->numRows) continue;
     296
     297        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     298            // check mask and value for this pixel
     299            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     300            if (isnan(vPix[row][col])) continue;
     301
     302            float xDiff = col - xCM;
     303            if (fabs(xDiff) > radius) continue;
     304
     305            // coordinate of mirror pixel
     306            int xFlip = xCM - xDiff;
     307            if (xFlip < 0) continue;
     308            if (xFlip >= source->pixels->numCols) continue;
     309
     310            // check mask and value for mirror pixel
     311            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     312            if (isnan(vPix[yFlip][xFlip])) continue;
     313
     314            // radius is just a function of (xDiff, yDiff)
     315            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     316            if (r2 > R2) continue;
     317
     318            float fDiff1 = vPix[row][col] - sky;
     319            float fDiff2 = vPix[yFlip][xFlip] - sky;
     320            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     321
     322            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
     323            float r = sqrt(r2);
     324            float rf = r * pDiff;
     325            float rh = sqrt(r) * pDiff;
     326            float rs = 0.5 * (fDiff1 + fDiff2);
     327
     328            float rfa = r * fDiff1;
     329            float rsa = fDiff1;
     330
     331            RF  += rf;
     332            RH  += rh;
     333            RS  += rs;
     334
     335            RFa  += rfa;
     336            RSa  += rsa;
     337        }
     338    }
     339
     340    source->moments->Mrf = RF/RS;
     341    source->moments->Mrh = RH/RS;
     342
     343    float R1 = RFa / RSa;
     344    if ((fabs(Xo - TEST_X1) < 3) && (fabs(Yo - TEST_Y1) < 3)) {
     345        fprintf (stderr, "R1: %f vs %f\n", R1, source->moments->Mrf);
     346    }
     347    if ((fabs(Xo - TEST_X2) < 3) && (fabs(Yo - TEST_Y2) < 3)) {
     348        fprintf (stderr, "R2: %f vs %f\n", R1, source->moments->Mrf);
     349    }
     350
     351    // fprintf (stderr, "Rad: %f vs %f\n", R1, source->moments->Mrf);
    275352
    276353    // if Mrf (first radial moment) is very small, we are getting into low-significance
    … …  
    294371    float SumOuter = 0.0;
    295372
     373    // calculate the Kron flux, and related fluxes (symmetrical averaging)
    296374    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    297 
     375       
    298376        float yDiff = row - yCM;
    299377        if (fabs(yDiff) > radKouter) continue;
     378       
     379        // coordinate of mirror pixel
     380        int yFlip = yCM - yDiff;
     381        if (yFlip < 0) continue;
     382        if (yFlip >= source->pixels->numRows) continue;
     383       
     384        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     385            // check mask and value for this pixel
     386            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     387            if (isnan(vPix[row][col])) continue;
     388           
     389            float xDiff = col - xCM;
     390            if (fabs(xDiff) > radKouter) continue;
     391           
     392            // coordinate of mirror pixel
     393            int xFlip = xCM - xDiff;
     394            if (xFlip < 0) continue;
     395            if (xFlip >= source->pixels->numCols) continue;
     396           
     397            // check mask and value for mirror pixel
     398            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     399            if (isnan(vPix[yFlip][xFlip])) continue;
     400           
     401            // radKron is just a function of (xDiff, yDiff)
     402            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     403
     404            float fDiff1 = vPix[row][col] - sky;
     405            float fDiff2 = vPix[yFlip][xFlip] - sky;
     406            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     407            // float pDiff = vPix[row][col] - sky;
     408            float wDiff = vWgt[row][col];
     409                                   
     410            // skip pixels below specified significance level.  this is allowed, but should be
     411            // avoided -- the over-weights the wings of bright stars compared to those of faint
     412            // stars.
     413            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     414           
     415# define WEIGHTED 0
     416# if (WEIGHTED)
     417            float z = r2 * rsigma2 / 4.0;
     418            assert (z >= 0.0);
     419            float weight  = exp(-z);
     420# else
     421            float weight  = 1.0;
     422# endif
     423
     424            float r  = sqrt(r2);
     425            if (r < radKron) {
     426                Sum += pDiff*weight;
     427                Var += wDiff*weight;
     428                nKronPix ++;
     429                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
     430            }
     431
     432            // use sigma (fixed by psf) not a radKron based value
     433            if (r < sigma) {
     434                SumCore += pDiff;
     435                VarCore += wDiff;
     436                nCorePix ++;
     437            }
     438
     439            if ((r > radKinner) && (r < radKron)) {
     440                SumInner += pDiff;
     441                nInner ++;
     442            }
     443            if ((r > radKron)  && (r < radKouter)) {
     444                SumOuter += pDiff;
     445                nOuter ++;
     446            }
     447        }
     448    }
     449    // *** should I rescale these fluxes by pi R^2 / nNpix?
     450    // XXX source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
     451    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
     452    // XXX source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
     453    // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
     454    // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
     455    // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
     456
     457    source->moments->KronCore    = SumCore;
     458    source->moments->KronCoreErr = sqrt(VarCore);
     459    source->moments->KronFlux    = Sum;
     460    source->moments->KronFluxErr = sqrt(Var);
     461    source->moments->KronFinner = SumInner;
     462    source->moments->KronFouter = SumOuter;
     463
     464    // XXX not sure I should save this here...
     465    source->moments->KronFluxPSF    = source->moments->KronFlux;
     466    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
     467    source->moments->KronRadiusPSF  = source->moments->Mrf;
     468
     469    psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  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",
     470             source->moments->Mrf,   source->moments->KronFlux,
     471             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
     472             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
     473             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     474
     475    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",
     476             source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
     477
     478    return(true);
     479}
     480
     481bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
     482
     483    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
     484    // Sum = SUM (z - sky)
     485    // X1  = SUM (x - xc)*(z - sky)
     486    // .. etc
     487
     488    float sky = 0.0;
     489
     490    float peakPixel = -PS_MAX_F32;
     491    psS32 numPixels = 0;
     492    float Sum = 0.0;
     493    float Var = 0.0;
     494    float X1 = 0.0;
     495    float Y1 = 0.0;
     496    float R2 = PS_SQR(radius);
     497    float minSN2 = PS_SQR(minSN);
     498    float rsigma2 = 0.5 / PS_SQR(sigma);
     499
     500    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
     501    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
     502
     503    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     504    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     505    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     506    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     507    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
     508
     509    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
     510    // not depend on the fractional pixel location of the source.  However, the aperture
     511    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
     512    // position of the expected centroid
     513
     514    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     515
     516        float yDiff = row + 0.5 - yPeak;
     517        if (fabs(yDiff) > radius) continue;
    300518
    301519        float *vPix = source->pixels->data.F32[row];
    302520        float *vWgt = source->variance->data.F32[row];
    303521
    304         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    305         // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
     522        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     523        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
    306524
    307525        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    … …  
    315533            if (isnan(*vPix)) continue;
    316534
    317             float xDiff = col - xCM;
    318             if (fabs(xDiff) > radKouter) continue;
    319 
    320             // radKron is just a function of (xDiff, yDiff)
    321             float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    322 
    323             float pDiff = *vPix - sky;
    324             float wDiff = *vWgt;
    325 
    326             // skip pixels below specified significance level.  this is allowed, but should be
    327             // avoided -- the over-weights the wings of bright stars compared to those of faint
    328             // stars.
    329             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    330 
    331 # define WEIGHTED 1
    332 # if (WEIGHTED)
    333             float z = r2 * rsigma2 / 4.0;
    334             assert (z >= 0.0);
    335             float weight  = exp(-z);
    336 # endif
    337 
    338             float r  = sqrt(r2);
    339             if (r < radKron) {
    340 # if (WEIGHTED)
    341                 Sum += pDiff*weight;
    342                 Var += wDiff*weight;
    343 # else
    344                 Sum += pDiff;
    345                 Var += wDiff;
    346 # endif
    347                 nKronPix ++;
    348                 // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
    349             }
    350 
    351             // use sigma (fixed by psf) not a radKron based value
    352             if (r < sigma) {
    353                 SumCore += pDiff;
    354                 VarCore += wDiff;
    355                 nCorePix ++;
    356             }
    357 
    358             if ((r > radKinner) && (r < radKron)) {
    359                 SumInner += pDiff;
    360                 nInner ++;
    361             }
    362             if ((r > radKron)  && (r < radKouter)) {
    363                 SumOuter += pDiff;
    364                 nOuter ++;
    365             }
    366         }
    367     }
    368     // *** should I rescale these fluxes by pi R^2 / nNpix?
    369     source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
    370     source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
    371     source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
    372     source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
    373     source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
    374     source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
    375 
    376     // XXX not sure I should save this here...
    377     source->moments->KronFluxPSF    = source->moments->KronFlux;
    378     source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
    379     source->moments->KronRadiusPSF  = source->moments->Mrf;
    380 
    381     psTrace ("psModules.objects", 4, "Mrf: %f  KronFlux: %f  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",
    382              source->moments->Mrf,   source->moments->KronFlux,
    383              source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
    384              source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
    385              source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
    386 
    387     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",
    388              source->peak->xf, source->peak->yf, source->peak->rawFlux, sqrt(source->peak->detValue), source->moments->Mx,   source->moments->My, Sum, source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy, sky, source->moments->nPixels);
    389 
    390     return(true);
    391 }
    392 
    393 bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
    394 
    395     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    396     // Sum = SUM (z - sky)
    397     // X1  = SUM (x - xc)*(z - sky)
    398     // .. etc
    399 
    400     float sky = 0.0;
    401 
    402     float peakPixel = -PS_MAX_F32;
    403     psS32 numPixels = 0;
    404     float Sum = 0.0;
    405     float Var = 0.0;
    406     float X1 = 0.0;
    407     float Y1 = 0.0;
    408     float R2 = PS_SQR(radius);
    409     float minSN2 = PS_SQR(minSN);
    410     float rsigma2 = 0.5 / PS_SQR(sigma);
    411 
    412     float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
    413     float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
    414 
    415     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
    416     // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
    417     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    418     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    419     // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    420 
    421     // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
    422     // not depend on the fractional pixel location of the source.  However, the aperture
    423     // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
    424     // position of the expected centroid
    425 
    426     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    427 
    428         float yDiff = row + 0.5 - yPeak;
    429         if (fabs(yDiff) > radius) continue;
    430 
    431         float *vPix = source->pixels->data.F32[row];
    432         float *vWgt = source->variance->data.F32[row];
    433 
    434         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    435         // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
    436 
    437         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    438             if (vMsk) {
    439                 if (*vMsk & maskVal) {
    440                     vMsk++;
    441                     continue;
    442                 }
    443                 vMsk++;
    444             }
    445             if (isnan(*vPix)) continue;
    446 
    447535            float xDiff = col + 0.5 - xPeak;
    448536            if (fabs(xDiff) > radius) continue;
Note: See TracChangeset for help on using the changeset viewer.