IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 18, 2006, 8:37:38 AM (20 years ago)
Author:
magnier
Message:

allowing for ROBUST stats in clipping, better trace levels

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psMinimizePolyFit.c

    r9540 r9625  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2006-10-13 22:04:58 $
     12 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-10-18 18:37:38 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4242
    4343# define USE_GAUSS_JORDAN 1
     44# define USE_ROBUST_STATS_FOR_CLIPPING 1
    4445
    4546#define PS_VECTOR_GEN_CHEBY_INDEX(VEC, SIZE, TYPE) \
    … …  
    811812    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
    812813    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     814    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    813815    psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
    814816    psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
    … …  
    859861
    860862        stats  = psVectorStats(stats, resid, NULL, mask, maskValue);
    861         psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
    862         psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
    863         psF32 minClipValue = -minClipSigma*stats->sampleStdev;
    864         psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
     863        if (stats == NULL) {
     864            psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector.  Returning NULL.\n");
     865            psFree(resid)
     866            psFree(fit)
     867            return(NULL);
     868        }
     869        # if (USE_ROBUST_STATS_FOR_CLIPPING)
     870            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
     871        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
     872        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     873        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     874        psF32 clipMedian = stats->robustMedian;
     875        # else
     876
     877            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
     878        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
     879        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     880        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     881        psF32 clipMedian = stats->sampleMedian;
     882        # endif
    865883
    866884        // set mask if pts are not valid
    … …  
    872890            }
    873891
    874             if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
    875                     (resid->data.F64[i] - stats->sampleMedian < minClipValue)) {
     892            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
    876893                if (f->type.type == PS_TYPE_F64) {
    877894                    psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
    … …  
    12461263    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
    12471264    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     1265    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    12481266    psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
    12491267    psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
    … …  
    12521270        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n", N);
    12531271        psS32 Nkeep = 0;
    1254         if (psTraceGetLevel("psLib.math") >= 6) {
     1272        if (psTraceGetLevel("psLib.math") >= 7) {
    12551273            if (mask != NULL) {
    12561274                for (psS32 i = 0 ; i < mask->n ; i++) {
    1257                     psTrace("psLib.math", 6,  "mask[%d] is %d\n", i, mask->data.U8[i]);
     1275                    psTrace("psLib.math", 7,  "mask[%d] is %d\n", i, mask->data.U8[i]);
    12581276                }
    12591277            }
    … …  
    12821300        }
    12831301
    1284         if (psTraceGetLevel("psLib.math") >= 6) {
     1302        if (psTraceGetLevel("psLib.math") >= 7) {
    12851303            if (mask != NULL) {
    12861304                for (psS32 i = 0 ; i < mask->n ; i++) {
    12871305                    if (!((mask != NULL) && (mask->data.U8[i] & maskValue))) {
    1288                         psTrace("psLib.math", 6,  "(f, fit)[%d] is (%f, %f).  resid is (%f)\n",
     1306                        psTrace("psLib.math", 7,  "(f, fit)[%d] is (%f, %f).  resid is (%f)\n",
    12891307                                i, f->data.F32[i], fit->data.F32[i], resid->data.F64[i]);
    12901308                    }
    … …  
    13001318            return(NULL);
    13011319        }
    1302         psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
    1303         psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
    1304         psF32 minClipValue = -minClipSigma*stats->sampleStdev;
    1305         psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
     1320        # if (USE_ROBUST_STATS_FOR_CLIPPING)
     1321            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
     1322        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
     1323        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     1324        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     1325        psF32 clipMedian = stats->robustMedian;
     1326        # else
     1327
     1328            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
     1329        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
     1330        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     1331        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     1332        psF32 clipMedian = stats->sampleMedian;
     1333        # endif
    13061334
    13071335        // set mask if pts are not valid
    … …  
    13131341            }
    13141342
    1315             if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
    1316                     (resid->data.F64[i] - stats->sampleMedian < minClipValue)) {
     1343            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
    13171344                if (fit->type.type == PS_TYPE_F64) {
    13181345                    psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
    … …  
    13311358        }
    13321359
    1333         psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
     1360        psTrace("psLib.math", 4, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
    13341361        psFree(fit);
    13351362    }
    … …  
    17661793    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
    17671794    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     1795    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    17681796    psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
    17691797    psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
    … …  
    18201848        }
    18211849
    1822         psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
    1823         psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
    1824         psF32 minClipValue = -minClipSigma*stats->sampleStdev;
    1825         psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
     1850        # if (USE_ROBUST_STATS_FOR_CLIPPING)
     1851            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
     1852        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
     1853        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     1854        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     1855        psF32 clipMedian = stats->robustMedian;
     1856        # else
     1857
     1858            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
     1859        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
     1860        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     1861        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     1862        psF32 clipMedian = stats->sampleMedian;
     1863        # endif
    18261864
    18271865        // set mask if pts are not valid
    … …  
    18331871            }
    18341872
    1835             if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
    1836                     (resid->data.F64[i] - stats->sampleMedian < minClipValue))  {
     1873            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue))  {
    18371874                if (f->type.type == PS_TYPE_F64) {
    18381875                    psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
    … …  
    23342371    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
    23352372    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     2373    stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    23362374    psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter);
    23372375    psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
    … …  
    23882426            return(NULL);
    23892427        }
    2390         psTrace("psLib.math", 6, "Median is %f\n", stats->sampleMedian);
    2391         psTrace("psLib.math", 6, "Stdev is %f\n", stats->sampleStdev);
    2392         psF32 minClipValue = -minClipSigma*stats->sampleStdev;
    2393         psF32 maxClipValue = +maxClipSigma*stats->sampleStdev;
     2428        # if (USE_ROBUST_STATS_FOR_CLIPPING)
     2429            psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian);
     2430        psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev);
     2431        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     2432        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     2433        psF32 clipMedian = stats->robustMedian;
     2434        # else
     2435
     2436            psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian);
     2437        psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev);
     2438        psF32 minClipValue = -minClipSigma*stats->robustStdev;
     2439        psF32 maxClipValue = +maxClipSigma*stats->robustStdev;
     2440        psF32 clipMedian = stats->sampleMedian;
     2441        # endif
    23942442
    23952443        // set mask if pts are not valid
    … …  
    24012449            }
    24022450
    2403             if ((resid->data.F64[i] - stats->sampleMedian > maxClipValue) ||
    2404                     (resid->data.F64[i] - stats->sampleMedian < minClipValue)) {
     2451            if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) {
    24052452                if (f->type.type == PS_TYPE_F64) {
    24062453                    psTrace("psLib.math", 6, "Masking element %d (%f).  resid->data.F64[%d] is %f\n",
Note: See TracChangeset for help on using the changeset viewer.