import torch


def compute_jacobian_on_surface(u, v, forward_transform, eps=0.01):
    """
    Computes the differentials: 

    [dX/dv, dY/dv, dX/dv, dX/du, dY/du, dX/du] 

    for the given projection function) using central differences. u and v are an orthogonal coordinate system on the surface and X, Y, Z are 3D Cartesian coordinates..

    Returns (u.shape[0], u.shape[1], 2, 2)
    """

    # Compute dX/du, dY/du, dZ, du
    x0, y0, z0 = forward_transform(u - eps, v)
    x1, y1, z1 = forward_transform(u + eps, v)
    dx_du = (x1 - x0) / (2 * eps)
    dy_du = (y1 - y0) / (2 * eps)
    dz_du = (z1 - z0) / (2 * eps)

    # Compute dX/dv, dY/dv, dZ/dv
    x2, y2, z2 = forward_transform(u, v - eps)
    x3, y3, z3 = forward_transform(u, v + eps)
    dx_dv = (x3 - x2) / (2 * eps)
    dy_dv = (y3 - y2) / (2 * eps)
    dz_dv = (z3 - z2) / (2 * eps)

    return torch.stack((torch.stack(
        (dx_du, dy_du, dz_du), -1), torch.stack((dx_dv, dy_dv, dz_dv), -1)),
                       -1)