IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26505


Ignore:
Timestamp:
Jan 4, 2010, 10:06:09 AM (16 years ago)
Author:
eugene
Message:

add skyErr, a way of increasing the constant portion of the variance in the kernel analysis: var = (readNoise2 + skyErr + gain*flux + sysErr*flux2); clean up testing verbosity; temporary dump residuals and quit in visualization code

Location:
branches/eam_branches/20091201/psModules/src/imcombine
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

    r26502 r26505  
    16951695        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
    16961696        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
    1697     }
    1698 
    1699     pmSubtractionVisualShowFit();
    1700     pmSubtractionVisualPlotFit(kernels);
     1697
     1698        pmSubtractionVisualShowFit(norm);
     1699        pmSubtractionVisualPlotFit(kernels);
     1700    }
    17011701
    17021702    psFree(residual);

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

    r26502 r26505  
    177177        }
    178178    }
    179 #if 1
     179#if 0
    180180    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, moment);
    181181#endif
    … …  
    211211    }
    212212
    213 #if 1
     213#if 0
    214214    sum = 0.0;   // Sum of kernel component
    215215    min = FLT_MAX;
    … …  
    409409                preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
    410410
    411                 pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false);
     411                // XXX do we use Alard-Lupton normalization (last param true) or not?
     412                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
    412413
    413414                // XXXX test demo that deconvolved kernel is valid

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c

    r26469 r26505  
    7575                                 float stampSpacing, // Spacing between stamps
    7676                                 float sysError,     // Relative systematic error in images
     77                                 float skyError,     // Relative systematic error in images
    7778                                 int size,         // Kernel half-size
    7879                                 int footprint,     // Convolution footprint for stamps
    … …  
    8586
    8687    *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2,
    87                                       size, footprint, stampSpacing, sysError, mode);
     88                                      size, footprint, stampSpacing, sysError, skyError, mode);
    8889    if (!*stamps) {
    8990        psError(psErrorCodeLast(), false, "Unable to find stamps.");
    … …  
    109110                                  int stride, // Size for convolution patches
    110111                                  float sysError,           // Systematic error in images
     112                                  float skyError,           // Systematic error in images
    111113                                  float kernelError, // Systematic error in kernel
    112114                                  psImageMaskType maskVal, // Value to mask for input
    … …  
    159161        PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false);
    160162    }
     163    if (isfinite(sysError)) {
     164        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false);
     165    }
    161166    if (isfinite(kernelError)) {
    162167        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false);
    … …  
    276281    }
    277282
    278     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, kernelError,
     283    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, NAN, NAN, kernelError,
    279284                               maskVal, maskBad, maskPoor, poorFrac, badFrac, mode)) {
    280285        psFree(kernels);
    … …  
    342347                        int inner, int ringsOrder, int binning, float penalty,
    343348                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    344                         int iter, float rej, float sysError, float kernelError, psImageMaskType maskVal,
     349                        int iter, float rej, float sysError, float skyError, float kernelError, psImageMaskType maskVal,
    345350                        psImageMaskType maskBad, psImageMaskType maskPoor, float poorFrac,
    346351                        float badFrac, pmSubtractionMode subMode)
    347352{
    348     if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, kernelError,
     353    if (!subtractionMatchCheck(conv1, conv2, ro1, ro2, stride, sysError, skyError, kernelError,
    349354                               maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) {
    350355        return false;
    … …  
    476481            if (stampsName && strlen(stampsName) > 0) {
    477482                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, size,
    478                                                         footprint, stampSpacing, sysError, subMode);
     483                                                        footprint, stampSpacing, sysError, skyError, subMode);
    479484            } else if (sources) {
    480485                stamps = pmSubtractionStampsSetFromSources(sources, ro1->image, subMask, region, size,
    481                                                            footprint, stampSpacing, sysError, subMode);
     486                                                           footprint, stampSpacing, sysError, skyError, subMode);
    482487            }
    483488
    … …  
    485490            // doesn't matter.
    486491            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
    487                                       stampSpacing, sysError, size, footprint, subMode)) {
     492                                      stampSpacing, sysError, skyError, size, footprint, subMode)) {
    488493                goto MATCH_ERROR;
    489494            }
    … …  
    557562
    558563                if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
    559                                           stampThresh1, stampThresh2, stampSpacing, sysError,
     564                                          stampThresh1, stampThresh2, stampSpacing, sysError, skyError,
    560565                                          size, footprint, subMode)) {
    561566                    goto MATCH_ERROR;

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h

    r26035 r26505  
    4040                        float rej,      ///< Rejection threshold
    4141                        float sysError, ///< Relative systematic error in images
     42                        float skyError, ///< Relative systematic error in images
    4243                        float kernelError, ///< Relative systematic error in kernel
    4344                        psImageMaskType maskVal, ///< Value to mask for input

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c

    r26502 r26505  
    181181
    182182pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
    183                                                     int footprint, float spacing, float sysErr)
     183                                                    int footprint, float spacing, float sysErr, float skyErr)
    184184{
    185185    pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
    … …  
    224224    list->footprint = footprint;
    225225    list->sysErr = sysErr;
     226    list->skyErr = skyErr;
    226227
    227228    return list;
    … …  
    326327                                                const psImage *image2, const psImage *subMask,
    327328                                                const psRegion *region, float thresh1, float thresh2,
    328                                                 int size, int footprint, float spacing, float sysErr,
     329                                                int size, int footprint, float spacing, float sysErr, float skyErr,
    329330                                                pmSubtractionMode mode)
    330331{
    … …  
    383384
    384385    if (!stamps) {
    385         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr);
     386        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
    386387    }
    387388
    … …  
    441442                    xStamp = xList->data.F32[j];
    442443                    yStamp = yList->data.F32[j];
    443                     fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
     444                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
    444445                }
    445446            } else {
    … …  
    496497                                               const psImage *image, const psImage *subMask,
    497498                                               const psRegion *region, int size, int footprint,
    498                                                float spacing, float sysErr, pmSubtractionMode mode)
     499                                               float spacing, float sysErr, float skyErr, pmSubtractionMode mode)
    499500
    500501{
    … …  
    517518    pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
    518519                                                                 region, footprint, spacing,
    519                                                                  sysErr); // Stamp list
     520                                                                 sysErr, skyErr); // Stamp list
    520521    int numStamps = stamps->num;        // Number of stamps
    521522
    … …  
    548549        }
    549550
    550         fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
     551        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
    551552
    552553        bool found = false;
    … …  
    742743            return false;
    743744        }
    744         fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
     745        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
    745746
    746747        // Catch memory leaks --- these should have been freed and NULLed before
    … …  
    765766            psImage *varSub = psImageSubset(variance, region); // Subimage with stamp
    766767            psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp
    767             if (isfinite(stamps->sysErr) && stamps->sysErr > 0) {
     768
     769            if ((isfinite(stamps->skyErr) && (stamps->skyErr > 0)) ||
     770                (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) {
    768771                float sysErr = 0.25 * PS_SQR(stamps->sysErr); // Systematic error
     772                float skyErr = stamps->skyErr;
    769773                psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images
    770774                for (int y = -size; y <= size; y++) {
    771775                    for (int x = -size; x <= size; x++) {
    772776                        float additional = image1->kernel[y][x] + image2->kernel[y][x];
    773                         weight->kernel[y][x] = 1.0 / (var->kernel[y][x] + sysErr * PS_SQR(additional));
     777                        weight->kernel[y][x] = 1.0 / (skyErr + var->kernel[y][x] + sysErr * PS_SQR(additional));
    774778                    }
    775779                }
    … …  
    795799pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
    796800                                                          const psImage *subMask, const psRegion *region,
    797                                                           int size, int footprint, float spacing, float sysErr,
     801                                                          int size, int footprint, float spacing, float sysErr, float skyErr,
    798802                                                          pmSubtractionMode mode)
    799803{
    … …  
    819823            y->data.F32[numOut] = source->peak->yf;
    820824        }
    821         fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
     825        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
    822826        numOut++;
    823827    }
    … …  
    826830
    827831    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
    828                                                             footprint, spacing, sysErr, mode); // Stamps
     832                                                            footprint, spacing, sysErr, skyErr, mode); // Stamps
    829833    psFree(x);
    830834    psFree(y);
    … …  
    840844pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
    841845                                                       const psImage *subMask, const psRegion *region,
    842                                                        int size, int footprint, float spacing, float sysErr,
     846                                                       int size, int footprint, float spacing, float sysErr, float skyErr,
    843847                                                       pmSubtractionMode mode)
    844848{
    … …  
    858862
    859863    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
    860                                                             spacing, sysErr, mode);
     864                                                            spacing, sysErr, skyErr, mode);
    861865    psFree(data);
    862866

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h

    r26342 r26505  
    2626    psKernel *window;                   ///< window function generated from ensemble of stamps
    2727    float sysErr;                       ///< Systematic error
     28    float skyErr;                       ///< Systematic error
    2829} pmSubtractionStampList;
    2930
    … …  
    3435                                                    int footprint, // Half-size of stamps
    3536                                                    float spacing, // Rough average spacing between stamps
    36                                                     float sysErr  // Relative systematic error or NAN
     37                                                    float sysErr,  // Relative systematic error or NAN
     38                                                    float skyErr  // Relative systematic error or NAN
    3739    );
    3840
    … …  
    9597                                                float spacing, ///< Rough spacing for stamps
    9698                                                float sysErr,  ///< Relative systematic error in images
     99                                                float skyErr,  ///< Relative systematic error in images
    97100                                                pmSubtractionMode mode ///< Mode for subtraction
    98101    );
    … …  
    108111                                               float spacing, ///< Rough spacing for stamps
    109112                                               float sysErr,  ///< Systematic error in images
     113                                               float skyErr,  ///< Systematic error in images
    110114                                               pmSubtractionMode mode ///< Mode for subtraction
    111115    );
    … …  
    121125    float spacing,                      ///< Rough spacing for stamps
    122126    float sysErr,                       ///< Systematic error in images
     127    float skyErr,                       ///< Systematic error in images
    123128    pmSubtractionMode mode              ///< Mode for subtraction
    124129    );
    … …  
    134139    float spacing,                      ///< Rough spacing for stamps
    135140    float sysErr,                       ///< Systematic error in images
     141    float skyErr,                       ///< Systematic error in images
    136142    pmSubtractionMode mode              ///< Mode for subtraction
    137143    );

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c

    r26502 r26505  
    535535}
    536536
    537 bool pmSubtractionVisualShowFit() {
     537bool pmSubtractionVisualShowFit(double norm) {
     538
     539    // XXX a dumb test : dump the residual image and exit
     540    {
     541        psMetadata *header = psMetadataAlloc();
     542        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
     543        psFits *fits = psFitsOpen("resid.fits", "w");
     544        psFitsWriteImage(fits, header, residualImage, 0, NULL);
     545        psFitsClose(fits);
     546        exit (0);
     547    }
    538548
    539549    // if (!pmVisualIsVisual()) return true;

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.h

    r26472 r26505  
    99bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps);
    1010bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index);
    11 bool pmSubtractionVisualShowFit();
     11bool pmSubtractionVisualShowFit(double norm);
    1212bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
    1313bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);
Note: See TracChangeset for help on using the changeset viewer.