IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33964


Ignore:
Timestamp:
May 30, 2012, 6:06:49 PM (14 years ago)
Author:
bills
Message:

Limit the integration window in the kron radius calculation based on the surface
brightness of the source. Add configuration variable KRON_SB_MIN_FACTOR to control
this limit.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/ippconfig/recipes/psphot.config

    r33937 r33964  
    192192KRON_APPLY_WEIGHT                   BOOL  FALSE
    193193KRON_APPLY_WINDOW                   BOOL  FALSE
    194 KRON_SMOOTH                         BOOL  FALSE
     194KRON_SMOOTH                         BOOL  TRUE
    195195KRON_SMOOTH_SIGMA                   F32   1.7
    196196KRON_SMOOTH_NSIGMA                  S32   2
     197KRON_SB_MIN_FACTOR                  F32   2
    197198
    198199# Extended source fit parameters
    … …  
    364365  PSPHOT.STACK.MATCH.PSF.SOURCE       STR   AUTO # which inputs to convolve? (RAW, CNV, AUTO)
    365366  PSPHOT.STACK.TARGET.PSF.AUTO        BOOL  F    # automatically determine target PSF size?
    366   PSPHOT.STACK.USE.RAW                BOOL  F    # perform photometry & morphology analysis on the raw image
     367  PSPHOT.STACK.USE.RAW                BOOL  T    # perform photometry & morphology analysis on the raw image
     368  KRON_SMOOTH                         BOOL  T    # do not smooth before kron radius measurement
    367369  @PSPHOT.STACK.TARGET.PSF.FWHM       F32   6.0 8.0 # FWHM of target PSF (if NOT AUTO sized; pixels)
    368370  RADIAL_APERTURES                    BOOL  T    # calculate flux in circular radial apertures?
    … …  
    450452
    451453  EXT_FIT_MAX_RADIUS                  F32   50.0
    452   PSF_MODEL                           STR   PS_MODEL_PS1_V1
     454  PSF_MODEL                           STR   PS_MODEL_P1_V1
    453455  EXT_MODEL                           STR   PS_MODEL_SERSIC
    454456END

  • trunk/psphot/src/psphotKronIterate.c

    r33959 r33964  
    108108        KRON_SMOOTH_NSIGMA = 2;
    109109    }
     110    /*
     111     *  Parameter for calculating maximum integration radius based on source's surface
     112     *  brightness
     113     *  Given minimum surface brightness SBmin and a flux the maximum radius is found from
     114     *
     115     *  SBmin = source->flux / (pi * Rmax**2)
     116     *  Rmax = sqrt (source->flux ) / sqrt (SBmin * pi)
     117     *
     118     *  Now what do we use for SBmin?
     119     *
     120     *  SBmin = ( some flux ) / (some area)
     121     *  some flux ~ flux of Ns sigma PSF source
     122     *  some area ~ K times area of a PSF
     123     * 
     124     * flux of Ns sigma source ~ Ns * SKY_STDEV * PSF_EFFECTIVE_AREA
     125     * PSF_EFFECTIVE_AREA = 4 pi sigma_PSF^2
     126     * (the 4 accounts for the fact that 1 sigma is not the total area, it is
     127     * actually a larger region).
     128     *
     129     *   SBmin = Ns * SKY_STDEV * 4 * pi * sigma_PSF^2 / (K * pi * sigma_PSF^2)
     130     *         = Ns * SKY_STDEV * 4 / K
     131     *
     132     * We combine the two parameters Ns and K, and the constant 4 into a single recipe value
     133     * KRON_SB_MIN_FACTOR with Ns = 5 and K = 2 the corresponding value
     134     * for KRON_SB_MIN_FACTOR is 5 * 4 / 2 = 10
     135     */
     136    float KRON_SB_MIN_FACTOR = psMetadataLookupF32 (&status, recipe, "KRON_SB_MIN_FACTOR");
     137    if (!status) {
     138        KRON_SB_MIN_FACTOR = 10;
     139    }
     140    float SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_DEV");
     141    float KRON_SB_MIN_DIVISOR = sqrt ( M_PI * KRON_SB_MIN_FACTOR * SKY_STDEV );
    110142
    111143    // bit-masks to test for good/bad pixels
    … …  
    193225            psArrayAdd(job->args, 1, kronWindow);
    194226            psArrayAdd(job->args, 1, cells->data[j]); // sources
     227            psArrayAdd(job->args, 1, smoothedImage);
    195228            PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
    196229            PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
    … …  
    201234            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
    202235            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
    203             psArrayAdd(job->args, 1, smoothedImage);
     236            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
    204237
    205238// set this to 0 to run without threading
    … …  
    236269    psFree (cellGroups);
    237270    psFree (kronWindow);
     271    if (KRON_SMOOTH) {
     272        for (int i = 0; i < sources->n; i++) {
     273            pmSource *source = sources->data[i];
     274            psFree(source->tmpPtr);
     275            source->tmpPtr = NULL;
     276        }
     277    }
    238278    psFree (smoothedImage);
    239279
    … …  
    247287    psImage *kronWindow             = job->args->data[1];
    248288    psArray *sources                = job->args->data[2];
    249     psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
    250     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
    251     float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
    252     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
    253     int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
    254     bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[8],S32);
    255     bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[9],S32);
    256     float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[10],F32);
    257     psImage *smoothedImage          = job->args->data[11];
    258 
    259     // psImage *smoothedPixels = NULL;
     289    psImage *smoothedImage          = job->args->data[3];
     290    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
     291    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     292    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
     293    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
     294    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
     295    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
     296    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
     297    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
     298    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
     299
    260300    for (int j = 0; j < KRON_ITERATIONS; j++) {
    261301        for (int i = 0; i < sources->n; i++) {
    … …  
    276316                smoothedPixels = (psImage *) source->tmpPtr;
    277317                if (smoothedPixels) {
    278                     // psFree(source->tmpPtr);
    279                     // smoothedPixels = psImageSubset(smoothedImage, source->region);
    280318                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0);
    281                     // source->tmpPtr = smoothedPixels;
    282319                }
    283320                reSubtract = true;
    284321            }
    285322
    286             // use S/N to control max window size
    287             // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
    288 
    289             // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux)
    290             // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2)
    291             // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0)
    292             // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
    293 
    294             // iterate to the window radius
    295             // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
    296 
    297             // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
    298             // if we lack the skyRadius (eg MATCHED sources), go to the default value
     323            // On first iteration set window radius to sky radius (if valid) on second iteration
     324            // use a factor times the previous radial moment value up to a maximum value that
     325            // depends on the surface brightness of the source
    299326            float maxWindow;
    300327            if (j == 0) {
    301328                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
    302329            } else {
    303                 maxWindow = isfinite(source->moments->Mrf) ?  6.0*source->moments->Mrf : RADIUS;
     330                if (KRON_SB_MIN_DIVISOR) {
     331                    if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) {
     332                        // Limit window radius based on surface brightness
     333                        float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;
     334
     335                        if (source->moments->Mrf > 0) {
     336                            maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
     337                        } else {
     338                            maxWindow = Rmax;
     339                        }
     340                    } else {
     341                        maxWindow = RADIUS;
     342                    }
     343                } else {
     344                    // old code
     345                    maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
     346                }
    304347            }
    305348            float windowRadius = PS_MAX(RADIUS, maxWindow);
    … …  
    317360
    318361            // clear the window function for this source based on the moments
    319             // Note this function also applies cuts on the source and returns false if it
     362            // Note: this function also applies cuts on the source and returns false if it
    320363            // does not meet the requirements for measuring the Kron Radius or Magnitude.
    321             // Not that it performs that function even if KRON_APPLY_WINDOW is false
     364            // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
    322365            if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
    323366
    … …  
    335378                if (smoothedPixels) {
    336379                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
    337                     if (j + 1 == KRON_ITERATIONS) {
    338                         // We're done with the smoothedPixels
    339                         psFree(source->tmpPtr);
    340                         source->tmpPtr = NULL;
    341                         smoothedPixels = NULL;
    342                     }
    343380                }
    344381            }
    345382        }
    346383    }
    347     // psFree(smoothedPixels);
     384
    348385    return true;
    349386}

  • trunk/psphot/src/psphotSetThreads.c

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