def plane_fit(points):
    """
    p, n = plane_fit(points)

    Given an array, points, of shape (d,...)
    representing points in d-dimensional space,
    fit an d-dimensional plane to the points.
    Return a point, p, on the plane (the point-cloud centroid),
    and the normal, n.
    """
    import numpy as np
    from numpy.linalg import svd
    points = np.reshape(points, (np.shape(points)[0], -1))
    assert points.shape[0] <= points.shape[1
        ], 'There are only {} points in {} dimensions.'.format(points.shape
        [1], points.shape[0])
    ctr = points.mean(axis=1)
    x = points - ctr[:, (np.newaxis)]
    M = np.dot(x, x.T)
    return ctr, svd(M)[0][:, (-1)]
