def random_sample_point(
    state: state_def.State,
    obj: bpy.types.Object,
    face_mask: np.ndarray,
    plane: tuple[str, int],
) -> Vector:
    """
    Given a plane, return a random point on the plane.
    """

    if obj.type != "MESH":
        raise ValueError(f"Unexpected {obj.type=}")

    plane_mask = state.planes.tagged_plane_mask(obj, face_mask, plane)
    if not np.any(plane_mask):
        logging.warning(
            f"No faces in object {obj.name} are coplanar with plane {plane}."
        )

    # Create a bmesh from the object mesh
    bm = bmesh.new()
    bm.from_mesh(obj.data)
    bm.faces.ensure_lookup_table()

    faces = [bm.faces[i] for i in np.where(plane_mask)[0]]

    # Calculate the area for each face and create a cumulative distribution
    areas = np.array([f.calc_area() for f in faces])
    cumulative_areas = np.cumsum(areas)
    total_area = cumulative_areas[-1]

    # Generate a random number and find the corresponding face
    random_area_point = np.random.rand() * total_area
    face_index = np.searchsorted(cumulative_areas, random_area_point)
    selected_face = faces[face_index]

    verts = [v.co for v in selected_face.verts]

    # Use barycentric coordinates to sample a random point in the triangle
    # Random weights for each vertex
    weights = np.random.rand(3)
    weights /= np.sum(weights)
    random_point_local = (
        weights[0] * verts[0] + weights[1] * verts[1] + weights[2] * verts[2]
    )
    random_point_global = obj.matrix_world @ Vector(random_point_local)

    bm.free()

    return random_point_global

    Given a plane, return a random point on the plane.
    
infinigen/core/constraints/example_solver/geometry/stability.py