IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 24, 2011, 10:09:38 PM (15 years ago)
Author:
eugene
Message:

determine the min kron radius = radius for bright PSF stars; mask interpolation used index instead of pixel coordinates in psImageInterpolate; add minKronRadius to pmSource; change errMag to psfMagErr; in growth curve, avoid perfect int radii (they can be inconsistent on + and - due to float precision; add pmGrowthCurveForSources; psfMagErr is always calculated from psfModel; apply ApTrend to all source psf mags and fluxes (not just psf sources); always use the psfModel for PSF_QF,_PERFECT; use BILINEAR to interpolate for aperture mags (BIQUAD implementation is wrong; apply growth curve to apMag & apFlux for all sources; save GrowthCurve with PSF model and read from PSF model

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r31327 r31361  
    6666// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
    6767
    68 bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
     68bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
    6969{
    7070    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    7474
    7575    // this function assumes the sky has been well-subtracted for the image
    76     psF32 sky = 0.0;
     76    float sky = 0.0;
    7777
    7878    if (source->moments == NULL) {
    … …  
    8080    }
    8181
    82     psF32 Sum = 0.0;
    83     psF32 Var = 0.0;
    84     psF32 SumCore = 0.0;
    85     psF32 VarCore = 0.0;
    86     psF32 R2 = PS_SQR(radius);
    87     psF32 minSN2 = PS_SQR(minSN);
    88     psF32 rsigma2 = 0.5 / PS_SQR(sigma);
     82    float Sum = 0.0;
     83    float Var = 0.0;
     84    float SumCore = 0.0;
     85    float VarCore = 0.0;
     86    float R2 = PS_SQR(radius);
     87    float minSN2 = PS_SQR(minSN);
     88    float rsigma2 = 0.5 / PS_SQR(sigma);
    8989
    9090    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    109109    // Xn  = SUM (x - xc)^n * (z - sky)
    110110
    111     psF32 RFW = 0.0;
    112     psF32 RHW = 0.0;
    113 
    114     psF32 RF = 0.0;
    115     psF32 RH = 0.0;
    116     psF32 RS = 0.0;
    117     psF32 XX = 0.0;
    118     psF32 XY = 0.0;
    119     psF32 YY = 0.0;
    120     psF32 XXX = 0.0;
    121     psF32 XXY = 0.0;
    122     psF32 XYY = 0.0;
    123     psF32 YYY = 0.0;
    124     psF32 XXXX = 0.0;
    125     psF32 XXXY = 0.0;
    126     psF32 XXYY = 0.0;
    127     psF32 XYYY = 0.0;
    128     psF32 YYYY = 0.0;
     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;
     117    float XX = 0.0;
     118    float XY = 0.0;
     119    float YY = 0.0;
     120    float XXX = 0.0;
     121    float XXY = 0.0;
     122    float XYY = 0.0;
     123    float YYY = 0.0;
     124    float XXXX = 0.0;
     125    float XXXY = 0.0;
     126    float XXYY = 0.0;
     127    float XYYY = 0.0;
     128    float YYYY = 0.0;
    129129
    130130    Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
    … …  
    140140    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
    141141    // xCM, yCM from pixel coords to pixel index here.
    142     psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
    143     psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     142    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     143    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    144144
    145145    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    146146
    147         psF32 yDiff = row - yCM;
     147        float yDiff = row - yCM;
    148148        if (fabs(yDiff) > radius) continue;
    149149
    150         psF32 *vPix = source->pixels->data.F32[row];
    151         psF32 *vWgt = source->variance->data.F32[row];
     150        float *vPix = source->pixels->data.F32[row];
     151        float *vWgt = source->variance->data.F32[row];
    152152
    153153        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    … …  
    164164            if (isnan(*vPix)) continue;
    165165
    166             psF32 xDiff = col - xCM;
     166            float xDiff = col - xCM;
    167167            if (fabs(xDiff) > radius) continue;
    168168
    169169            // radius is just a function of (xDiff, yDiff)
    170             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     170            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    171171            if (r2 > R2) continue;
    172172
    173             psF32 fDiff = *vPix - sky;
    174             psF32 pDiff = fDiff;
    175             psF32 wDiff = *vWgt;
     173            float fDiff = *vPix - sky;
     174            float pDiff = fDiff;
     175            float wDiff = *vWgt;
    176176
    177177            // skip pixels below specified significance level.  this is allowed, but should be
    … …  
    184184            // weighting over weights the sky for faint sources
    185185            if (sigma > 0.0) {
    186                 psF32 z = r2 * rsigma2;
     186                float z = r2 * rsigma2;
    187187                assert (z >= 0.0);
    188                 psF32 weight  = exp(-z);
     188                float weight  = exp(-z);
    189189
    190190                wDiff *= weight;
    … …  
    195195
    196196            // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
    197             psF32 r = sqrt(r2);
    198             psF32 rf = r * fDiff;
    199             psF32 rh = sqrt(r) * fDiff;
    200             psF32 rs = fDiff;
    201 
    202             psF32 rfw = r * pDiff;
    203             psF32 rhw = sqrt(r) * pDiff;
    204 
    205             psF32 x = xDiff * pDiff;
    206             psF32 y = yDiff * pDiff;
    207 
    208             psF32 xx = xDiff * x;
    209             psF32 xy = xDiff * y;
    210             psF32 yy = yDiff * y;
    211 
    212             psF32 xxx = xDiff * xx / r;
    213             psF32 xxy = xDiff * xy / r;
    214             psF32 xyy = xDiff * yy / r;
    215             psF32 yyy = yDiff * yy / r;
    216 
    217             psF32 xxxx = xDiff * xxx / r2;
    218             psF32 xxxy = xDiff * xxy / r2;
    219             psF32 xxyy = xDiff * xyy / r2;
    220             psF32 xyyy = xDiff * yyy / r2;
    221             psF32 yyyy = yDiff * yyy / r2;
     197            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;
     204
     205            float x = xDiff * pDiff;
     206            float y = yDiff * pDiff;
     207
     208            float xx = xDiff * x;
     209            float xy = xDiff * y;
     210            float yy = yDiff * y;
     211
     212            float xxx = xDiff * xx / r;
     213            float xxy = xDiff * xy / r;
     214            float xyy = xDiff * yy / r;
     215            float yyy = yDiff * yy / r;
     216
     217            float xxxx = xDiff * xxx / r2;
     218            float xxxy = xDiff * xxy / r2;
     219            float xxyy = xDiff * xyy / r2;
     220            float xyyy = xDiff * yyy / r2;
     221            float yyyy = yDiff * yyy / r2;
    222222
    223223            RF  += rf;
    … …  
    245245    }
    246246
    247     // if Mrf (first radial moment) is << sigma, we are getting into low-significance
    248     // territory.  saturate at 0.75*sigma.  conversely, if Mrf is > radius, we are clearly
    249     // making an error.  saturate at radius.
    250     source->moments->Mrf = MIN(radius, MAX(0.75*sigma, RF/RS));
    251     source->moments->Mrh = MIN(radius, MAX(0.75*sigma, RH/RS));
     247    source->moments->Mrf = RF/RS;
     248    source->moments->Mrh = RH/RS;
    252249
    253250    source->moments->Mxx = XX/Sum;
    … …  
    266263    source->moments->Myyyy = YYYY/Sum;
    267264
    268     // XXX TEST:
    269     // source->moments->KronFinner = RFW/Sum;
    270     // source->moments->KronFouter = sigma;
    271 
    272     // Calculate the Kron magnitude (make this block optional?)
    273     float radKinner = 1.0*source->moments->Mrf;
    274     float radKron   = 2.5*source->moments->Mrf;
    275     float radKouter = 4.0*source->moments->Mrf;
     265    // if Mrf (first radial moment) is very small, we are getting into low-significance
     266    // territory.  saturate at minKronRadius.  conversely, if Mrf is > radius, we are clearly
     267    // making an error.  saturate at radius.
     268    float kronRefRadius = MIN(radius, MAX(minKronRadius, source->moments->Mrf));
     269
     270    float radKinner = 1.0*kronRefRadius;
     271    float radKron   = 2.5*kronRefRadius;
     272    float radKouter = 4.0*kronRefRadius;
    276273
    277274    int nKronPix = 0;
    … …  
    285282    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    286283
    287         psF32 yDiff = row - yCM;
     284        float yDiff = row - yCM;
    288285        if (fabs(yDiff) > radKouter) continue;
    289286
    290         psF32 *vPix = source->pixels->data.F32[row];
    291         psF32 *vWgt = source->variance->data.F32[row];
     287        float *vPix = source->pixels->data.F32[row];
     288        float *vWgt = source->variance->data.F32[row];
    292289
    293290        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    … …  
    304301            if (isnan(*vPix)) continue;
    305302
    306             psF32 xDiff = col - xCM;
     303            float xDiff = col - xCM;
    307304            if (fabs(xDiff) > radKouter) continue;
    308305
    309306            // radKron is just a function of (xDiff, yDiff)
    310             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    311 
    312             psF32 pDiff = *vPix - sky;
    313             psF32 wDiff = *vWgt;
     307            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     308
     309            float pDiff = *vPix - sky;
     310            float wDiff = *vWgt;
    314311
    315312            // skip pixels below specified significance level.  this is allowed, but should be
    … …  
    318315            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    319316
    320             psF32 r  = sqrt(r2);
     317            float r  = sqrt(r2);
    321318            if (r < radKron) {
    322319                Sum += pDiff;
    … …  
    363360}
    364361
    365 bool pmSourceMomentsGetCentroid(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
     362bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
    366363
    367364    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    … …  
    370367    // .. etc
    371368
    372     psF32 sky = 0.0;
    373 
    374     psF32 peakPixel = -PS_MAX_F32;
     369    float sky = 0.0;
     370
     371    float peakPixel = -PS_MAX_F32;
    375372    psS32 numPixels = 0;
    376     psF32 Sum = 0.0;
    377     psF32 Var = 0.0;
    378     psF32 X1 = 0.0;
    379     psF32 Y1 = 0.0;
    380     psF32 R2 = PS_SQR(radius);
    381     psF32 minSN2 = PS_SQR(minSN);
    382     psF32 rsigma2 = 0.5 / PS_SQR(sigma);
     373    float Sum = 0.0;
     374    float Var = 0.0;
     375    float X1 = 0.0;
     376    float Y1 = 0.0;
     377    float R2 = PS_SQR(radius);
     378    float minSN2 = PS_SQR(minSN);
     379    float rsigma2 = 0.5 / PS_SQR(sigma);
    383380
    384381    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
    … …  
    386383
    387384    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
    388     // psF32 weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     385    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
    389386    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    390387    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    … …  
    398395    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    399396
    400         psF32 yDiff = row + 0.5 - yPeak;
     397        float yDiff = row + 0.5 - yPeak;
    401398        if (fabs(yDiff) > radius) continue;
    402399
    403         psF32 *vPix = source->pixels->data.F32[row];
    404         psF32 *vWgt = source->variance->data.F32[row];
     400        float *vPix = source->pixels->data.F32[row];
     401        float *vWgt = source->variance->data.F32[row];
    405402
    406403        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    … …  
    417414            if (isnan(*vPix)) continue;
    418415
    419             psF32 xDiff = col + 0.5 - xPeak;
     416            float xDiff = col + 0.5 - xPeak;
    420417            if (fabs(xDiff) > radius) continue;
    421418
    422419            // radius is just a function of (xDiff, yDiff)
    423             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     420            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    424421            if (r2 > R2) continue;
    425422
    426             psF32 pDiff = *vPix - sky;
    427             psF32 wDiff = *vWgt;
     423            float pDiff = *vPix - sky;
     424            float wDiff = *vWgt;
    428425
    429426            // skip pixels below specified significance level.  for a PSFs, this
    … …  
    438435            // weighting over weights the sky for faint sources
    439436            if (sigma > 0.0) {
    440                 psF32 z  = r2*rsigma2;
     437                float z  = r2*rsigma2;
    441438                assert (z >= 0.0);
    442                 psF32 weight  = exp(-z);
     439                float weight  = exp(-z);
    443440
    444441                wDiff *= weight;
    … …  
    449446            Sum += pDiff;
    450447
    451             psF32 xWght = xDiff * pDiff;
    452             psF32 yWght = yDiff * pDiff;
     448            float xWght = xDiff * pDiff;
     449            float yWght = yDiff * pDiff;
    453450
    454451            X1  += xWght;
    … …  
    507504    return true;
    508505}
     506
     507float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
     508
     509    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
     510
     511    for (int i = 0; i < sources->n; i++) {
     512        pmSource *src = sources->data[i]; // Source of interest
     513        if (!src || !src->moments) {
     514            continue;
     515        }
     516
     517        if (src->mode & PM_SOURCE_MODE_BLEND) {
     518            continue;
     519        }
     520
     521        if (!src->moments->nPixels) continue;
     522
     523        if (src->moments->SN < PSF_SN_LIM) continue;
     524
     525        // XXX put in Mxx,Myy cut based on clump location
     526
     527        psVectorAppend(radii, src->moments->Mrf);
     528    }
     529
     530    // find the peak in this image
     531    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     532
     533    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
     534        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
     535        psFree(stats);
     536        return NAN;
     537    }
     538
     539    float minRadius = stats->sampleMedian;
     540
     541    psFree(radii);
     542    psFree(stats);
     543    return minRadius;
     544}
     545
Note: See TracChangeset for help on using the changeset viewer.