Index: /trunk/ippconfig/recipes/psphot.config =================================================================== --- /trunk/ippconfig/recipes/psphot.config (revision 33963) +++ /trunk/ippconfig/recipes/psphot.config (revision 33964) @@ -192,7 +192,8 @@ KRON_APPLY_WEIGHT BOOL FALSE KRON_APPLY_WINDOW BOOL FALSE -KRON_SMOOTH BOOL FALSE +KRON_SMOOTH BOOL TRUE KRON_SMOOTH_SIGMA F32 1.7 KRON_SMOOTH_NSIGMA S32 2 +KRON_SB_MIN_FACTOR F32 2 # Extended source fit parameters @@ -364,5 +365,6 @@ PSPHOT.STACK.MATCH.PSF.SOURCE STR AUTO # which inputs to convolve? (RAW, CNV, AUTO) PSPHOT.STACK.TARGET.PSF.AUTO BOOL F # automatically determine target PSF size? - PSPHOT.STACK.USE.RAW BOOL F # perform photometry & morphology analysis on the raw image + PSPHOT.STACK.USE.RAW BOOL T # perform photometry & morphology analysis on the raw image + KRON_SMOOTH BOOL T # do not smooth before kron radius measurement @PSPHOT.STACK.TARGET.PSF.FWHM F32 6.0 8.0 # FWHM of target PSF (if NOT AUTO sized; pixels) RADIAL_APERTURES BOOL T # calculate flux in circular radial apertures? @@ -450,5 +452,5 @@ EXT_FIT_MAX_RADIUS F32 50.0 - PSF_MODEL STR PS_MODEL_PS1_V1 + PSF_MODEL STR PS_MODEL_P1_V1 EXT_MODEL STR PS_MODEL_SERSIC END Index: /trunk/psphot/src/psphotKronIterate.c =================================================================== --- /trunk/psphot/src/psphotKronIterate.c (revision 33963) +++ /trunk/psphot/src/psphotKronIterate.c (revision 33964) @@ -108,4 +108,36 @@ KRON_SMOOTH_NSIGMA = 2; } + /* + * Parameter for calculating maximum integration radius based on source's surface + * brightness + * Given minimum surface brightness SBmin and a flux the maximum radius is found from + * + * SBmin = source->flux / (pi * Rmax**2) + * Rmax = sqrt (source->flux ) / sqrt (SBmin * pi) + * + * Now what do we use for SBmin? + * + * SBmin = ( some flux ) / (some area) + * some flux ~ flux of Ns sigma PSF source + * some area ~ K times area of a PSF + * + * flux of Ns sigma source ~ Ns * SKY_STDEV * PSF_EFFECTIVE_AREA + * PSF_EFFECTIVE_AREA = 4 pi sigma_PSF^2 + * (the 4 accounts for the fact that 1 sigma is not the total area, it is + * actually a larger region). + * + * SBmin = Ns * SKY_STDEV * 4 * pi * sigma_PSF^2 / (K * pi * sigma_PSF^2) + * = Ns * SKY_STDEV * 4 / K + * + * We combine the two parameters Ns and K, and the constant 4 into a single recipe value + * KRON_SB_MIN_FACTOR with Ns = 5 and K = 2 the corresponding value + * for KRON_SB_MIN_FACTOR is 5 * 4 / 2 = 10 + */ + float KRON_SB_MIN_FACTOR = psMetadataLookupF32 (&status, recipe, "KRON_SB_MIN_FACTOR"); + if (!status) { + KRON_SB_MIN_FACTOR = 10; + } + float SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_DEV"); + float KRON_SB_MIN_DIVISOR = sqrt ( M_PI * KRON_SB_MIN_FACTOR * SKY_STDEV ); // bit-masks to test for good/bad pixels @@ -193,4 +225,5 @@ psArrayAdd(job->args, 1, kronWindow); psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, smoothedImage); PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); @@ -201,5 +234,5 @@ PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32); PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA, PS_TYPE_F32); - psArrayAdd(job->args, 1, smoothedImage); + PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32); // set this to 0 to run without threading @@ -236,4 +269,11 @@ psFree (cellGroups); psFree (kronWindow); + if (KRON_SMOOTH) { + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + psFree(source->tmpPtr); + source->tmpPtr = NULL; + } + } psFree (smoothedImage); @@ -247,15 +287,15 @@ psImage *kronWindow = job->args->data[1]; psArray *sources = job->args->data[2]; - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); - float RADIUS = PS_SCALAR_VALUE(job->args->data[5],F32); - float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); - int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[7],S32); - bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[8],S32); - bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[9],S32); - float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[10],F32); - psImage *smoothedImage = job->args->data[11]; - - // psImage *smoothedPixels = NULL; + psImage *smoothedImage = job->args->data[3]; + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + float RADIUS = PS_SCALAR_VALUE(job->args->data[6],F32); + float MIN_KRON_RADIUS = PS_SCALAR_VALUE(job->args->data[7],F32); + int KRON_ITERATIONS = PS_SCALAR_VALUE(job->args->data[8],S32); + bool KRON_APPLY_WEIGHT = PS_SCALAR_VALUE(job->args->data[9],S32); + bool KRON_APPLY_WINDOW = PS_SCALAR_VALUE(job->args->data[10],S32); + float KRON_SMOOTH_SIGMA = PS_SCALAR_VALUE(job->args->data[11],F32); + float KRON_SB_MIN_DIVISOR = PS_SCALAR_VALUE(job->args->data[12],F32); + for (int j = 0; j < KRON_ITERATIONS; j++) { for (int i = 0; i < sources->n; i++) { @@ -276,30 +316,33 @@ smoothedPixels = (psImage *) source->tmpPtr; if (smoothedPixels) { - // psFree(source->tmpPtr); - // smoothedPixels = psImageSubset(smoothedImage, source->region); pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0); - // source->tmpPtr = smoothedPixels; } reSubtract = true; } - // use S/N to control max window size - // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr; - - // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux) - // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2) - // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0) - // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS; - - // iterate to the window radius - // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow); - - // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS - // if we lack the skyRadius (eg MATCHED sources), go to the default value + // On first iteration set window radius to sky radius (if valid) on second iteration + // use a factor times the previous radial moment value up to a maximum value that + // depends on the surface brightness of the source float maxWindow; if (j == 0) { maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS; } else { - maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; + if (KRON_SB_MIN_DIVISOR) { + if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) { + // Limit window radius based on surface brightness + float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR; + + if (source->moments->Mrf > 0) { + maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax); + } else { + maxWindow = Rmax; + } + } else { + maxWindow = RADIUS; + } + } else { + // old code + maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS; + } } float windowRadius = PS_MAX(RADIUS, maxWindow); @@ -317,7 +360,7 @@ // clear the window function for this source based on the moments - // Note this function also applies cuts on the source and returns false if it + // Note: this function also applies cuts on the source and returns false if it // does not meet the requirements for measuring the Kron Radius or Magnitude. - // Not that it performs that function even if KRON_APPLY_WINDOW is false + // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) { @@ -335,15 +378,9 @@ if (smoothedPixels) { pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0); - if (j + 1 == KRON_ITERATIONS) { - // We're done with the smoothedPixels - psFree(source->tmpPtr); - source->tmpPtr = NULL; - smoothedPixels = NULL; - } } } } } - // psFree(smoothedPixels); + return true; } Index: /trunk/psphot/src/psphotSetThreads.c =================================================================== --- /trunk/psphot/src/psphotSetThreads.c (revision 33963) +++ /trunk/psphot/src/psphotSetThreads.c (revision 33964) @@ -30,5 +30,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 12); + task = psThreadTaskAlloc("PSPHOT_KRON_ITERATE", 13); task->function = &psphotKronIterate_Threaded; psThreadTaskAdd(task);