IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 20, 2010, 1:14:11 PM (16 years ago)
Author:
eugene
Message:

Many changes:

  • psphot-related issues

1) added a new feature to psLib/src/fft/psImageFFT to allow for a pre-FFTed kernel generated for a specific image dimension which can then be applied to an arbitrary number of images for convolution -- this avoids the kernel FFT for every image convolution.
2) updated the psMinimizeLMM to include 2 levels of tolerance: minTol and maxTol: if the minimization reaches minTol, it stops. if it only reaches maxTol within maxIter, it stops and is considered successful, if it fails to reach maxTol, then it fails. This allows us to accept fits that are actually acceptable (within error) even if they are not as close as ideally possible.
3) add new stat: psfWeightNotPoor (vs psfWeightNotBad) : the first gives the fraction of psf-weighted unmasked pixels considering any mask bits (except the internal 'mark' bit), while the former considers only 'bad' mask bits -- these are written to QF_PSF and QF_PSF_PERFECT in the CMF files.
4) define user-set parameters for max and min valid flux in the linear fit analysis: Note this was the cause of the non-negative fluxes in forced photometry -- the min limit was hard-wired to 0.0.
5) significance image is now constructed as (image + (image/1000)2) / sqrt(variance) so that bright sources have well-defined peaks (other wise, they become somewhat flat-topped as image = sqrt(variance), leading to ill-defined peaks).
6) modification of the visualization functions to accept facility and level values akin to psTrace
7) modification of the source fitting APIs to allow thread- and source- independent fit options (iterations, tolerances, etc)
8) use Kron magnitude as test for source size (CR vs EXT) instead of PSF-based aperture
9) set a min systematic error in the aperture mags when used for size classification significance
10) threaded the psf-convolved model (PCM) fitting process (extended source fits)
11) for extended sources, adjust the radius based on the footprint and re-calculate moments for guess; save the psf-based moments for output in psf table
12) for Sersic fitting, choose the best index with a grid search using only a few iterations; iterate fully on the selected value.
13) same for Sersic fitting using the PCM fitting process
14) fixed a bug in which the PSF candidate stars which failed in the psf fitting were not correctly unmarked as being PSF stars
15) define galaxy fit radius based on sky stdev (global, not local)
16) move the PCM code to psModules
17) save and report the raw aperture mag in addition to the curve-of-growth corrected one (PS1_V3)
18) save and report the Kron parameters and higher-order moments
19) some psModule header reorganization for code clarity
20) fixed a bug in which the diff stats were counting 'marked' pixels (outside area of interest) as 'masked'.
21) save the radial profile aperture sizes in the headers
22) better guess for Sersic parameters based on moments (requires setting the functional form so that the scale length is right)

  • ppSub-related:

1) ensure masked pixels are NANed in output diff image
2) add code to flag detections if they have bright positive neighbors (+ output of these in PS1_DV2)
3) define separate penalties for each image (based on their fwhm values) (this requires the penalties to be calculated later in the code).
4) define separate apertures for each image for flux normalization
5) choose aperture based on curve-of-growth (was based on fixed fraction of full aperture flux, and thus noisy)
6) some fine tuning of the penalty factor (this still seems arbitrary, and results are somewhat sensitive to the right value)

Location:
trunk/psModules
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/imcombine/pmSubtractionKernels.c

    r27335 r29004  
    2626    psFree(kernels->vStop);
    2727    psFree(kernels->preCalc);
    28     psFree(kernels->penalties);
     28    psFree(kernels->penalties1);
     29    psFree(kernels->penalties2);
    2930    psFree(kernels->solution1);
    3031    psFree(kernels->solution2);
    … …  
    140141    kernels->v = psVectorRealloc(kernels->v, start + numNew);
    141142    kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew);
    142     kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     143
     144    kernels->penalties1 = psVectorRealloc(kernels->penalties1, start + numNew);
     145    kernels->penalties2 = psVectorRealloc(kernels->penalties2, start + numNew);
     146
    143147    kernels->inner = start;
    144148    kernels->num += numNew;
    … …  
    156160            kernels->v->data.S32[index] = v;
    157161            kernels->preCalc->data[index] = NULL;
    158             kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
    159             psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
     162
     163            // XXX this needs to be changed to use the *convolved* second moment
     164            kernels->penalties1->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
     165            psAssert (isfinite(kernels->penalties1->data.F32[index]), "invalid penalty");
     166
     167            kernels->penalties2->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
     168            psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
     169
    160170            psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
    161171        }
    … …  
    166176    kernels->v->n = start + numNew;
    167177    kernels->preCalc->n = start + numNew;
    168     kernels->penalties->n = start + numNew;
     178
     179    kernels->penalties1->n = start + numNew;
     180    kernels->penalties2->n = start + numNew;
    169181
    170182    return true;
    171183}
    172184
    173 bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
    174                                          int index, int size, int uOrder, int vOrder, float fwhm,
    175                                          bool AlardLuptonStyle, bool forceZeroNull)
     185static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
     186                                                int index, int uOrder, int vOrder, float fwhm,
     187                                                bool AlardLuptonStyle, bool forceZeroNull)
    176188{
    177189    // we have 4 cases here:
    … …  
    182194
    183195    // Calculate moments
    184     double penalty = 0.0;                   // Moment, for penalty
    185196    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
    186197    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    187     for (int v = -size; v <= size; v++) {
    188         for (int u = -size; u <= size; u++) {
     198
     199    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     200        for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
    189201            double value = preCalc->kernel->kernel[v][u];
    190202            double value2 = PS_SQR(value);
    191203            sum += value;
    192204            sum2 += value2;
    193             penalty += value2 * PS_SQR((PS_SQR(u) + PS_SQR(v)));
    194205            min = PS_MIN(value, min);
    195206            max = PS_MAX(value, max);
    … …  
    198209
    199210#if 0
    200     fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
     211    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max);
    201212#endif
    202213
    … …  
    207218        if (uOrder % 2 == 0 && vOrder % 2 == 0) {
    208219            // Even functions: normalise to unit sum and subtract null pixel so that sum is zero
    209             scale2D = 1.0 / fabs(sum);
     220            // Re-normalize
     221            // scale2D  = 1.0 / fabs(sum);
     222            scale2D  = 1.0 / sqrt(sum2);
    210223            zeroNull = true;
    211224        } else {
    … …  
    239252
    240253    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
    241     penalty *= 1.0 / sum2;
    242254
    243255    if (zeroNull) {
    244         preCalc->kernel->kernel[0][0] -= 1.0;
    245     }
    246 
    247 #if 0
     256        // preCalc->kernel->kernel[0][0] -= 1.0;
     257        preCalc->kernel->kernel[0][0] -= sum / sqrt (sum2);
     258    }
     259
     260#if 1
    248261    {
    249         double sum = 0.0;   // Sum of kernel component
     262        double Sum = 0.0;   // Sum of kernel component
     263        double Sum2 = 0.0;   // Sum of kernel component
    250264        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
    251         for (int v = -size; v <= size; v++) {
    252             for (int u = -size; u <= size; u++) {
    253                 sum += preCalc->kernel->kernel[v][u];
     265        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
     266            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
     267                double value = preCalc->kernel->kernel[v][u];
     268                Sum += value;
     269                Sum2 += PS_SQR(value);
    254270                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
    255271                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
    256272            }
    257273        }
    258         fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
     274        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
    259275    }
    260276#endif
    … …  
    267283    }
    268284    kernels->preCalc->data[index] = preCalc;
    269     kernels->penalties->data.F32[index] = kernels->penalty * penalty;
    270     psAssert (isfinite(kernels->penalties->data.F32[index]), "invalid penalty");
    271     psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
     285    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
    272286
    273287    return true;
    … …  
    321335
    322336                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    323                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
    324                 // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false, false);
     337                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
     338                // pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], false, false);
    325339            }
    326340        }
    … …  
    379393            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    380394                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    381                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     395                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    382396            }
    383397        }
    … …  
    385399            // XXX modify size for hermitians to account for sqrt(2) in Hermitian definition (relative to ISIS Gaussian)
    386400            pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS_RADIAL, order, order, size, sigma / sqrt(2.0)); // structure to hold precalculated values
    387             pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, order, order, fwhms->data.F32[i], true, true);
     401            pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, order, order, fwhms->data.F32[i], true, true);
    388402        }
    389403    }
    … …  
    437451            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
    438452                pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
    439                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     453                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    440454            }
    441455        }
    … …  
    506520
    507521                // XXX do we use Alard-Lupton normalization (last param true) or not?
    508                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true, false);
     522                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, uOrder, vOrder, fwhms->data.F32[i], true, false);
    509523
    510524                // XXXX test demo that deconvolved kernel is valid
    … …  
    572586    kernels->preCalc = psArrayAlloc(numBasisFunctions);
    573587    kernels->penalty = penalty;
    574     kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     588    kernels->penalties1 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     589    psVectorInit(kernels->penalties1, NAN);
     590    kernels->penalties2 = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     591    psVectorInit(kernels->penalties2, NAN);
     592    kernels->havePenalties = false;
    575593    kernels->size = size;
    576594    kernels->inner = 0;
    … …  
    771789
    772790    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
    773     psVectorInit(kernels->penalties, 0.0);
     791    psVectorInit(kernels->penalties1, 0.0);
     792    psVectorInit(kernels->penalties2, 0.0);
    774793
    775794    return kernels;
    … …  
    866885
    867886    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
    868     psVectorInit(kernels->penalties, 0.0);
     887    psVectorInit(kernels->penalties1, 0.0);
     888    psVectorInit(kernels->penalties2, 0.0);
    869889
    870890    return kernels;
    … …  
    10401060                kernels->u->data.S32[index] = uOrder;
    10411061                kernels->v->data.S32[index] = vOrder;
    1042                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
    1043                 if (!isfinite(kernels->penalties->data.F32[index])) {
     1062
     1063                // XXX convert to use the convolved 2nd moment
     1064                kernels->penalties1->data.F32[index] = kernels->penalty * fabsf(moment);
     1065                if (!isfinite(kernels->penalties1->data.F32[index])) {
     1066                    psAbort ("invalid penalty");
     1067                }
     1068                kernels->penalties2->data.F32[index] = kernels->penalty * fabsf(moment);
     1069                if (!isfinite(kernels->penalties2->data.F32[index])) {
    10441070                    psAbort ("invalid penalty");
    10451071                }
    … …  
    12471273    out->preCalc = psMemIncrRefCounter(in->preCalc);
    12481274    out->penalty = in->penalty;
    1249     out->penalties = psMemIncrRefCounter(in->penalties);
     1275    out->penalties1 = psMemIncrRefCounter(in->penalties1);
     1276    out->penalties2 = psMemIncrRefCounter(in->penalties2);
    12501277    out->uStop = psMemIncrRefCounter(in->uStop);
    12511278    out->vStop = psMemIncrRefCounter(in->vStop);
Note: See TracChangeset for help on using the changeset viewer.