def affine_matrix_from_points(v0, v1, shear=True, scale=True, usesvd=True):
v0 = numpy.array(v0, dtype=numpy.float64, copy=True)<EOL>v1 = numpy.array(v1, dtype=numpy.float64, copy=True)<EOL>ndims = v0.shape[<NUM_LIT:0>]<EOL>if ndims < <NUM_LIT:2> or v0.shape[<NUM_LIT:1>] < ndims or v0.shape != v1.shape:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>t0 = -numpy.mean(v0, axis=<NUM_LIT:1>)<EOL>M0 = numpy.identity(ndims+<NUM_LIT:1>)<EOL>M0[:ndims, ndims] = t0<EOL>v0 += t0.reshape(ndims, <NUM_LIT:1>)<EOL>t1 = -numpy.mean(v1, axis=<NUM_LIT:1>)<EOL>M1 = numpy.identity(ndims+<NUM_LIT:1>)<EOL>M1[:ndims, ndims] = t1<EOL>v1 += t1.reshape(ndims, <NUM_LIT:1>)<EOL>if shear:<EOL><INDENT>A = numpy.concatenate((v0, v1), axis=<NUM_LIT:0>)<EOL>u, s, vh = numpy.linalg.svd(A.T)<EOL>vh = vh[:ndims].T<EOL>B = vh[:ndims]<EOL>C = vh[ndims:<NUM_LIT:2>*ndims]<EOL>t = numpy.dot(C, numpy.linalg.pinv(B))<EOL>t = numpy.concatenate((t, numpy.zeros((ndims, <NUM_LIT:1>))), axis=<NUM_LIT:1>)<EOL>M = numpy.vstack((t, ((<NUM_LIT:0.0>,)*ndims) + (<NUM_LIT:1.0>,)))<EOL><DEDENT>elif usesvd or ndims != <NUM_LIT:3>:<EOL><INDENT>u, s, vh = numpy.linalg.svd(numpy.dot(v1, v0.T))<EOL>R = numpy.dot(u, vh)<EOL>if numpy.linalg.det(R) < <NUM_LIT:0.0>:<EOL><INDENT>R -= numpy.outer(u[:, ndims-<NUM_LIT:1>], vh[ndims-<NUM_LIT:1>, :]*<NUM_LIT>)<EOL>s[-<NUM_LIT:1>] *= -<NUM_LIT:1.0><EOL><DEDENT>M = numpy.identity(ndims+<NUM_LIT:1>)<EOL>M[:ndims, :ndims] = R<EOL><DEDENT>else:<EOL><INDENT>xx, yy, zz = numpy.sum(v0 * v1, axis=<NUM_LIT:1>)<EOL>xy, yz, zx = numpy.sum(v0 * numpy.roll(v1, -<NUM_LIT:1>, axis=<NUM_LIT:0>), axis=<NUM_LIT:1>)<EOL>xz, yx, zy = numpy.sum(v0 * numpy.roll(v1, -<NUM_LIT:2>, axis=<NUM_LIT:0>), axis=<NUM_LIT:1>)<EOL>N = [[xx+yy+zz, <NUM_LIT:0.0>,      <NUM_LIT:0.0>,      <NUM_LIT:0.0>],<EOL>[yz-zy,    xx-yy-zz, <NUM_LIT:0.0>,      <NUM_LIT:0.0>],<EOL>[zx-xz,    xy+yx,    yy-xx-zz, <NUM_LIT:0.0>],<EOL>[xy-yx,    zx+xz,    yz+zy,    zz-xx-yy]]<EOL>w, V = numpy.linalg.eigh(N)<EOL>q = V[:, numpy.argmax(w)]<EOL>q /= vector_norm(q)  <EOL>M = quaternion_matrix(q)<EOL><DEDENT>if scale and not shear:<EOL><INDENT>v0 *= v0<EOL>v1 *= v1<EOL>M[:ndims, :ndims] *= math.sqrt(numpy.sum(v1) / numpy.sum(v0))<EOL><DEDENT>M = numpy.dot(numpy.linalg.inv(M1), numpy.dot(M, M0))<EOL>M /= M[ndims, ndims]<EOL>return M<EOL>
Return affine transform matrix to register two point sets.

    v0 and v1 are shape (ndims, *) arrays of at least ndims non-homogeneous
    coordinates, where ndims is the dimensionality of the coordinate space.

    If shear is False, a similarity transformation matrix is returned.
    If also scale is False, a rigid/Euclidean transformation matrix
    is returned.

    By default the algorithm by Hartley and Zissermann [15] is used.
    If usesvd is True, similarity and Euclidean transformation matrices
    are calculated by minimizing the weighted sum of squared deviations
    (RMSD) according to the algorithm by Kabsch [8].
    Otherwise, and if ndims is 3, the quaternion based algorithm by Horn [9]
    is used, which is slower when using this Python implementation.

    The returned matrix performs rotation, translation and uniform scaling
    (if specified).

    >>> v0 = [[0, 1031, 1031, 0], [0, 0, 1600, 1600]]
    >>> v1 = [[675, 826, 826, 677], [55, 52, 281, 277]]
    >>> affine_matrix_from_points(v0, v1)
    array([[   0.14549,    0.00062,  675.50008],
           [   0.00048,    0.14094,   53.24971],
           [   0.     ,    0.     ,    1.     ]])
    >>> T = translation_matrix(numpy.random.random(3)-0.5)
    >>> R = random_rotation_matrix(numpy.random.random(3))
    >>> S = scale_matrix(random.random())
    >>> M = concatenate_matrices(T, R, S)
    >>> v0 = (numpy.random.rand(4, 100) - 0.5) * 20
    >>> v0[3] = 1
    >>> v1 = numpy.dot(M, v0)
    >>> v0[:3] += numpy.random.normal(0, 1e-8, 300).reshape(3, -1)
    >>> M = affine_matrix_from_points(v0[:3], v1[:3])
    >>> numpy.allclose(v1, numpy.dot(M, v0))
    True

    More examples in superimposition_matrix()