IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 32216


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

Location:
branches/eam_branches/ipp-20110710
Files:
2 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;

  • branches/eam_branches/ipp-20110710/psphot/src/psphotKronIterate.c

    r32205 r32216  
    1010{
    1111    bool status = true;
     12
     13    // return true;
    1214
    1315    // select the appropriate recipe information
    … …  
    127129    psphotSaveImage (NULL, kronWindow, "kron.window.v0.fits");
    128130
    129     for (int i = 0; i < sources->n; i++) {
    130 
    131         pmSource *source = sources->data[i];
    132         if (!source->peak) continue; // XXX how can we have a peak-less source?
    133 
    134         // allocate space for moments
    135         if (!source->moments) continue;
    136 
    137         // replace object in image
    138         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    139             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     131    for (int j = 0; j < 5; j++) {
     132        for (int i = 0; i < sources->n; i++) {
     133
     134            pmSource *source = sources->data[i];
     135            if (!source->peak) continue; // XXX how can we have a peak-less source?
     136
     137            // allocate space for moments
     138            if (!source->moments) continue;
     139
     140            // replace object in image
     141            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     142                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     143            }
     144
     145            // iterate to the window radius
     146            float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
     147
     148            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
     149            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
     150
     151            // clear the window function for this source based on the moments
     152            psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
     153            // psphotKronWindowSetSource (source, kronWindow, false, false);
     154            // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
     155
     156            // 165, 539;
     157            if ((fabs(source->peak->xf - 165) < 3) && (fabs(source->peak->yf - 539) < 3)) {
     158                fprintf (stderr, "test obj\n");
     159            }
     160
     161            // this function populates moments->Mrf,KronFlux,KronFluxErr
     162            psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
     163            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     164
     165            // set a window function for each source based on the moments
     166            psphotKronWindowSetSource (source, kronWindow, true, true);
     167            // psphotKronWindowSetSource (source, kronWindow, false, true);
     168
     169            // test source fluxes
     170            pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
     171            float kmag = -2.5*log10(source->moments->KronFlux);
     172# define TEST_X1 167
     173# define TEST_Y1 299
     174# define TEST_X2 180
     175# define TEST_Y2 300
     176            if ((fabs(source->peak->xf - TEST_X1) < 3) && (fabs(source->peak->yf - TEST_Y1) < 3)) {
     177                fprintf (stderr, "R1: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
     178            }
     179            if ((fabs(source->peak->xf - TEST_X2) < 3) && (fabs(source->peak->yf - TEST_Y2) < 3)) {
     180                fprintf (stderr, "R2: %f vs %f  (%f) (%f)\n", source->moments->KronRadiusPSF, source->moments->Mrf, kmag, windowRadius);
     181            }
     182
     183            // re-subtract the object, leave local sky
     184            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    140185        }
    141 
    142         // iterate to the window radius
    143         float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
    144 
    145         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    146         pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    147 
    148         // clear the window function for this source based on the moments
    149         psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
    150         // psphotKronWindowSetSource (source, kronWindow, false, false);
    151         // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
    152 
    153         // this function populates moments->Mrf,KronFlux,KronFluxErr
    154         psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    155         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    156 
    157         // set a window function for each source based on the moments
    158         psphotKronWindowSetSource (source, kronWindow, true, true);
    159         // psphotKronWindowSetSource (source, kronWindow, false, true);
    160 
    161         // test source fluxes
    162         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    163         float kmag = -2.5*log10(source->moments->KronFlux);
    164         if (source->psfMag - kmag > 0.25) {
    165             // fprintf (stderr, "continue\n");
    166         }
    167 
    168         // re-subtract the object, leave local sky
    169         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    170     }
    171 
    172     psphotSaveImage (NULL, kronWindow, "kron.window.v1.fits");
    173 
    174     for (int i = 0; i < sources->n; i++) {
    175 
    176         pmSource *source = sources->data[i];
    177         if (!source->peak) continue; // XXX how can we have a peak-less source?
    178 
    179         // allocate space for moments
    180         if (!source->moments) continue;
    181 
    182         // replace object in image
    183         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    184             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    185         }
    186 
    187         // iterate to the window radius
    188         float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), EXT_FIT_MAX_RADIUS);
    189 
    190         // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    191         pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    192 
    193         // clear the window function for this source based on the moments
    194         psphotKronWindowSetSource (source, kronWindow, (i > 0), false);
    195         // psphotKronWindowSetSource (source, kronWindow, false, false);
    196         // psphotVisualRangeImage (kapa, kronWindow, "kronwin", 1, 0.0, 1.0);
    197 
    198         // this function populates moments->Mrf,KronFlux,KronFluxErr
    199         psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    200         psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    201 
    202         // set a window function for each source based on the moments
    203         psphotKronWindowSetSource (source, kronWindow, true, true);
    204         // psphotKronWindowSetSource (source, kronWindow, false, true);
    205 
    206         // test source fluxes
    207         pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius);
    208         float kmag = -2.5*log10(source->moments->KronFlux);
    209         if (source->psfMag - kmag > 0.25) {
    210             // fprintf (stderr, "continue\n");
    211         }
    212 
    213         // re-subtract the object, leave local sky
    214         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    215     }
    216 
    217     psphotSaveImage (NULL, kronWindow, "kron.window.v2.fits");
     186        char name[64];
     187        sprintf (name, "kron.window.v%d.fits", j+1);
     188        psphotSaveImage (NULL, kronWindow, name);
     189    }
    218190    psFree (kronWindow);
    219191
    … …  
    231203
    232204    psF32 R2 = PS_SQR(radius);
     205    float rsigma2 = 0.5 / PS_SQR(radius/2.0);
    233206
    234207    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    261234    int Ywo = source->pixels->row0;
    262235
     236    psF32 **vPix = source->pixels->data.F32;
     237    psF32 **vWin = kronWindow->data.F32;
     238    psF32 **vWgt = source->variance->data.F32;
     239   
     240    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     241
    263242    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    264243
    … …  
    266245        if (fabs(yDiff) > radius) continue;
    267246
    268         psF32 *vPix = source->pixels->data.F32[row];
    269         psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
    270 
    271         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    272 
    273         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWin++) {
    274             if (vMsk) {
    275                 if (*vMsk & maskVal) {
    276                     vMsk++;
    277                     continue;
    278                 }
    279                 vMsk++;
    280             }
    281             if (isnan(*vPix)) continue;
     247        // coordinate of mirror pixel
     248        int yFlip = yCM - yDiff;
     249        if (yFlip < 0) continue;
     250        if (yFlip >= source->pixels->numRows) continue;
     251
     252        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     253            // check mask and value for this pixel
     254            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     255            if (isnan(vPix[row][col])) continue;
    282256
    283257            psF32 xDiff = col - xCM;
    284258            if (fabs(xDiff) > radius) continue;
     259
     260            // coordinate of mirror pixel
     261            int xFlip = xCM - xDiff;
     262            if (xFlip < 0) continue;
     263            if (xFlip >= source->pixels->numCols) continue;
     264
     265            // check mask and value for mirror pixel
     266            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     267            if (isnan(vPix[yFlip][xFlip])) continue;
    285268
    286269            // radius is just a function of (xDiff, yDiff)
    … …  
    289272
    290273            // flux * window
    291             float weight  = *vWin;
     274            float z = r2 * rsigma2;
     275            assert (z >= 0.0);
     276
    292277            // float weight  = 1.0;
    293             psF32 pDiff = *vPix * weight;
     278            float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
     279            float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
     280            // float weight1  = vWin[row+Ywo][col+Xwo];
     281            // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
     282
     283            float fDiff1 = vPix[row][col]*weight1;
     284            float fDiff2 = vPix[yFlip][xFlip]*weight2;
     285
     286            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    294287
    295288            // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    296289            psF32 rf = pDiff * sqrt(r2);
    297             psF32 rs = pDiff;
     290            psF32 rs = 0.5 * (fDiff1 + fDiff2);
    298291
    299292            RF  += rf;
    … …  
    322315        if (fabs(yDiff) > radKron) continue;
    323316
    324         psF32 *vPix = source->pixels->data.F32[row];
    325         psF32 *vWgt = source->variance->data.F32[row];
    326         psF32 *vWin = &kronWindow->data.F32[row + Ywo][Xwo];
    327 
    328         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    329 
    330         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++, vWin++) {
    331             if (vMsk) {
    332                 if (*vMsk & maskVal) {
    333                     vMsk++;
    334                     continue;
    335                 }
    336                 vMsk++;
    337             }
    338             if (isnan(*vPix)) continue;
     317        // coordinate of mirror pixel
     318        int yFlip = yCM - yDiff;
     319        if (yFlip < 0) continue;
     320        if (yFlip >= source->pixels->numRows) continue;
     321
     322        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     323            // check mask and value for this pixel
     324            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     325            if (isnan(vPix[row][col])) continue;
    339326
    340327            psF32 xDiff = col - xCM;
    341328            if (fabs(xDiff) > radKron) continue;
     329
     330            // coordinate of mirror pixel
     331            int xFlip = xCM - xDiff;
     332            if (xFlip < 0) continue;
     333            if (xFlip >= source->pixels->numCols) continue;
     334
     335            // check mask and value for mirror pixel
     336            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     337            if (isnan(vPix[yFlip][xFlip])) continue;
    342338
    343339            // radKron is just a function of (xDiff, yDiff)
    … …  
    345341            if (r2 > radKron2) continue;
    346342
    347             float weight  = *vWin;
     343            // float z = r2 * rsigma2;
     344            // assert (z >= 0.0);
     345
    348346            // float weight  = 1.0;
    349             psF32 pDiff = *vPix * weight;
    350             psF32 wDiff = *vWgt * weight;
     347            // float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
     348            // float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
     349            float weight1  = vWin[row+Ywo][col+Xwo];
     350            float weight2  = vWin[yFlip+Ywo][xFlip+Xwo];
     351
     352            float fDiff1 = vPix[row][col]*weight1;
     353            float fDiff2 = vPix[yFlip][xFlip]*weight2;
     354
     355            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     356            psF32 wDiff = vWgt[row][col] * weight1;
    351357
    352358            Sum += pDiff;
    353359            Var += wDiff;
    354             Win += weight;
     360            Win += weight1;
    355361            nKronPix ++;
    356362        }
    … …  
    387393    float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    388394
     395    // float kratio = source->moments->KronFinner / source->moments->KronFlux;
     396
     397    float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
    389398    // float scale = useKronRadius ? 2.0 * source->moments->Mrf / Mmajor : 2.0;
    390     float scale = 3.0 * source->moments->Mrf / Mmajor;
    391     // float scale = 2.0;
     399    // float scale = (kratio > 0.4) ? 9.0 * source->moments->Mrf / Mmajor : 3.0 * source->moments->Mrf / Mmajor;
    392400
    393401    float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
Note: See TracChangeset for help on using the changeset viewer.