Calculating / Estimating the normal at any location on a face

Question

In a recent answer I employed Object.ray_cast which returns the location and normal of a ray hitting a face. The normal is needed to calculate the rays reflection, ray.reflect(norm). The normal returned for any hit on a face is the face's normal, not as I was hoping the normal for that location on the face.

In image camera is projecting a grid of rays onto sphere and reflecting them back to create the plane.

The result of using face normal can be clearly seen in reflected mesh. Mimics the rings / loops of the sphere.

To get around this, as an estimate a barycentric interpolation using the triangle constructed from the centre of the face, and the vertices of the outer edge that the hit location falls within. Using the normals at these three points, the location of hit on face can be plugged into mathutils.geometry.barycentric_transform(...) to estimate the normal at the point of hit.

def calc_normal(mesh, face_index, point):
    face = mesh.polygons[face_index]
    o = face.center
    p = point - o
    # find the edge

    for ek in face.edge_keys:
        v0, v1 = (mesh.vertices[i] for i in ek)
        if (v0.co - o).angle(p) <= (v0.co - o).angle(v1.co - o):
            break

    return barycentric_transform(point, o, v0.co, v1.co,
                                face.normal, v0.normal, v1.normal)



def obj_raycast(obj, ray_origin, ray_target, matrix=None):
    # returns global coords
    if matrix is None:
        matrix = obj.matrix_world

    # get the ray relative to the object
    matrix_inv = matrix.inverted()
    ray_origin_obj = matrix_inv * ray_origin
    ray_target_obj = matrix_inv * ray_target
    ray_direction_obj = ray_target_obj - ray_origin_obj

    # cast the ray
    success, location, normal, face_index = obj.ray_cast(ray_origin_obj, ray_direction_obj)

    if success:
        n = calc_normal(obj.data, face_index, location)
        return matrix * location, (matrix * (location + n) - matrix * location).normalized(), face_index
    else:
        return None, None, None

Which produces this "ok" result.

Are there other ways to calculate / estimate the surface normal at any given point. Perhaps using BVHTree Utilities ray_cast