def cartesian_to_barycentric_3D(tri, xy):
    '''
    cartesian_to_barycentric_3D(tri,xy) is identical to cartesian_to_barycentric_2D(tri,xy) except
    it works on 3D data. Note that if tri is a 3 x 3 x n, a 3 x n x 3 or an n x 3 x 3 matrix, the
    first dimension must always be the triangle vertices and the second 3-sized dimension must be
    the (x,y,z) coordinates.
    '''
    xy = np.asarray(xy)
    tri = np.asarray(tri)
    if len(xy.shape) == 1:
        return cartesian_to_barycentric_3D(np.transpose(np.asarray([tri]), (1,2,0)),
                                           np.asarray([xy]).T)[:,0]
    xy = xy if xy.shape[0] == 3 else xy.T
    if tri.shape[0] == 3:
        tri = tri if tri.shape[1] == 3 else np.transpose(tri, (0,2,1))
    elif tri.shape[1] == 3:
        tri = tri.T if tri.shape[0] == 3 else np.transpose(tri, (1,2,0))
    elif tri.shape[2] == 3:
        tri = np.transpose(tri, (2,1,0) if tri.shape[1] == 3 else (2,0,1))
    if tri.shape[0] != 3 or tri.shape[1] != 3:
        raise ValueError('Triangle array did not have dimensions of sizes 3 and 3')
    if xy.shape[0] != 3:
        raise ValueError('coordinate matrix did not have a dimension of size 3')
    if tri.shape[2] != xy.shape[1]:
        raise ValueError('number of triangles and coordinates must match')
    # The algorithm here is borrowed from this stack-exchange post:
    # http://gamedev.stackexchange.com/questions/23743
    # in which it is attributed to Christer Ericson's book Real-Time Collision Detection.
    v0 = tri[1] - tri[0]
    v1 = tri[2] - tri[0]
    v2 = xy - tri[0]
    d00 = np.sum(v0 * v0, axis=0)
    d01 = np.sum(v0 * v1, axis=0)
    d11 = np.sum(v1 * v1, axis=0)
    d20 = np.sum(v2 * v0, axis=0)
    d21 = np.sum(v2 * v1, axis=0)
    den = d00*d11 - d01*d01
    zero = np.isclose(den, 0)
    unit = 1 - zero
    den += zero
    l2 = unit * (d11 * d20 - d01 * d21) / den
    l3 = unit * (d00 * d21 - d01 * d20) / den
    return np.asarray([1.0 - l2 - l3, l2])
def cartesian_to_barycentric_3D ( tri , xy ) : xy = np . asarray ( xy ) tri = np . asarray ( tri ) if len ( xy . shape ) == 1 : return cartesian_to_barycentric_3D ( np . transpose ( np . asarray ( [ tri ] ) , ( 1 , 2 , 0 ) ) , np . asarray ( [ xy ] ) . T ) [ : , 0 ] xy = xy if xy . shape [ 0 ] == 3 else xy . T if tri . shape [ 0 ] == 3 : tri = tri if tri . shape [ 1 ] == 3 else np . transpose ( tri , ( 0 , 2 , 1 ) ) elif tri . shape [ 1 ] == 3 : tri = tri . T if tri . shape [ 0 ] == 3 else np . transpose ( tri , ( 1 , 2 , 0 ) ) elif tri . shape [ 2 ] == 3 : tri = np . transpose ( tri , ( 2 , 1 , 0 ) if tri . shape [ 1 ] == 3 else ( 2 , 0 , 1 ) ) if tri . shape [ 0 ] != 3 or tri . shape [ 1 ] != 3 : raise ValueError ( 'Triangle array did not have dimensions of sizes 3 and 3' ) if xy . shape [ 0 ] != 3 : raise ValueError ( 'coordinate matrix did not have a dimension of size 3' ) if tri . shape [ 2 ] != xy . shape [ 1 ] : raise ValueError ( 'number of triangles and coordinates must match' ) v0 = tri [ 1 ] - tri [ 0 ] v1 = tri [ 2 ] - tri [ 0 ] v2 = xy - tri [ 0 ] d00 = np . sum ( v0 * v0 , axis = 0 ) d01 = np . sum ( v0 * v1 , axis = 0 ) d11 = np . sum ( v1 * v1 , axis = 0 ) d20 = np . sum ( v2 * v0 , axis = 0 ) d21 = np . sum ( v2 * v1 , axis = 0 ) den = d00 * d11 - d01 * d01 zero = np . isclose ( den , 0 ) unit = 1 - zero den += zero l2 = unit * ( d11 * d20 - d01 * d21 ) / den l3 = unit * ( d00 * d21 - d01 * d20 ) / den return np . asarray ( [ 1.0 - l2 - l3 , l2 ] )
cartesian_to_barycentric_3D(tri,xy) is identical to cartesian_to_barycentric_2D(tri,xy) except
    it works on 3D data. Note that if tri is a 3 x 3 x n, a 3 x n x 3 or an n x 3 x 3 matrix, the
    first dimension must always be the triangle vertices and the second 3-sized dimension must be
    the (x,y,z) coordinates.