IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 6, 2012, 6:32:19 PM (14 years ago)
Author:
watersc1
Message:

merge from trunk. Preliminary versions of stacktool/warptool updates and regenerate_background.pl script. The warp code is finished and tested, and I still need to get the stack version resolved.

Location:
branches/czw_branch/20120906/psphot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20120906/psphot

  • branches/czw_branch/20120906/psphot/src

  • branches/czw_branch/20120906/psphot/src/psphotKronIterate.c

    r34404 r34772  
    11# include "psphotInternal.h"
    22
    3 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
    4 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);
    5 
    6 
    7 bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass)
    8 {
     3bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal, psImage *smoothedPixels);
     4bool psphotKronFluxMeasure (pmSource *source, psImageMaskType maskVal);
     5
     6
     7bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule, int pass) {
    98    bool status = true;
    109
    … …  
    9897    }
    9998
    100     bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
    101     if (!status) {
    102         KRON_APPLY_WEIGHT = true;
    103     }
    104 
    105     bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
    106     if (!status) {
    107         KRON_APPLY_WINDOW = false;
    108     }
    10999    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
    110100    if (!status) {
    … …  
    170160    }
    171161
    172     // generate the window image: multiply the flux by this to downweight neighbors
    173     // XXX: we don't need this image if we aren't going to apply the window
    174     psImage *kronWindow = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
    175     psImageInit (kronWindow, 1.0);
    176 
    177     // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
    178     for (int i = 0; i < sources->n; i++) {
    179 
    180         pmSource *source = sources->data[i];
    181 
    182         // set a window function for each source based on the moments
    183         // (this skips really bad sources (no peak, no moments, DEFECT)
    184         psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
     162    // pass 3 we only measure fluxes for matched sources so we don't need smooth pixels or multiple iterations
     163    if (pass == 3) {
     164        KRON_SMOOTH = false;
     165        KRON_ITERATIONS = 1;
    185166    }
    186167
    … …  
    235216
    236217            psArrayAdd(job->args, 1, readout);
    237             psArrayAdd(job->args, 1, kronWindow);
    238218            psArrayAdd(job->args, 1, cells->data[j]); // sources
    239219            psArrayAdd(job->args, 1, smoothedImage);
    … …  
    243223            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
    244224            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
    245             PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
    246             PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
    247225            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
    248226            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
    … …  
    281259    }
    282260    psFree (cellGroups);
    283     psFree (kronWindow);
    284261    if (KRON_SMOOTH) {
    285262        for (int i = 0; i < sources->n; i++) {
    … …  
    298275
    299276    pmReadout *readout              = job->args->data[0];
    300     psImage *kronWindow             = job->args->data[1];
    301     psArray *sources                = job->args->data[2];
    302     psImage *smoothedImage          = job->args->data[3];
    303     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    304     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    305     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
    306     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
    307     int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
    308     bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
    309     bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
    310     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
    311     float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
    312     int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
    313 #ifndef REVERT_ON_BAD_MEASUREMENT
    314     (void) pass;
    315 #endif
     277    psArray *sources                = job->args->data[1];
     278    psImage *smoothedImage          = job->args->data[2];
     279    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
     280    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     281    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
     282    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
     283    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
     284    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[8],F32);
     285    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[9],F32);
     286    int pass                        = PS_SCALAR_VALUE(job->args->data[10],S32);
     287
     288    bool measureRadius = true;
     289    if (pass == 3) {
     290        measureRadius = false;
     291    }
    316292
    317293    for (int j = 0; j < KRON_ITERATIONS; j++) {
    … …  
    320296            pmSource *source = sources->data[i];
    321297            if (!source->peak) continue; // XXX how can we have a peak-less source?
    322 
    323             // check status of this source's moments
    324             if (!source->moments) continue;
    325             if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
    326             if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     298            if (!source->moments) continue;
     299            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     300            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     301            // skip saturated stars modeled with a radial profile
     302            if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
     303
     304            // in pass 3 we only measure the flux for matched sources
     305            if ((pass == 3) && !(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue;
     306
     307# if (0)
     308# define TEST_X 653
     309# define TEST_Y 466
     310            if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     311                fprintf (stderr, "test object\n");
     312            }
     313# undef TEST_X
     314# undef TEST_Y
     315# endif
     316
     317            if (measureRadius) {
     318                // check status of this source's moments
     319                // XXX: I don't think that we have to apply these restrictions since we dropped the window function
     320                if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
     321                if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
     322            }
     323
    327324            if (!isfinite(source->moments->Mrf)) {
    328325                // Once we save a bad Mrf measurement we give up on this source
    329                 // checking here allows us to avoid adding and subtracting the model
     326                // XXX: is this the right thing to do?
    330327                continue;
    331328            }
    332 
    333             // skip saturated stars modeled with a radial profile
    334             if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
    335329
    336330            // replace object in image
    … …  
    348342            // On first iteration set window radius to sky radius (if valid). We also use this on subsequent
    349343            // iterations if we cannot find a better limit
    350             float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     344            float maxWindow;
     345            if (measureRadius) {
     346                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     347            } else {
     348                maxWindow = source->moments->Mrf;
     349            }
    351350            if (j > 0) {
    352351                // on subsequent iterations we use a factor times the previous radial moment value
    … …  
    370369            float windowRadius = PS_MAX(RADIUS, maxWindow);
    371370
    372 #ifdef REVERT_ON_BAD_MEASURMENT
    373             // save previous measurements. We might revert back to them if this round fails
    374             float MrfPrior = source->moments->Mrf;
    375             float KronFluxPrior = source->moments->KronFlux;
    376             float KronFluxErrPrior = source->moments->KronFluxErr;
    377 #endif
    378 
    379371            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    380372            bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    381             psAssert (source->pixels, "WTF?");
     373            psAssert (source->pixels, "redefine pixels failed?");
    382374            if (extend && smoothedPixels) {
    383375                psFree(source->tmpPtr);
    384376                smoothedPixels = psImageSubset(smoothedImage, source->region);
    385                 psAssert (smoothedPixels, "WTF?");
     377                psAssert (smoothedPixels, "redefine smoothed pixels failed?");
    386378                source->tmpPtr = (psPtr) smoothedPixels ;
    387379            }
    388380
    389 
    390             // clear the window function for this source based on the moments
    391             // Note: this function also applies cuts on the source and returns false if it
    392             // does not meet the requirements for measuring the Kron Radius or Magnitude.
    393             // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
    394             if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
    395 
    396                 // this function populates moments->Mrf,KronFlux,KronFluxErr
    397                 psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
    398                 psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    399 
    400                 // set a window function for each source based on the moments
    401                 psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
    402             }
    403 
    404 #ifdef REVERT_ON_BAD_MEASUREMENT
    405             // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
    406             // in pass 1 keep that measurement
    407             bool reverted = false;
    408             if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
    409                 source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
    410                 source->moments->KronFlux = KronFluxPrior;
    411                 source->moments->KronFluxErr = KronFluxErrPrior;
    412                 reverted = true;
    413             }
    414 #endif
     381            bool measureFlux = true;
     382            if (measureRadius) {
     383                measureFlux = psphotKronRadiusMeasure (source, windowRadius, MIN_KRON_RADIUS, maskVal, smoothedPixels);
     384            }
     385
     386            if (measureFlux) {
     387                // Make sure the sources images are large enough for the measured Kron Radius
     388                bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y,
     389                    2.5 * source->moments->Mrf + 2);
     390
     391                psAssert (source->pixels, "redefine pixels failed?");
     392                if (extend && smoothedPixels) {
     393                    psFree(source->tmpPtr);
     394                    smoothedPixels = psImageSubset(smoothedImage, source->region);
     395                    psAssert (smoothedPixels, "redefine smoothed pixels failed?");
     396                    source->tmpPtr = (psPtr) smoothedPixels ;
     397                }
     398
     399                // this function populates moments->Mrf,KronFlux,KronFluxErr, KronFinner, and KronFouter
     400                psphotKronFluxMeasure (source, maskVal);
     401            }
     402
     403            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     404
    415405#ifdef DUMP_KRS
    416 #ifndef REVERT_ON_BAD_MEASUREMENT
    417             bool reverted = false;
    418 #endif
    419             fprintf(dumpFile, "%7d %1d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, reverted, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
     406            fprintf(dumpFile, "%7d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %2d\n", source->id, source->moments->Mrf, MrfPrior, maxWindow, windowRadius, source->peak->xf, source->peak->yf, source->imageID);
    420407#endif
    421408
    … …  
    433420}
    434421
    435 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
    436     bool applyWeight, psImage *smoothedPixels) {
     422bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal,
     423    psImage *smoothedPixels) {
    437424
    438425    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    442429
    443430    psF32 R2 = PS_SQR(radius);
    444     float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
    445431
    446432    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    455441    // Xn  = SUM (x - xc)^n * (z - sky)
    456442
    457 
    458443    psF32 RF = 0.0;
    459444    psF32 RS = 0.0;
    … …  
    471456    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    472457
    473     int Xwo = source->pixels->col0;
    474     int Ywo = source->pixels->row0;
    475 
    476458    psF32 **vPix;
    477459    if (smoothedPixels) {
    … …  
    480462        vPix = source->pixels->data.F32;
    481463    }
    482     psF32 **vWin = kronWindow->data.F32;
     464   
     465    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     466
     467    {
     468        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     469
     470            psF32 yDiff = row - yCM;
     471            if (fabs(yDiff) > radius) continue;
     472
     473            // coordinate of mirror pixel
     474            int yFlip = yCM - yDiff;
     475            if (yFlip < 0) continue;
     476            if (yFlip >= source->pixels->numRows) continue;
     477
     478            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     479                // check mask and value for this pixel
     480                if (vMsk && (vMsk[row][col] & maskVal)) continue;
     481                if (isnan(vPix[row][col])) continue;
     482
     483                psF32 xDiff = col - xCM;
     484                if (fabs(xDiff) > radius) continue;
     485
     486                // coordinate of mirror pixel
     487                int xFlip = xCM - xDiff;
     488                if (xFlip < 0) continue;
     489                if (xFlip >= source->pixels->numCols) continue;
     490
     491                // check mask and value for mirror pixel
     492                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     493                if (isnan(vPix[yFlip][xFlip])) continue;
     494
     495                // radius is just a function of (xDiff, yDiff)
     496                psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     497                if (r2 > R2) continue;
     498
     499                float fDiff1 = vPix[row][col];
     500                float fDiff2 = vPix[yFlip][xFlip];
     501
     502                float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
     503
     504                // Kron Flux uses the 1st radial moment
     505                psF32 rf = pDiff * sqrt(r2);
     506                psF32 rs = 0.5 * (fDiff1 + fDiff2);
     507
     508                RF  += rf;
     509                RS  += rs;
     510            }
     511        }
     512
     513        float MrfTry = RF/RS;
     514        if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
     515            // We did not get a good measurement
     516            source->moments->Mrf = NAN;
     517            source->moments->KronFlux  = NAN;
     518            source->moments->KronFluxErr  = NAN;
     519            return false;
     520        }
     521
     522        float Mrf = MAX(minKronRadius, MrfTry);
     523        // Saturate the 1st radial moment
     524        if (sqrt(source->peak->detValue) < 10.0) {
     525            Mrf = MIN (radius, Mrf);
     526        }
     527        source->moments->Mrf = Mrf;
     528    }
     529    return true;
     530}
     531
     532bool psphotKronFluxMeasure(pmSource *source, psImageMaskType maskVal) {
     533
     534    PS_ASSERT_PTR_NON_NULL(source, false);
     535    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     536    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     537
     538    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
     539    // image coordinates.  the source->pixels image has an offset relative to its parent of
     540    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
     541    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
     542    // minimal round-off error in the sums.  since these values are subtracted just to minimize
     543    // the dynamic range and are added back below, the exact value does not matter. these are
     544    // (int) so they can be used in the image index below.
     545
     546    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
     547    // Xn  = SUM (x - xc)^n * (z - sky)
     548
     549    // the peak position is less accurate but less subject to extreme deviations
     550    float dX = source->moments->Mx - source->peak->xf;
     551    float dY = source->moments->My - source->peak->yf;
     552    float dR = hypot(dX, dY);
     553    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
     554    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
     555
     556    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     557    // xCM, yCM from pixel coords to pixel index here.
     558    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     559    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     560
     561    psF32 **vPix = source->pixels->data.F32;
    483562    psF32 **vWgt = source->variance->data.F32;
    484563   
    485564    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    486565
     566    source->moments->KronFlux  = NAN;
     567    source->moments->KronFluxErr  = NAN;
     568    source->moments->KronFinner = NAN;
     569    source->moments->KronFouter = NAN;
     570
     571    // Calculate the Kron fluxes
     572    float radKinner = 1.0 * source->moments->Mrf;
     573    float radKron   = 2.5 * source->moments->Mrf;
     574    float radKouter = 4.0 * source->moments->Mrf;
     575
     576    float limitRadius = MIN (radKouter, source->windowRadius);
     577    if (radKouter > source->windowRadius) {
     578        // This happens but the measurement isn't important enough to allocate the extra pixels
     579        // psWarning ("outer kron radius: %f is larger than windowRadius: %f for %d\n",
     580        //     radKouter, source->windowRadius, source->id);
     581        limitRadius = MIN (radKron, source->windowRadius);
     582    }
     583    if (radKron > source->windowRadius) {
     584        // caller should have prevented this from happening
     585        psWarning ("kron radius: %f is larger than windowRadius: %f for %d\n",
     586            radKron, source->windowRadius, source->id);
     587        return false;
     588    }
     589
     590    float Sum = 0.0;
     591    float Var = 0.0;
     592    float SumInner = 0.0;
     593    float SumOuter = 0.0;
     594
     595    // set vPix to the source pixels (it may have been set to the smoothed image above)
     596    vPix = source->pixels->data.F32;
     597
    487598    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    488599
    489600        psF32 yDiff = row - yCM;
    490         if (fabs(yDiff) > radius) continue;
    491 
    492         // coordinate of mirror pixel
    493         int yFlip = yCM - yDiff;
    494         if (yFlip < 0) continue;
    495         if (yFlip >= source->pixels->numRows) continue;
     601        if (fabs(yDiff) > limitRadius) continue;
    496602
    497603        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    … …  
    501607
    502608            psF32 xDiff = col - xCM;
    503             if (fabs(xDiff) > radius) continue;
    504 
    505             // coordinate of mirror pixel
    506             int xFlip = xCM - xDiff;
    507             if (xFlip < 0) continue;
    508             if (xFlip >= source->pixels->numCols) continue;
    509 
    510             // check mask and value for mirror pixel
    511             if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
    512             if (isnan(vPix[yFlip][xFlip])) continue;
    513 
    514             // radius is just a function of (xDiff, yDiff)
     609            if (fabs(xDiff) > limitRadius) continue;
     610
    515611            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    516             if (r2 > R2) continue;
    517 
    518             // flux * window
    519             float z = r2 * rsigma2;
    520             assert (z >= 0.0);
    521 
    522             // weight by window image and wide Gaussian
    523             float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
    524             float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
    525 
    526             float fDiff1 = vPix[row][col]*weight1;
    527             float fDiff2 = vPix[yFlip][xFlip]*weight2;
    528 
    529             float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    530 
    531             // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
    532             psF32 rf = pDiff * sqrt(r2);
    533             psF32 rs = 0.5 * (fDiff1 + fDiff2);
    534 
    535             RF  += rf;
    536             RS  += rs;
     612            psF32 r = sqrt(r2);
     613
     614            float pDiff = vPix[row][col];
     615            psF32 wDiff = vWgt[row][col];
     616
     617            if (r > radKinner && r < radKron) {
     618                SumInner += pDiff;
     619            }
     620            if (r < radKron) {
     621                Sum += pDiff;
     622                Var += wDiff;
     623            }
     624            if (r > radKron && r < radKouter) {
     625                SumOuter += pDiff;
     626            }
    537627        }
    538628    }
    539629
    540     float MrfTry = RF/RS;
    541     if (RF <= 0. || RS <= 0 || !isfinite(MrfTry)) {
    542         // We did not get a good measurement
    543         source->moments->Mrf = NAN;
    544         source->moments->KronFlux  = NAN;
    545         source->moments->KronFluxErr  = NAN;
    546         return false;
    547     }
    548 
    549     float Mrf = MAX(minKronRadius, MrfTry);
    550     // Saturate the 1st radial moment
    551     if (sqrt(source->peak->detValue) < 10.0) {
    552         Mrf = MIN (radius, Mrf);
    553     }
    554 
    555     // Calculate the Kron magnitude (make this block optional?)
    556     float radKron  = 2.5*Mrf;
    557     float radKron2 = radKron*radKron;
    558 
    559     int nKronPix = 0;
    560     float Sum = 0.0;
    561     float Var = 0.0;
    562     float Win = 0.0;
    563 
    564     // set vPix to the source pixels (it may have been set to the
    565     // smoothed image above)
    566     vPix = source->pixels->data.F32;
    567 
    568 
    569     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    570 
    571         psF32 yDiff = row - yCM;
    572         if (fabs(yDiff) > radKron) continue;
    573 
    574         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    575             // check mask and value for this pixel
    576             if (vMsk && (vMsk[row][col] & maskVal)) continue;
    577             if (isnan(vPix[row][col])) continue;
    578 
    579             psF32 xDiff = col - xCM;
    580             if (fabs(xDiff) > radKron) continue;
    581 
    582             // radKron is just a function of (xDiff, yDiff)
    583             psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    584             if (r2 > radKron2) continue;
    585 
    586             float weight1  = vWin[row+Ywo][col+Xwo];
    587             float fDiff1 = vPix[row][col]*weight1;
    588 
    589             float pDiff = fDiff1;
    590             psF32 wDiff = vWgt[row][col] * weight1;
    591 
    592             Sum += pDiff;
    593             Var += wDiff;
    594             Win += weight1;
    595             nKronPix ++;
    596         }
    597     }
    598 
    599     source->moments->Mrf         = Mrf;
    600630    source->moments->KronFlux    = Sum;
    601631    source->moments->KronFluxErr = sqrt(Var);
     632    source->moments->KronFinner  = SumInner;
     633
     634    // only save radKouter if the radius is inside the integration radius limit
     635    if (radKouter <= limitRadius) {
     636        source->moments->KronFouter = SumOuter;
     637    }
    602638
    603639    return true;
    604640}
    605 
    606 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
    607 
    608     if (!source) return false;
    609     if (!source->peak) return false; // XXX how can we have a peak-less source?
    610     if (!source->moments) return false;
    611     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    612     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    613     if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
    614     if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    615     psAssert(kronWindow, "need a window");
    616 
    617     // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
    618     // We should give the this module a chance to measure a good value.
    619     // However experiments show that it hardly ever succeeds in getting a better value
    620     if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
    621 
    622     if (!applyWindow) {
    623         return true;
    624     }
    625 
    626     // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has
    627     // these values (and peak of 1.0)
    628 
    629     int Nx = kronWindow->numCols;
    630     int Ny = kronWindow->numRows;
    631 
    632     float Xo = source->moments->Mx;
    633     float Yo = source->moments->My;
    634 
    635     psEllipseMoments moments;
    636     moments.x2 = source->moments->Mxx;
    637     moments.y2 = source->moments->Myy;
    638     moments.xy = source->moments->Mxy;
    639 
    640     psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
    641     if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
    642         // Shall we log a proper warning? This happens often with matched sources (forced photometry)
    643         // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
    644         return false;
    645     }
    646 
    647     // Why this factor of 0.5 ?
    648     float scale = 0.5 * source->moments->Mrf / axes.major;
    649     axes.major *= scale;
    650     axes.minor *= scale;
    651 
    652     psEllipseShape shape = psEllipseAxesToShape(axes);
    653     if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
    654         // Shall we log a proper warning? This happens often with matched sources (forced photometry)
    655         // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
    656         return false;
    657     }
    658 
    659     float Smajor = axes.major;
    660 
    661     int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    662     int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1);
    663     int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1);
    664     int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
    665 
    666     float rMxx = 0.5 / PS_SQR(shape.sx);
    667     float rMyy = 0.5 / PS_SQR(shape.sy);
    668     float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
    669                                     // implementation. XXX: Is this correct?
    670 
    671     for (int iy = minY; iy < maxY; iy++) {
    672         for (int ix = minX; ix < maxX; ix++) {
    673 
    674             float dX = (ix + 0.5 - Xo);
    675             float dY = (iy + 0.5 - Yo);
    676 
    677             float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
    678 
    679             float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
    680 
    681             kronWindow->data.F32[iy][ix] *= f;
    682         }
    683     }
    684 
    685     return true;
    686 }
Note: See TracChangeset for help on using the changeset viewer.