IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 6, 2011, 1:02:53 PM (15 years ago)
Author:
eugene
Message:
  • add concept of saddlePoints to peaks (not actually used in the end)
  • add tmp flags to mark sources for analysis or not in psphotStack
  • autocode the pmSourceIO_CMF_PS1_* functions
  • use 1D gauss approx for convolution in PCM fitting
  • added pmSourceExtFitPars (not actually used)
  • in model guess, use 1st radial moments to define size (if it exists)
  • include PSF_INST_MAG, AP_MAG, KRON_MAG in xfit output
  • fix the position for extended source fits (avoid instability)
  • Sersic-like models (incl. Exp and Dev) use Reff, not sigma; conversion tools need to respect this
  • only use a single pass on the centroid (unwindowed, but limited to 1.5 sigma radius) - this avoids moving the centroid because of nearby neighbors
  • use symmetrical averaging (geometric mean) to calculated 1st radial moment (and avoid neighbor biases), do not use symm. averaging for the flux
  • fix the integration of the sersic, pgauss, and related model functions.
  • fix the central pixel to have the full flux for sersic-like models (interpolated value)
Location:
trunk/psModules/src/objects
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        55*.la
        66*.lo
         7pmSourceIO_CMF_PS1_V1.c
         8pmSourceIO_CMF_PS1_V2.c
         9pmSourceIO_CMF_PS1_V3.c

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

    r31670 r32347  
    9696    // (int) so they can be used in the image index below.
    9797
    98     // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to
    99     // get an unbiased (but probably noisy) centroid
    100     if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
    101         return false;
    102     }
    103     // second pass applies the Gaussian window and uses the centroid from the first pass
    104     if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
     98    // XXX // do 2 passes : the first pass should use a somewhat smaller radius and no sigma window to
     99    // XXX // get an unbiased (but probably noisy) centroid
     100    // XXX if (!pmSourceMomentsGetCentroid (source, 0.75*radius, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
     101    // XXX      return false;
     102    // XXX }
     103    // XXX // second pass applies the Gaussian window and uses the centroid from the first pass
     104    // XXX if (!pmSourceMomentsGetCentroid (source, radius, sigma, minSN, maskVal, source->moments->Mx, source->moments->My)) {
     105    // XXX      return false;
     106    // XXX }
     107
     108    // If we use a large radius for the centroid, it will be biased by any neighbors.  The flux
     109    // of any object drops pretty quickly outside 1-2 sigmas.  (The exception is bright
     110    // saturated stars, for which we need to use a very large radius here)
     111    if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
    105112        return false;
    106113    }
    … …  
    108115    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
    109116    // Xn  = SUM (x - xc)^n * (z - sky)
    110 
    111     float RFW = 0.0;
    112     float RHW = 0.0;
    113 
    114     float RF = 0.0;
    115     float RH = 0.0;
    116     float RS = 0.0;
    117117    float XX = 0.0;
    118118    float XY = 0.0;
    … …  
    143143    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    144144
     145    // calculate the higher-order moments using Xo,Yo
    145146    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    146147
    … …  
    194195            Sum += pDiff;
    195196
    196             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
    197197            float r = sqrt(r2);
    198             float rf = r * fDiff;
    199             float rh = sqrt(r) * fDiff;
    200             float rs = fDiff;
    201 
    202             float rfw = r * pDiff;
    203             float rhw = sqrt(r) * pDiff;
    204198
    205199            float x = xDiff * pDiff;
    … …  
    221215            float yyyy = yDiff * yyy / r2;
    222216
    223             RF  += rf;
    224             RH  += rh;
    225             RS  += rs;
    226 
    227             RFW  += rfw;
    228             RHW  += rhw;
    229 
    230217            XX  += xx;
    231218            XY  += xy;
    … …  
    242229            XYYY  += xyyy;
    243230            YYYY  += yyyy;
     231
     232            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
     233            // XXX float r = sqrt(r2);
     234            // XXX float rf = r * fDiff;
     235            // XXX float rh = sqrt(r) * fDiff;
     236            // XXX float rs = fDiff;
     237            // XXX
     238            // XXX float rfw = r * pDiff;
     239            // XXX float rhw = sqrt(r) * pDiff;
     240            // XXX
     241            // XXX RF  += rf;
     242            // XXX RH  += rh;
     243            // XXX RS  += rs;
     244            // XXX
     245            // XXX RFW  += rfw;
     246            // XXX RHW  += rhw;
    244247        }
    245248    }
    246 
    247     source->moments->Mrf = RF/RS;
    248     source->moments->Mrh = RH/RS;
    249 
    250249    source->moments->Mxx = XX/Sum;
    251250    source->moments->Mxy = XY/Sum;
    … …  
    262261    source->moments->Mxyyy = XYYY/Sum;
    263262    source->moments->Myyyy = YYYY/Sum;
     263
     264    // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging
     265
     266    float **vPix = source->pixels->data.F32;
     267    float **vWgt = source->variance->data.F32;
     268    psImageMaskType  **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     269
     270    float RF = 0.0;
     271    float RH = 0.0;
     272    float RS = 0.0;
     273
     274    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     275
     276        float yDiff = row - yCM;
     277        if (fabs(yDiff) > radius) continue;
     278
     279        // coordinate of mirror pixel
     280        int yFlip = yCM - yDiff;
     281        if (yFlip < 0) continue;
     282        if (yFlip >= source->pixels->numRows) continue;
     283
     284        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     285            // check mask and value for this pixel
     286            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     287            if (isnan(vPix[row][col])) continue;
     288
     289            float xDiff = col - xCM;
     290            if (fabs(xDiff) > radius) continue;
     291
     292            // coordinate of mirror pixel
     293            int xFlip = xCM - xDiff;
     294            if (xFlip < 0) continue;
     295            if (xFlip >= source->pixels->numCols) continue;
     296
     297            // check mask and value for mirror pixel
     298            if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     299            if (isnan(vPix[yFlip][xFlip])) continue;
     300
     301            // radius is just a function of (xDiff, yDiff)
     302            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     303            if (r2 > R2) continue;
     304
     305            float fDiff1 = vPix[row][col] - sky;
     306            float fDiff2 = vPix[yFlip][xFlip] - sky;
     307            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     308
     309            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
     310            float r = sqrt(r2);
     311            float rf = r * pDiff;
     312            float rh = sqrt(r) * pDiff;
     313            float rs = 0.5 * (fDiff1 + fDiff2);
     314
     315            RF  += rf;
     316            RH  += rh;
     317            RS  += rs;
     318        }
     319    }
     320
     321    source->moments->Mrf = RF/RS;
     322    source->moments->Mrh = RH/RS;
    264323
    265324    // if Mrf (first radial moment) is very small, we are getting into low-significance
    … …  
    270329        kronRefRadius = MIN(radius, kronRefRadius);
    271330    }
     331
     332    // *** now calculate the kron flux values using the 1st radial moment
    272333
    273334    float radKinner = 1.0*kronRefRadius;
    … …  
    283344    float SumOuter = 0.0;
    284345
     346    // calculate the Kron flux, and related fluxes (NO symmetrical averaging)
    285347    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    286 
     348       
    287349        float yDiff = row - yCM;
    288350        if (fabs(yDiff) > radKouter) continue;
     351       
     352        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     353            // check mask and value for this pixel
     354            if (vMsk && (vMsk[row][col] & maskVal)) continue;
     355            if (isnan(vPix[row][col])) continue;
     356           
     357            float xDiff = col - xCM;
     358            if (fabs(xDiff) > radKouter) continue;
     359           
     360            // radKron is just a function of (xDiff, yDiff)
     361            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     362
     363            float fDiff1 = vPix[row][col] - sky;
     364            float pDiff = fDiff1;
     365            float wDiff = vWgt[row][col];
     366                                   
     367            // skip pixels below specified significance level.  this is allowed, but should be
     368            // avoided -- the over-weights the wings of bright stars compared to those of faint
     369            // stars.
     370            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
     371           
     372            float r  = sqrt(r2);
     373            if (r < radKron) {
     374                Sum += pDiff;
     375                Var += wDiff;
     376                nKronPix ++;
     377                // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
     378            }
     379
     380            // use sigma (fixed by psf) not a radKron based value
     381            if (r < sigma) {
     382                SumCore += pDiff;
     383                VarCore += wDiff;
     384                nCorePix ++;
     385            }
     386
     387            if ((r > radKinner) && (r < radKron)) {
     388                SumInner += pDiff;
     389                nInner ++;
     390            }
     391            if ((r > radKron)  && (r < radKouter)) {
     392                SumOuter += pDiff;
     393                nOuter ++;
     394            }
     395        }
     396    }
     397    // *** should I rescale these fluxes by pi R^2 / nNpix?
     398    // XXX source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
     399    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
     400    // XXX source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
     401    // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
     402    // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
     403    // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
     404
     405    source->moments->KronCore    = SumCore;
     406    source->moments->KronCoreErr = sqrt(VarCore);
     407    source->moments->KronFlux    = Sum;
     408    source->moments->KronFluxErr = sqrt(Var);
     409    source->moments->KronFinner = SumInner;
     410    source->moments->KronFouter = SumOuter;
     411
     412    // XXX not sure I should save this here...
     413    source->moments->KronFluxPSF    = source->moments->KronFlux;
     414    source->moments->KronFluxPSFErr = source->moments->KronFluxErr;
     415    source->moments->KronRadiusPSF  = source->moments->Mrf;
     416
     417    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",
     418             source->moments->Mrf,   source->moments->KronFlux,
     419             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
     420             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
     421             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     422
     423    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",
     424             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);
     425
     426    return(true);
     427}
     428
     429bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
     430
     431    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
     432    // Sum = SUM (z - sky)
     433    // X1  = SUM (x - xc)*(z - sky)
     434    // .. etc
     435
     436    float sky = 0.0;
     437
     438    float peakPixel = -PS_MAX_F32;
     439    psS32 numPixels = 0;
     440    float Sum = 0.0;
     441    float Var = 0.0;
     442    float X1 = 0.0;
     443    float Y1 = 0.0;
     444    float R2 = PS_SQR(radius);
     445    float minSN2 = PS_SQR(minSN);
     446    float rsigma2 = 0.5 / PS_SQR(sigma);
     447
     448    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
     449    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
     450
     451    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     452    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     453    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     454    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     455    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
     456
     457    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
     458    // not depend on the fractional pixel location of the source.  However, the aperture
     459    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
     460    // position of the expected centroid
     461
     462    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     463
     464        float yDiff = row + 0.5 - yPeak;
     465        if (fabs(yDiff) > radius) continue;
    289466
    290467        float *vPix = source->pixels->data.F32[row];
    291468        float *vWgt = source->variance->data.F32[row];
    292469
    293         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    294         // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
     470        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     471        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
    295472
    296473        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    … …  
    304481            if (isnan(*vPix)) continue;
    305482
    306             float xDiff = col - xCM;
    307             if (fabs(xDiff) > radKouter) continue;
    308 
    309             // radKron is just a function of (xDiff, yDiff)
    310             float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    311 
    312             float pDiff = *vPix - sky;
    313             float wDiff = *vWgt;
    314 
    315             // skip pixels below specified significance level.  this is allowed, but should be
    316             // avoided -- the over-weights the wings of bright stars compared to those of faint
    317             // stars.
    318             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    319 
    320             float r  = sqrt(r2);
    321             if (r < radKron) {
    322                 Sum += pDiff;
    323                 Var += wDiff;
    324                 nKronPix ++;
    325                 // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
    326             }
    327 
    328             // use sigma (fixed by psf) not a radKron based value
    329             if (r < sigma) {
    330                 SumCore += pDiff;
    331                 VarCore += wDiff;
    332                 nCorePix ++;
    333             }
    334 
    335             if ((r > radKinner) && (r < radKron)) {
    336                 SumInner += pDiff;
    337                 nInner ++;
    338             }
    339             if ((r > radKron)  && (r < radKouter)) {
    340                 SumOuter += pDiff;
    341                 nOuter ++;
    342             }
    343         }
    344     }
    345     // *** should I rescale these fluxes by pi R^2 / nNpix?
    346     source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
    347     source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
    348     source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
    349     source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
    350     source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
    351     source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
    352 
    353     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",
    354              source->moments->Mrf,   source->moments->KronFlux,
    355              source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
    356              source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
    357              source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
    358 
    359     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",
    360              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);
    361 
    362     return(true);
    363 }
    364 
    365 bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
    366 
    367     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    368     // Sum = SUM (z - sky)
    369     // X1  = SUM (x - xc)*(z - sky)
    370     // .. etc
    371 
    372     float sky = 0.0;
    373 
    374     float peakPixel = -PS_MAX_F32;
    375     psS32 numPixels = 0;
    376     float Sum = 0.0;
    377     float Var = 0.0;
    378     float X1 = 0.0;
    379     float Y1 = 0.0;
    380     float R2 = PS_SQR(radius);
    381     float minSN2 = PS_SQR(minSN);
    382     float rsigma2 = 0.5 / PS_SQR(sigma);
    383 
    384     float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
    385     float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
    386 
    387     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
    388     // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
    389     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    390     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    391     // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    392 
    393     // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
    394     // not depend on the fractional pixel location of the source.  However, the aperture
    395     // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
    396     // position of the expected centroid
    397 
    398     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    399 
    400         float yDiff = row + 0.5 - yPeak;
    401         if (fabs(yDiff) > radius) continue;
    402 
    403         float *vPix = source->pixels->data.F32[row];
    404         float *vWgt = source->variance->data.F32[row];
    405 
    406         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    407         // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
    408 
    409         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    410             if (vMsk) {
    411                 if (*vMsk & maskVal) {
    412                     vMsk++;
    413                     continue;
    414                 }
    415                 vMsk++;
    416             }
    417             if (isnan(*vPix)) continue;
    418 
    419483            float xDiff = col + 0.5 - xPeak;
    420484            if (fabs(xDiff) > radius) continue;
Note: See TracChangeset for help on using the changeset viewer.