IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30383


Ignore:
Timestamp:
Jan 26, 2011, 5:21:42 PM (15 years ago)
Author:
eugene
Message:

report the diff chisq as well as the residual chisq; truncate SVD at dynamic range of 1e10; fix the scaling process

Location:
branches/eam_branches/ipp-20101205/psModules/src/imcombine
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c

    r30333 r30383  
    10721072        // SINGLE solution
    10731073# if (1)
    1074         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1074        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
    10751075        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
    10761076# endif
    … …  
    11761176        // DUAL solution
    11771177# if (1)
    1178         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
     1178        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
    11791179        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
    11801180# endif
    … …  
    13061306
    13071307// given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq
    1308 bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) {
     1308bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window) {
    13091309
    13101310# ifndef USE_WEIGHT
    … …  
    13161316
    13171317    int npix = 0;
    1318     float chisq = 0;
     1318    float chisqR = 0;
     1319    float chisqD = 0;
    13191320
    13201321    // get the chisq
    13211322    for (int y = residual->yMin; y <= residual->yMax; y++) {
    13221323        for (int x = residual->xMin; x <= residual->xMax; x++) {
    1323             float value = PS_SQR(residual->kernel[y][x]);
     1324            float valueR = PS_SQR(residual->kernel[y][x]);
    13241325            if (weight) {
    1325                 value *= weight->kernel[y][x];
     1326                valueR *= weight->kernel[y][x];
    13261327            }
    1327             // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
    1328             if (false && window) {
    1329                 value *= window->kernel[y][x];
     1328            // XXX NOTE: do NOT apply the window to the chisq portions of the calculation (that would bias the chisq)
     1329            chisqR += valueR;
     1330
     1331            float valueD = PS_SQR(difference->kernel[y][x]);
     1332            if (weight) {
     1333                valueD *= weight->kernel[y][x];
    13301334            }
    1331             chisq += value;
     1335            chisqD += valueD;
    13321336            npix ++;
    13331337        }
    13341338    }
    1335     psVectorAppend(chisqVector, chisq / npix);
     1339    psVectorAppend(chisqRVector, chisqR / npix);
     1340    psVectorAppend(chisqDVector, chisqD / npix);
    13361341
    13371342    float value1 = 0;
    … …  
    14301435    int Nelem = fluxesVector->n - 1;
    14311436    bool valid = true;
    1432     valid &= isfinite(chisqVector->data.F32[Nelem]);
     1437    valid &= isfinite(chisqRVector->data.F32[Nelem]);
    14331438    valid &= isfinite(fluxesVector->data.F32[Nelem]);
    14341439    valid &= isfinite(momentVector->data.F32[Nelem]);
    … …  
    14531458    // XXX need to save these somewhere
    14541459    psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    1455     psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1460    psVector *chisqD = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     1461    psVector *chisqR = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    14561462    psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
    14571463    psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK);
    … …  
    14641470    // storage for the image (convolved2 is not used in SINGLE mode)
    14651471    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     1472    psKernel *difference = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
    14661473    psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
    14671474    psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
    … …  
    14741481            psVectorAppend(moments, NAN);
    14751482            psVectorAppend(fluxes, NAN);
    1476             psVectorAppend(chisq, NAN);
     1483            psVectorAppend(chisqD, NAN);
     1484            psVectorAppend(chisqR, NAN);
    14771485            psVectorAppend(stampMask, 0x01);
    14781486            continue;
    … …  
    14871495        // Calculate residuals
    14881496        psImageInit(residual->image, 0.0);
     1497        psImageInit(difference->image, 0.0);
    14891498
    14901499        psKernel *weight = NULL;
    … …  
    15371546            }
    15381547
    1539             // generate the residual image
     1548            // Generate the difference, residual, and convolved source images.  Note the we
     1549            // accumulate the convolution of (A-B), so we need to replace it to generate the
     1550            // images of the convolved source image.
    15401551            for (int y = - footprint; y <= footprint; y++) {
    15411552                for (int x = - footprint; x <= footprint; x++) {
     1553                    difference->kernel[y][x] = target->kernel[y][x] - source->kernel[y][x] * norm - background;
     1554                    residual->kernel[y][x] = difference->kernel[y][x] - convolved1->kernel[y][x];
    15421555                    convolved1->kernel[y][x] += source->kernel[y][x] * norm;
    1543                     residual->kernel[y][x] = target->kernel[y][x] - convolved1->kernel[y][x] - background;
    15441556                }
    15451557            }
     1558
    15461559            // XXX if we want to have a weight and window, we'll need to pass through to here
    1547             pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, NULL, residual, weight, window);
     1560            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, NULL, difference, residual, weight, window);
    15481561
    15491562        } else {
    … …  
    15741587            }
    15751588
     1589            // Generate the difference, residual, and convolved source images.  Note the we
     1590            // accumulate the convolutions of (A-B), so we need to replace (A or B) to generate
     1591            // the images of the convolved source images.
    15761592            for (int y = - footprint; y <= footprint; y++) {
    15771593                for (int x = - footprint; x <= footprint; x++) {
     1594                    difference->kernel[y][x] = image2->kernel[y][x] - image1->kernel[y][x] * norm - background;
     1595                    residual->kernel[y][x] = difference->kernel[y][x] + convolved2->kernel[y][x] - convolved1->kernel[y][x];
    15781596                    convolved1->kernel[y][x] += image1->kernel[y][x] * norm;
    15791597                    convolved2->kernel[y][x] += image2->kernel[y][x];
    1580                     residual->kernel[y][x] = convolved2->kernel[y][x] - convolved1->kernel[y][x] - background;
    15811598                }
    15821599            }
    … …  
    15891606            }
    15901607
    1591             pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, convolved2, residual, weight, window);
     1608            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, convolved2, difference, residual, weight, window);
    15921609        }
    15931610    }
    … …  
    15951612    // find the mean chisq and mean moment
    15961613    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    1597     psVectorStats (stats, chisq, NULL, stampMask, 0xff);
    1598     float chisqValue = stats->sampleMean;
     1614    psVectorStats (stats, chisqD, NULL, stampMask, 0xff);
     1615    float chisqDValue = stats->sampleMean;
     1616
     1617    psStatsInit(stats);
     1618    psVectorStats (stats, chisqR, NULL, stampMask, 0xff);
     1619    float chisqRValue = stats->sampleMean;
    15991620
    16001621    psStatsInit(stats);
    … …  
    16371658    // penalized by increasing the score somewhat.  the 0.01 value is not well-chosen.
    16381659    float orderFactor = 0.01 * kernels->spatialOrder;
    1639     float score = 2.0 * chisqValue / (sumKernel1 + sumKernel2) + orderFactor;
    1640     psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqValue, momentValue, sumKernel1, sumKernel2, score);
     1660    float score = 2.0 * chisqRValue / (sumKernel1 + sumKernel2) + orderFactor;
     1661    psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, chisqD: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqRValue, chisqDValue, momentValue, sumKernel1, sumKernel2, score);
    16411662
    16421663    // save this result if it is the first or the best (skip if bestMatch is NULL)
    … …  
    16631684            match->nGood        = nGood;
    16641685            match->fluxes       = psMemIncrRefCounter(fluxes);
    1665             match->chisq        = psMemIncrRefCounter(chisq);
     1686            match->chisq        = psMemIncrRefCounter(chisqR);
    16661687            match->moments      = psMemIncrRefCounter(moments);
    16671688            match->stampMask    = psMemIncrRefCounter(stampMask);
    … …  
    16691690    }
    16701691
    1671     pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments);
     1692    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisqR, moments);
    16721693
    16731694    psFree(stats);
    1674     psFree(chisq);
     1695    psFree(chisqR);
     1696    psFree(chisqD);
    16751697    psFree(fluxes);
    16761698    psFree(moments);

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h

    r30322 r30383  
    8282bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window);
    8383
    84 bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window);
     84bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window);
    8585
    8686bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

    r30362 r30383  
    476476                        const psArray *sources, const char *stampsName,
    477477                        pmSubtractionKernelsType type, int size, int spatialOrder,
    478                         const psVector *isisWidths, const psVector *isisOrders,
     478                        psVector *isisWidths, const psVector *isisOrders,
    479479                        int inner, int ringsOrder, int binning, float penalty,
    480480                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
    … …  
    695695            // check on the kernel scaling -- if the kron-based radial moments are very different, adjust to match them
    696696            {
    697                 float fwhm1;
    698                 float fwhm2;
    699                 pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
    700                 psAssert(isfinite(fwhm1), "fwhm 1 not set");
    701                 psAssert(isfinite(fwhm2), "fwhm 2 not set");
     697                // float fwhm1;
     698                // float fwhm2;
     699                // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     700                // psAssert(isfinite(fwhm1), "fwhm 1 not set");
     701                // psAssert(isfinite(fwhm2), "fwhm 2 not set");
    702702
    703703                // XXX this is BAD: depends on the relationship below:
    … …  
    706706                float radMoment1 = stamps->normWindow1 / 2.75;
    707707                float radMoment2 = stamps->normWindow2 / 2.75;
    708 
    709                 float maxFWHM = PS_MAX(fwhm1, fwhm2);
    710                 float maxRadial = PS_MAX(radMoment1, radMoment2);
     708                pmSubtractionParamsScale(NULL, NULL, isisWidths, radMoment1, radMoment2);
     709
     710                // float maxFWHM = PS_MAX(fwhm1, fwhm2);
     711                // float maxRadial = PS_MAX(radMoment1, radMoment2);
    711712               
    712713                // if (fabs(2.0*(maxFWHM - maxRadial)/(maxFWHM + maxRadial)) > 0.25) {
    713                 if (1) {
    714 
    715                     float scale = maxRadial / maxFWHM;
    716                     psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale);
    717                    
    718                     for (int i = 0; i < isisWidths->n; i++) {
    719                         isisWidths->data.F32[i] *= scale;
    720                     }
    721                 }
     714                // if (1) {
     715                //
     716                //     float scale = maxRadial / maxFWHM;
     717                //     psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale);
     718                //    
     719                //     for (int i = 0; i < isisWidths->n; i++) {
     720                //      isisWidths->data.F32[i] *= scale;
     721                //     }
     722                // }
    722723            }
    723724
    … …  
    12521253}
    12531254
    1254 
    1255 bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
    1256                               float scaleRef, float scaleMin, float scaleMax)
     1255static float scaleRefOption = NAN;
     1256static float scaleMinOption = NAN;
     1257static float scaleMaxOption = NAN;
     1258static bool  scaleOption = false;
     1259
     1260bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax) {
     1261
     1262    if (scale) {
     1263        PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
     1264        PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
     1265        PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
     1266        PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
     1267    }
     1268
     1269    scaleRefOption = scaleRef;
     1270    scaleMinOption = scaleMin;
     1271    scaleMaxOption = scaleMax;
     1272    scaleOption = scale;
     1273   
     1274    return true;
     1275}
     1276
     1277bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2)
    12571278{
    1258     PS_ASSERT_PTR_NON_NULL(kernelSize, false);
    1259     PS_ASSERT_PTR_NON_NULL(stampSize, false);
     1279    // PS_ASSERT_PTR_NON_NULL(kernelSize, false);
     1280    // PS_ASSERT_PTR_NON_NULL(stampSize, false);
    12601281    PS_ASSERT_VECTOR_NON_NULL(widths, false);
    12611282    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
    1262     PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
    1263     PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
    1264     PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
    1265     PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
    1266 
    1267     float fwhm1;
    1268     float fwhm2;
    1269 
    1270     pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
    1271     psAssert(isfinite(fwhm1), "fwhm 1 not set");
    1272     psAssert(isfinite(fwhm2), "fwhm 2 not set");
     1283
     1284    if (!scaleOption) return true;
     1285
     1286    // pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     1287    // psAssert(isfinite(fwhm1), "fwhm 1 not set");
     1288    // psAssert(isfinite(fwhm2), "fwhm 2 not set");
    12731289   
    12741290    // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
    1275     float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
    1276 
    1277     if (isfinite(scaleMin) && scale < scaleMin) {
    1278         scale = scaleMin;
    1279     }
    1280     if (isfinite(scaleMax) && scale > scaleMax) {
    1281         scale = scaleMax;
     1291    float scale = PS_MAX(fwhm1, fwhm2) / scaleRefOption;      // Scaling factor
     1292
     1293    if (isfinite(scaleMinOption) && scale < scaleMinOption) {
     1294        scale = scaleMinOption;
     1295    }
     1296    if (isfinite(scaleMaxOption) && scale > scaleMaxOption) {
     1297        scale = scaleMaxOption;
    12821298    }
    12831299
    … …  
    12851301        widths->data.F32[i] *= scale;
    12861302    }
    1287     *kernelSize = *kernelSize * scale + 0.5;
    1288     *stampSize = *stampSize * scale + 0.5;
    1289 
    1290     psLogMsg("psModules.imcombine", PS_LOG_INFO,
    1291              "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     1303    if (kernelSize) {
     1304        *kernelSize = *kernelSize * scale + 0.5;
     1305    }
     1306    if (stampSize) {
     1307        *stampSize = *stampSize * scale + 0.5;
     1308    }
     1309
     1310    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Scaling kernel parameters by %f", scale);
     1311    if (kernelSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified kernel size %d", *kernelSize);
     1312    if (stampSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified stamp size %d", *stampSize);
    12921313
    12931314    return true;

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h

    r30288 r30383  
    2727                        int size,       ///< Kernel half-size
    2828                        int order,      ///< Spatial polynomial order
    29                         const psVector *widths, ///< ISIS Gaussian widths
     29                        psVector *widths, ///< ISIS Gaussian widths
    3030                        const psVector *orders, ///< ISIS Polynomial orders
    3131                        int inner,      ///< Inner radius for various kernel types
    … …  
    101101
    102102/// Scale subtraction parameters according to the FWHMs of the inputs
    103 bool pmSubtractionParamsScale(
    104     int *kernelSize,                    ///< Half-size of the kernel
    105     int *stampSize,                     ///< Half-size of the stamp (footprint)
    106     psVector *widths,                   ///< ISIS widths
    107     float scaleRef,                     ///< Reference width for scaling
    108     float scaleMin,                     ///< Minimum scaling ratio, or NAN
    109     float scaleMax                      ///< Maximum scaling ratio, or NAN
    110     );
     103// bool pmSubtractionParamsScale(
     104//     int *kernelSize,                    ///< Half-size of the kernel
     105//     int *stampSize,                     ///< Half-size of the stamp (footprint)
     106//     psVector *widths,                   ///< ISIS widths
     107//     float scaleRef,                     ///< Reference width for scaling
     108//     float scaleMin,                     ///< Minimum scaling ratio, or NAN
     109//     float scaleMax                      ///< Maximum scaling ratio, or NAN
     110//     );
     111
     112bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2);
     113
     114bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax);
    111115
    112116bool pmSubtractionMatchAttempt(
Note: See TracChangeset for help on using the changeset viewer.