IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42824


Ignore:
Timestamp:
May 8, 2025, 4:46:11 PM (12 months ago)
Author:
eugene
Message:

add IRLS versions of polynomial fits

Location:
trunk/psLib/src/math
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r41896 r42824  
    6969if ((ORIG != NULL) && (ORIG->type.type != PS_TYPE_F64)) { psFree(TEMP); }
    7070
     71psVector *psVector_GetModifiedErrors_Caucy
     72( const psVector *f,
     73  const psVector *fEval,
     74  const psVector *fErr,
     75  const psVector *mask,
     76  psVectorMaskType maskValue);
     77
    7178/*****************************************************************************/
    7279/* TYPE DEFINITIONS                                                          */
    … …  
    9097returned as a psVector sums.
    9198*****************************************************************************/
    92 static psVector *BuildSums1D(
    93     psVector* sums,
    94     psF64 x,
    95     psS32 nTerm)
     99static psVector *BuildSums1D
     100( psVector* sums,
     101  psF64 x,
     102  psS32 nTerm)
    96103{
    97104    psS32 nSum = 0;
    … …  
    710717    switch (poly->type) {
    711718    case PS_POLYNOMIAL_ORD:
    712       if ((f64->n < 10000)&&(poly->nX > 1)) {
     719      // if x is NULL, the domain is the index of the vector, in which case we do not
     720      // need to rescale (range is guaranteed to be < 10000)
     721      if ((x != NULL) && (f64->n < 10000) && (poly->nX > 1)) {
    713722        scale = true;
    714723
    … …  
    752761        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
    753762        psFree(z64); // Done with this.
    754       }
    755       else {
     763      } else {
    756764        result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64);
    757765      }
    … …  
    10301038}
    10311039
     1040// These should probably be tunable:
     1041# define FIT_TOLERANCE 1e-4
     1042# define FLT_TOLERANCE 1e-6
     1043# define WEIGHT_THRESHOLD 0.3
     1044
     1045// This function accepts F32 and F64 input vectors.
     1046bool psVectorIRLSFitPolynomial1D(
     1047    psPolynomial1D *poly,
     1048    psStats *stats,
     1049    const psVector *mask,
     1050    psVectorMaskType maskValue,
     1051    const psVector *f,
     1052    const psVector *fErr,
     1053    const psVector *xIn)
     1054{
     1055    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
     1056    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
     1057    PS_ASSERT_POLY_NON_NULL(poly, false);
     1058    PS_ASSERT_VECTOR_NON_NULL(f, false);
     1059    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
     1060    if (mask != NULL) {
     1061        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
     1062        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
     1063    }
     1064    if (fErr != NULL) {
     1065        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
     1066        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
     1067    }
     1068    if (xIn != NULL) {
     1069        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
     1070        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
     1071    }
     1072
     1073    // Internal pointers for possibly NULL vectors.
     1074    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     1075
     1076    // initial fit with nominal errors
     1077    if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x)) {
     1078        psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     1079        return false;
     1080    }
     1081
     1082    // use polyOld to save the last fit
     1083    psPolynomial1D *polyOld = NULL;
     1084
     1085    // use clipIter as max number of iterations
     1086    bool converged = false;
     1087    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
     1088        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial1D()\n", N);
     1089
     1090        // evaluate the fit at the input positions
     1091        psVector *fEval = psPolynomial1DEvalVector (poly, x);
     1092
     1093        // calculate modified errors based on the deviation from the fit
     1094        psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
     1095        psFree (fEval);
     1096
     1097        // save the last fit (recycle the structure once allocated)
     1098        polyOld = psPolynomial1DCopy (polyOld, poly);
     1099
     1100        // calculate a new fit with modified errors:
     1101        if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, modErr, x)) {
     1102            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     1103            psFree(x);
     1104            psFree(modErr);
     1105            return false;
     1106        }
     1107
     1108        // has the solution converged?
     1109        converged = true;
     1110        for (int ix = 0; ix <= poly->nX; ix++) {
     1111          if ((fabs(poly->coeff[ix] - polyOld->coeff[ix]) > FIT_TOLERANCE * fabs(poly->coeff[ix])) &&
     1112              (fabs(poly->coeff[ix] - polyOld->coeff[ix]) > FLT_TOLERANCE))
     1113            converged = false;
     1114        }
     1115
     1116# if (0)       
     1117        // XXX test:
     1118        FILE *ftest = fopen ("irls.wt.dat", "w");
     1119        for (int i = 0; i < modErr->n; i++) {
     1120            if (modErr->type.type == PS_TYPE_F64) {
     1121                fprintf (ftest, "%d %f\n", i, modErr->data.F64[i]);
     1122            } else {
     1123                fprintf (ftest, "%d %f\n", i, modErr->data.F32[i]);
     1124            }
     1125        }
     1126        fclose (ftest);
     1127# endif
     1128        psFree (modErr);
     1129    }
     1130
     1131    // Free local temporary variables
     1132    psFree(x);
     1133    psFree(polyOld);
     1134
     1135    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     1136    return true;
     1137}
    10321138
    10331139/******************************************************************************
    … …  
    16131719}
    16141720
     1721// This function accepts F32 and F64 input vectors.
     1722bool psVectorIRLSFitPolynomial2D(
     1723    psPolynomial2D *poly,
     1724    psStats *stats,
     1725    const psVector *mask,
     1726    psVectorMaskType maskValue,
     1727    const psVector *f,
     1728    const psVector *fErr,
     1729    const psVector *xIn,
     1730    const psVector *yIn)
     1731{
     1732    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
     1733
     1734    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
     1735    PS_ASSERT_POLY_NON_NULL(poly, false);
     1736    PS_ASSERT_VECTOR_NON_NULL(f, false);
     1737    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
     1738    if (mask != NULL) {
     1739        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
     1740        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
     1741    }
     1742    if (fErr != NULL) {
     1743        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
     1744        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
     1745    }
     1746    if (xIn != NULL) {
     1747        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
     1748        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
     1749    }
     1750    if (yIn != NULL) {
     1751        PS_ASSERT_VECTORS_SIZE_EQUAL(yIn, f, false);
     1752        PS_ASSERT_VECTOR_TYPE(yIn, f->type.type, false);
     1753    }
     1754
     1755    // Internal pointers for possibly NULL vectors. 
     1756    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     1757    psVector *y = (yIn != NULL) ? psMemIncrRefCounter((psVector *) yIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     1758
     1759    // initial fit with nominal errors
     1760    if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y)) {
     1761        psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     1762        psFree(x);
     1763        psFree(y);
     1764        return false;
     1765    }
     1766
     1767    // use polyOld to save the last fit
     1768    psPolynomial2D *polyOld = NULL;
     1769
     1770    // use clipIter as max number of iterations
     1771    bool converged = false;
     1772    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
     1773        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial2D()\n", N);
     1774
     1775        // evaluate the fit at the input positions
     1776        psVector *fEval = psPolynomial2DEvalVector (poly, x, y);
     1777
     1778        // calculate modified errors based on the deviation from the fit
     1779        psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
     1780        psFree (fEval);
     1781
     1782        // save the last fit (recycle the structure once allocated)
     1783        polyOld = psPolynomial2DCopy (polyOld, poly);
     1784
     1785        // calculate a new fit with modified errors:
     1786        if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, modErr, x, y)) {
     1787            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     1788            psFree(x);
     1789            psFree(y);
     1790            psFree(modErr);
     1791            return false;
     1792        }
     1793
     1794        // has the solution converged?
     1795        converged = true;
     1796        for (int ix = 0; ix <= poly->nX; ix++) {
     1797            for (int iy = 0; iy <= poly->nY; iy++) {
     1798                if ((fabs(poly->coeff[ix][iy] - polyOld->coeff[ix][iy]) > FIT_TOLERANCE * fabs(poly->coeff[ix][iy])) &&
     1799                    (fabs(poly->coeff[ix][iy] - polyOld->coeff[ix][iy]) > FLT_TOLERANCE))
     1800                    converged = false;
     1801            }
     1802        }
     1803        psFree (modErr);
     1804    }
     1805
     1806    // Free local temporary variables
     1807    psFree(x);
     1808    psFree(y);
     1809    psFree(polyOld);
     1810
     1811    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     1812    return true;
     1813}
    16151814
    16161815/******************************************************************************
    … …  
    20162215    // Free local temporary variables
    20172216    psFree(resid);
     2217
     2218    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     2219    return true;
     2220}
     2221
     2222// This function accepts F32 and F64 input vectors.
     2223bool psVectorIRLSFitPolynomial3D(
     2224    psPolynomial3D *poly,
     2225    psStats *stats,
     2226    const psVector *mask,
     2227    psVectorMaskType maskValue,
     2228    const psVector *f,
     2229    const psVector *fErr,
     2230    const psVector *xIn,
     2231    const psVector *yIn,
     2232    const psVector *zIn)
     2233{
     2234    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
     2235
     2236    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
     2237    PS_ASSERT_POLY_NON_NULL(poly, false);
     2238    PS_ASSERT_VECTOR_NON_NULL(f, false);
     2239    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
     2240    if (mask != NULL) {
     2241        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
     2242        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
     2243    }
     2244    if (fErr != NULL) {
     2245        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
     2246        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
     2247    }
     2248    if (xIn != NULL) {
     2249        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
     2250        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
     2251    }
     2252    if (yIn != NULL) {
     2253        PS_ASSERT_VECTORS_SIZE_EQUAL(yIn, f, false);
     2254        PS_ASSERT_VECTOR_TYPE(yIn, f->type.type, false);
     2255    }
     2256    if (zIn != NULL) {
     2257        PS_ASSERT_VECTORS_SIZE_EQUAL(zIn, f, false);
     2258        PS_ASSERT_VECTOR_TYPE(zIn, f->type.type, false);
     2259    }
     2260
     2261    // Internal pointers for possibly NULL vectors. 
     2262    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2263    psVector *y = (yIn != NULL) ? psMemIncrRefCounter((psVector *) yIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2264    psVector *z = (zIn != NULL) ? psMemIncrRefCounter((psVector *) zIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2265
     2266    // initial fit with nominal errors
     2267    if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z)) {
     2268        psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     2269        psFree(x);
     2270        psFree(y);
     2271        psFree(z);
     2272        return false;
     2273    }
     2274
     2275    // use polyOld to save the last fit
     2276    psPolynomial3D *polyOld = NULL;
     2277
     2278    // use clipIter as max number of iterations
     2279    bool converged = false;
     2280    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
     2281        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial3D()\n", N);
     2282
     2283        // evaluate the fit at the input positions
     2284        psVector *fEval = psPolynomial3DEvalVector (poly, x, y, z);
     2285
     2286        // calculate modified errors based on the deviation from the fit
     2287        psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
     2288        psFree (fEval);
     2289
     2290        // save the last fit (recycle the structure once allocated)
     2291        polyOld = psPolynomial3DCopy (polyOld, poly);
     2292
     2293        // calculate a new fit with modified errors:
     2294        if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, modErr, x, y, z)) {
     2295            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     2296            psFree(x);
     2297            psFree(y);
     2298            psFree(z);
     2299            psFree(modErr);
     2300            return false;
     2301        }
     2302
     2303        // has the solution converged?
     2304        converged = true;
     2305        for (int ix = 0; ix <= poly->nX; ix++) {
     2306            for (int iy = 0; iy <= poly->nY; iy++) {
     2307                for (int iz = 0; iz <= poly->nZ; iz++) {
     2308                    if ((fabs(poly->coeff[ix][iy][iz] - polyOld->coeff[ix][iy][iz]) > FIT_TOLERANCE * fabs(poly->coeff[ix][iy][iz])) &&
     2309                        (fabs(poly->coeff[ix][iy][iz] - polyOld->coeff[ix][iy][iz]) > FLT_TOLERANCE))
     2310                        converged = false;
     2311                }
     2312            }
     2313        }
     2314        psFree (modErr);
     2315    }
     2316
     2317    // Free local temporary variables
     2318    psFree(x);
     2319    psFree(y);
     2320    psFree(z);
     2321    psFree(polyOld);
    20182322
    20192323    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
    … …  
    24552759    return true;
    24562760}
     2761
     2762// This function accepts F32 and F64 input vectors.
     2763bool psVectorIRLSFitPolynomial4D(
     2764    psPolynomial4D *poly,
     2765    psStats *stats,
     2766    const psVector *mask,
     2767    psVectorMaskType maskValue,
     2768    const psVector *f,
     2769    const psVector *fErr,
     2770    const psVector *xIn,
     2771    const psVector *yIn,
     2772    const psVector *zIn,
     2773    const psVector *tIn)
     2774{
     2775    psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__);
     2776
     2777    PS_ASSERT (poly->type == PS_POLYNOMIAL_ORD, false); // XXX for now, only allow ORD
     2778    PS_ASSERT_POLY_NON_NULL(poly, false);
     2779    PS_ASSERT_VECTOR_NON_NULL(f, false);
     2780    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false);
     2781    if (mask != NULL) {
     2782        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false);
     2783        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false);
     2784    }
     2785    if (fErr != NULL) {
     2786        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false);
     2787        PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false);
     2788    }
     2789    if (xIn != NULL) {
     2790        PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false);
     2791        PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false);
     2792    }
     2793    if (yIn != NULL) {
     2794        PS_ASSERT_VECTORS_SIZE_EQUAL(yIn, f, false);
     2795        PS_ASSERT_VECTOR_TYPE(yIn, f->type.type, false);
     2796    }
     2797    if (zIn != NULL) {
     2798        PS_ASSERT_VECTORS_SIZE_EQUAL(zIn, f, false);
     2799        PS_ASSERT_VECTOR_TYPE(zIn, f->type.type, false);
     2800    }
     2801    if (tIn != NULL) {
     2802        PS_ASSERT_VECTORS_SIZE_EQUAL(tIn, f, false);
     2803        PS_ASSERT_VECTOR_TYPE(tIn, f->type.type, false);
     2804    }
     2805
     2806    // Internal pointers for possibly NULL vectors. 
     2807    psVector *x = (xIn != NULL) ? psMemIncrRefCounter((psVector *) xIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2808    psVector *y = (yIn != NULL) ? psMemIncrRefCounter((psVector *) yIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2809    psVector *z = (zIn != NULL) ? psMemIncrRefCounter((psVector *) zIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2810    psVector *t = (tIn != NULL) ? psMemIncrRefCounter((psVector *) tIn) : psVectorCreate(NULL, 0, f->n, 1, f->type.type);
     2811
     2812    // initial fit with nominal errors
     2813    if (!psVectorFitPolynomial4D(poly, mask, maskValue, f, fErr, x, y, z, t)) {
     2814        psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     2815        psFree(x);
     2816        psFree(y);
     2817        psFree(z);
     2818        psFree(t);
     2819        return false;
     2820    }
     2821
     2822    // use polyOld to save the last fit
     2823    psPolynomial4D *polyOld = NULL;
     2824
     2825    // use clipIter as max number of iterations
     2826    bool converged = false;
     2827    for (psS32 N = 0; !converged && (N < stats->clipIter); N++) {
     2828        psTrace("psLib.math", 6, "Loop iteration %d.  Calling psVectorFitPolynomial4D()\n", N);
     2829
     2830        // evaluate the fit at the input positions
     2831        psVector *fEval = psPolynomial4DEvalVector (poly, x, y, z, t);
     2832
     2833        // calculate modified errors based on the deviation from the fit
     2834        psVector *modErr = psVector_GetModifiedErrors_Caucy (f, fEval, fErr, mask, maskValue);
     2835        psFree (fEval);
     2836
     2837        // save the last fit (recycle the structure once allocated)
     2838        polyOld = psPolynomial4DCopy (polyOld, poly);
     2839
     2840        // calculate a new fit with modified errors:
     2841        if (!psVectorFitPolynomial4D(poly, mask, maskValue, f, modErr, x, y, z, t)) {
     2842            psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial.  Returning false.\n");
     2843            psFree(x);
     2844            psFree(y);
     2845            psFree(z);
     2846            psFree(t);
     2847            psFree(modErr);
     2848            return false;
     2849        }
     2850
     2851        // has the solution converged?
     2852        converged = true;
     2853        for (int ix = 0; ix <= poly->nX; ix++) {
     2854            for (int iy = 0; iy <= poly->nY; iy++) {
     2855                for (int iz = 0; iz <= poly->nZ; iz++) {
     2856                    for (int it = 0; it <= poly->nT; it++) {
     2857                        if ((fabs(poly->coeff[ix][iy][iz][it] - polyOld->coeff[ix][iy][iz][it]) > FIT_TOLERANCE * fabs(poly->coeff[ix][iy][iz][it])) &&
     2858                            (fabs(poly->coeff[ix][iy][iz][it] - polyOld->coeff[ix][iy][iz][it]) > FLT_TOLERANCE))
     2859                            converged = false;
     2860                    }
     2861                }
     2862            }
     2863        }
     2864        psFree (modErr);
     2865    }
     2866
     2867    // Free local temporary variables
     2868    psFree(x);
     2869    psFree(y);
     2870    psFree(z);
     2871    psFree(t);
     2872    psFree(polyOld);
     2873
     2874    psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
     2875    return true;
     2876}
     2877
     2878// ######################## utilities ###################
     2879
     2880// Used by IRLS fitting
     2881// This function assumes the input vectors (f, fEval, fErr) all have the same type
     2882// This requirement is already enforced in the calling function (psVectorIRLSFitPolynomial1D)
     2883psVector *psVector_GetModifiedErrors_Caucy (
     2884    const psVector *f,
     2885    const psVector *fEval,
     2886    const psVector *fErr,
     2887    const psVector *mask,
     2888    psVectorMaskType maskValue) {
     2889
     2890  psVector *mErr = psVectorAlloc (f->n, f->type.type);
     2891
     2892  switch (f->type.type) {
     2893    case PS_TYPE_F64:
     2894      {
     2895          psF64 *fPtr     = f->data.F64;
     2896          psF64 *fErrPtr  = fErr->data.F64;
     2897          psF64 *mErrPtr  = mErr->data.F64;
     2898          psF64 *fEvalPtr = fEval->data.F64;
     2899          for (int i = 0; i < f->n; i++, fPtr++, fEvalPtr++, fErrPtr++, mErrPtr++) {
     2900              if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValue)) continue;
     2901              double dF = (*fPtr - *fEvalPtr) / 2.385;
     2902              double dV = PS_SQR(*fErrPtr) + PS_SQR(dF);
     2903              *mErrPtr = sqrt(dV);
     2904          }
     2905      }
     2906      break;
     2907    case PS_TYPE_F32:
     2908      {
     2909          psF32 *fPtr     = f->data.F32;
     2910          psF32 *fErrPtr  = fErr->data.F32;
     2911          psF32 *mErrPtr  = mErr->data.F32;
     2912          psF32 *fEvalPtr = fEval->data.F32;
     2913          for (int i = 0; i < f->n; i++, fPtr++, fEvalPtr++, fErrPtr++, mErrPtr++) {
     2914              if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValue)) continue;
     2915              double dF = (*fPtr - *fEvalPtr) / 2.385;
     2916              double dV = PS_SQR(*fErrPtr) + PS_SQR(dF);
     2917              *mErrPtr = sqrt(dV);
     2918          }
     2919      }
     2920      break;
     2921    default:
     2922        psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n");
     2923        return NULL;
     2924  }
     2925  return mErr;
     2926}
     2927

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

    r21183 r42824  
    134134);
    135135
     136bool psVectorIRLSFitPolynomial1D(
     137    psPolynomial1D *poly,
     138    psStats *stats,
     139    const psVector *mask,
     140    psVectorMaskType maskValue,
     141    const psVector *f,
     142    const psVector *fErr,
     143    const psVector *x
     144);
     145
     146bool psVectorIRLSFitPolynomial2D(
     147    psPolynomial2D *poly,
     148    psStats *stats,
     149    const psVector *mask,
     150    psVectorMaskType maskValue,
     151    const psVector *f,
     152    const psVector *fErr,
     153    const psVector *x,
     154    const psVector *y
     155);
     156
     157bool psVectorIRLSFitPolynomial3D(
     158    psPolynomial3D *poly,
     159    psStats *stats,
     160    const psVector *mask,
     161    psVectorMaskType maskValue,
     162    const psVector *f,
     163    const psVector *fErr,
     164    const psVector *x,
     165    const psVector *y,
     166    const psVector *z
     167);
     168
     169bool psVectorIRLSFitPolynomial4D(
     170    psPolynomial4D *poly,
     171    psStats *stats,
     172    const psVector *mask,
     173    psVectorMaskType maskValue,
     174    const psVector *f,
     175    const psVector *fErr,
     176    const psVector *x,
     177    const psVector *y,
     178    const psVector *z,
     179    const psVector *t
     180);
     181
    136182/// @}
    137183#endif // #ifndef PS_MINIMIZE_POLYFIT_H
Note: See TracChangeset for help on using the changeset viewer.