// Symmetric point-to-plane alignment.
static bool align_symm(const ::std::vector<PtPair> &pairs,
    float scale, point &centroid1, point &centroid2, xform &alignxf)
{
    size_t npairs = pairs.size();
    double A[6][6] = { { 0 } }, b[6] = { 0 };

    for (size_t i = 0; i < npairs; i++) {
        dvec3 p1 = dvec3(scale * (pairs[i].p1 - centroid1));
        dvec3 p2 = dvec3(scale * (pairs[i].p2 - centroid2));
        dvec3 n = dvec3(pairs[i].n1 + pairs[i].n2);
        dvec3 p = p1 + p2;
        dvec3 c = p CROSS n;
        dvec3 d = p1 - p2;

        double x[6] = { c[0], c[1], c[2], n[0], n[1], n[2] };
        double dn = d DOT n;

        for (int j = 0; j < 6; j++) {
            b[j] += dn * x[j];
            for (int k = j; k < 6; k++)
                A[j][k] += x[j] * x[k];
        }
    }

    // Make matrix symmetric
    for (int j = 1; j < 6; j++)
        for (int k = 0; k < j; k++)
            A[j][k] = A[k][j];

    // Eigen-decomposition and inverse
    double eval[6], einv[6];
    eigdc<double,6>(A, eval);
    for (int i = 0; i < 6; i++) {
        if (eval[i] < 1e-6)
            return false;
        else
            einv[i] = 1.0 / eval[i];
    }

    // Solve system
    eigmult<double,6>(A, einv, b);

    // Extract rotation and translation
    dvec3 rot(b[0], b[1], b[2]), trans(b[3], b[4], b[5]);
    double rotangle = atan(len(rot));
    trans *= 1.0 / scale;
    trans *= cos(rotangle);

    xform R = xform::rot(rotangle, rot);
    alignxf = xform::trans(centroid1) *
              R * xform::trans(trans) * R *
              xform::trans(-centroid2);

    return true;
}