IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 1, 2012, 11:36:26 AM (14 years ago)
Author:
bills
Message:

Kron Radius and magnitude measurement improvements

  1. Reimplent kron window function using shape and ellipse functions to avoid wildly out of range values
  2. Don't attempt to measure sources with bad or unmeasured moments
  3. two iterations using 6 * Mrf on the second iteration
  4. New recipe values to control make weighting, window application, and smoothing optional.

Make put minimum sky slope value in the recipe SKY_SLOPE_MIN

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotKronIterate.c

    r33690 r33837  
    11# include "psphotInternal.h"
    2 # ifndef ROUND
    3 # define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
    4 # endif
    5 
    6 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert);
    7 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal);
     2
     3bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
     4bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, bool smooth);
    85
    96bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
    107{
    118    bool status = true;
    12 
    13     // return true;
    149
    1510    fprintf (stdout, "\n");
    … …  
    8378    if (!status) {
    8479        MIN_KRON_RADIUS = 0.25*RADIUS;
     80    }
     81
     82    int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS");
     83    if (!status) {
     84        KRON_ITERATIONS = 1;
     85    }
     86
     87    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
     88    if (!status) {
     89        KRON_APPLY_WEIGHT = true;
     90    }
     91
     92    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
     93    if (!status) {
     94        KRON_APPLY_WINDOW = false;
     95    }
     96    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
     97    if (!status) {
     98        KRON_SMOOTH = false;
    8599    }
    86100
    … …  
    114128        // set a window function for each source based on the moments
    115129        // (this skips really bad sources (no peak, no moments, DEFECT)
    116         psphotKronWindowSetSource (source, kronWindow, false, true);
     130        psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
    117131    }
    118132
    … …  
    140154            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,             PS_TYPE_F32);
    141155            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
     156            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
     157            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
     158            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
     159            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_SMOOTH, PS_TYPE_S32);
    142160
    143161// set this to 0 to run without threading
    144 # if (1)
     162# if (0)
    145163            if (!psThreadJobAddPending(job)) {
    146164                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    … …  
    188206    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
    189207    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
    190 
    191     // XXX TEST : set iteration to 1
    192     for (int j = 0; j < 1; j++) {
     208    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
     209    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[8],S32);
     210    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[9],S32);
     211    bool KRON_SMOOTH                = PS_SCALAR_VALUE(job->args->data[10],S32);
     212
     213    for (int j = 0; j < KRON_ITERATIONS; j++) {
    193214        for (int i = 0; i < sources->n; i++) {
    194215
    … …  
    196217            if (!source->peak) continue; // XXX how can we have a peak-less source?
    197218
    198             // allocate space for moments
     219            // check status of this source's moments
    199220            if (!source->moments) continue;
     221            if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
     222            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
    200223
    201224            // replace object in image
    … …  
    219242            // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
    220243            // if we lack the skyRadius (eg MATCHED sources), go to the default value
    221             float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     244            float maxWindow;
     245            if (j == 0) {
     246                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
     247            } else {
     248                maxWindow = isfinite(source->moments->Mrf) ?  6.0*source->moments->Mrf : RADIUS;
     249            }
    222250            float windowRadius = PS_MAX(RADIUS, maxWindow);
    223251
    224 #ifdef notdef
    225             fprintf(stderr, "Redefining pixels for source: %d %4d %4d new radius: %f\n",
    226                                         i, source->peak->x, source->peak->y, windowRadius+2);
    227 #endif
    228252            // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
    229253            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
    … …  
    231255
    232256            // clear the window function for this source based on the moments
    233             psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
    234 
    235             // this function populates moments->Mrf,KronFlux,KronFluxErr
    236             psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
    237             psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    238 
    239             // set a window function for each source based on the moments
    240             psphotKronWindowSetSource (source, kronWindow, true, true);
     257            if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
     258
     259                // this function populates moments->Mrf,KronFlux,KronFluxErr
     260                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, KRON_SMOOTH);
     261                psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
     262
     263                // set a window function for each source based on the moments
     264                psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
     265            }
    241266
    242267            // if we subtracted it above, re-subtract the object, leave local sky
    … …  
    249274}
    250275
    251 bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) {
     276bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
     277bool applyWeight, bool smooth) {
    252278
    253279    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    257283
    258284    psF32 R2 = PS_SQR(radius);
    259     float rsigma2 = 0.5 / PS_SQR(radius/2.0);
     285    float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
    260286
    261287    // a note about coordinates: coordinates of objects throughout psphot refer to the primary
    … …  
    288314    int Ywo = source->pixels->row0;
    289315
    290     psF32 **vPix = source->pixels->data.F32;
     316    psF32 **vPix;
     317    psImage *smoothedImage = NULL;;
     318    if (smooth) {
     319        smoothedImage = psImageCopy(NULL, source->pixels, PS_TYPE_F32);
     320        psImageSmooth(smoothedImage, 1.7, 2);
     321        vPix = smoothedImage->data.F32;
     322    } else {
     323        vPix = source->pixels->data.F32;
     324    }
    291325    psF32 **vWin = kronWindow->data.F32;
    292326    psF32 **vWgt = source->variance->data.F32;
    … …  
    361395    float Var = 0.0;
    362396    float Win = 0.0;
     397
     398    // set vPix to the source pixels (it may have been set to the
     399    // smoothed image above)
     400    vPix = source->pixels->data.F32;
    363401
    364402    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    … …  
    395433    source->moments->KronFlux    = Sum;
    396434    source->moments->KronFluxErr = sqrt(Var);
     435    if (smoothedImage) {
     436        psFree(smoothedImage);
     437    }
    397438
    398439    return true;
    399440}
    400441
    401 bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) {
     442bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
    402443
    403444    if (!source) return false;
    … …  
    406447    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
    407448    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
     449    if (!source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) return false;
     450    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
    408451    psAssert(kronWindow, "need a window");
     452
     453    if (source->moments->Mrf < 0) return false;
     454
     455    if (!applyWindow) {
     456        return true;
     457    }
    409458
    410459    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has
    … …  
    417466    float Yo = source->moments->My;
    418467
    419     float Mxx = source->moments->Mxx;
    420     float Mxy = source->moments->Mxy;
    421     float Myy = source->moments->Myy;
    422 
    423     float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    424     float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
    425 
    426     // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width
    427     float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
    428 
    429     float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
    430     float Sxy = Mxy / (scale * Mmajor * Mminor);
    431     float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2
    432 
    433     float Smajor = sqrt(Mmajor);
     468    psEllipseMoments moments;
     469    moments.x2 = source->moments->Mxx;
     470    moments.y2 = source->moments->Myy;
     471    moments.xy = source->moments->Mxy;
     472
     473    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
     474    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
     475        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
     476        // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
     477        return false;
     478    }
     479
     480    float scale = fabs(0.5 * source->moments->Mrf / axes.major);
     481    axes.major *= scale;
     482    axes.minor *= scale;
     483
     484    psEllipseShape shape = psEllipseAxesToShape(axes);
     485    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
     486        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
     487        // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
     488        return false;
     489    }
     490
     491    float Sxx = PS_SQR(shape.sx);
     492    float Syy = PS_SQR(shape.sy);
     493    float Sxy = -1. * shape.sxy;
     494    float Smajor = axes.major;
    434495
    435496    int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
    … …  
    440501    float rMxx = 0.5 / Sxx;
    441502    float rMyy = 0.5 / Syy;
    442    
     503
    443504    for (int iy = minY; iy < maxY; iy++) {
    444505        for (int ix = minX; ix < maxX; ix++) {
    … …  
    451512            float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
    452513
    453             kronWindow->data.F32[iy][ix] *= f;
     514            kronWindow->data.F32[iy][ix] *= f;
    454515        }
    455516    }
     517
    456518    return true;
    457519}
Note: See TracChangeset for help on using the changeset viewer.