import numpy as np
from scipy.spatial import Delaunay


points = np.array(
    [
        [0.0, 0.0],
        [2.0, 0.0],
        [2.0, 1.2],
        [0.0, 1.0],
        [0.9, 0.4],
        [1.3, 0.8],
    ]
)

triangulation = Delaunay(points)

print(f"input points: {len(points)}")
print(f"triangles: {len(triangulation.simplices)}")
print("simplices:")
for simplex_index, simplex in enumerate(triangulation.simplices):
    vertices = ", ".join(str(int(vertex)) for vertex in simplex)
    print(f"  {simplex_index}: [{vertices}]")

first_simplex = triangulation.simplices[0]
first_coordinates = [
    [float(f"{points[vertex, 0]:.1f}"), float(f"{points[vertex, 1]:.1f}")]
    for vertex in first_simplex
]
print("first triangle coordinates:")
for coordinate in first_coordinates:
    print(f"  {coordinate}")

query_points = np.array(
    [
        [0.2, 0.2],
        [1.1, 0.6],
        [2.3, 0.5],
    ]
)

print("point lookup:")
for point, simplex_index in zip(query_points, triangulation.find_simplex(query_points)):
    x, y = point
    if simplex_index == -1:
        print(f"  ({x:.1f}, {y:.1f}) -> outside triangulation")
        continue
    vertices = triangulation.simplices[simplex_index].tolist()
    print(f"  ({x:.1f}, {y:.1f}) -> simplex {int(simplex_index)} vertices {vertices}")