IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34559 for trunk


Ignore:
Timestamp:
Oct 24, 2012, 1:50:31 PM (14 years ago)
Author:
bills
Message:

Changes to Kron flux measurements

  1. Insure that the source's window radius > Kron radius. (It wasn't in some cases previously)
  2. Remove the window and weight functions as they are no longer used
  3. Re-measure KronFinner
  4. Re-measure KronFouter if the sources window radius is large enough, otherwise set to NAN
Location:
trunk/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotKronIterate.c

    r34542 r34559  
    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, bool measureRadius);
    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
    … …  
    170169    }
    171170
    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);
    185     }
    186 
     171    // pass 3 we only measure fluxes for matched sources so we don't need smooth pixels or multiple iterations
    187172    if (pass == 3) {
    188173        KRON_SMOOTH = false;
    … …  
    240225
    241226            psArrayAdd(job->args, 1, readout);
    242             psArrayAdd(job->args, 1, kronWindow);
    243227            psArrayAdd(job->args, 1, cells->data[j]); // sources
    244228            psArrayAdd(job->args, 1, smoothedImage);
    … …  
    248232            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
    249233            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
    250             PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
    251             PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
    252234            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
    253235            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
    … …  
    286268    }
    287269    psFree (cellGroups);
    288     psFree (kronWindow);
    289270    if (KRON_SMOOTH) {
    290271        for (int i = 0; i < sources->n; i++) {
    … …  
    303284
    304285    pmReadout *readout              = job->args->data[0];
    305     psImage *kronWindow             = job->args->data[1];
    306     psArray *sources                = job->args->data[2];
    307     psImage *smoothedImage          = job->args->data[3];
    308     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    309     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
    310     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
    311     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
    312     int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
    313     bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
    314     bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
    315     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
    316     float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
    317     int pass                        = PS_SCALAR_VALUE(job->args->data[13],S32);
     286    psArray *sources                = job->args->data[1];
     287    psImage *smoothedImage          = job->args->data[2];
     288    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
     289    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     290    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
     291    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
     292    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
     293    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[8],F32);
     294    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[9],F32);
     295    int pass                        = PS_SCALAR_VALUE(job->args->data[10],S32);
    318296
    319297    bool measureRadius = true;
    … …  
    327305            pmSource *source = sources->data[i];
    328306            if (!source->peak) continue; // XXX how can we have a peak-less source?
    329 
    330             if (pass == 3) {
    331                 // in pass 3 we only measure the flux for matched sources
    332                 if (!(source->mode2 & PM_SOURCE_MODE2_MATCHED)) {
    333                     continue;
    334                 }
    335             }
     307            if (!source->moments) continue;
     308            if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
     309            if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     310            // skip saturated stars modeled with a radial profile
     311            if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
     312
     313            // in pass 3 we only measure the flux for matched sources
     314            if ((pass == 3) && !(source->mode2 & PM_SOURCE_MODE2_MATCHED)) continue;
    336315
    337316# if (0)
    338317# define TEST_X 653
    339318# define TEST_Y 466
    340         if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
    341             fprintf (stderr, "test object\n");
    342         }
     319            if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
     320                fprintf (stderr, "test object\n");
     321            }
    343322# undef TEST_X
    344323# undef TEST_Y
    … …  
    347326            if (measureRadius) {
    348327                // check status of this source's moments
    349                 if (!source->moments) continue;
     328                // XXX: I don't think that we have to apply these restrictions since we dropped the window function
    350329                if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
    351330                if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
    … …  
    354333            if (!isfinite(source->moments->Mrf)) {
    355334                // Once we save a bad Mrf measurement we give up on this source
    356                 // checking here allows us to avoid adding and subtracting the model
     335                // XXX: is this the right thing to do?
    357336                continue;
    358337            }
    359 
    360             // skip saturated stars modeled with a radial profile
    361             if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
    362338
    363339            // replace object in image
    … …  
    402378            float windowRadius = PS_MAX(RADIUS, maxWindow);
    403379
    404 #ifdef REVERT_ON_BAD_MEASURMENT
    405             // save previous measurements. We might revert back to them if this round fails
    406             float MrfPrior = source->moments->Mrf;
    407             float KronFluxPrior = source->moments->KronFlux;
    408             float KronFluxErrPrior = source->moments->KronFluxErr;
    409 #endif
    410 
    411380            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    412381            bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    413             psAssert (source->pixels, "WTF?");
     382            psAssert (source->pixels, "redefine pixels failed?");
    414383            if (extend && smoothedPixels) {
    415384                psFree(source->tmpPtr);
    416385                smoothedPixels = psImageSubset(smoothedImage, source->region);
    417                 psAssert (smoothedPixels, "WTF?");
     386                psAssert (smoothedPixels, "redefine smoothed pixels failed?");
    418387                source->tmpPtr = (psPtr) smoothedPixels ;
    419388            }
    420389
    421 
    422             // clear the window function for this source based on the moments
    423             // Note: this function also applies cuts on the source and returns false if it
    424             // does not meet the requirements for measuring the Kron Radius or Magnitude.
    425             // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
    426             if (!measureRadius || psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
    427 
    428                 // this function populates moments->Mrf,KronFlux,KronFluxErr
    429                 psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal,
    430                     KRON_APPLY_WEIGHT, smoothedPixels, measureRadius);
    431 
    432                 psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    433 
    434                 // set a window function for each source based on the moments
    435                 psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
    436             }
    437 
    438 #ifdef REVERT_ON_BAD_MEASUREMENT
    439             // on pass 2 if we get an invalid measurement on a pass 1 source where we had a good one previously
    440             // in pass 1 keep that measurement
    441             bool reverted = false;
    442             if (pass > 1 && !isfinite(source->moments->Mrf) && (source->mode2 & PM_SOURCE_MODE2_PASS1_SRC)) {
    443                 source->moments->Mrf = MrfPrior;    // This is finite otherwise we wouldn't have gotten here
    444                 source->moments->KronFlux = KronFluxPrior;
    445                 source->moments->KronFluxErr = KronFluxErrPrior;
    446                 reverted = true;
    447             }
    448 #endif
     390            bool measureFlux = true;
     391            if (measureRadius) {
     392                measureFlux = psphotKronRadiusMeasure (source, windowRadius, MIN_KRON_RADIUS, maskVal, smoothedPixels);
     393            }
     394
     395            if (measureFlux) {
     396                // Make sure the sources images are large enough for the measured Kron Radius
     397                bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y,
     398                    2.5 * source->moments->Mrf + 2);
     399
     400                psAssert (source->pixels, "redefine pixels failed?");
     401                if (extend && smoothedPixels) {
     402                    psFree(source->tmpPtr);
     403                    smoothedPixels = psImageSubset(smoothedImage, source->region);
     404                    psAssert (smoothedPixels, "redefine smoothed pixels failed?");
     405                    source->tmpPtr = (psPtr) smoothedPixels ;
     406                }
     407
     408                // this function populates moments->Mrf,KronFlux,KronFluxErr, KronFinner, and KronFouter
     409                psphotKronFluxMeasure (source, maskVal);
     410            }
     411
     412            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     413
    449414#ifdef DUMP_KRS
    450 #ifndef REVERT_ON_BAD_MEASUREMENT
    451             bool reverted = false;
    452 #endif
    453             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);
     415            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);
    454416#endif
    455417
    … …  
    467429}
    468430
    469 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
    470     bool applyWeight, psImage *smoothedPixels, bool measureRadius) {
     431bool psphotKronRadiusMeasure (pmSource *source, float radius, float minKronRadius, psImageMaskType maskVal,
     432    psImage *smoothedPixels) {
    471433
    472434    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    476438
    477439    psF32 R2 = PS_SQR(radius);
    478     float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
    479440
    480441    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    489450    // Xn  = SUM (x - xc)^n * (z - sky)
    490451
    491 
    492452    psF32 RF = 0.0;
    493453    psF32 RS = 0.0;
    … …  
    505465    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    506466
    507     int Xwo = source->pixels->col0;
    508     int Ywo = source->pixels->row0;
    509 
    510467    psF32 **vPix;
    511468    if (smoothedPixels) {
    … …  
    514471        vPix = source->pixels->data.F32;
    515472    }
    516     psF32 **vWin = kronWindow->data.F32;
    517     psF32 **vWgt = source->variance->data.F32;
    518473   
    519474    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    520475
    521     if (measureRadius) {
     476    {
    522477        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    523478
    … …  
    551506                if (r2 > R2) continue;
    552507
    553                 // flux * window
    554                 float z = r2 * rsigma2;
    555                 assert (z >= 0.0);
    556 
    557                 // weight by window image and wide Gaussian
    558                 float weight1  = vWin[row+Ywo][col+Xwo]*exp(-z);
    559                 float weight2  = vWin[yFlip+Ywo][xFlip+Xwo]*exp(-z);
    560 
    561                 float fDiff1 = vPix[row][col]*weight1;
    562                 float fDiff2 = vPix[yFlip][xFlip]*weight2;
     508                float fDiff1 = vPix[row][col];
     509                float fDiff2 = vPix[yFlip][xFlip];
    563510
    564511                float pDiff = (fDiff1 > 0.0) ? sqrt(fabs(fDiff1*fDiff2)) : -sqrt(fabs(fDiff1*fDiff2));
    565512
    566                 // Kron Flux uses the 1st radial moment (maybe Gaussian windowed?)
     513                // Kron Flux uses the 1st radial moment
    567514                psF32 rf = pDiff * sqrt(r2);
    568515                psF32 rs = 0.5 * (fDiff1 + fDiff2);
    … …  
    589536        source->moments->Mrf = Mrf;
    590537    }
    591 
    592     // Calculate the Kron magnitude (make this block optional?)
    593     float radKron  = 2.5 * source->moments->Mrf;
    594     float radKron2 = radKron*radKron;
    595 
    596     int nKronPix = 0;
     538    return true;
     539}
     540
     541bool psphotKronFluxMeasure(pmSource *source, psImageMaskType maskVal) {
     542
     543    PS_ASSERT_PTR_NON_NULL(source, false);
     544    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     545    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     546
     547    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
     548    // image coordinates.  the source->pixels image has an offset relative to its parent of
     549    // col0,row0: a pixel (x,y) in the primary image has coordinates of (x-col0, y-row0) in
     550    // this subimage.  we subtract off the peak coordinates, adjusted to this subimage, to have
     551    // minimal round-off error in the sums.  since these values are subtracted just to minimize
     552    // the dynamic range and are added back below, the exact value does not matter. these are
     553    // (int) so they can be used in the image index below.
     554
     555    // Now calculate higher-order moments, using the above-calculated first moments to adjust coordinates
     556    // Xn  = SUM (x - xc)^n * (z - sky)
     557
     558    // the peak position is less accurate but less subject to extreme deviations
     559    float dX = source->moments->Mx - source->peak->xf;
     560    float dY = source->moments->My - source->peak->yf;
     561    float dR = hypot(dX, dY);
     562    float Xo = (dR < 2.0) ? source->moments->Mx : source->peak->xf;
     563    float Yo = (dR < 2.0) ? source->moments->My : source->peak->yf;
     564
     565    // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     566    // xCM, yCM from pixel coords to pixel index here.
     567    psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     568    psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     569
     570    psF32 **vPix = source->pixels->data.F32;
     571    psF32 **vWgt = source->variance->data.F32;
     572   
     573    psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     574
     575    source->moments->KronFlux  = NAN;
     576    source->moments->KronFluxErr  = NAN;
     577    source->moments->KronFinner = NAN;
     578    source->moments->KronFouter = NAN;
     579
     580    // Calculate the Kron fluxes
     581    float radKinner = 1.0 * source->moments->Mrf;
     582    float radKron   = 2.5 * source->moments->Mrf;
     583    float radKouter = 4.0 * source->moments->Mrf;
     584
     585    float limitRadius = MIN (radKouter, source->windowRadius);
     586    if (radKouter > source->windowRadius) {
     587        // This happens but the measurement isn't important enough to allocate the extra pixels
     588        // psWarning ("outer kron radius: %f is larger than windowRadius: %f for %d\n",
     589        //     radKouter, source->windowRadius, source->id);
     590        limitRadius = MIN (radKron, source->windowRadius);
     591    }
     592    if (radKron > source->windowRadius) {
     593        // caller should have prevented this from happening
     594        psWarning ("kron radius: %f is larger than windowRadius: %f for %d\n",
     595            radKron, source->windowRadius, source->id);
     596        return false;
     597    }
     598
    597599    float Sum = 0.0;
    598600    float Var = 0.0;
    599     float Win = 0.0;
    600 
    601     // set vPix to the source pixels (it may have been set to the
    602     // smoothed image above)
     601    float SumInner = 0.0;
     602    float SumOuter = 0.0;
     603
     604    // set vPix to the source pixels (it may have been set to the smoothed image above)
    603605    vPix = source->pixels->data.F32;
    604606
    … …  
    606608
    607609        psF32 yDiff = row - yCM;
    608         if (fabs(yDiff) > radKron) continue;
     610        if (fabs(yDiff) > limitRadius) continue;
    609611
    610612        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    … …  
    614616
    615617            psF32 xDiff = col - xCM;
    616             if (fabs(xDiff) > radKron) continue;
    617 
    618             // radKron is just a function of (xDiff, yDiff)
     618            if (fabs(xDiff) > limitRadius) continue;
     619
    619620            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
    620             if (r2 > radKron2) continue;
    621 
    622             float weight1  = vWin[row+Ywo][col+Xwo];
    623             float fDiff1 = vPix[row][col]*weight1;
    624 
    625             float pDiff = fDiff1;
    626             psF32 wDiff = vWgt[row][col] * weight1;
    627 
    628             Sum += pDiff;
    629             Var += wDiff;
    630             Win += weight1;
    631             nKronPix ++;
     621            psF32 r = sqrt(r2);
     622
     623            float pDiff = vPix[row][col];
     624            psF32 wDiff = vWgt[row][col];
     625
     626            if (r > radKinner && r < radKron) {
     627                SumInner += pDiff;
     628            }
     629            if (r < radKron) {
     630                Sum += pDiff;
     631                Var += wDiff;
     632            }
     633            if (r > radKron && r < radKouter) {
     634                SumOuter += pDiff;
     635            }
    632636        }
    633637    }
    … …  
    635639    source->moments->KronFlux    = Sum;
    636640    source->moments->KronFluxErr = sqrt(Var);
     641    source->moments->KronFinner  = SumInner;
     642
     643    // only save radKouter if the radius is inside the integration radius limit
     644    if (radKouter <= limitRadius) {
     645        source->moments->KronFouter = SumOuter;
     646    }
    637647
    638648    return true;
    639649}
    640 
    641 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
    642 
    643     if (!source) return false;
    644     if (!source->peak) return false; // XXX how can we have a peak-less source?
    645     if (!source->moments) return false;
    646     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    647     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
    648     if (source->mode2 & PM_SOURCE_MODE2_MATCHED) return true;
    649 
    650     if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
    651     if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    652     psAssert(kronWindow, "need a window");
    653 
    654     // XXX: If we are not applying the window then we don't need to check for valid Mrf here.
    655     // We should give the this module a chance to measure a good value.
    656     // However experiments show that it hardly ever succeeds in getting a better value
    657     if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
    658 
    659     if (!applyWindow) {
    660         return true;
    661     }
    662 
    663     // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has
    664     // these values (and peak of 1.0)
    665 
    666     int Nx = kronWindow->numCols;
    667     int Ny = kronWindow->numRows;
    668 
    669     float Xo = source->moments->Mx;
    670     float Yo = source->moments->My;
    671 
    672     psEllipseMoments moments;
    673     moments.x2 = source->moments->Mxx;
    674     moments.y2 = source->moments->Myy;
    675     moments.xy = source->moments->Mxy;
    676 
    677     psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
    678     if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
    679         // Shall we log a proper warning? This happens often with matched sources (forced photometry)
    680         // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
    681         return false;
    682     }
    683 
    684     // Why this factor of 0.5 ?
    685     float scale = 0.5 * source->moments->Mrf / axes.major;
    686     axes.major *= scale;
    687     axes.minor *= scale;
    688 
    689     psEllipseShape shape = psEllipseAxesToShape(axes);
    690     if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
    691         // Shall we log a proper warning? This happens often with matched sources (forced photometry)
    692         // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
    693         return false;
    694     }
    695 
    696     float Smajor = axes.major;
    697 
    698     int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    699     int maxX = PS_MIN(PS_MAX(Xo + 5*Smajor, 0), Nx - 1);
    700     int minY = PS_MIN(PS_MAX(Yo - 5*Smajor, 0), Ny - 1);
    701     int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
    702 
    703     float rMxx = 0.5 / PS_SQR(shape.sx);
    704     float rMyy = 0.5 / PS_SQR(shape.sy);
    705     float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
    706                                     // implementation. XXX: Is this correct?
    707 
    708     for (int iy = minY; iy < maxY; iy++) {
    709         for (int ix = minX; ix < maxX; ix++) {
    710 
    711             float dX = (ix + 0.5 - Xo);
    712             float dY = (iy + 0.5 - Yo);
    713 
    714             float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
    715 
    716             float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
    717 
    718             kronWindow->data.F32[iy][ix] *= f;
    719         }
    720     }
    721 
    722     return true;
    723 }

  • trunk/psphot/src/psphotSetThreads.c

    r34317 r34559  
    3030    psFree(task);
    3131
    32     task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 14);
     32    task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 11);
    3333    task->function = &psphotKronIterate_Threaded;
    3434    psThreadTaskAdd(task);
Note: See TracChangeset for help on using the changeset viewer.