IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 17, 2014, 12:30:45 PM (12 years ago)
Author:
eugene
Message:

merge changes (from past YEAR) into trunk

Location:
branches/eam_branches/ipp-ops-20130712/psModules
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-ops-20130712/psModules

    • Property svn:mergeinfo deleted

  • branches/eam_branches/ipp-ops-20130712/psModules/src/objects

    • Property svn:ignore

      •  

        old new  
        1212pmSourceIO_CMF_PS1_V1.v1.c
        1313pmSourceIO_CMF_PS1_V4.c
         14pmSourceIO_CMF_PS1_V5.c
        1415pmSourceIO_CMF_PS1_SV1.c
        1516pmSourceIO_CMF_PS1_SV2.c
         17pmSourceIO_CMF_PS1_SV3.c
        1618pmSourceIO_CMF_PS1_DV1.c
        1719pmSourceIO_CMF_PS1_DV2.c
        1820pmSourceIO_CMF_PS1_DV3.c
        19 
         21pmSourceIO_CMF_PS1_DV4.c

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

    r35560 r37066  
    3535#include "pmMoments.h"
    3636#include "pmModelFuncs.h"
     37#include "pmModelClass.h"
    3738#include "pmModel.h"
    3839#include "pmModelUtils.h"
    39 #include "pmModelClass.h"
    4040#include "pmSourceMasks.h"
    4141#include "pmSourceExtendedPars.h"
    4242#include "pmSourceDiffStats.h"
    4343#include "pmSourceSatstar.h"
     44#include "pmSourceLensing.h"
    4445#include "pmSource.h"
    4546
    … …  
    6566void pmSourceMomentsSetVerbose(bool state){ beVerbose = state; }
    6667
     68bool pmSourceMomentsHighOrder    (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal);
     69bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal);
     70bool pmSourceMomentsKronFluxes   (pmSource *source, float sigma,  float minSN, psImageMaskType maskVal);
     71
    6772// if mode & EXTERNAL or mode2 & MATCHED, do not re-calculate the centroid (use peak as centroid)
    68 
    6973bool pmSourceMoments(pmSource *source, float radius, float sigma, float minSN, float minKronRadius, psImageMaskType maskVal)
    7074{
    … …  
    7478    PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
    7579
    76     // this function assumes the sky has been well-subtracted for the image
    77     float sky = 0.0;
    78 
    7980    if (source->moments == NULL) {
    8081      source->moments = pmMomentsAlloc();
    8182    }
    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);
    9083
    9184    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    110103    // of any object drops pretty quickly outside 1-2 sigmas.  (The exception is bright
    111104    // saturated stars, for which we need to use a very large radius here)
     105    // NOTE: if (mode & EXTERNAL) or (mode2 & MATCHED), do not re-calculate the centroid (use peak as centroid)
     106    // (we still call this function because it sets moments->Sum,SN,Peak,nPixels
    112107    if (!pmSourceMomentsGetCentroid (source, 1.5*sigma, 0.0, minSN, maskVal, source->peak->xf, source->peak->yf)) {
    113108        return false;
    114109    }
    115110
     111    pmSourceMomentsHighOrder (source, radius, sigma, minSN, maskVal);
     112
     113    // now calculate the 1st radial moment (for kron flux) using symmetrical averaging
     114    pmSourceMomentsRadialMoment (source, radius, minKronRadius, maskVal);
     115
     116    // now calculate the kron flux values using the 1st radial moment
     117    pmSourceMomentsKronFluxes (source, sigma, minSN, maskVal);
     118
     119    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",
     120             source->moments->Mrf,   source->moments->KronFlux,
     121             source->moments->Mxx,   source->moments->Mxy,   source->moments->Myy,
     122             source->moments->Mxxx,  source->moments->Mxxy,  source->moments->Mxyy,  source->moments->Myyy,
     123             source->moments->Mxxxx, source->moments->Mxxxy, source->moments->Mxxyy, source->moments->Mxyyy, source->moments->Myyyy);
     124
     125    psTrace ("psModules.objects", 3, "peak %f %f (%f = %f) Mx: %f  My: %f  Sum: %f  Mxx: %f  Mxy: %f  Myy: %f  Npix: %d\n",
     126             source->peak->xf, source->peak->yf,
     127             source->peak->rawFlux, sqrt(source->peak->detValue),
     128             source->moments->Mx, source->moments->My,
     129             source->moments->Sum,
     130             source->moments->Mxx, source->moments->Mxy, source->moments->Myy,
     131             source->moments->nPixels);
     132
     133    return(true);
     134}
     135
     136bool pmSourceMomentsGetCentroid(pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal, float xGuess, float yGuess) {
     137
     138    // First Pass: calculate the first moments (these are subtracted from the coordinates below)
     139    // Sum = SUM (z - sky)
     140    // X1  = SUM (x - xc)*(z - sky)
     141    // .. etc
     142
     143    float sky = 0.0;
     144
     145    float peakPixel = -PS_MAX_F32;
     146    psS32 numPixels = 0;
     147    float Sum = 0.0;
     148    float Var = 0.0;
     149    float X1 = 0.0;
     150    float Y1 = 0.0;
     151    float R2 = PS_SQR(radius);
     152    float minSN2 = PS_SQR(minSN);
     153    float rsigma2 = 0.5 / PS_SQR(sigma);
     154
     155    float xPeak = xGuess - source->pixels->col0; // coord of peak in subimage
     156    float yPeak = yGuess - source->pixels->row0; // coord of peak in subimage
     157
     158    // we are guaranteed to have a valid pixel and variance at this location (right? right?)
     159    // float weightNorm = source->pixels->data.F32[yPeak][xPeak] / sqrt (source->variance->data.F32[yPeak][xPeak]);
     160    // psAssert (isfinite(source->pixels->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     161    // psAssert (isfinite(source->variance->data.F32[yPeak][xPeak]), "peak must be on valid pixel");
     162    // psAssert (source->variance->data.F32[yPeak][xPeak] > 0, "peak must be on valid pixel");
     163
     164    // the moments [Sum(x*f) / Sum(f)] are calculated in pixel index values, and should
     165    // not depend on the fractional pixel location of the source.  However, the aperture
     166    // (radius) and the Gaussian window (sigma) depend subtly on the fractional pixel
     167    // position of the expected centroid
     168
     169    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     170
     171        float yDiff = row + 0.5 - yPeak;
     172        if (fabs(yDiff) > radius) continue;
     173
     174        float *vPix = source->pixels->data.F32[row];
     175        float *vWgt = source->variance ? source->variance->data.F32[row] : source->pixels->data.F32[row];
     176
     177        psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[row];
     178        // psImageMaskType *vMsk = (source->maskView == NULL) ? NULL : source->maskView->data.PS_TYPE_IMAGE_MASK_DATA[row];
     179
     180        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
     181            if (vMsk) {
     182                if (*vMsk & maskVal) {
     183                    vMsk++;
     184                    continue;
     185                }
     186                vMsk++;
     187            }
     188            if (isnan(*vPix)) continue;
     189
     190            float xDiff = col + 0.5 - xPeak;
     191            if (fabs(xDiff) > radius) continue;
     192
     193            // radius is just a function of (xDiff, yDiff)
     194            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     195            if (r2 > R2) continue;
     196
     197            float pDiff = *vPix - sky;
     198            float wDiff = *vWgt;
     199
     200            // skip pixels below specified significance level.  for a PSFs, this
     201            // over-weights the wings of bright stars compared to those of faint stars.
     202            // for the estimator used for extended source analysis (where the window
     203            // function is allowed to be arbitrarily large), we need to clip to avoid
     204            // negative second moments.
     205            if (PS_SQR(pDiff) < minSN2*wDiff) continue; //
     206            if ((minSN > 0.0) && (pDiff < 0)) continue; //
     207
     208            // Apply a Gaussian window function.  Be careful with the window function.  S/N
     209            // weighting over weights the sky for faint sources
     210            if (sigma > 0.0) {
     211                float z  = r2*rsigma2;
     212                assert (z >= 0.0);
     213                float weight  = exp(-z);
     214
     215                wDiff *= weight;
     216                pDiff *= weight;
     217            }
     218
     219            Var += wDiff;
     220            Sum += pDiff;
     221
     222            float xWght = xDiff * pDiff;
     223            float yWght = yDiff * pDiff;
     224
     225            X1  += xWght;
     226            Y1  += yWght;
     227
     228            peakPixel = PS_MAX (*vPix, peakPixel);
     229            numPixels++;
     230        }
     231    }
     232
     233    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
     234    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
     235
     236    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     237    if ((numPixels < minPixels) || (Sum <= 0)) {
     238        psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
     239        return (false);
     240    }
     241
     242    // calculate the first moment.
     243    float Mx = X1/Sum;
     244    float My = Y1/Sum;
     245    if ((fabs(Mx) > radius) || (fabs(My) > radius)) {
     246        psTrace ("psModules.objects", 3, "extreme centroid swing; invalid peak %d, %d\n", source->peak->x, source->peak->y);
     247        return (false);
     248    }
     249    if ((fabs(Mx) > 2.0) || (fabs(My) > 2.0)) {
     250        psTrace ("psModules.objects", 3, " big centroid swing; ok peak? %d, %d\n", source->peak->x, source->peak->y);
     251    }
     252
     253    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);
     254
     255    // add back offset of peak in primary image
     256    // also offset from pixel index to pixel coordinate
     257    // (the calculation above uses pixel index instead of coordinate)
     258    // 0.5 PIX: moments are calculated using the pixel index and converted here to pixel coords
     259
     260    // we only update the centroid if the position is not supplied from elsewhere
     261    bool skipCentroid = false;
     262    skipCentroid |= (source->mode  & PM_SOURCE_MODE_EXTERNAL); // skip externally supplied positions
     263    skipCentroid |= (source->mode2 & PM_SOURCE_MODE2_MATCHED); // skip sources defined by other image positions
     264
     265    if (skipCentroid) {
     266        source->moments->Mx = source->peak->xf;
     267        source->moments->My = source->peak->yf;
     268    } else {
     269        source->moments->Mx = Mx + xGuess;
     270        source->moments->My = My + yGuess;
     271    }
     272
     273    source->moments->Sum = Sum;
     274    source->moments->SN  = Sum / sqrt(Var);
     275    source->moments->Peak = peakPixel;
     276    source->moments->nPixels = numPixels;
     277
     278    return true;
     279}
     280
     281float pmSourceMinKronRadius(psArray *sources, float PSF_SN_LIM) {
     282
     283    psVector *radii = psVectorAllocEmpty(100, PS_TYPE_F32);
     284
     285    for (int i = 0; i < sources->n; i++) {
     286        pmSource *src = sources->data[i]; // Source of interest
     287        if (!src || !src->moments) {
     288            continue;
     289        }
     290
     291        if (src->mode & PM_SOURCE_MODE_BLEND) {
     292            continue;
     293        }
     294
     295        if (!src->moments->nPixels) continue;
     296
     297        if (src->moments->SN < PSF_SN_LIM) continue;
     298
     299        // XXX put in Mxx,Myy cut based on clump location
     300
     301        psVectorAppend(radii, src->moments->Mrf);
     302    }
     303
     304    // find the peak in this image
     305    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     306
     307    if (!psVectorStats (stats, radii, NULL, NULL, 0)) {
     308        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
     309        psFree(stats);
     310        return NAN;
     311    }
     312
     313    float minRadius = stats->sampleMedian;
     314
     315    psFree(radii);
     316    psFree(stats);
     317    return minRadius;
     318}
     319
     320bool pmSourceMomentsHighOrder (pmSource *source, float radius, float sigma, float minSN, psImageMaskType maskVal) {
     321
     322    // this function assumes the sky has been well-subtracted for the image
     323    float Sum = 0.0;
     324    float R2 = PS_SQR(radius);
     325    float minSN2 = PS_SQR(minSN);
     326    float rsigma2 = 0.5 / PS_SQR(sigma);
     327
    116328    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
    117     // Xn  = SUM (x - xc)^n * (z - sky)
     329    // Xn  = SUM (x - xc)^n * (z - sky) -- note that sky is 0.0 by definition here
    118330    float XX = 0.0;
    119331    float XY = 0.0;
    … …  
    129341    float YYYY = 0.0;
    130342
    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;
    138343    float Xo = source->moments->Mx;
    139344    float Yo = source->moments->My;
    … …  
    154359
    155360        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];
    157361
    158362        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    … …  
    173377            if (r2 > R2) continue;
    174378
    175             float fDiff = *vPix - sky;
     379            float fDiff = *vPix;
    176380            float pDiff = fDiff;
    177381            float wDiff = *vWgt;
    … …  
    181385            // stars.
    182386            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    183             if ((minSN > 0.0) && (pDiff < 0)) continue; //
     387            if ((minSN > 0.0) && (pDiff < 0)) continue;
    184388
    185389            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    186             // weighting over weights the sky for faint sources
     390            // weighting over-weights the sky for faint sources
    187391            if (sigma > 0.0) {
    188392                float z = r2 * rsigma2;
    … …  
    230434            XYYY  += xyyy;
    231435            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;
    248436        }
    249437    }
     438    // NOT needed : source->moments->wSum = Sum;
     439
    250440    source->moments->Mxx = XX/Sum;
    251441    source->moments->Mxy = XY/Sum;
    … …  
    263453    source->moments->Myyyy = YYYY/Sum;
    264454
    265     // *** now calculate the 1st radial moment (for kron flux) -- symmetrical averaging
     455    return true;
     456}
     457
     458bool pmSourceMomentsRadialMoment (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal) {
     459
    266460
    267461    float **vPix = source->pixels->data.F32;
    268     float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
    269462    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    270463
    … …  
    272465    float RH = 0.0;
    273466    float RS = 0.0;
     467
     468    // centroid around which to calculate the moments
     469    float Xo = source->moments->Mx;
     470    float Yo = source->moments->My;
     471
     472    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     473    // xCM, yCM from pixel coords to pixel index here.
     474    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     475    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     476
     477    float R2 = PS_SQR(radius);
    274478
    275479    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    304508            if (r2 > R2) continue;
    305509
    306             float fDiff1 = vPix[row][col] - sky;
    307             float fDiff2 = vPix[yFlip][xFlip] - sky;
     510            float fDiff1 = vPix[row][col];
     511            float fDiff2 = vPix[yFlip][xFlip];
    308512            float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    309513
    … …  
    329533        kronRefRadius = MIN(radius, kronRefRadius);
    330534    }
    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;
     535
     536    // if source is externally supplied and it already has a finite Mrf do not change it
     537    if (! ((source->mode & PM_SOURCE_MODE_EXTERNAL) && isfinite(source->moments->Mrf))) {
     538        source->moments->Mrf = kronRefRadius;
     539    }
     540
     541    return true;
     542}
     543
     544bool pmSourceMomentsKronFluxes (pmSource *source, float sigma, float minSN, psImageMaskType maskVal) {
     545
     546    float **vPix = source->pixels->data.F32;
     547    float **vWgt = source->variance ? source->variance->data.F32 : source->pixels->data.F32;
     548    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     549
     550    float radKinner = 1.0*source->moments->Mrf;
     551    float radKron   = 2.5*source->moments->Mrf;
     552    float radKouter = 4.0*source->moments->Mrf;
    338553
    339554    int nKronPix = 0;
    … …  
    341556    int nInner = 0;
    342557    int nOuter = 0;
    343     Sum = Var = 0.0;
     558   
     559    float Sum = 0.0;
     560    float Var = 0.0;
     561    float SumCore = 0.0;
     562    float VarCore = 0.0;
    344563    float SumInner = 0.0;
    345564    float SumOuter = 0.0;
     565
     566    // centroid around which to calculate the moments
     567    float Xo = source->moments->Mx;
     568    float Yo = source->moments->My;
     569
     570    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     571    // xCM, yCM from pixel coords to pixel index here.
     572    float xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     573    float yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     574
     575    float minSN2 = PS_SQR(minSN);
    346576
    347577    // calculate the Kron flux, and related fluxes (NO symmetrical averaging)
    … …  
    362592            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    363593
    364             float fDiff1 = vPix[row][col] - sky;
     594            float fDiff1 = vPix[row][col];
    365595            float pDiff = fDiff1;
    366596            float wDiff = vWgt[row][col];
    … …  
    376606                Var += wDiff;
    377607                nKronPix ++;
    378                 // if (beVerbose) fprintf (stderr, "mome: %d %d  %f  %f  %f\n", col, row, sky, *vPix, Sum);
    379608            }
    380609
    … …  
    397626    }
    398627    // *** 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;
     628    // XXX source->moments->KronCore    = SumCore       * M_PI *  PS_SQR(sigma) / nCorePix;
     629    // XXX source->moments->KronCoreErr = sqrt(VarCore) * M_PI *  PS_SQR(sigma) / nCorePix;
     630    // XXX source->moments->KronFlux    = Sum           * M_PI * PS_SQR(radKron) / nKronPix;
     631    // XXX source->moments->KronFluxErr = sqrt(Var)     * M_PI * PS_SQR(radKron) / nKronPix;
     632    // XXX source->moments->KronFinner  = SumInner      * M_PI * (PS_SQR(radKron)   - PS_SQR(radKinner)) / nInner;
     633    // XXX source->moments->KronFouter  = SumOuter      * M_PI * (PS_SQR(radKouter) -   PS_SQR(radKron)) / nOuter;
    405634
    406635    source->moments->KronCore    = SumCore;
    … …  
    408637    source->moments->KronFlux    = Sum;
    409638    source->moments->KronFluxErr = sqrt(Var);
    410     source->moments->KronFinner = SumInner;
    411     source->moments->KronFouter = SumOuter;
     639    source->moments->KronFinner  = SumInner;
     640    source->moments->KronFouter  = SumOuter;
    412641
    413642    // XXX not sure I should save this here...
    … …  
    416645    source->moments->KronRadiusPSF  = source->moments->Mrf;
    417646
    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 
    572647    return true;
    573648}
    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.