src/ptafunc1.c
l_int32
ptaGetQuadraticLSF(PTA        *pta,
                   l_float32  *pa,
                   l_float32  *pb,
                   l_float32  *pc,
                   NUMA      **pnafit)
{
l_int32     n, i, ret;
l_float32   x, y, sx, sy, sx2, sx3, sx4, sxy, sx2y;
l_float32  *xa, *ya;
l_float32  *f[3];
l_float32   g[3];
NUMA       *nafit;

    PROCNAME("ptaGetQuadraticLSF");

    if (!pa && !pb && !pc && !pnafit)
        return ERROR_INT("no output requested", procName, 1);
    if (pa) *pa = 0.0;
    if (pb) *pb = 0.0;
    if (pc) *pc = 0.0;
    if (pnafit) *pnafit = NULL;
    if (!pta)
        return ERROR_INT("pta not defined", procName, 1);

    if ((n = ptaGetCount(pta)) < 3)
        return ERROR_INT("less than 3 pts not found", procName, 1);
    xa = pta->x;  /* not a copy */
    ya = pta->y;  /* not a copy */

    sx = sy = sx2 = sx3 = sx4 = sxy = sx2y = 0.;
    for (i = 0; i < n; i++) {
        x = xa[i];
        y = ya[i];
        sx += x;
        sy += y;
        sx2 += x * x;
        sx3 += x * x * x;
        sx4 += x * x * x * x;
        sxy += x * y;
        sx2y += x * x * y;
    }

    for (i = 0; i < 3; i++)
        f[i] = (l_float32 *)CALLOC(3, sizeof(l_float32));
    f[0][0] = sx4;
    f[0][1] = sx3;
    f[0][2] = sx2;
    f[1][0] = sx3;
    f[1][1] = sx2;
    f[1][2] = sx;
    f[2][0] = sx2;
    f[2][1] = sx;
    f[2][2] = n;
    g[0] = sx2y;
    g[1] = sxy;
    g[2] = sy;

        /* Solve for the unknowns, also putting f-inverse into f */
    ret = gaussjordan(f, g, 3);
    for (i = 0; i < 3; i++)
        FREE(f[i]);
    if (ret)
        return ERROR_INT("quadratic solution failed", procName, 1);

    if (pa) *pa = g[0];
    if (pb) *pb = g[1];
    if (pc) *pc = g[2];
    if (pnafit) {
        nafit = numaCreate(n);
        *pnafit = nafit;
        for (i = 0; i < n; i++) {
            x = xa[i];
            y = g[0] * x * x + g[1] * x + g[2];
            numaAddNumber(nafit, y);
        }
    }

    return 0;
}
/*!
 *  ptaGetQuadraticLSF()
 *
 *      Input:  pta
 *              &a  (<optional return> coeff a of LSF: y = ax^2 + bx + c)
 *              &b  (<optional return> coeff b of LSF: y = ax^2 + bx + c)
 *              &c  (<optional return> coeff c of LSF: y = ax^2 + bx + c)
 *              &nafit (<optional return> numa of least square fit)
 *      Return: 0 if OK, 1 on error
 *
 *  Notes:
 *      (1) This does a quadratic least square fit to the set of points
 *          in @pta.  That is, it finds coefficients a, b and c that minimize:
 *
 *              sum (yi - a*xi*xi -b*xi -c)^2
 *               i
 *
 *          The method is simple: differentiate this expression w/rt
 *          a, b and c, and solve the resulting three equations for these
 *          coefficients in terms of various sums over the input data (xi, yi).
 *          The three equations are in the form:
 *             f[0][0]a + f[0][1]b + f[0][2]c = g[0]
 *             f[1][0]a + f[1][1]b + f[1][2]c = g[1]
 *             f[2][0]a + f[2][1]b + f[2][2]c = g[2]
 *      (2) If @nafit is defined, this returns an array of fitted values,
 *          corresponding to the two implicit Numa arrays (nax and nay) in pta.
 *          Thus, just as you can plot the data in pta as nay vs. nax,
 *          you can plot the linear least square fit as nafit vs. nax.
 */


(1) This does a quadratic least square fit to the set of points
in @pta.  That is, it finds coefficients a, b and c that minimize.



The method is simple: differentiate this expression w/rt
a, b and c, and solve the resulting three equations for these
coefficients in terms of various sums over the input data (xi, yi).
( 1 ) This does a quadratic least square fit to the set of points in @pta . That is it finds coefficients a b and c that minimize . The method is simple : differentiate this expression w / rt a b and c and solve the resulting three equations for these coefficients in terms of various sums over the input data ( xi yi ) .
l_int32
ptaGetQuadraticLSF(PTA        *pta,
                   l_float32  *pa,
                   l_float32  *pb,
                   l_float32  *pc,
                   NUMA      **pnafit)
{
l_int32     n, i, ret;
l_float32   x, y, sx, sy, sx2, sx3, sx4, sxy, sx2y;
l_float32  *xa, *ya;
l_float32  *f[3];
l_float32   g[3];
NUMA       *nafit;
    PROCNAME("ptaGetQuadraticLSF");
    if (!pa && !pb && !pc && !pnafit)
        return ERROR_INT("no output requested", procName, 1);
    if (pa) *pa = 0.0;
    if (pb) *pb = 0.0;
    if (pc) *pc = 0.0;
    if (pnafit) *pnafit = NULL;
    if (!pta)
        return ERROR_INT("pta not defined", procName, 1);
    if ((n = ptaGetCount(pta)) < 3)
        return ERROR_INT("less than 3 pts not found", procName, 1);
    xa = pta->x;  
    ya = pta->y;  
    sx = sy = sx2 = sx3 = sx4 = sxy = sx2y = 0.;
    for (i = 0; i < n; i++) {
        x = xa[i];
        y = ya[i];
        sx += x;
        sy += y;
        sx2 += x * x;
        sx3 += x * x * x;
        sx4 += x * x * x * x;
        sxy += x * y;
        sx2y += x * x * y;
    }
    for (i = 0; i < 3; i++)
        f[i] = (l_float32 *)CALLOC(3, sizeof(l_float32));
    f[0][0] = sx4;
    f[0][1] = sx3;
    f[0][2] = sx2;
    f[1][0] = sx3;
    f[1][1] = sx2;
    f[1][2] = sx;
    f[2][0] = sx2;
    f[2][1] = sx;
    f[2][2] = n;
    g[0] = sx2y;
    g[1] = sxy;
    g[2] = sy;
    ret = gaussjordan(f, g, 3);
    for (i = 0; i < 3; i++)
        FREE(f[i]);
    if (ret)
        return ERROR_INT("quadratic solution failed", procName, 1);
    if (pa) *pa = g[0];
    if (pb) *pb = g[1];
    if (pc) *pc = g[2];
    if (pnafit) {
        nafit = numaCreate(n);
        *pnafit = nafit;
        for (i = 0; i < n; i++) {
            x = xa[i];
            y = g[0] * x * x + g[1] * x + g[2];
            numaAddNumber(nafit, y);
        }
    }
    return 0;
}
l_int32 ptaGetQuadraticLSF ( PTA * pta , l_float32 * pa , l_float32 * pb , l_float32 * pc , NUMA * * pnafit ) { l_int32 n , i , ret ; l_float32 x , y , sx , sy , sx2 , sx3 , sx4 , sxy , sx2y ; l_float32 * xa , * ya ; l_float32 * f [ 3 ] ; l_float32 g [ 3 ] ; NUMA * nafit ; PROCNAME ( " " ) ; if ( ! pa && ! pb && ! pc && ! pnafit ) return ERROR_INT ( " " , procName , 1 ) ; if ( pa ) * pa = 0.0 ; if ( pb ) * pb = 0.0 ; if ( pc ) * pc = 0.0 ; if ( pnafit ) * pnafit = NULL ; if ( ! pta ) return ERROR_INT ( " " , procName , 1 ) ; if ( ( n = ptaGetCount ( pta ) ) < 3 ) return ERROR_INT ( " " , procName , 1 ) ; xa = pta -> x ; ya = pta -> y ; sx = sy = sx2 = sx3 = sx4 = sxy = sx2y = 0. ; for ( i = 0 ; i < n ; i ++ ) { x = xa [ i ] ; y = ya [ i ] ; sx += x ; sy += y ; sx2 += x * x ; sx3 += x * x * x ; sx4 += x * x * x * x ; sxy += x * y ; sx2y += x * x * y ; } for ( i = 0 ; i < 3 ; i ++ ) f [ i ] = ( l_float32 * ) CALLOC ( 3 , sizeof ( l_float32 ) ) ; f [ 0 ] [ 0 ] = sx4 ; f [ 0 ] [ 1 ] = sx3 ; f [ 0 ] [ 2 ] = sx2 ; f [ 1 ] [ 0 ] = sx3 ; f [ 1 ] [ 1 ] = sx2 ; f [ 1 ] [ 2 ] = sx ; f [ 2 ] [ 0 ] = sx2 ; f [ 2 ] [ 1 ] = sx ; f [ 2 ] [ 2 ] = n ; g [ 0 ] = sx2y ; g [ 1 ] = sxy ; g [ 2 ] = sy ; ret = gaussjordan ( f , g , 3 ) ; for ( i = 0 ; i < 3 ; i ++ ) FREE ( f [ i ] ) ; if ( ret ) return ERROR_INT ( " " , procName , 1 ) ; if ( pa ) * pa = g [ 0 ] ; if ( pb ) * pb = g [ 1 ] ; if ( pc ) * pc = g [ 2 ] ; if ( pnafit ) { nafit = numaCreate ( n ) ; * pnafit = nafit ; for ( i = 0 ; i < n ; i ++ ) { x = xa [ i ] ; y = g [ 0 ] * x * x + g [ 1 ] * x + g [ 2 ] ; numaAddNumber ( nafit , y ) ; } } return 0 ; }
ptaGetQuadraticLSF()
      Input:  pta
              &a  (<optional return> coeff a of LSF: y = ax^2 + bx + c)
              &b  (<optional return> coeff b of LSF: y = ax^2 + bx + c)
              &c  (<optional return> coeff c of LSF: y = ax^2 + bx + c)
              &nafit (<optional return> numa of least square fit)
      Return: 0 if OK, 1 on error
ptaGetQuadraticLSF () Input : pta &a ( <optional return > coeff a of LSF : y = ax^2 + bx + c ) &b ( <optional return > coeff b of LSF : y = ax^2 + bx + c ) &c ( <optional return > coeff c of LSF : y = ax^2 + bx + c ) &nafit ( <optional return > numa of least square fit ) Return : 0 if OK 1 on error
/* not a copy */ /* not a copy */ /* Solve for the unknowns, also putting f-inverse into f */