IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 10, 2013, 2:55:11 PM (12 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20130904

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psModules

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

    r35560 r36375  
    6565void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
    6666
     67bool pmSourceMomentsHighOrder    (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal);
     68bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
     69bool pmSourceMomentsKronFluxes   (pmSource *source, float sigma,  float minSN, psImageMaskType maskVal);
     70
    6771// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
    68 
    6972bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
    7073{
    … …  
    7477    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    7578
    76     // this function assumes the sky has been well-subtracted for the image
    77     float sky = 0.0;
    78 
    7979    if (source->moments == NULL) {
    8080      source->moments = pmMomentsAlloc();
    8181    }
    82 
    83     float Sum = 0.0;
    84     float Var = 0.0;
    85     float SumCore = 0.0;
    86     float VarCore = 0.0;
    87     float R2 = PS_SQR(radius);
    88     float minSN2 = PS_SQR(minSN);
    89     float rsigma2 = 0.5 / PS_SQR(sigma);
    9082
    9183    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    110102    // of any object drops pretty quickly outside 1-2 sigmas.  (The exception is bright
    111103    // saturated stars, for which we need to use a very large radius here)
     104    // NOTE: if (mode & EXTERNAL) or (mode2 & MATCHED), do not re-calculate the centroid (use peak as centroid)
     105    // (we still call this function because it sets moments->Sum,SN,Peak,nPixels
    112106    if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
    113107        return false;
    114108    }
    115109
     110    pmSourceMomentsHighOrder (source, radius, sigma, minSN, maskVal);
     111
     112    // now calculate the 1st radial moment (for kron flux) using symmetrical averaging
     113    pmSourceMomentsRadialMoment (source, radius, minKronRadius, maskVal);
     114
     115    // now calculate the kron flux values using the 1st radial moment
     116    pmSourceMomentsKronFluxes (source, sigma, minSN, maskVal);
     117
     118    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",
     119             source->moments->Mrf,   source->moments->KronFlux,
     120             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
     121             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
     122             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     123
     124    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  Npix: %d\n",
     125             source->peak->xf, source->peak->yf,
     126             source->peak->rawFlux, sqrt(source->peak->detValue),
     127             source->moments->Mx, source->moments->My,
     128             source->moments->Sum,
     129             source->moments->Mxx, source->moments->Mxy, source->moments->Myy,
     130             source->moments->nPixels);
     131
     132    return(true);
     133}
     134
     135bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
     136
     137    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
     138    // Sum = SUM (z - sky)
     139    // X1  = SUM (x - xc)*(z - sky)
     140    // .. etc
     141
     142    float sky = 0.0;
     143
     144    float peakPixel = -PS_MAX_F32;
     145    psS32 numPixels = 0;
     146    float Sum = 0.0;
     147    float Var = 0.0;
     148    float X1 = 0.0;
     149    float Y1 = 0.0;
     150    float R2 = PS_SQR(radius);
     151    float minSN2 = PS_SQR(minSN);
     152    float rsigma2 = 0.5 / PS_SQR(sigma);
     153
     154    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
     155    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
     156
     157    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     158    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     159    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     160    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     161    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
     162
     163    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
     164    // not depend on the fractional pixel location of the source.  However, the aperture
     165    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
     166    // position of the expected centroid
     167
     168    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     169
     170        float yDiff = row + 0.5 - yPeak;
     171        if (fabs(yDiff) > radius) continue;
     172
     173        float *vPix = source->pixels->data.F32[row];
     174        float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
     175
     176        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     177        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
     178
     179        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     180            if (vMsk) {
     181                if (*vMsk & maskVal) {
     182                    vMsk++;
     183                    continue;
     184                }
     185                vMsk++;
     186            }
     187            if (isnan(*vPix)) continue;
     188
     189            float xDiff = col + 0.5 - xPeak;
     190            if (fabs(xDiff) > radius) continue;
     191
     192            // radius is just a function of (xDiff, yDiff)
     193            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     194            if (r2 > R2) continue;
     195
     196            float pDiff = *vPix - sky;
     197            float wDiff = *vWgt;
     198
     199            // skip pixels below specified significance level.  for a PSFs, this
     200            // over-weights the wings of bright stars compared to those of faint stars.
     201            // for the estimator used for extended source analysis (where the window
     202            // function is allowed to be arbitrarily large), we need to clip to avoid
     203            // negative second moments.
     204            if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
     205            if ((minSN > 0.0) && (pDiff < 0)) continue; //
     206
     207            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     208            // weighting over weights the sky for faint sources
     209            if (sigma > 0.0) {
     210                float z  = r2*rsigma2;
     211                assert (z >= 0.0);
     212                float weight  = exp(-z);
     213
     214                wDiff *= weight;
     215                pDiff *= weight;
     216            }
     217
     218            Var += wDiff;
     219            Sum += pDiff;
     220
     221            float xWght = xDiff * pDiff;
     222            float yWght = yDiff * pDiff;
     223
     224            X1  += xWght;
     225            Y1  += yWght;
     226
     227            peakPixel = PS_MAX (*vPix, peakPixel);
     228            numPixels++;
     229        }
     230    }
     231
     232    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
     233    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
     234
     235    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     236    if ((numPixels < minPixels) || (Sum <= 0)) {
     237        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
     238        return (false);
     239    }
     240
     241    // calculate the first moment.
     242    float Mx = X1/Sum;
     243    float My = Y1/Sum;
     244    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
     245        psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
     246        return (false);
     247    }
     248    if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) {
     249        psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y);
     250    }
     251
     252    psTrace ("psModules.objects", 5, "id: %d, sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels);
     253
     254    // add back offset of peak in primary image
     255    // also offset from pixel index to pixel coordinate
     256    // (the calculation above uses pixel index instead of coordinate)
     257    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
     258
     259    // we only update the centroid if the position is not supplied from elsewhere
     260    bool skipCentroid = false;
     261    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
     262    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
     263
     264    if (skipCentroid) {
     265        source->moments->Mx = source->peak->xf;
     266        source->moments->My = source->peak->yf;
     267    } else {
     268        source->moments->Mx = Mx + xGuess;
     269        source->moments->My = My + yGuess;
     270    }
     271
     272    source->moments->Sum = Sum;
     273    source->moments->SN  = Sum / sqrt(Var);
     274    source->moments->Peak = peakPixel;
     275    source->moments->nPixels = numPixels;
     276
     277    return true;
     278}
     279
     280float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
     281
     282    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
     283
     284    for (int i = 0; i < sources->n; i++) {
     285        pmSource *src = sources->data[i]; // Source of interest
     286        if (!src || !src->moments) {
     287            continue;
     288        }
     289
     290        if (src->mode & PM_SOURCE_MODE_BLEND) {
     291            continue;
     292        }
     293
     294        if (!src->moments->nPixels) continue;
     295
     296        if (src->moments->SN < PSF_SN_LIM) continue;
     297
     298        // XXX put in Mxx,Myy cut based on clump location
     299
     300        psVectorAppend(radii, src->moments->Mrf);
     301    }
     302
     303    // find the peak in this image
     304    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     305
     306    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
     307        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
     308        psFree(stats);
     309        return NAN;
     310    }
     311
     312    float minRadius = stats->sampleMedian;
     313
     314    psFree(radii);
     315    psFree(stats);
     316    return minRadius;
     317}
     318
     319bool pmSourceMomentsHighOrder (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal) {
     320
     321    // this function assumes the sky has been well-subtracted for the image
     322    float Sum = 0.0;
     323    float R2 = PS_SQR(radius);
     324    float minSN2 = PS_SQR(minSN);
     325    float rsigma2 = 0.5 / PS_SQR(sigma);
     326
    116327    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
    117     // Xn  = SUM (x - xc)^n * (z - sky)
     328    // Xn  = SUM (x - xc)^n * (z - sky) -- note that sky is 0.0 by definition here
    118329    float XX = 0.0;
    119330    float XY = 0.0;
    … …  
    129340    float YYYY = 0.0;
    130341
    131     Sum = 0.0;  // the second pass may include slightly different pixels, re-determine Sum
    132 
    133     // float dX = source->moments->Mx - source->peak->xf;
    134     // float dY = source->moments->My - source->peak->yf;
    135     // float dR = hypot(dX, dY);
    136     // float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
    137     // float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
    138342    float Xo = source->moments->Mx;
    139343    float Yo = source->moments->My;
    … …  
    154358
    155359        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    156         // psImageMaskType  *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
    157360
    158361        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    … …  
    173376            if (r2 > R2) continue;
    174377
    175             float fDiff = *vPix - sky;
     378            float fDiff = *vPix;
    176379            float pDiff = fDiff;
    177380            float wDiff = *vWgt;
    … …  
    181384            // stars.
    182385            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    183             if ((minSN > 0.0) && (pDiff < 0)) continue; //
     386            if ((minSN > 0.0) && (pDiff < 0)) continue;
    184387
    185388            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    186             // weighting over weights the sky for faint sources
     389            // weighting over-weights the sky for faint sources
    187390            if (sigma > 0.0) {
    188391                float z = r2 * rsigma2;
    … …  
    230433            XYYY  += xyyy;
    231434            YYYY  += yyyy;
    232 
    233             // Kron Flux uses the 1st radial moment (NOT Gaussian windowed?)
    234             // XXX float r = sqrt(r2);
    235             // XXX float rf = r * fDiff;
    236             // XXX float rh = sqrt(r) * fDiff;
    237             // XXX float rs = fDiff;
    238             // XXX
    239             // XXX float rfw = r * pDiff;
    240             // XXX float rhw = sqrt(r) * pDiff;
    241             // XXX
    242             // XXX RF  += rf;
    243             // XXX RH  += rh;
    244             // XXX RS  += rs;
    245             // XXX
    246             // XXX RFW  += rfw;
    247             // XXX RHW  += rhw;
    248435        }
    249436    }
    … …  
    263450    source->moments->Myyyy = YYYY/Sum;
    264451
    265     // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging
     452    return true;
     453}
     454
     455bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
     456
    266457
    267458    float **vPix = source->pixels->data.F32;
    268     float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
    269459    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    270460
    … …  
    272462    float RH = 0.0;
    273463    float RS = 0.0;
     464
     465    // centroid around which to calculate the moments
     466    float Xo = source->moments->Mx;
     467    float Yo = source->moments->My;
     468
     469    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     470    // xCM, yCM from pixel coords to pixel index here.
     471    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     472    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     473
     474    float R2 = PS_SQR(radius);
    274475
    275476    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    304505            if (r2 > R2) continue;
    305506
    306             float fDiff1 = vPix[row][col] - sky;
    307             float fDiff2 = vPix[yFlip][xFlip] - sky;
     507            float fDiff1 = vPix[row][col];
     508            float fDiff2 = vPix[yFlip][xFlip];
    308509            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    309510
    … …  
    329530        kronRefRadius = MIN(radius, kronRefRadius);
    330531    }
    331     source->moments->Mrf = kronRefRadius;
    332 
    333     // *** now calculate the kron flux values using the 1st radial moment
    334 
    335     float radKinner = 1.0*kronRefRadius;
    336     float radKron   = 2.5*kronRefRadius;
    337     float radKouter = 4.0*kronRefRadius;
     532
     533    // if source is externally supplied and it already has a finite Mrf do not change it
     534    if (! ((source->mode & PM_SOURCE_MODE_EXTERNAL) && isfinite(source->moments->Mrf))) {
     535        source->moments->Mrf = kronRefRadius;
     536    }
     537
     538    return true;
     539}
     540
     541bool pmSourceMomentsKronFluxes (pmSource *source, float sigma, float minSN, psImageMaskType maskVal) {
     542
     543    float **vPix = source->pixels->data.F32;
     544    float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
     545    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     546
     547    float radKinner = 1.0*source->moments->Mrf;
     548    float radKron   = 2.5*source->moments->Mrf;
     549    float radKouter = 4.0*source->moments->Mrf;
    338550
    339551    int nKronPix = 0;
    … …  
    341553    int nInner = 0;
    342554    int nOuter = 0;
    343     Sum = Var = 0.0;
     555   
     556    float Sum = 0.0;
     557    float Var = 0.0;
     558    float SumCore = 0.0;
     559    float VarCore = 0.0;
    344560    float SumInner = 0.0;
    345561    float SumOuter = 0.0;
     562
     563    // centroid around which to calculate the moments
     564    float Xo = source->moments->Mx;
     565    float Yo = source->moments->My;
     566
     567    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     568    // xCM, yCM from pixel coords to pixel index here.
     569    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     570    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     571
     572    float minSN2 = PS_SQR(minSN);
    346573
    347574    // calculate the Kron flux, and related fluxes (NO symmetrical averaging)
    … …  
    362589            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    363590
    364             float fDiff1 = vPix[row][col] - sky;
     591            float fDiff1 = vPix[row][col];
    365592            float pDiff = fDiff1;
    366593            float wDiff = vWgt[row][col];
    … …  
    376603                Var += wDiff;
    377604                nKronPix ++;
    378                 // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
    379605            }
    380606
    … …  
    397623    }
    398624    // *** should I rescale these fluxes by pi R^2 / nNpix?
    399     // XXX source->moments->KronCore    = SumCore       * M_PI * PS_SQR(sigma) / nCorePix;
    400     // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI * PS_SQR(sigma) / nCorePix;
    401     // XXX source->moments->KronFlux    = Sum       * M_PI * PS_SQR(radKron) / nKronPix;
    402     // XXX source->moments->KronFluxErr = sqrt(Var) * M_PI * PS_SQR(radKron) / nKronPix;
    403     // XXX source->moments->KronFinner = SumInner * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
    404     // XXX source->moments->KronFouter = SumOuter * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
     625    // XXX source->moments->KronCore    = SumCore       * M_PI *  PS_SQR(sigma) / nCorePix;
     626    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI *  PS_SQR(sigma) / nCorePix;
     627    // XXX source->moments->KronFlux    = Sum           * M_PI * PS_SQR(radKron) / nKronPix;
     628    // XXX source->moments->KronFluxErr = sqrt(Var)     * M_PI * PS_SQR(radKron) / nKronPix;
     629    // XXX source->moments->KronFinner  = SumInner      * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
     630    // XXX source->moments->KronFouter  = SumOuter      * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
    405631
    406632    source->moments->KronCore    = SumCore;
    … …  
    408634    source->moments->KronFlux    = Sum;
    409635    source->moments->KronFluxErr = sqrt(Var);
    410     source->moments->KronFinner = SumInner;
    411     source->moments->KronFouter = SumOuter;
     636    source->moments->KronFinner  = SumInner;
     637    source->moments->KronFouter  = SumOuter;
    412638
    413639    // XXX not sure I should save this here...
    … …  
    416642    source->moments->KronRadiusPSF  = source->moments->Mrf;
    417643
    418     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",
    419              source->moments->Mrf,   source->moments->KronFlux,
    420              source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
    421              source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
    422              source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
    423 
    424     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",
    425              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);
    426 
    427     return(true);
    428 }
    429 
    430 bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
    431 
    432     // First Pass: calculate the first moments (these are subtracted from the coordinates below)
    433     // Sum = SUM (z - sky)
    434     // X1  = SUM (x - xc)*(z - sky)
    435     // .. etc
    436 
    437     float sky = 0.0;
    438 
    439     float peakPixel = -PS_MAX_F32;
    440     psS32 numPixels = 0;
    441     float Sum = 0.0;
    442     float Var = 0.0;
    443     float X1 = 0.0;
    444     float Y1 = 0.0;
    445     float R2 = PS_SQR(radius);
    446     float minSN2 = PS_SQR(minSN);
    447     float rsigma2 = 0.5 / PS_SQR(sigma);
    448 
    449     float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
    450     float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
    451 
    452     // we are guaranteed to have a valid pixel and variance at this location (right? right?)
    453     // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
    454     // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    455     // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
    456     // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
    457 
    458     // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
    459     // not depend on the fractional pixel location of the source.  However, the aperture
    460     // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
    461     // position of the expected centroid
    462 
    463     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    464 
    465         float yDiff = row + 0.5 - yPeak;
    466         if (fabs(yDiff) > radius) continue;
    467 
    468         float *vPix = source->pixels->data.F32[row];
    469         float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
    470 
    471         psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
    472         // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
    473 
    474         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    475             if (vMsk) {
    476                 if (*vMsk & maskVal) {
    477                     vMsk++;
    478                     continue;
    479                 }
    480                 vMsk++;
    481             }
    482             if (isnan(*vPix)) continue;
    483 
    484             float xDiff = col + 0.5 - xPeak;
    485             if (fabs(xDiff) > radius) continue;
    486 
    487             // radius is just a function of (xDiff, yDiff)
    488             float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    489             if (r2 > R2) continue;
    490 
    491             float pDiff = *vPix - sky;
    492             float wDiff = *vWgt;
    493 
    494             // skip pixels below specified significance level.  for a PSFs, this
    495             // over-weights the wings of bright stars compared to those of faint stars.
    496             // for the estimator used for extended source analysis (where the window
    497             // function is allowed to be arbitrarily large), we need to clip to avoid
    498             // negative second moments.
    499             if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
    500             if ((minSN > 0.0) && (pDiff < 0)) continue; //
    501 
    502             // Apply a Gaussian window function.  Be careful with the window function.  S/N
    503             // weighting over weights the sky for faint sources
    504             if (sigma > 0.0) {
    505                 float z  = r2*rsigma2;
    506                 assert (z >= 0.0);
    507                 float weight  = exp(-z);
    508 
    509                 wDiff *= weight;
    510                 pDiff *= weight;
    511             }
    512 
    513             Var += wDiff;
    514             Sum += pDiff;
    515 
    516             float xWght = xDiff * pDiff;
    517             float yWght = yDiff * pDiff;
    518 
    519             X1  += xWght;
    520             Y1  += yWght;
    521 
    522             peakPixel = PS_MAX (*vPix, peakPixel);
    523             numPixels++;
    524         }
    525     }
    526 
    527     // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
    528     int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
    529 
    530     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    531     if ((numPixels < minPixels) || (Sum <= 0)) {
    532         psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
    533         return (false);
    534     }
    535 
    536     // calculate the first moment.
    537     float Mx = X1/Sum;
    538     float My = Y1/Sum;
    539     if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
    540         psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
    541         return (false);
    542     }
    543     if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) {
    544         psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y);
    545     }
    546 
    547     psTrace ("psModules.objects", 5, "id: %d, sky: %f  Mx: %f  My: %f  Sum: %f  X1: %f  Y1: %f  Npix: %d\n", source->id, sky, Mx, My, Sum, X1, Y1, numPixels);
    548 
    549     // add back offset of peak in primary image
    550     // also offset from pixel index to pixel coordinate
    551     // (the calculation above uses pixel index instead of coordinate)
    552     // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
    553 
    554     // we only update the centroid if the position is not supplied from elsewhere
    555     bool skipCentroid = false;
    556     skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
    557     skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
    558 
    559     if (skipCentroid) {
    560         source->moments->Mx = source->peak->xf;
    561         source->moments->My = source->peak->yf;
    562     } else {
    563         source->moments->Mx = Mx + xGuess;
    564         source->moments->My = My + yGuess;
    565     }
    566 
    567     source->moments->Sum = Sum;
    568     source->moments->SN  = Sum / sqrt(Var);
    569     source->moments->Peak = peakPixel;
    570     source->moments->nPixels = numPixels;
    571 
    572644    return true;
    573645}
    574 
    575 float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
    576 
    577     psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
    578 
    579     for (int i = 0; i < sources->n; i++) {
    580         pmSource *src = sources->data[i]; // Source of interest
    581         if (!src || !src->moments) {
    582             continue;
    583         }
    584 
    585         if (src->mode & PM_SOURCE_MODE_BLEND) {
    586             continue;
    587         }
    588 
    589         if (!src->moments->nPixels) continue;
    590 
    591         if (src->moments->SN < PSF_SN_LIM) continue;
    592 
    593         // XXX put in Mxx,Myy cut based on clump location
    594 
    595         psVectorAppend(radii, src->moments->Mrf);
    596     }
    597 
    598     // find the peak in this image
    599     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    600 
    601     if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
    602         psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
    603         psFree(stats);
    604         return NAN;
    605     }
    606 
    607     float minRadius = stats->sampleMedian;
    608 
    609     psFree(radii);
    610     psFree(stats);
    611     return minRadius;
    612 }
    613 
Note: See TracChangeset for help on using the changeset viewer.