Function AlignmentMat

jni/stasm/misc.cpp:522–565

// www.isbe.man.ac.uk/~bim/Mods/app_models.pdf.

const MAT AlignmentMat(
    const Shape&  shape,       // in
    const Shape&  anchorshape, // in
    const double* weights)     // in: if NULL (default) all points equally weighted
{
    double W = 0;
    double sx = 0, sy = 0, sx1 = 0, sy1 = 0;
    double sxx_syy = 0, sxx1_syy1 = 0, sxy1_syx1 = 0;

    for (int i = 0; i < shape.rows; i++)
    {
        const double x  = shape(i, IX);
        const double y  = shape(i, IY);
        const double x1 = anchorshape(i, IX);
        const double y1 = anchorshape(i, IY);
        if (PointUsed(x, y) && PointUsed(x1, y1))
        {
            const double w = (weights? weights[i]: 1.);
            W   += w;
            sx  += w * x;
            sy  += w * y;
            sx1 += w * x1;
            sy1 += w * y1;
            sxy1_syx1 += w * (x * y1 - y * x1);
            sxx1_syy1 += w * (x * x1 + y * y1);
            sxx_syy   += w * (x * x  + y * y);
        }
    }
    MAT soln_data = (MAT(4,4) <<  sxx_syy,       0,  sx,  sy,
                                        0, sxx_syy, -sy,  sx,
                                       sx,     -sy,   W,   0,
                                       sy,      sx,   0,   W );

    VEC vec_data = (MAT(4,1) << sxx1_syy1,
                                sxy1_syx1,
                                      sx1,
                                      sy1 );

    const VEC soln(LinSolveLu(soln_data, vec_data));

    return (MAT(3, 3) << soln(0), -soln(1), soln(2),    //  a -b tx
                         soln(1),  soln(0), soln(3),    //  b  a ty
                               0,        0,       1 );  //  0  0  1
}

Shape ShiftShape( // add xshift and yshift to shape coords, skipping unused points
    const Shape& shape,   // in