From e40fff5de5106ef51d87a0e0d4ab8e3ae4e3c5f2 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 19:27:04 +0530
Subject: [PATCH 01/29] Implement spherical triangulation

---
 .../spherical_delaunay_stereographic.jl       | 168 ++++++++++++++++++
 1 file changed, 168 insertions(+)
 create mode 100644 docs/src/experiments/spherical_delaunay_stereographic.jl

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
new file mode 100644
index 000000000..ee195ffcd
--- /dev/null
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -0,0 +1,168 @@
+#=
+# Spherical Delaunay triangulation using stereographic projections
+
+This is the approach which d3-geo-voronoi and friends use.  The alternative is STRIPACK which computes in 3D on the sphere.
+The 3D approach is basically that the 3d convex hull of a set of points on the sphere is its Delaunay triangulation.
+=#
+
+import GeometryOps as GO, GeoInterface as GI
+import Proj # for easy stereographic projection - TODO implement in Julia
+import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
+import CoordinateTransformations, Rotations
+
+using Downloads # does what it says on the tin
+using JSON3 # to load data
+using CairoMakie # for plotting
+import Makie: Point3d
+
+points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+
+pivot_ind = findfirst(isfinite, points)
+pivot_point = points[pivot_ind]
+necessary_rotation = Rotations.RotY(pi) * Rotations.RotY(-pivot_point[2]) * Rotations.RotZ(-pivot_point[1])
+
+net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
+
+stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+
+# Iterate through the points and find infinite/invalid points
+max_radius = 1
+really_far_idxs = Int[]
+for (i, point) in enumerate(stereographic_points)
+    m = hypot(point[1], point[2])
+    if !isfinite(m) || m > 1e32
+        push!(really_far_idxs, i)
+    elseif m > max_radius
+        max_radius = m
+    end
+end
+far_value = 1e6 * sqrt(max_radius)
+if !isempty(really_far_idxs)
+    stereographic_points[really_far_idxs] .= Point2(far_value, 0)
+end
+
+# Add infinite horizon points
+triangulation_points = copy(stereographic_points)
+push!(triangulation_points, Point2(0, far_value))
+push!(triangulation_points, Point2(-far_value, 0))
+push!(triangulation_points, Point2(0, -far_value))
+
+triangulation = DelTri.triangulate(triangulation_points)
+
+triangulation.points[1:end-3] .= stereographic_points
+triangulation.points[end-2:end] .= (stereographic_points[pivot_ind],)
+
+f, a, p = triplot(triangulation; axis = (; type = Axis3,));
+m = p.plots[1].plots[1][1][]
+
+geo_mesh_points = vcat(points, repeat([pivot_point], 3))
+cartesian_mesh_points = UnitCartesianFromGeographic().(geo_mesh_points)
+f, a, p = mesh(cartesian_mesh_points, getfield(getfield(m, :simplices), :faces))
+wireframe(p.mesh[])
+CairoMakie.GeometryBasics.Mesh(UnitCartesianFromGeographic().(vcat(points, repeat([pivot_point], 3))), splat(CairoMakie.GeometryBasics.TriangleFace).(collect(triangulation.triangles))) |> wireframe
+
+geo_triangulation = deepcopy(triangulation)
+geo_triangulation.points .= geo_mesh_points
+
+itp = NaturalNeighbours.interpolate(geo_triangulation, last.(geo_mesh_points); derivatives = true)
+lons = LinRange(-180, 180, 360)
+lats = LinRange(-90, 90, 180)
+_x = vec([x for x in lons, _ in lats])
+_y = vec([y for _ in lons, y in lats])
+
+sibson_1_vals = itp(_x, _y; method=Sibson(1))
+
+
+struct StereographicFromCartesian <: CoordinateTransformations.Transformation
+end
+
+function (::StereographicFromCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    return Point2(x/(1-z), y/(1-z))
+end
+
+struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
+end
+
+function (::UnitCartesianFromGeographic)(geographic_point)
+    # Longitude is directly translatable to a spherical coordinate
+    # θ (azimuth)
+    θ = deg2rad(GI.x(geographic_point))
+    # The polar angle is 90 degrees minus the latitude
+    # ϕ (polar angle)
+    ϕ = deg2rad(90 - GI.y(geographic_point))
+    # Since this is the unit sphere, the radius is assumed to be 1,
+    # and we don't need to multiply by it.
+    return Point3(
+        sin(ϕ) * cos(θ),
+        sin(ϕ) * sin(θ),
+        cos(ϕ)
+    )
+end
+
+struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation 
+end
+
+function (::GeographicFromUnitCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    return Point2(
+        atan(y, x),
+        atan(hypot(x, y), z),
+    )
+end
+
+transformed_points = GO.transform(points) do point
+    # first, spherical to Cartesian
+    longitude, latitude = deg2rad(GI.x(point)), deg2rad(GI.y(point))
+    radius = 1 # we will operate on the unit sphere
+    xyz = Point3d(
+        radius * sin(latitude) * cos(longitude),
+        radius * sin(latitude) * sin(longitude),
+        radius * cos(latitude)
+    )
+    # then, rotate Cartesian so the pivot point is at the south pole
+    rotated_xyz = Rotations.AngleAxis
+end
+
+
+
+
+function randsphericalangles(n)
+    θ = 2π .* rand(n)
+    ϕ = acos.(2 .* rand(n) .- 1)
+    return Point2.(θ, ϕ)
+end
+
+function randsphere(n)
+    θϕ = randsphericalangles(n)
+    return Point3.(
+        sin.(last.(θϕs)) .* cos.(first.(θϕs)),
+        sin.(last.(θϕs)) .* sin.(first.(θϕs)),
+        cos.(last.(θϕs))
+    )
+end
+
+θϕs = randsphericalangles(50)
+θs, ϕs = first.(θϕs), last.(θϕs)
+pts = Point3.(
+    sin.(ϕs) .* cos.(θs),
+    sin.(ϕs) .* sin.(θs),
+    cos.(ϕs)
+)
+
+f, a, p = scatter(pts; color = [fill(:blue, 49); :red])
+
+function Makie.rotate!(t::Makie.Transformable, rotation::Rotations.Rotation)
+    quat = Rotations.QuatRotation(rotation)
+    Makie.rotate!(Makie.Absolute, t, Makie.Quaternion(quat.q.v1, quat.q.v2, quat.q.v3, quat.q.s))
+end
+
+
+rotate!(p, Rotations.RotX(π/2))
+rotate!(p, Rotations.RotX(0))
+
+pivot = θϕs[end]
+rotate!(p, )
+

From ef028e459a99346522355130f5a75b4c8a7a20d0 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 21:53:38 +0530
Subject: [PATCH 02/29] Use DelaunayTriangulation's API to get better results

---
 .../spherical_delaunay_stereographic.jl       | 73 ++++++++++++++++---
 1 file changed, 64 insertions(+), 9 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index ee195ffcd..89c57e85c 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -15,15 +15,28 @@ using JSON3 # to load data
 using CairoMakie # for plotting
 import Makie: Point3d
 
+# These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
 
 pivot_ind = findfirst(isfinite, points)
+
+# debug
+point_colors = fill(:blue, length(points))
+point_colors[pivot_ind] = :red
+# end debug
+
 pivot_point = points[pivot_ind]
-necessary_rotation = Rotations.RotY(pi) * Rotations.RotY(-pivot_point[2]) * Rotations.RotZ(-pivot_point[1])
+necessary_rotation = Rotations.RotY(-π) * Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
 
 net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
 
+scatter(map(net_transformation_to_corrected_cartesian, points); color = point_colors)
+#
 stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+scatter(stereographic_points; color = point_colors)
+
+triangulation_points = copy(stereographic_points)
+
 
 # Iterate through the points and find infinite/invalid points
 max_radius = 1
@@ -36,18 +49,48 @@ for (i, point) in enumerate(stereographic_points)
         max_radius = m
     end
 end
+@show max_radius length(really_far_idxs)
 far_value = 1e6 * sqrt(max_radius)
+
 if !isempty(really_far_idxs)
-    stereographic_points[really_far_idxs] .= Point2(far_value, 0)
+    triangulation_points[really_far_idxs] .= Point2(far_value, 0)
 end
 
-# Add infinite horizon points
+boundary_points = reverse([
+    (-far_value, -far_value),
+    (-far_value, far_value),
+    (far_value, far_value),
+    (far_value, -far_value),
+    (-far_value, -far_value),
+])
 triangulation_points = copy(stereographic_points)
-push!(triangulation_points, Point2(0, far_value))
-push!(triangulation_points, Point2(-far_value, 0))
-push!(triangulation_points, Point2(0, -far_value))
 
-triangulation = DelTri.triangulate(triangulation_points)
+boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+
+triangulation = DelTri.triangulate(pts; boundary_nodes)
+triplot(triangulation)
+#  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
+
+# First, get all the "solid" faces, ie, faces not attached to boundary nodes
+faces = map(DelTri.each_solid_triangle(triangulation)) do face_tuple
+    CairoMakie.GeometryBasics.TriangleFace(map(face_tuple) do i; i in DelTri.get_boundary_nodes(triangulation) ? pivot_ind : i end)
+end
+for ghost_face in DelTri.each_ghost_triangle(triangulation)
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
+end
+
+wireframe(CairoMakie.GeometryBasics.normal_mesh(UnitCartesianFromGeographic().(points), faces))
+
+_f, _a, _tri_p = triplot(triangulation; axis = (; type = Axis3));
+msh = _tri_p.plots[1].plots[1][1][]
+
+msh3d = GeometryBasics.Mesh(Makie.to_ndim(Point3d, GeometryBasics.coordinates(msh), 0), msh.faces)
+f, a, p = wireframe(msh3d; axis = (; type = LScene));
+p.transformation.transform_func[] = Makie.PointTrans{3}() do point
+    CartesianFromStereographic()(point)
+end
+
+f
 
 triangulation.points[1:end-3] .= stereographic_points
 triangulation.points[end-2:end] .= (stereographic_points[pivot_ind],)
@@ -82,6 +125,15 @@ function (::StereographicFromCartesian)(xyz::AbstractVector)
     return Point2(x/(1-z), y/(1-z))
 end
 
+struct CartesianFromStereographic <: CoordinateTransformations.Transformation
+end
+
+function (::CartesianFromStereographic)(stereographic_point)
+    X, Y = stereographic_point
+    x2y2_1 = X^2 + Y^2 + 1
+    return Point3(2X/x2y2_1, 2Y/x2y2_1, (x2y2_1 - 2)/x2y2_1)
+end
+
 struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
 end
 
@@ -163,6 +215,9 @@ end
 rotate!(p, Rotations.RotX(π/2))
 rotate!(p, Rotations.RotX(0))
 
-pivot = θϕs[end]
-rotate!(p, )
+pivot_point = θϕs[end]
+necessary_rotation = Rotations.RotY(pi) * Rotations.RotY(-pivot_point[2]) * Rotations.RotZ(-pivot_point[1])
+
+rotate!(p, necessary_rotation)
 
+f
\ No newline at end of file

From 27a71248c69d62b4224916b77b10c8ab58c1a776 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 22:39:49 +0530
Subject: [PATCH 03/29] North pole is at infinity if I use wikipedia (which I
 do)

---
 docs/src/experiments/spherical_delaunay_stereographic.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 89c57e85c..abb0a8e4e 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -26,7 +26,7 @@ point_colors[pivot_ind] = :red
 # end debug
 
 pivot_point = points[pivot_ind]
-necessary_rotation = Rotations.RotY(-π) * Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
+necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
 
 net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
 
@@ -106,6 +106,8 @@ CairoMakie.GeometryBasics.Mesh(UnitCartesianFromGeographic().(vcat(points, repea
 
 geo_triangulation = deepcopy(triangulation)
 geo_triangulation.points .= geo_mesh_points
+# try NaturalNeighbours, fail miserably
+using NaturalNeighbours
 
 itp = NaturalNeighbours.interpolate(geo_triangulation, last.(geo_mesh_points); derivatives = true)
 lons = LinRange(-180, 180, 360)
@@ -115,6 +117,7 @@ _y = vec([y for _ in lons, y in lats])
 
 sibson_1_vals = itp(_x, _y; method=Sibson(1))
 
+# necessary coordinate transformations
 
 struct StereographicFromCartesian <: CoordinateTransformations.Transformation
 end
@@ -122,6 +125,8 @@ end
 function (::StereographicFromCartesian)(xyz::AbstractVector)
     @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
     x, y, z = xyz
+    # The Wikipedia definition has the north pole at infinity,
+    # this implementation has the south pole at infinity.
     return Point2(x/(1-z), y/(1-z))
 end
 

From 0476e3e001ff8209264477cc312e4a55ab32687b Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 22:42:05 +0530
Subject: [PATCH 04/29] visualize the transform etc for better understanding

---
 .../spherical_delaunay_stereographic.jl       | 59 +++++++++++--------
 1 file changed, 33 insertions(+), 26 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index abb0a8e4e..9f9a95767 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -65,51 +65,58 @@ boundary_points = reverse([
 ])
 triangulation_points = copy(stereographic_points)
 
-boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+# boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
 
-triangulation = DelTri.triangulate(pts; boundary_nodes)
+# triangulation = DelTri.triangulate(triangulation_points; boundary_nodes)
+triangulation = DelTri.triangulate(triangulation_points)
 triplot(triangulation)
 #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
 
 # First, get all the "solid" faces, ie, faces not attached to boundary nodes
-faces = map(DelTri.each_solid_triangle(triangulation)) do face_tuple
-    CairoMakie.GeometryBasics.TriangleFace(map(face_tuple) do i; i in DelTri.get_boundary_nodes(triangulation) ? pivot_ind : i end)
-end
+boundary_faces = Iterators.filter(Base.Fix1(DelTri.is_boundary_triangle, triangulation), DelTri.each_solid_triangle(triangulation)) |> collect
+faces = map(CairoMakie.GeometryBasics.TriangleFace, DelTri.each_solid_triangle(triangulation))
+
 for ghost_face in DelTri.each_ghost_triangle(triangulation)
     push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
 end
 
-wireframe(CairoMakie.GeometryBasics.normal_mesh(UnitCartesianFromGeographic().(points), faces))
+wireframe(CairoMakie.GeometryBasics.normal_mesh(map(UnitCartesianFromGeographic(), points), faces))
 
-_f, _a, _tri_p = triplot(triangulation; axis = (; type = Axis3));
-msh = _tri_p.plots[1].plots[1][1][]
+grid_lons = LinRange(-180, 180, 10)
+grid_lats = LinRange(-90, 90, 10)
+lon_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((grid_lons[j], -90)), Point{2, Float64}((grid_lons[j], 90))]) for j in 1:10] |> vec
+lat_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((-180, grid_lats[i])), Point{2, Float64}((180, grid_lats[i]))]) for i in 1:10] |> vec
+
+sampled_grid = GO.segmentize(lat_grid; max_distance = 0.01)
+
+geographic_grid = GO.transform(UnitCartesianFromGeographic(), sampled_grid) .|> x -> GI.LineString{true, false}(x.geom)
+stereographic_grid = GO.transform(sampled_grid) do point
+    (StereographicFromCartesian() ∘ UnitCartesianFromGeographic())(point)
+end .|> x -> GI.LineString{false, false}(x.geom)
 
-msh3d = GeometryBasics.Mesh(Makie.to_ndim(Point3d, GeometryBasics.coordinates(msh), 0), msh.faces)
-f, a, p = wireframe(msh3d; axis = (; type = LScene));
-p.transformation.transform_func[] = Makie.PointTrans{3}() do point
-    CartesianFromStereographic()(point)
-end
 
-f
 
-triangulation.points[1:end-3] .= stereographic_points
-triangulation.points[end-2:end] .= (stereographic_points[pivot_ind],)
+fig = Figure(); ax = LScene(fig[1, 1])
+for (i, line) in enumerate(geographic_grid)
+    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
+end
+fig
+
+#
+fig = Figure(); ax = LScene(fig[1, 1])
+for (i, line) in enumerate(stereographic_grid)
+    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
+end
+fig
 
-f, a, p = triplot(triangulation; axis = (; type = Axis3,));
-m = p.plots[1].plots[1][1][]
+all(reduce(vcat, faces) |> sort! |> unique! .== 1:length(points))
 
-geo_mesh_points = vcat(points, repeat([pivot_point], 3))
-cartesian_mesh_points = UnitCartesianFromGeographic().(geo_mesh_points)
-f, a, p = mesh(cartesian_mesh_points, getfield(getfield(m, :simplices), :faces))
-wireframe(p.mesh[])
-CairoMakie.GeometryBasics.Mesh(UnitCartesianFromGeographic().(vcat(points, repeat([pivot_point], 3))), splat(CairoMakie.GeometryBasics.TriangleFace).(collect(triangulation.triangles))) |> wireframe
+faces[findall(x -> 1 in x, faces)]
 
-geo_triangulation = deepcopy(triangulation)
-geo_triangulation.points .= geo_mesh_points
 # try NaturalNeighbours, fail miserably
 using NaturalNeighbours
 
-itp = NaturalNeighbours.interpolate(geo_triangulation, last.(geo_mesh_points); derivatives = true)
+itp = NaturalNeighbours.interpolate(triangulation, last.(points); derivatives = true)
 lons = LinRange(-180, 180, 360)
 lats = LinRange(-90, 90, 180)
 _x = vec([x for x in lons, _ in lats])

From c6da37fa3936d77ebed34cb33ce8fdfad2b29c96 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 22:49:33 +0530
Subject: [PATCH 05/29] Working spherical triangulation function

---
 .../spherical_delaunay_stereographic.jl       | 61 ++++++++++++++++++-
 1 file changed, 59 insertions(+), 2 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 9f9a95767..5d4512db9 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -15,9 +15,67 @@ using JSON3 # to load data
 using CairoMakie # for plotting
 import Makie: Point3d
 
+function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+    # @assert GI.crstrait(first(input_points)) isa GI.AbstractGeographicCRSTrait
+    points = GO.tuples(input_points)
+
+    pivot_ind = findfirst(x -> all(isfinite, x), points)
+    
+    pivot_point = points[pivot_ind]
+    necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
+
+    net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
+    stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+    triangulation_points = copy(stereographic_points)
+    
+    # Iterate through the points and find infinite/invalid points
+    max_radius = 1
+    really_far_idxs = Int[]
+    for (i, point) in enumerate(stereographic_points)
+        m = hypot(point[1], point[2])
+        if !isfinite(m) || m > 1e32
+            push!(really_far_idxs, i)
+        elseif m > max_radius
+            max_radius = m
+        end
+    end
+    @debug max_radius length(really_far_idxs)
+    far_value = 1e6 * sqrt(max_radius)
+
+    if !isempty(really_far_idxs)
+        triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
+    end
+
+    boundary_points = reverse([
+        (-far_value, -far_value),
+        (-far_value, far_value),
+        (far_value, far_value),
+        (far_value, -far_value),
+        (-far_value, -far_value),
+    ])
+
+    # boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+
+    # triangulation = DelTri.triangulate(triangulation_points; boundary_nodes)
+    triangulation = DelTri.triangulate(triangulation_points)
+    #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
+
+    # First, get all the "solid" faces, ie, faces not attached to boundary nodes
+    boundary_faces = Iterators.filter(Base.Fix1(DelTri.is_boundary_triangle, triangulation), DelTri.each_solid_triangle(triangulation)) |> collect
+    faces = map(facetype, DelTri.each_solid_triangle(triangulation))
+
+    for ghost_face in DelTri.each_ghost_triangle(triangulation)
+        push!(faces, facetype(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
+    end
+
+    return faces
+end
+
 # These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
 
+delaunay_triangulate_spherical(points)
+
 pivot_ind = findfirst(isfinite, points)
 
 # debug
@@ -53,7 +111,7 @@ end
 far_value = 1e6 * sqrt(max_radius)
 
 if !isempty(really_far_idxs)
-    triangulation_points[really_far_idxs] .= Point2(far_value, 0)
+    triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
 end
 
 boundary_points = reverse([
@@ -63,7 +121,6 @@ boundary_points = reverse([
     (far_value, -far_value),
     (-far_value, -far_value),
 ])
-triangulation_points = copy(stereographic_points)
 
 # boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
 

From 1bcba2ccb4e34d5ed86a8bde065c719664ca550c Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 23:27:19 +0530
Subject: [PATCH 06/29] Plot correctly

---
 .../spherical_delaunay_stereographic.jl            | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 5d4512db9..516bcf54c 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -5,14 +5,14 @@ This is the approach which d3-geo-voronoi and friends use.  The alternative is S
 The 3D approach is basically that the 3d convex hull of a set of points on the sphere is its Delaunay triangulation.
 =#
 
-import GeometryOps as GO, GeoInterface as GI
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
 import Proj # for easy stereographic projection - TODO implement in Julia
 import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
 import CoordinateTransformations, Rotations
 
 using Downloads # does what it says on the tin
 using JSON3 # to load data
-using CairoMakie # for plotting
+using CairoMakie, GeoMakie # for plotting
 import Makie: Point3d
 
 function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
@@ -73,9 +73,13 @@ end
 
 # These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+faces = delaunay_triangulate_spherical(points)
 
-delaunay_triangulate_spherical(points)
+# This is the super-cool scrollable 3D globe (though it's a bit deformed... :D)
+f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
 
+# We can also replicate the observable notebook almost exactly (just missing ExactPredicates):
+f, a, p = Makie.mesh(points, faces; axis = (; type = GeoAxis, dest = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"), color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
 pivot_ind = findfirst(isfinite, points)
 
 # debug
@@ -85,14 +89,14 @@ point_colors[pivot_ind] = :red
 
 pivot_point = points[pivot_ind]
 necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
-
+#
 net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
 
 scatter(map(net_transformation_to_corrected_cartesian, points); color = point_colors)
 #
 stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
 scatter(stereographic_points; color = point_colors)
-
+#
 triangulation_points = copy(stereographic_points)
 
 

From 733faacfc94113db9eabba03a2643287d39b14c9 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 21 Jul 2024 23:27:34 +0530
Subject: [PATCH 07/29] First shot at hacking NaturalNeighbours interpolation

---
 .../spherical_delaunay_stereographic.jl       | 36 +++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 516bcf54c..857f49807 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -80,6 +80,42 @@ f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color =
 
 # We can also replicate the observable notebook almost exactly (just missing ExactPredicates):
 f, a, p = Makie.mesh(points, faces; axis = (; type = GeoAxis, dest = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"), color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
+# Whoops, this doesn't look so good!  Let's try to do this more "manually" instead.
+# We'll use NaturalNeighbours.jl for this, but we first have to reconstruct the triangulation
+# with the same faces, but in a Bertin projection...
+using NaturalNeighbours
+lonlat2bertin = Proj.Transformation(GFT.EPSG(4326), GFT.ProjString("+proj=bertin1953 +type=crs +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"); always_xy = true)
+lons = LinRange(-180, 180, 300)
+lats = LinRange(-90, 90, 150)
+bertin_points = lonlat2bertin.(lons, lats')
+
+projected_points = GO.reproject(GO.tuples(points), source_crs = GFT.EPSG(4326), target_crs = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09")
+
+ch = DelTri.convex_hull(projected_points) # assumes each point is in the triangulation
+boundary_nodes = DelTri.get_vertices(ch) 
+bertin_boundary_poly = GI.Polygon([GI.LineString(DelTri.get_points(ch)[DelTri.get_vertices(ch)])])
+
+tri = DelTri.Triangulation(projected_points, faces, boundary_nodes)
+itp = NaturalNeighbours.interpolate(tri, last.(points); derivatives = true)
+
+mat = [
+    if GO.contains(bertin_boundary_poly, (x, y))
+        itp(x, y; method = Nearest())
+    else
+        NaN
+    end
+    for (x, y) in bertin_points
+]
+# TODO: this currently doesn't work, because some points are not inside a triangle and so cannot be located.
+# Options are:
+# 1. Reject all points outside the convex hull of the projected points.  (Tried but failed)
+# 2. ...
+
+
+
+
+
+# Now, we proceed with the script implementation which has quite a bit of debug information + plots
 pivot_ind = findfirst(isfinite, points)
 
 # debug

From 6979ce421decb5dacf03f27b5df0cda05e877791 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 22 Jul 2024 00:20:50 +0530
Subject: [PATCH 08/29] Update
 docs/src/experiments/spherical_delaunay_stereographic.jl

Co-authored-by: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com>
---
 docs/src/experiments/spherical_delaunay_stereographic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 857f49807..5da83bb68 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -65,7 +65,7 @@ function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.Geom
     faces = map(facetype, DelTri.each_solid_triangle(triangulation))
 
     for ghost_face in DelTri.each_ghost_triangle(triangulation)
-        push!(faces, facetype(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
+        push!(faces, facetype(map(ghost_face) do i; DelaunayTriangulation.is_ghost_vertex(i) ? pivot_ind : i end))
     end
 
     return faces

From 5f1c6bb8ff1b7896b54d80080be9f63d6ac13311 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 22 Jul 2024 08:40:23 +0530
Subject: [PATCH 09/29] Remove degenerate faces

---
 .../spherical_delaunay_stereographic.jl       | 41 ++++++++++++++-----
 1 file changed, 30 insertions(+), 11 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 5da83bb68..d542f661f 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -61,11 +61,20 @@ function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.Geom
     #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
 
     # First, get all the "solid" faces, ie, faces not attached to boundary nodes
-    boundary_faces = Iterators.filter(Base.Fix1(DelTri.is_boundary_triangle, triangulation), DelTri.each_solid_triangle(triangulation)) |> collect
-    faces = map(facetype, DelTri.each_solid_triangle(triangulation))
-
+    original_triangles = collect(DelTri.each_solid_triangle(triangulation))
+    boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
+    faces = map(facetype, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
+    
+    for boundary_face in view(original_triangles, boundary_face_inds)
+        push!(faces, facetype(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
+    end
+    
     for ghost_face in DelTri.each_ghost_triangle(triangulation)
-        push!(faces, facetype(map(ghost_face) do i; DelaunayTriangulation.is_ghost_vertex(i) ? pivot_ind : i end))
+        push!(faces, facetype(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
+    end
+    # Remove degenerate triangles
+    filter!(faces) do face
+        !(face[1] == face[2] || face[2] == face[3] || face[1] == face[3])
     end
 
     return faces
@@ -100,7 +109,7 @@ itp = NaturalNeighbours.interpolate(tri, last.(points); derivatives = true)
 
 mat = [
     if GO.contains(bertin_boundary_poly, (x, y))
-        itp(x, y; method = Nearest())
+        itp(x, y; method = Nearest(), project = false)
     else
         NaN
     end
@@ -162,19 +171,29 @@ boundary_points = reverse([
     (-far_value, -far_value),
 ])
 
-# boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
 
-# triangulation = DelTri.triangulate(triangulation_points; boundary_nodes)
-triangulation = DelTri.triangulate(triangulation_points)
+triangulation = DelTri.triangulate(pts; boundary_nodes)
+# triangulation = DelTri.triangulate(triangulation_points)
 triplot(triangulation)
+DelTri.validate_triangulation(triangulation)
 #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
 
 # First, get all the "solid" faces, ie, faces not attached to boundary nodes
-boundary_faces = Iterators.filter(Base.Fix1(DelTri.is_boundary_triangle, triangulation), DelTri.each_solid_triangle(triangulation)) |> collect
-faces = map(CairoMakie.GeometryBasics.TriangleFace, DelTri.each_solid_triangle(triangulation))
+original_triangles = collect(DelTri.each_solid_triangle(triangulation))
+boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
+faces = map(CairoMakie.GeometryBasics.TriangleFace, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
+
+for boundary_face in view(original_triangles, boundary_face_inds)
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
+end
 
 for ghost_face in DelTri.each_ghost_triangle(triangulation)
-    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; i ≤ 0 ? pivot_ind : i end))
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
+end
+# Remove degenerate triangles
+filter!(faces) do face
+    !(face[1] == face[2] || face[2] == face[3] || face[1] == face[3])
 end
 
 wireframe(CairoMakie.GeometryBasics.normal_mesh(map(UnitCartesianFromGeographic(), points), faces))

From 72c395c18f208ca8f05be200963aa38fc567c118 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Tue, 23 Jul 2024 10:20:17 +0530
Subject: [PATCH 10/29] Use a bit more DelTri API

---
 docs/src/experiments/spherical_delaunay_stereographic.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index d542f661f..177d1c1e8 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -74,7 +74,7 @@ function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.Geom
     end
     # Remove degenerate triangles
     filter!(faces) do face
-        !(face[1] == face[2] || face[2] == face[3] || face[1] == face[3])
+        !(DelTri.geti(face) == DelTri.getj(face) || DelTri.getj(face) == DelTri.getk(face) || DelTri.geti(face) == DelTri.getk(face))
     end
 
     return faces
@@ -84,7 +84,9 @@ end
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
 faces = delaunay_triangulate_spherical(points)
 
-# This is the super-cool scrollable 3D globe (though it's a bit deformed... :D)
+
+
+# This is the super-cool scrollable 3D globe (though it's a bit deformed...)
 f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
 
 # We can also replicate the observable notebook almost exactly (just missing ExactPredicates):

From f1ffce7b677eee21a87003314caaa12075ba75db Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 26 Jul 2024 10:45:16 +0530
Subject: [PATCH 11/29] Add a Quickhull.jl method for 3d triangulation

---
 .../spherical_delaunay_stereographic.jl       | 23 ++++++++++++++++++-
 1 file changed, 22 insertions(+), 1 deletion(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 177d1c1e8..01fff728c 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -15,10 +15,18 @@ using JSON3 # to load data
 using CairoMakie, GeoMakie # for plotting
 import Makie: Point3d
 
-function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+abstract type SphericalTriangulationAlgorithm end
+struct StereographicDelaunayTriangulation <: SphericalTriangulationAlgorithm end
+struct SphericalConvexHull <: SphericalTriangulationAlgorithm end 
+
+spherical_triangulation(input_points; kwargs...) = spherical_triangulation(StereographicDelaunayTriangulation(), input_points; kwargs...)
+
+function spherical_triangulation(::StereographicDelaunayTriangulation, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
     # @assert GI.crstrait(first(input_points)) isa GI.AbstractGeographicCRSTrait
     points = GO.tuples(input_points)
+    # @assert points isa Vector{GI.Point}
 
+    # In 
     pivot_ind = findfirst(x -> all(isfinite, x), points)
     
     pivot_point = points[pivot_ind]
@@ -80,6 +88,19 @@ function delaunay_triangulate_spherical(input_points; facetype = CairoMakie.Geom
     return faces
 end
 
+function spherical_triangulation(::SphericalConvexHull, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+    points = GO.tuples(input_points) # we have to decompose the points into tuples, so they work with Quickhull.jl
+    # @assert points isa Vector{GI.Point}
+    cartesian_points = map(UnitCartesianFromGeographic(), points)
+    # The Delaunay triangulation of points on a sphere is simply the convex hull of those points in 3D Cartesian space.
+    # We can use e.g Quickhull.jl to get us such a convex hull.
+    hull = Quickhull.quickhull(cartesian_points)
+    # We return only the faces from these triangulation methods, so we simply map
+    # the facetype to the returned values from `Quickhull.facets`.
+    return map(facetype, Quickhull.facets(hull))
+end
+
+
 # These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
 faces = delaunay_triangulate_spherical(points)

From 41a3e2e398c51f036c7fa714ecab65c60698663d Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 26 Jul 2024 10:45:28 +0530
Subject: [PATCH 12/29] Move all coordinate transforms on top

---
 .../spherical_delaunay_stereographic.jl       | 109 +++++++++---------
 1 file changed, 56 insertions(+), 53 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 01fff728c..fde12bed9 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -101,6 +101,62 @@ function spherical_triangulation(::SphericalConvexHull, input_points; facetype =
 end
 
 
+# necessary coordinate transformations
+
+struct StereographicFromCartesian <: CoordinateTransformations.Transformation
+end
+
+function (::StereographicFromCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    # The Wikipedia definition has the north pole at infinity,
+    # this implementation has the south pole at infinity.
+    return Point2(x/(1-z), y/(1-z))
+end
+
+struct CartesianFromStereographic <: CoordinateTransformations.Transformation
+end
+
+function (::CartesianFromStereographic)(stereographic_point)
+    X, Y = stereographic_point
+    x2y2_1 = X^2 + Y^2 + 1
+    return Point3(2X/x2y2_1, 2Y/x2y2_1, (x2y2_1 - 2)/x2y2_1)
+end
+
+struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
+end
+
+function (::UnitCartesianFromGeographic)(geographic_point)
+    # Longitude is directly translatable to a spherical coordinate
+    # θ (azimuth)
+    θ = deg2rad(GI.x(geographic_point))
+    # The polar angle is 90 degrees minus the latitude
+    # ϕ (polar angle)
+    ϕ = deg2rad(90 - GI.y(geographic_point))
+    # Since this is the unit sphere, the radius is assumed to be 1,
+    # and we don't need to multiply by it.
+    return Point3(
+        sin(ϕ) * cos(θ),
+        sin(ϕ) * sin(θ),
+        cos(ϕ)
+    )
+end
+
+struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation 
+end
+
+function (::GeographicFromUnitCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    return Point2(
+        atan(y, x),
+        atan(hypot(x, y), z),
+    )
+end
+
+
+
+
 # These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
 faces = delaunay_triangulate_spherical(points)
@@ -263,59 +319,6 @@ _y = vec([y for _ in lons, y in lats])
 
 sibson_1_vals = itp(_x, _y; method=Sibson(1))
 
-# necessary coordinate transformations
-
-struct StereographicFromCartesian <: CoordinateTransformations.Transformation
-end
-
-function (::StereographicFromCartesian)(xyz::AbstractVector)
-    @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
-    x, y, z = xyz
-    # The Wikipedia definition has the north pole at infinity,
-    # this implementation has the south pole at infinity.
-    return Point2(x/(1-z), y/(1-z))
-end
-
-struct CartesianFromStereographic <: CoordinateTransformations.Transformation
-end
-
-function (::CartesianFromStereographic)(stereographic_point)
-    X, Y = stereographic_point
-    x2y2_1 = X^2 + Y^2 + 1
-    return Point3(2X/x2y2_1, 2Y/x2y2_1, (x2y2_1 - 2)/x2y2_1)
-end
-
-struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
-end
-
-function (::UnitCartesianFromGeographic)(geographic_point)
-    # Longitude is directly translatable to a spherical coordinate
-    # θ (azimuth)
-    θ = deg2rad(GI.x(geographic_point))
-    # The polar angle is 90 degrees minus the latitude
-    # ϕ (polar angle)
-    ϕ = deg2rad(90 - GI.y(geographic_point))
-    # Since this is the unit sphere, the radius is assumed to be 1,
-    # and we don't need to multiply by it.
-    return Point3(
-        sin(ϕ) * cos(θ),
-        sin(ϕ) * sin(θ),
-        cos(ϕ)
-    )
-end
-
-struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation 
-end
-
-function (::GeographicFromUnitCartesian)(xyz::AbstractVector)
-    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
-    x, y, z = xyz
-    return Point2(
-        atan(y, x),
-        atan(hypot(x, y), z),
-    )
-end
-
 transformed_points = GO.transform(points) do point
     # first, spherical to Cartesian
     longitude, latitude = deg2rad(GI.x(point)), deg2rad(GI.y(point))

From 893b647546f92552b3766c6ff3c72078e12c5697 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 26 Jul 2024 11:05:36 +0530
Subject: [PATCH 13/29] add another minimal line + test for Quickhull

---
 docs/src/experiments/spherical_delaunay_stereographic.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index fde12bed9..72d331e8c 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -160,8 +160,9 @@ end
 # These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
 faces = delaunay_triangulate_spherical(points)
-
-
+# or
+# faces = Quickhull.quickhull(map(UnitCartesianFromGeographic(), points)) |> Quickhull.faces |> collect
+# @assert isempty(setdiff(unique!(sort!(reduce(vcat, faces))), 1:length(points)))
 
 # This is the super-cool scrollable 3D globe (though it's a bit deformed...)
 f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)

From ecd799628b710fa6189126eca97fe0968d17c3cd Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 26 Jul 2024 11:15:40 +0530
Subject: [PATCH 14/29] Add license and more description

---
 .../spherical_delaunay_stereographic.jl       | 22 +++++++++++++++++--
 1 file changed, 20 insertions(+), 2 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index 72d331e8c..f13dbb33b 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -1,8 +1,26 @@
 #=
 # Spherical Delaunay triangulation using stereographic projections
 
-This is the approach which d3-geo-voronoi and friends use.  The alternative is STRIPACK which computes in 3D on the sphere.
-The 3D approach is basically that the 3d convex hull of a set of points on the sphere is its Delaunay triangulation.
+This file encodes two approaches to spherical Delaunay triangulation.
+
+1. `StereographicDelaunayTriangulation` projects the coordinates to a rotated stereographic projection, then computes the Delaunay triangulation on the plane.
+    This approach was taken from d3-geo-voronoi.
+    ```
+    Copyright 2018-2021 Philippe Rivière
+
+    Permission to use, copy, modify, and/or distribute this software for any purpose
+    with or without fee is hereby granted, provided that the above copyright notice
+    and this permission notice appear in all copies.
+
+    THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+    REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+    FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+    INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+    OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+    TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+    THIS SOFTWARE.
+    ```
+2. `SphericalConvexHull` computes the convex hull of the points in 3D Cartesian space, which is by definition the Delaunay triangulation.  This triangulation is currently done using the unreleased Quickhull.jl.
 =#
 
 import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT

From 7e831077759de53f30a92839b49f4b4cfd52cf1b Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 28 Jul 2024 22:18:53 +0530
Subject: [PATCH 15/29] First buggy implementation of Laplace interpolation

---
 .../spherical_delaunay_stereographic.jl       |   2 +-
 ...herical_natural_neighbour_interpolation.jl | 137 ++++++++++++++++++
 2 files changed, 138 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/experiments/spherical_natural_neighbour_interpolation.jl

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index f13dbb33b..d70c58161 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -177,7 +177,7 @@ end
 
 # These points are known to be good points, i.e., lon, lat, alt
 points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
-faces = delaunay_triangulate_spherical(points)
+faces = spherical_triangulation(points)
 # or
 # faces = Quickhull.quickhull(map(UnitCartesianFromGeographic(), points)) |> Quickhull.faces |> collect
 # @assert isempty(setdiff(unique!(sort!(reduce(vcat, faces))), 1:length(points)))
diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
new file mode 100644
index 000000000..636583068
--- /dev/null
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -0,0 +1,137 @@
+#=
+# Spherical natural neighbour interpolatoin
+=#
+
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+import Proj # for easy stereographic projection - TODO implement in Julia
+import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
+import CoordinateTransformations, Rotations
+
+using Downloads # does what it says on the tin
+using JSON3 # to load data
+using CairoMakie, GeoMakie # for plotting
+import Makie: Point3d
+
+# include(joinpath(@__DIR__, "spherical_delaunay_stereographic.jl"))
+
+using LinearAlgebra
+
+struct SphericalCap{T}
+    point::Point3{T}
+    radius::T
+end
+
+function SphericalCap(point::Point3{T1}, radius::T2) where {T1, T2}
+    return SphericalCap{promote_type(T1, T2)}(point, radius)
+end
+
+
+function circumcenter_on_unit_sphere(a, b, c)
+    LinearAlgebra.normalize(a × b + b × c + c × a)
+end
+
+spherical_distance(x::Point3, y::Point3) = acos(clamp(x ⋅ y, -1.0, 1.0))
+
+"Get the circumcenter of the triangle (a, b, c) on the unit sphere.  Returns a normalized 3-vector."
+function SphericalCap(a::Point3, b::Point3, c::Point3)
+    circumcenter = circumcenter_on_unit_sphere(a, b, c)
+    circumradius = spherical_distance(a, circumcenter)
+    return SphericalCap(circumcenter, circumradius)
+end
+
+function bowyer_watson_envelope!(applicable_cap_indices, query_point, points, faces, caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces)))
+    # brute force for now, but try the jump and search algorithm later
+    # can use e.g GeometryBasics.decompose(Point3{Float64}, GeometryBasics.Sphere(Point3(0.0), 1.0), 5) 
+    # to get starting points, or similar
+    empty!(applicable_cap_indices)
+    for (i, cap) in enumerate(caps)
+        if cap.radius > spherical_distance(query_point, cap.point)
+            push!(applicable_cap_indices, i)
+        end
+    end
+    # Now that we have the face indices, we need to get the applicable points
+    applicable_points = Int64[]
+    for i in applicable_cap_indices
+        current_face = faces[i]
+        edge_reoccurs = false
+        for current_edge in ((current_face[1], current_face[2]), (current_face[2], current_face[3]), (current_face[3], current_face[1]))
+            for j in applicable_cap_indices
+                if j == i
+                    continue # can't compare a triangle to itself
+                end
+                face_to_compare = faces[j]
+                for edge_to_compare in ((face_to_compare[1], face_to_compare[2]), (face_to_compare[2], face_to_compare[3]), (face_to_compare[3], face_to_compare[1]))
+                    if edge_to_compare == current_edge || reverse(edge_to_compare) == current_edge
+                        edge_reoccurs = true
+                        break
+                    end
+                end
+            end
+            if !edge_reoccurs # edge is unique
+                push!(applicable_points, current_edge[1])
+                push!(applicable_points, current_edge[2])
+            end
+        end
+    end
+    return unique!(applicable_points)
+end
+
+
+
+
+
+
+
+# These points are known to be good points, i.e., lon, lat, alt
+geographic_points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+z_values = last.(geographic_points)
+faces = spherical_triangulation(geographic_points)
+# correct the faces, since the order seems to be off
+faces = reverse.(faces)
+
+unique!(sort!(reduce(vcat, faces))) # so how am I getting this index?
+
+cartesian_points = UnitCartesianFromGeographic().(geographic_points)
+
+caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces))
+
+lons = -180.0:180.0
+lats = -90.0:90.0
+
+eval_laplace_coords(cartesian_points, z_values, faces, Point3(1.0, 0.0, 0.0))
+
+values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do point
+    eval_laplace_coords(cartesian_points, z_values, faces, point)
+end
+
+heatmap(lons, lats, values)
+
+# diagnostics
+# f, a, p = scatter(reduce(vcat, (view(cartesian_points, face) for face in view(faces, neighbour_inds))))
+# scatter!(query_point; color = :red, markersize = 40)
+import NaturalNeighbours: previndex_circular, nextindex_circular
+function laplace_ratio(points, envelope, i #= current vertex index =#, interpolation_point)
+    u = envelope[i]
+    prev_u = envelope[previndex_circular(envelope, i)]
+    next_u = envelope[nextindex_circular(envelope, i)]
+    p = points[u]
+    q, r = points[prev_u], points[next_u]
+    g1 = circumcenter_on_unit_sphere(q, p, interpolation_point)
+    g2 = circumcenter_on_unit_sphere(p, r, interpolation_point)
+    ℓ = spherical_distance(g1, g2)
+    d = spherical_distance(p, interpolation_point)
+    w = ℓ / d
+    return w, u, prev_u, next_u
+end
+
+function eval_laplace_coords(points, zs, faces, interpolation_point; envelope = Int64[])
+    envelope = bowyer_watson_envelope!(envelope, interpolation_point, points, faces, caps)
+    weighted_sum = 0.0
+    weight = 0.0
+    for i in eachindex(envelope)
+        w, u, prev_u, next_u = laplace_ratio(points, envelope, i, interpolation_point)
+        weighted_sum += w * zs[i]
+        weight += w
+    end
+    return weighted_sum / weight
+end
\ No newline at end of file

From 937e66fd0828fc552aa7d12aa14d21acd993929c Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 28 Jul 2024 22:20:02 +0530
Subject: [PATCH 16/29] Fix incorrect z indexing

---
 .../experiments/spherical_natural_neighbour_interpolation.jl    | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 636583068..5b51abf9b 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -130,7 +130,7 @@ function eval_laplace_coords(points, zs, faces, interpolation_point; envelope =
     weight = 0.0
     for i in eachindex(envelope)
         w, u, prev_u, next_u = laplace_ratio(points, envelope, i, interpolation_point)
-        weighted_sum += w * zs[i]
+        weighted_sum += w * zs[envelope[i]]
         weight += w
     end
     return weighted_sum / weight

From 0a7a4e609ce9f46f01db539083a95a07b6fe22b1 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 28 Jul 2024 22:35:06 +0530
Subject: [PATCH 17/29] More correctness

---
 ...herical_natural_neighbour_interpolation.jl | 88 ++++++++++++-------
 1 file changed, 55 insertions(+), 33 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 5b51abf9b..9b2af4a31 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -39,7 +39,7 @@ function SphericalCap(a::Point3, b::Point3, c::Point3)
     return SphericalCap(circumcenter, circumradius)
 end
 
-function bowyer_watson_envelope!(applicable_cap_indices, query_point, points, faces, caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces)))
+function bowyer_watson_envelope!(applicable_points, query_point, points, faces, caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces)); applicable_cap_indices = Int64[])
     # brute force for now, but try the jump and search algorithm later
     # can use e.g GeometryBasics.decompose(Point3{Float64}, GeometryBasics.Sphere(Point3(0.0), 1.0), 5) 
     # to get starting points, or similar
@@ -77,6 +77,48 @@ function bowyer_watson_envelope!(applicable_cap_indices, query_point, points, fa
 end
 
 
+import NaturalNeighbours: previndex_circular, nextindex_circular
+function laplace_ratio(points, envelope, i #= current vertex index =#, interpolation_point)
+    u = envelope[i]
+    prev_u = envelope[previndex_circular(envelope, i)]
+    next_u = envelope[nextindex_circular(envelope, i)]
+    p = points[u]
+    q, r = points[prev_u], points[next_u]
+    g1 = circumcenter_on_unit_sphere(q, p, interpolation_point)
+    g2 = circumcenter_on_unit_sphere(p, r, interpolation_point)
+    ℓ = spherical_distance(g1, g2)
+    d = spherical_distance(p, interpolation_point)
+    w = ℓ / d
+    return w, u, prev_u, next_u
+end
+
+struct NaturalCoordinates{F, I}
+    coordinates::Vector{F}
+    indices::Vector{I}
+    interpolation_point::Point3{Float64}
+end
+
+function laplace_nearest_neighbour_coords(points, faces, interpolation_point; envelope = Int64[], cap_idxs = Int64[])
+    envelope = bowyer_watson_envelope!(envelope, interpolation_point, points, faces, caps; applicable_cap_indices = cap_idxs)
+    weighted_sum = 0.0
+    weight = 0.0
+    coords = NaturalCoordinates(Float64[], Int64[], interpolation_point)
+    for i in eachindex(envelope)
+        w, u, prev_u, next_u = laplace_ratio(points, envelope, i, interpolation_point)
+        push!(coords.coordinates, w)
+        push!(coords.indices, u)
+    end
+    coords.coordinates ./= sum(coords.coordinates)
+    return coords
+end
+
+function eval_laplace_coordinates(points, faces, zs, interpolation_point)
+    coords = laplace_nearest_neighbour_coords(points, faces, interpolation_point)
+    if isempty(coords.coordinates)
+        return NaN
+    end
+    return sum((coord * zs[idx] for (coord, idx) in zip(coords.coordinates, coords.indices)))
+end
 
 
 
@@ -95,43 +137,23 @@ cartesian_points = UnitCartesianFromGeographic().(geographic_points)
 
 caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces))
 
-lons = -180.0:180.0
-lats = -90.0:90.0
+lons = -180.0:0.5:180.0
+lats = -90.0:0.5:90.0
 
-eval_laplace_coords(cartesian_points, z_values, faces, Point3(1.0, 0.0, 0.0))
+eval_laplace_coordinates(cartesian_points, faces, z_values, Point3(1.0, 0.0, 0.0))
 
 values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do point
-    eval_laplace_coords(cartesian_points, z_values, faces, point)
+    eval_laplace_coordinates(cartesian_points, faces, z_values, point)
 end
 
-heatmap(lons, lats, values)
+heatmap(lons, lats, values; axis = (; aspect = DataAspect()))
 
+f = Figure();
+a = LScene(f[1, 1])
+p = meshimage!(a, lons, lats, rotl90(values))
+p.transformation.transform_func[] = Makie.PointTrans{3}(UnitCartesianFromGeographic())
+scatter!(a, cartesian_points)
+f # not entirely sure what's going on here
 # diagnostics
 # f, a, p = scatter(reduce(vcat, (view(cartesian_points, face) for face in view(faces, neighbour_inds))))
-# scatter!(query_point; color = :red, markersize = 40)
-import NaturalNeighbours: previndex_circular, nextindex_circular
-function laplace_ratio(points, envelope, i #= current vertex index =#, interpolation_point)
-    u = envelope[i]
-    prev_u = envelope[previndex_circular(envelope, i)]
-    next_u = envelope[nextindex_circular(envelope, i)]
-    p = points[u]
-    q, r = points[prev_u], points[next_u]
-    g1 = circumcenter_on_unit_sphere(q, p, interpolation_point)
-    g2 = circumcenter_on_unit_sphere(p, r, interpolation_point)
-    ℓ = spherical_distance(g1, g2)
-    d = spherical_distance(p, interpolation_point)
-    w = ℓ / d
-    return w, u, prev_u, next_u
-end
-
-function eval_laplace_coords(points, zs, faces, interpolation_point; envelope = Int64[])
-    envelope = bowyer_watson_envelope!(envelope, interpolation_point, points, faces, caps)
-    weighted_sum = 0.0
-    weight = 0.0
-    for i in eachindex(envelope)
-        w, u, prev_u, next_u = laplace_ratio(points, envelope, i, interpolation_point)
-        weighted_sum += w * zs[envelope[i]]
-        weight += w
-    end
-    return weighted_sum / weight
-end
\ No newline at end of file
+# scatter!(query_point; color = :red, markersize = 40)
\ No newline at end of file

From 37190f7bd7e63743c0bf20adec68462585ee3f21 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 28 Jul 2024 23:28:53 +0530
Subject: [PATCH 18/29] Experiment: why is my Bowyer-Watson wrong

---
 ...herical_natural_neighbour_interpolation.jl | 39 +++++++++++++++----
 1 file changed, 32 insertions(+), 7 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 9b2af4a31..26b048d51 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -15,6 +15,7 @@ import Makie: Point3d
 # include(joinpath(@__DIR__, "spherical_delaunay_stereographic.jl"))
 
 using LinearAlgebra
+using GeometryBasics
 
 struct SphericalCap{T}
     point::Point3{T}
@@ -45,12 +46,12 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
     # to get starting points, or similar
     empty!(applicable_cap_indices)
     for (i, cap) in enumerate(caps)
-        if cap.radius > spherical_distance(query_point, cap.point)
+        if cap.radius ≥ spherical_distance(query_point, cap.point)
             push!(applicable_cap_indices, i)
         end
     end
     # Now that we have the face indices, we need to get the applicable points
-    applicable_points = Int64[]
+    empty!(applicable_points)
     for i in applicable_cap_indices
         current_face = faces[i]
         edge_reoccurs = false
@@ -73,6 +74,15 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
             end
         end
     end
+    # Start at point 1, find the first occurrence of point 1 in the applicable_points list.
+    # This is the last point of the edge coming from point 1.
+    # Now, swap the element before that with point 3.  Then continue on doing this.
+    # for (i, point_idx) in enumerate(applicable_points)
+    #     if i % 2 == 0
+    #         continue
+    #     end
+    #     applicable_points[i] = findfirst(==(applicable_points[i+1]), points)
+    # end
     return unique!(applicable_points)
 end
 
@@ -100,8 +110,6 @@ end
 
 function laplace_nearest_neighbour_coords(points, faces, interpolation_point; envelope = Int64[], cap_idxs = Int64[])
     envelope = bowyer_watson_envelope!(envelope, interpolation_point, points, faces, caps; applicable_cap_indices = cap_idxs)
-    weighted_sum = 0.0
-    weight = 0.0
     coords = NaturalCoordinates(Float64[], Int64[], interpolation_point)
     for i in eachindex(envelope)
         w, u, prev_u, next_u = laplace_ratio(points, envelope, i, interpolation_point)
@@ -150,10 +158,27 @@ heatmap(lons, lats, values; axis = (; aspect = DataAspect()))
 
 f = Figure();
 a = LScene(f[1, 1])
-p = meshimage!(a, lons, lats, rotl90(values))
+p = meshimage!(a, lons, lats, rotl90(values); npoints = (720, 360))
 p.transformation.transform_func[] = Makie.PointTrans{3}(UnitCartesianFromGeographic())
-scatter!(a, cartesian_points)
+# scatter!(a, cartesian_points)
 f # not entirely sure what's going on here
 # diagnostics
 # f, a, p = scatter(reduce(vcat, (view(cartesian_points, face) for face in view(faces, neighbour_inds))))
-# scatter!(query_point; color = :red, markersize = 40)
\ No newline at end of file
+# scatter!(query_point; color = :red, markersize = 40)
+
+query_point = LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0))
+pt_inds = bowyer_watson_envelope!(Int64[], query_point, cartesian_points, faces, caps)
+
+f, a, p = scatter([query_point]; markersize = 30, color = :green, axis = (; type = LScene));
+scatter!(a, cartesian_points)
+scatter!(a, view(cartesian_points, pt_inds); color = eachindex(pt_inds), colormap = :turbo, markersize = 20)
+wireframe!(a, GeometryBasics.Mesh(cartesian_points, faces); alpha = 0.3)
+f
+
+function Makie.convert_arguments(::Type{Makie.Mesh}, cap::SphericalCap)
+    offset_point = LinearAlgebra.normalize(cap.point + LinearAlgebra.normalize(Point3(cos(cap.radius), sin(cap.radius), 0.0)))
+    points = [cap.point + Rotations.AngleAxis(θ, cap.point...) * offset_point for θ in LinRange(0, 2π, 20)]
+    push!(points, cap.point)
+    faces = [GeometryBasics.TriangleFace(i, i+1, 21) for i in 1:19]
+    return (GeometryBasics.normal_mesh(points, faces),)
+end

From 7a3667e8ce11a16e78f64002c7cec40edb87b9b0 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sun, 28 Jul 2024 23:53:15 +0530
Subject: [PATCH 19/29] Reset `edge_reoccurs`

---
 .../spherical_natural_neighbour_interpolation.jl      | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 26b048d51..9290d4b85 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -52,6 +52,9 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
     end
     # Now that we have the face indices, we need to get the applicable points
     empty!(applicable_points)
+    # for face_idx in applicable_cap_indices
+    #     append!(applicable_points, faces[face_idx])
+    # end
     for i in applicable_cap_indices
         current_face = faces[i]
         edge_reoccurs = false
@@ -72,6 +75,7 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
                 push!(applicable_points, current_edge[1])
                 push!(applicable_points, current_edge[2])
             end
+            edge_reoccurs = false
         end
     end
     # Start at point 1, find the first occurrence of point 1 in the applicable_points list.
@@ -120,10 +124,10 @@ function laplace_nearest_neighbour_coords(points, faces, interpolation_point; en
     return coords
 end
 
-function eval_laplace_coordinates(points, faces, zs, interpolation_point)
-    coords = laplace_nearest_neighbour_coords(points, faces, interpolation_point)
+function eval_laplace_coordinates(points, faces, zs, interpolation_point; envelope = Int64[], cap_idxs = Int64[])
+    coords = laplace_nearest_neighbour_coords(points, faces, interpolation_point; envelope, cap_idxs)
     if isempty(coords.coordinates)
-        return NaN
+        return Inf
     end
     return sum((coord * zs[idx] for (coord, idx) in zip(coords.coordinates, coords.indices)))
 end
@@ -173,6 +177,7 @@ f, a, p = scatter([query_point]; markersize = 30, color = :green, axis = (; type
 scatter!(a, cartesian_points)
 scatter!(a, view(cartesian_points, pt_inds); color = eachindex(pt_inds), colormap = :turbo, markersize = 20)
 wireframe!(a, GeometryBasics.Mesh(cartesian_points, faces); alpha = 0.3)
+
 f
 
 function Makie.convert_arguments(::Type{Makie.Mesh}, cap::SphericalCap)

From bd608e7fe47677cfcd7a6cdd316a8942faa851be Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 29 Jul 2024 00:47:32 +0530
Subject: [PATCH 20/29] Fix ambiguity

---
 docs/src/experiments/spherical_delaunay_stereographic.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
index d70c58161..89201453c 100644
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ b/docs/src/experiments/spherical_delaunay_stereographic.jl
@@ -106,7 +106,7 @@ function spherical_triangulation(::StereographicDelaunayTriangulation, input_poi
     return faces
 end
 
-function spherical_triangulation(::SphericalConvexHull, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+function spherical_triangulation(::SphericalConvexHull, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace{Int32})
     points = GO.tuples(input_points) # we have to decompose the points into tuples, so they work with Quickhull.jl
     # @assert points isa Vector{GI.Point}
     cartesian_points = map(UnitCartesianFromGeographic(), points)

From e638b1763854529902eebff8902e11768f196f38 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 29 Jul 2024 00:47:46 +0530
Subject: [PATCH 21/29] Correct the winding order

---
 ...herical_natural_neighbour_interpolation.jl | 28 +++++++++++++++++--
 1 file changed, 26 insertions(+), 2 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 9290d4b85..0085574ac 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -40,6 +40,23 @@ function SphericalCap(a::Point3, b::Point3, c::Point3)
     return SphericalCap(circumcenter, circumradius)
 end
 
+function _is_ccw_unit_sphere(v_0,v_c,v_i)
+    # checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW
+    return(LinearAlgebra.dot(LinearAlgebra.cross(v_c - v_0,v_i - v_c), v_i) < 0)
+end
+
+function angle_between(a, b, c)
+    ab = b - a
+    bc = c - b
+    norm_dot = (ab ⋅ bc) / (LinearAlgebra.norm(ab) * LinearAlgebra.norm(bc))
+    angle =  acos(clamp(norm_dot, -1.0, 1.0))
+    if _is_ccw_unit_sphere(a, b, c)
+        return angle
+    else
+        return 2π - angle
+    end
+end
+
 function bowyer_watson_envelope!(applicable_points, query_point, points, faces, caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces)); applicable_cap_indices = Int64[])
     # brute force for now, but try the jump and search algorithm later
     # can use e.g GeometryBasics.decompose(Point3{Float64}, GeometryBasics.Sphere(Point3(0.0), 1.0), 5) 
@@ -78,6 +95,11 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
             edge_reoccurs = false
         end
     end
+    unique!(applicable_points)
+    # Find a plane that can house most points OK
+    # since natural neighbours are local ish (probably wont be too far) 
+    # we can get away with this
+
     # Start at point 1, find the first occurrence of point 1 in the applicable_points list.
     # This is the last point of the edge coming from point 1.
     # Now, swap the element before that with point 3.  Then continue on doing this.
@@ -87,10 +109,12 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
     #     end
     #     applicable_points[i] = findfirst(==(applicable_points[i+1]), points)
     # end
-    return unique!(applicable_points)
+    three_d_points = view(points, applicable_points)
+    angles = [angle_between(three_d_points[1], query_point, point) for point in three_d_points]
+    pt_inds = sortperm(angles)
+    return applicable_points[pt_inds]
 end
 
-
 import NaturalNeighbours: previndex_circular, nextindex_circular
 function laplace_ratio(points, envelope, i #= current vertex index =#, interpolation_point)
     u = envelope[i]

From 4de0db76b14a3c6e03b4ebdc6e5986e230d3fd93 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 29 Jul 2024 08:07:17 +0530
Subject: [PATCH 22/29] Refactor files into a core spherical triangulation file
 + experiments

---
 docs/Project.toml                             |   1 +
 docs/src/experiments/spherical_delaunay.jl    | 215 ++++++++++
 .../spherical_delaunay_stereographic.jl       | 396 ------------------
 .../spherical_triangulation_experiments.jl    | 180 ++++++++
 4 files changed, 396 insertions(+), 396 deletions(-)
 create mode 100644 docs/src/experiments/spherical_delaunay.jl
 delete mode 100644 docs/src/experiments/spherical_delaunay_stereographic.jl
 create mode 100644 docs/src/experiments/spherical_triangulation_experiments.jl

diff --git a/docs/Project.toml b/docs/Project.toml
index fccffcc39..1cc496415 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -37,6 +37,7 @@ MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
+Quickhull = "baca2653-9c88-49d2-a488-12dbf25fc4fb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/docs/src/experiments/spherical_delaunay.jl b/docs/src/experiments/spherical_delaunay.jl
new file mode 100644
index 000000000..10e9ca964
--- /dev/null
+++ b/docs/src/experiments/spherical_delaunay.jl
@@ -0,0 +1,215 @@
+#=
+# Spherical Delaunay triangulation using stereographic projections
+
+This file encodes two approaches to spherical Delaunay triangulation.
+
+1. `StereographicDelaunayTriangulation` projects the coordinates to a rotated stereographic projection, then computes the Delaunay triangulation on the plane.
+    This approach was taken from d3-geo-voronoi.
+    ```
+    Copyright 2018-2021 Philippe Rivière
+
+    Permission to use, copy, modify, and/or distribute this software for any purpose
+    with or without fee is hereby granted, provided that the above copyright notice
+    and this permission notice appear in all copies.
+
+    THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+    REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+    FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+    INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+    OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+    TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+    THIS SOFTWARE.
+    ```
+2. `SphericalConvexHull` computes the convex hull of the points in 3D Cartesian space, which is by definition the Delaunay triangulation.  This triangulation is currently done using the unreleased Quickhull.jl.
+=#
+
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+import Proj # for easy stereographic projection - TODO implement in Julia
+import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
+import Quickhull # convex hulls in d+1 are Delaunay triangulations in dimension d
+import CoordinateTransformations, Rotations
+
+using Downloads # does what it says on the tin
+using JSON3 # to load data
+using CairoMakie, GeoMakie # for plotting
+import Makie: Point3d
+
+abstract type SphericalTriangulationAlgorithm end
+struct StereographicDelaunayTriangulation <: SphericalTriangulationAlgorithm end
+struct SphericalConvexHull <: SphericalTriangulationAlgorithm end 
+
+spherical_triangulation(input_points; kwargs...) = spherical_triangulation(StereographicDelaunayTriangulation(), input_points; kwargs...)
+
+function spherical_triangulation(::StereographicDelaunayTriangulation, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
+    # @assert GI.crstrait(first(input_points)) isa GI.AbstractGeographicCRSTrait
+    points = GO.tuples(input_points)
+    # @assert points isa Vector{GI.Point}
+
+    # In 
+    pivot_ind = findfirst(x -> all(isfinite, x), points)
+    
+    pivot_point = points[pivot_ind]
+    necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
+
+    net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
+    stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+    triangulation_points = copy(stereographic_points)
+    
+    # Iterate through the points and find infinite/invalid points
+    max_radius = 1
+    really_far_idxs = Int[]
+    for (i, point) in enumerate(stereographic_points)
+        m = hypot(point[1], point[2])
+        if !isfinite(m) || m > 1e32
+            push!(really_far_idxs, i)
+        elseif m > max_radius
+            max_radius = m
+        end
+    end
+    @debug max_radius length(really_far_idxs)
+    far_value = 1e6 * sqrt(max_radius)
+
+    if !isempty(really_far_idxs)
+        triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
+    end
+
+    boundary_points = reverse([
+        (-far_value, -far_value),
+        (-far_value, far_value),
+        (far_value, far_value),
+        (far_value, -far_value),
+        (-far_value, -far_value),
+    ])
+
+    # boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+
+    # triangulation = DelTri.triangulate(triangulation_points; boundary_nodes)
+    triangulation = DelTri.triangulate(triangulation_points)
+    #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
+
+    # First, get all the "solid" faces, ie, faces not attached to boundary nodes
+    original_triangles = collect(DelTri.each_solid_triangle(triangulation))
+    boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
+    faces = map(facetype, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
+    
+    for boundary_face in view(original_triangles, boundary_face_inds)
+        push!(faces, facetype(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
+    end
+    
+    for ghost_face in DelTri.each_ghost_triangle(triangulation)
+        push!(faces, facetype(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
+    end
+    # Remove degenerate triangles
+    filter!(faces) do face
+        !(DelTri.geti(face) == DelTri.getj(face) || DelTri.getj(face) == DelTri.getk(face) || DelTri.geti(face) == DelTri.getk(face))
+    end
+
+    return faces
+end
+
+function spherical_triangulation(::SphericalConvexHull, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace{Int32})
+    points = GO.tuples(input_points) # we have to decompose the points into tuples, so they work with Quickhull.jl
+    # @assert points isa Vector{GI.Point}
+    cartesian_points = map(UnitCartesianFromGeographic(), points)
+    # The Delaunay triangulation of points on a sphere is simply the convex hull of those points in 3D Cartesian space.
+    # We can use e.g Quickhull.jl to get us such a convex hull.
+    hull = Quickhull.quickhull(cartesian_points)
+    # We return only the faces from these triangulation methods, so we simply map
+    # the facetype to the returned values from `Quickhull.facets`.
+    return map(facetype, Quickhull.facets(hull))
+end
+
+
+# necessary coordinate transformations
+
+struct StereographicFromCartesian <: CoordinateTransformations.Transformation
+end
+
+function (::StereographicFromCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    # The Wikipedia definition has the north pole at infinity,
+    # this implementation has the south pole at infinity.
+    return Point2(x/(1-z), y/(1-z))
+end
+
+struct CartesianFromStereographic <: CoordinateTransformations.Transformation
+end
+
+function (::CartesianFromStereographic)(stereographic_point)
+    X, Y = stereographic_point
+    x2y2_1 = X^2 + Y^2 + 1
+    return Point3(2X/x2y2_1, 2Y/x2y2_1, (x2y2_1 - 2)/x2y2_1)
+end
+
+struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
+end
+
+function (::UnitCartesianFromGeographic)(geographic_point)
+    # Longitude is directly translatable to a spherical coordinate
+    # θ (azimuth)
+    θ = deg2rad(GI.x(geographic_point))
+    # The polar angle is 90 degrees minus the latitude
+    # ϕ (polar angle)
+    ϕ = deg2rad(90 - GI.y(geographic_point))
+    # Since this is the unit sphere, the radius is assumed to be 1,
+    # and we don't need to multiply by it.
+    return Point3(
+        sin(ϕ) * cos(θ),
+        sin(ϕ) * sin(θ),
+        cos(ϕ)
+    )
+end
+
+struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation 
+end
+
+function (::GeographicFromUnitCartesian)(xyz::AbstractVector)
+    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
+    x, y, z = xyz
+    return Point2(
+        atan(y, x),
+        atan(hypot(x, y), z),
+    )
+end
+
+function randsphericalangles(n)
+    θ = 2π .* rand(n)
+    ϕ = acos.(2 .* rand(n) .- 1)
+    return Point2.(θ, ϕ)
+end
+
+function randsphere(n)
+    θϕ = randsphericalangles(n)
+    return Point3.(
+        sin.(last.(θϕs)) .* cos.(first.(θϕs)),
+        sin.(last.(θϕs)) .* sin.(first.(θϕs)),
+        cos.(last.(θϕs))
+    )
+end
+
+# θϕs = randsphericalangles(50)
+# θs, ϕs = first.(θϕs), last.(θϕs)
+# pts = Point3.(
+#     sin.(ϕs) .* cos.(θs),
+#     sin.(ϕs) .* sin.(θs),
+#     cos.(ϕs)
+# )
+
+# f, a, p = scatter(pts; color = [fill(:blue, 49); :red])
+
+# function Makie.rotate!(t::Makie.Transformable, rotation::Rotations.Rotation)
+#     quat = Rotations.QuatRotation(rotation)
+#     Makie.rotate!(Makie.Absolute, t, Makie.Quaternion(quat.q.v1, quat.q.v2, quat.q.v3, quat.q.s))
+# end
+
+
+# rotate!(p, Rotations.RotX(π/2))
+# rotate!(p, Rotations.RotX(0))
+
+# pivot_point = θϕs[end]
+# necessary_rotation = Rotations.RotY(pi) * Rotations.RotY(-pivot_point[2]) * Rotations.RotZ(-pivot_point[1])
+
+# rotate!(p, necessary_rotation)
+
+# f
\ No newline at end of file
diff --git a/docs/src/experiments/spherical_delaunay_stereographic.jl b/docs/src/experiments/spherical_delaunay_stereographic.jl
deleted file mode 100644
index 89201453c..000000000
--- a/docs/src/experiments/spherical_delaunay_stereographic.jl
+++ /dev/null
@@ -1,396 +0,0 @@
-#=
-# Spherical Delaunay triangulation using stereographic projections
-
-This file encodes two approaches to spherical Delaunay triangulation.
-
-1. `StereographicDelaunayTriangulation` projects the coordinates to a rotated stereographic projection, then computes the Delaunay triangulation on the plane.
-    This approach was taken from d3-geo-voronoi.
-    ```
-    Copyright 2018-2021 Philippe Rivière
-
-    Permission to use, copy, modify, and/or distribute this software for any purpose
-    with or without fee is hereby granted, provided that the above copyright notice
-    and this permission notice appear in all copies.
-
-    THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
-    REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
-    FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
-    INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
-    OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
-    TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-    THIS SOFTWARE.
-    ```
-2. `SphericalConvexHull` computes the convex hull of the points in 3D Cartesian space, which is by definition the Delaunay triangulation.  This triangulation is currently done using the unreleased Quickhull.jl.
-=#
-
-import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
-import Proj # for easy stereographic projection - TODO implement in Julia
-import DelaunayTriangulation as DelTri # Delaunay triangulation on the 2d plane
-import CoordinateTransformations, Rotations
-
-using Downloads # does what it says on the tin
-using JSON3 # to load data
-using CairoMakie, GeoMakie # for plotting
-import Makie: Point3d
-
-abstract type SphericalTriangulationAlgorithm end
-struct StereographicDelaunayTriangulation <: SphericalTriangulationAlgorithm end
-struct SphericalConvexHull <: SphericalTriangulationAlgorithm end 
-
-spherical_triangulation(input_points; kwargs...) = spherical_triangulation(StereographicDelaunayTriangulation(), input_points; kwargs...)
-
-function spherical_triangulation(::StereographicDelaunayTriangulation, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace)
-    # @assert GI.crstrait(first(input_points)) isa GI.AbstractGeographicCRSTrait
-    points = GO.tuples(input_points)
-    # @assert points isa Vector{GI.Point}
-
-    # In 
-    pivot_ind = findfirst(x -> all(isfinite, x), points)
-    
-    pivot_point = points[pivot_ind]
-    necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
-
-    net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
-    stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
-    triangulation_points = copy(stereographic_points)
-    
-    # Iterate through the points and find infinite/invalid points
-    max_radius = 1
-    really_far_idxs = Int[]
-    for (i, point) in enumerate(stereographic_points)
-        m = hypot(point[1], point[2])
-        if !isfinite(m) || m > 1e32
-            push!(really_far_idxs, i)
-        elseif m > max_radius
-            max_radius = m
-        end
-    end
-    @debug max_radius length(really_far_idxs)
-    far_value = 1e6 * sqrt(max_radius)
-
-    if !isempty(really_far_idxs)
-        triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
-    end
-
-    boundary_points = reverse([
-        (-far_value, -far_value),
-        (-far_value, far_value),
-        (far_value, far_value),
-        (far_value, -far_value),
-        (-far_value, -far_value),
-    ])
-
-    # boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
-
-    # triangulation = DelTri.triangulate(triangulation_points; boundary_nodes)
-    triangulation = DelTri.triangulate(triangulation_points)
-    #  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
-
-    # First, get all the "solid" faces, ie, faces not attached to boundary nodes
-    original_triangles = collect(DelTri.each_solid_triangle(triangulation))
-    boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
-    faces = map(facetype, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
-    
-    for boundary_face in view(original_triangles, boundary_face_inds)
-        push!(faces, facetype(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
-    end
-    
-    for ghost_face in DelTri.each_ghost_triangle(triangulation)
-        push!(faces, facetype(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
-    end
-    # Remove degenerate triangles
-    filter!(faces) do face
-        !(DelTri.geti(face) == DelTri.getj(face) || DelTri.getj(face) == DelTri.getk(face) || DelTri.geti(face) == DelTri.getk(face))
-    end
-
-    return faces
-end
-
-function spherical_triangulation(::SphericalConvexHull, input_points; facetype = CairoMakie.GeometryBasics.TriangleFace{Int32})
-    points = GO.tuples(input_points) # we have to decompose the points into tuples, so they work with Quickhull.jl
-    # @assert points isa Vector{GI.Point}
-    cartesian_points = map(UnitCartesianFromGeographic(), points)
-    # The Delaunay triangulation of points on a sphere is simply the convex hull of those points in 3D Cartesian space.
-    # We can use e.g Quickhull.jl to get us such a convex hull.
-    hull = Quickhull.quickhull(cartesian_points)
-    # We return only the faces from these triangulation methods, so we simply map
-    # the facetype to the returned values from `Quickhull.facets`.
-    return map(facetype, Quickhull.facets(hull))
-end
-
-
-# necessary coordinate transformations
-
-struct StereographicFromCartesian <: CoordinateTransformations.Transformation
-end
-
-function (::StereographicFromCartesian)(xyz::AbstractVector)
-    @assert length(xyz) == 3 "StereographicFromCartesian expects a 3D Cartesian vector"
-    x, y, z = xyz
-    # The Wikipedia definition has the north pole at infinity,
-    # this implementation has the south pole at infinity.
-    return Point2(x/(1-z), y/(1-z))
-end
-
-struct CartesianFromStereographic <: CoordinateTransformations.Transformation
-end
-
-function (::CartesianFromStereographic)(stereographic_point)
-    X, Y = stereographic_point
-    x2y2_1 = X^2 + Y^2 + 1
-    return Point3(2X/x2y2_1, 2Y/x2y2_1, (x2y2_1 - 2)/x2y2_1)
-end
-
-struct UnitCartesianFromGeographic <: CoordinateTransformations.Transformation 
-end
-
-function (::UnitCartesianFromGeographic)(geographic_point)
-    # Longitude is directly translatable to a spherical coordinate
-    # θ (azimuth)
-    θ = deg2rad(GI.x(geographic_point))
-    # The polar angle is 90 degrees minus the latitude
-    # ϕ (polar angle)
-    ϕ = deg2rad(90 - GI.y(geographic_point))
-    # Since this is the unit sphere, the radius is assumed to be 1,
-    # and we don't need to multiply by it.
-    return Point3(
-        sin(ϕ) * cos(θ),
-        sin(ϕ) * sin(θ),
-        cos(ϕ)
-    )
-end
-
-struct GeographicFromUnitCartesian <: CoordinateTransformations.Transformation 
-end
-
-function (::GeographicFromUnitCartesian)(xyz::AbstractVector)
-    @assert length(xyz) == 3 "GeographicFromUnitCartesian expects a 3D Cartesian vector"
-    x, y, z = xyz
-    return Point2(
-        atan(y, x),
-        atan(hypot(x, y), z),
-    )
-end
-
-
-
-
-# These points are known to be good points, i.e., lon, lat, alt
-points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
-faces = spherical_triangulation(points)
-# or
-# faces = Quickhull.quickhull(map(UnitCartesianFromGeographic(), points)) |> Quickhull.faces |> collect
-# @assert isempty(setdiff(unique!(sort!(reduce(vcat, faces))), 1:length(points)))
-
-# This is the super-cool scrollable 3D globe (though it's a bit deformed...)
-f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
-
-# We can also replicate the observable notebook almost exactly (just missing ExactPredicates):
-f, a, p = Makie.mesh(points, faces; axis = (; type = GeoAxis, dest = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"), color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
-# Whoops, this doesn't look so good!  Let's try to do this more "manually" instead.
-# We'll use NaturalNeighbours.jl for this, but we first have to reconstruct the triangulation
-# with the same faces, but in a Bertin projection...
-using NaturalNeighbours
-lonlat2bertin = Proj.Transformation(GFT.EPSG(4326), GFT.ProjString("+proj=bertin1953 +type=crs +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"); always_xy = true)
-lons = LinRange(-180, 180, 300)
-lats = LinRange(-90, 90, 150)
-bertin_points = lonlat2bertin.(lons, lats')
-
-projected_points = GO.reproject(GO.tuples(points), source_crs = GFT.EPSG(4326), target_crs = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09")
-
-ch = DelTri.convex_hull(projected_points) # assumes each point is in the triangulation
-boundary_nodes = DelTri.get_vertices(ch) 
-bertin_boundary_poly = GI.Polygon([GI.LineString(DelTri.get_points(ch)[DelTri.get_vertices(ch)])])
-
-tri = DelTri.Triangulation(projected_points, faces, boundary_nodes)
-itp = NaturalNeighbours.interpolate(tri, last.(points); derivatives = true)
-
-mat = [
-    if GO.contains(bertin_boundary_poly, (x, y))
-        itp(x, y; method = Nearest(), project = false)
-    else
-        NaN
-    end
-    for (x, y) in bertin_points
-]
-# TODO: this currently doesn't work, because some points are not inside a triangle and so cannot be located.
-# Options are:
-# 1. Reject all points outside the convex hull of the projected points.  (Tried but failed)
-# 2. ...
-
-
-
-
-
-# Now, we proceed with the script implementation which has quite a bit of debug information + plots
-pivot_ind = findfirst(isfinite, points)
-
-# debug
-point_colors = fill(:blue, length(points))
-point_colors[pivot_ind] = :red
-# end debug
-
-pivot_point = points[pivot_ind]
-necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
-#
-net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
-
-scatter(map(net_transformation_to_corrected_cartesian, points); color = point_colors)
-#
-stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
-scatter(stereographic_points; color = point_colors)
-#
-triangulation_points = copy(stereographic_points)
-
-
-# Iterate through the points and find infinite/invalid points
-max_radius = 1
-really_far_idxs = Int[]
-for (i, point) in enumerate(stereographic_points)
-    m = hypot(point[1], point[2])
-    if !isfinite(m) || m > 1e32
-        push!(really_far_idxs, i)
-    elseif m > max_radius
-        max_radius = m
-    end
-end
-@show max_radius length(really_far_idxs)
-far_value = 1e6 * sqrt(max_radius)
-
-if !isempty(really_far_idxs)
-    triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
-end
-
-boundary_points = reverse([
-    (-far_value, -far_value),
-    (-far_value, far_value),
-    (far_value, far_value),
-    (far_value, -far_value),
-    (-far_value, -far_value),
-])
-
-boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
-
-triangulation = DelTri.triangulate(pts; boundary_nodes)
-# triangulation = DelTri.triangulate(triangulation_points)
-triplot(triangulation)
-DelTri.validate_triangulation(triangulation)
-#  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
-
-# First, get all the "solid" faces, ie, faces not attached to boundary nodes
-original_triangles = collect(DelTri.each_solid_triangle(triangulation))
-boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
-faces = map(CairoMakie.GeometryBasics.TriangleFace, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
-
-for boundary_face in view(original_triangles, boundary_face_inds)
-    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
-end
-
-for ghost_face in DelTri.each_ghost_triangle(triangulation)
-    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
-end
-# Remove degenerate triangles
-filter!(faces) do face
-    !(face[1] == face[2] || face[2] == face[3] || face[1] == face[3])
-end
-
-wireframe(CairoMakie.GeometryBasics.normal_mesh(map(UnitCartesianFromGeographic(), points), faces))
-
-grid_lons = LinRange(-180, 180, 10)
-grid_lats = LinRange(-90, 90, 10)
-lon_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((grid_lons[j], -90)), Point{2, Float64}((grid_lons[j], 90))]) for j in 1:10] |> vec
-lat_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((-180, grid_lats[i])), Point{2, Float64}((180, grid_lats[i]))]) for i in 1:10] |> vec
-
-sampled_grid = GO.segmentize(lat_grid; max_distance = 0.01)
-
-geographic_grid = GO.transform(UnitCartesianFromGeographic(), sampled_grid) .|> x -> GI.LineString{true, false}(x.geom)
-stereographic_grid = GO.transform(sampled_grid) do point
-    (StereographicFromCartesian() ∘ UnitCartesianFromGeographic())(point)
-end .|> x -> GI.LineString{false, false}(x.geom)
-
-
-
-fig = Figure(); ax = LScene(fig[1, 1])
-for (i, line) in enumerate(geographic_grid)
-    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
-end
-fig
-
-#
-fig = Figure(); ax = LScene(fig[1, 1])
-for (i, line) in enumerate(stereographic_grid)
-    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
-end
-fig
-
-all(reduce(vcat, faces) |> sort! |> unique! .== 1:length(points))
-
-faces[findall(x -> 1 in x, faces)]
-
-# try NaturalNeighbours, fail miserably
-using NaturalNeighbours
-
-itp = NaturalNeighbours.interpolate(triangulation, last.(points); derivatives = true)
-lons = LinRange(-180, 180, 360)
-lats = LinRange(-90, 90, 180)
-_x = vec([x for x in lons, _ in lats])
-_y = vec([y for _ in lons, y in lats])
-
-sibson_1_vals = itp(_x, _y; method=Sibson(1))
-
-transformed_points = GO.transform(points) do point
-    # first, spherical to Cartesian
-    longitude, latitude = deg2rad(GI.x(point)), deg2rad(GI.y(point))
-    radius = 1 # we will operate on the unit sphere
-    xyz = Point3d(
-        radius * sin(latitude) * cos(longitude),
-        radius * sin(latitude) * sin(longitude),
-        radius * cos(latitude)
-    )
-    # then, rotate Cartesian so the pivot point is at the south pole
-    rotated_xyz = Rotations.AngleAxis
-end
-
-
-
-
-function randsphericalangles(n)
-    θ = 2π .* rand(n)
-    ϕ = acos.(2 .* rand(n) .- 1)
-    return Point2.(θ, ϕ)
-end
-
-function randsphere(n)
-    θϕ = randsphericalangles(n)
-    return Point3.(
-        sin.(last.(θϕs)) .* cos.(first.(θϕs)),
-        sin.(last.(θϕs)) .* sin.(first.(θϕs)),
-        cos.(last.(θϕs))
-    )
-end
-
-θϕs = randsphericalangles(50)
-θs, ϕs = first.(θϕs), last.(θϕs)
-pts = Point3.(
-    sin.(ϕs) .* cos.(θs),
-    sin.(ϕs) .* sin.(θs),
-    cos.(ϕs)
-)
-
-f, a, p = scatter(pts; color = [fill(:blue, 49); :red])
-
-function Makie.rotate!(t::Makie.Transformable, rotation::Rotations.Rotation)
-    quat = Rotations.QuatRotation(rotation)
-    Makie.rotate!(Makie.Absolute, t, Makie.Quaternion(quat.q.v1, quat.q.v2, quat.q.v3, quat.q.s))
-end
-
-
-rotate!(p, Rotations.RotX(π/2))
-rotate!(p, Rotations.RotX(0))
-
-pivot_point = θϕs[end]
-necessary_rotation = Rotations.RotY(pi) * Rotations.RotY(-pivot_point[2]) * Rotations.RotZ(-pivot_point[1])
-
-rotate!(p, necessary_rotation)
-
-f
\ No newline at end of file
diff --git a/docs/src/experiments/spherical_triangulation_experiments.jl b/docs/src/experiments/spherical_triangulation_experiments.jl
new file mode 100644
index 000000000..adc835bca
--- /dev/null
+++ b/docs/src/experiments/spherical_triangulation_experiments.jl
@@ -0,0 +1,180 @@
+
+include(joinpath(@__DIR__, "spherical_delaunay.jl"))
+
+# These points are known to be good points, i.e., lon, lat, alt
+points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+faces = spherical_triangulation(points)
+# or
+# faces = Quickhull.quickhull(map(UnitCartesianFromGeographic(), points)) |> Quickhull.faces |> collect
+# @assert isempty(setdiff(unique!(sort!(reduce(vcat, faces))), 1:length(points)))
+
+# This is the super-cool scrollable 3D globe (though it's a bit deformed...)
+f, a, p = Makie.mesh(map(UnitCartesianFromGeographic(), points), faces; color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
+
+# We can also replicate the observable notebook almost exactly (just missing ExactPredicates):
+f, a, p = Makie.mesh(points, faces; axis = (; type = GeoAxis, dest = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"), color = last.(points), colormap = Reverse(:RdBu), colorrange = (-20, 40), shading = NoShading)
+# Whoops, this doesn't look so good!  Let's try to do this more "manually" instead.
+# We'll use NaturalNeighbours.jl for this, but we first have to reconstruct the triangulation
+# with the same faces, but in a Bertin projection...
+using NaturalNeighbours
+lonlat2bertin = Proj.Transformation(GFT.EPSG(4326), GFT.ProjString("+proj=bertin1953 +type=crs +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09"); always_xy = true)
+lons = LinRange(-180, 180, 300)
+lats = LinRange(-90, 90, 150)
+bertin_points = lonlat2bertin.(lons, lats')
+
+projected_points = GO.reproject(GO.tuples(points), source_crs = GFT.EPSG(4326), target_crs = "+proj=bertin1953 +lon_0=-16.5 +lat_0=-42 +x_0=7.93 +y_0=0.09")
+
+ch = DelTri.convex_hull(projected_points) # assumes each point is in the triangulation
+boundary_nodes = DelTri.get_vertices(ch) 
+bertin_boundary_poly = GI.Polygon([GI.LineString(DelTri.get_points(ch)[DelTri.get_vertices(ch)])])
+
+tri = DelTri.Triangulation(projected_points, faces, boundary_nodes)
+itp = NaturalNeighbours.interpolate(tri, last.(points); derivatives = true)
+
+mat = [
+    if GO.contains(bertin_boundary_poly, (x, y))
+        itp(x, y; method = Nearest(), project = false)
+    else
+        NaN
+    end
+    for (x, y) in bertin_points
+]
+# TODO: this currently doesn't work, because some points are not inside a triangle and so cannot be located.
+# Options are:
+# 1. Reject all points outside the convex hull of the projected points.  (Tried but failed)
+# 2. ...
+
+
+
+
+
+# Now, we proceed with the script implementation which has quite a bit of debug information + plots
+pivot_ind = findfirst(isfinite, points)
+
+# debug
+point_colors = fill(:blue, length(points))
+point_colors[pivot_ind] = :red
+# end debug
+
+pivot_point = points[pivot_ind]
+necessary_rotation = #=Rotations.RotY(-π) *=# Rotations.RotY(-deg2rad(90-pivot_point[2])) * Rotations.RotZ(-deg2rad(pivot_point[1]))
+#
+net_transformation_to_corrected_cartesian = CoordinateTransformations.LinearMap(necessary_rotation) ∘ UnitCartesianFromGeographic()
+
+scatter(map(net_transformation_to_corrected_cartesian, points); color = point_colors)
+#
+stereographic_points = (StereographicFromCartesian() ∘ net_transformation_to_corrected_cartesian).(points)
+scatter(stereographic_points; color = point_colors)
+#
+triangulation_points = copy(stereographic_points)
+
+
+# Iterate through the points and find infinite/invalid points
+max_radius = 1
+really_far_idxs = Int[]
+for (i, point) in enumerate(stereographic_points)
+    m = hypot(point[1], point[2])
+    if !isfinite(m) || m > 1e32
+        push!(really_far_idxs, i)
+    elseif m > max_radius
+        max_radius = m
+    end
+end
+@show max_radius length(really_far_idxs)
+far_value = 1e6 * sqrt(max_radius)
+
+if !isempty(really_far_idxs)
+    triangulation_points[really_far_idxs] .= (Point2(far_value, 0.0),)
+end
+
+boundary_points = reverse([
+    (-far_value, -far_value),
+    (-far_value, far_value),
+    (far_value, far_value),
+    (far_value, -far_value),
+    (-far_value, -far_value),
+])
+
+boundary_nodes, pts = DelTri.convert_boundary_points_to_indices(boundary_points; existing_points=triangulation_points)
+
+triangulation = DelTri.triangulate(pts; boundary_nodes)
+# triangulation = DelTri.triangulate(triangulation_points)
+triplot(triangulation)
+DelTri.validate_triangulation(triangulation)
+#  for diagnostics, try `fig, ax, sc = triplot(triangulation, show_constrained_edges=true, show_convex_hull=true)`
+
+# First, get all the "solid" faces, ie, faces not attached to boundary nodes
+original_triangles = collect(DelTri.each_solid_triangle(triangulation))
+boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
+faces = map(CairoMakie.GeometryBasics.TriangleFace, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
+
+for boundary_face in view(original_triangles, boundary_face_inds)
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
+end
+
+for ghost_face in DelTri.each_ghost_triangle(triangulation)
+    push!(faces, CairoMakie.GeometryBasics.TriangleFace(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
+end
+# Remove degenerate triangles
+filter!(faces) do face
+    !(face[1] == face[2] || face[2] == face[3] || face[1] == face[3])
+end
+
+wireframe(CairoMakie.GeometryBasics.normal_mesh(map(UnitCartesianFromGeographic(), points), faces))
+
+grid_lons = LinRange(-180, 180, 10)
+grid_lats = LinRange(-90, 90, 10)
+lon_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((grid_lons[j], -90)), Point{2, Float64}((grid_lons[j], 90))]) for j in 1:10] |> vec
+lat_grid = [CairoMakie.GeometryBasics.LineString([Point{2, Float64}((-180, grid_lats[i])), Point{2, Float64}((180, grid_lats[i]))]) for i in 1:10] |> vec
+
+sampled_grid = GO.segmentize(lat_grid; max_distance = 0.01)
+
+geographic_grid = GO.transform(UnitCartesianFromGeographic(), sampled_grid) .|> x -> GI.LineString{true, false}(x.geom)
+stereographic_grid = GO.transform(sampled_grid) do point
+    (StereographicFromCartesian() ∘ UnitCartesianFromGeographic())(point)
+end .|> x -> GI.LineString{false, false}(x.geom)
+
+
+
+fig = Figure(); ax = LScene(fig[1, 1])
+for (i, line) in enumerate(geographic_grid)
+    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
+end
+fig
+
+#
+fig = Figure(); ax = LScene(fig[1, 1])
+for (i, line) in enumerate(stereographic_grid)
+    lines!(ax, line; linewidth = 3, color = Makie.resample_cmap(:viridis, length(stereographic_grid))[i])
+end
+fig
+
+all(reduce(vcat, faces) |> sort! |> unique! .== 1:length(points))
+
+faces[findall(x -> 1 in x, faces)]
+
+# try NaturalNeighbours, fail miserably
+using NaturalNeighbours
+
+itp = NaturalNeighbours.interpolate(triangulation, last.(points); derivatives = true)
+lons = LinRange(-180, 180, 360)
+lats = LinRange(-90, 90, 180)
+_x = vec([x for x in lons, _ in lats])
+_y = vec([y for _ in lons, y in lats])
+
+sibson_1_vals = itp(_x, _y; method=Sibson(1))
+
+transformed_points = GO.transform(points) do point
+    # first, spherical to Cartesian
+    longitude, latitude = deg2rad(GI.x(point)), deg2rad(GI.y(point))
+    radius = 1 # we will operate on the unit sphere
+    xyz = Point3d(
+        radius * sin(latitude) * cos(longitude),
+        radius * sin(latitude) * sin(longitude),
+        radius * cos(latitude)
+    )
+    # then, rotate Cartesian so the pivot point is at the south pole
+    rotated_xyz = Rotations.AngleAxis
+end
+
+

From cf3465c80c571425257ed3cc9d2c6e748aa7cee0 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 29 Jul 2024 08:07:32 +0530
Subject: [PATCH 23/29] Make the natural neighbour interpolation file natively
 runnable

---
 .../spherical_natural_neighbour_interpolation.jl    | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 0085574ac..3ca4be716 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -12,7 +12,7 @@ using JSON3 # to load data
 using CairoMakie, GeoMakie # for plotting
 import Makie: Point3d
 
-# include(joinpath(@__DIR__, "spherical_delaunay_stereographic.jl"))
+include(joinpath(@__DIR__, "spherical_delaunay.jl"))
 
 using LinearAlgebra
 using GeometryBasics
@@ -167,6 +167,8 @@ faces = spherical_triangulation(geographic_points)
 # correct the faces, since the order seems to be off
 faces = reverse.(faces)
 
+faces = spherical_triangulation(SphericalConvexHull(), geographic_points)
+
 unique!(sort!(reduce(vcat, faces))) # so how am I getting this index?
 
 cartesian_points = UnitCartesianFromGeographic().(geographic_points)
@@ -194,16 +196,21 @@ f # not entirely sure what's going on here
 # f, a, p = scatter(reduce(vcat, (view(cartesian_points, face) for face in view(faces, neighbour_inds))))
 # scatter!(query_point; color = :red, markersize = 40)
 
-query_point = LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0))
+query_point = LinearAlgebra.normalize(Point3(1.0, 1.0, 0.5))
 pt_inds = bowyer_watson_envelope!(Int64[], query_point, cartesian_points, faces, caps)
+three_d_points = cartesian_points[pt_inds]
+angles = [angle_between(three_d_points[1], query_point, point) for point in three_d_points]
 
+# pushfirst!(angles, 0)
+pt_inds[sortperm(angles)]
 f, a, p = scatter([query_point]; markersize = 30, color = :green, axis = (; type = LScene));
 scatter!(a, cartesian_points)
-scatter!(a, view(cartesian_points, pt_inds); color = eachindex(pt_inds), colormap = :turbo, markersize = 20)
+scatter!(a, view(cartesian_points, pt_inds[sortperm(angles)]); color = eachindex(pt_inds), colormap = :RdBu, markersize = 20)
 wireframe!(a, GeometryBasics.Mesh(cartesian_points, faces); alpha = 0.3)
 
 f
 
+
 function Makie.convert_arguments(::Type{Makie.Mesh}, cap::SphericalCap)
     offset_point = LinearAlgebra.normalize(cap.point + LinearAlgebra.normalize(Point3(cos(cap.radius), sin(cap.radius), 0.0)))
     points = [cap.point + Rotations.AngleAxis(θ, cap.point...) * offset_point for θ in LinRange(0, 2π, 20)]

From 3f366343faea836b40631d6a890cf0b53f1f013b Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Tue, 30 Jul 2024 09:23:38 +0530
Subject: [PATCH 24/29] More fixes for point order from Bowyer Watson

---
 ...herical_natural_neighbour_interpolation.jl | 28 ++++++++-----------
 1 file changed, 11 insertions(+), 17 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 3ca4be716..49788b033 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -96,23 +96,16 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
         end
     end
     unique!(applicable_points)
-    # Find a plane that can house most points OK
-    # since natural neighbours are local ish (probably wont be too far) 
-    # we can get away with this
-
-    # Start at point 1, find the first occurrence of point 1 in the applicable_points list.
-    # This is the last point of the edge coming from point 1.
-    # Now, swap the element before that with point 3.  Then continue on doing this.
-    # for (i, point_idx) in enumerate(applicable_points)
-    #     if i % 2 == 0
-    #         continue
-    #     end
-    #     applicable_points[i] = findfirst(==(applicable_points[i+1]), points)
-    # end
+    # Pick a random point (in this case, the first point) and 
+    # compute the angle subtended by the first point, the query
+    # point, and each individual point.  Sorting these angles
+    # allows us to "wind" the polygon in a clockwise way.
     three_d_points = view(points, applicable_points)
     angles = [angle_between(three_d_points[1], query_point, point) for point in three_d_points]
     pt_inds = sortperm(angles)
-    return applicable_points[pt_inds]
+    permute!(applicable_points, pt_inds)
+    push!(applicable_points, applicable_points[begin])
+    return applicable_points
 end
 
 import NaturalNeighbours: previndex_circular, nextindex_circular
@@ -162,23 +155,24 @@ end
 
 # These points are known to be good points, i.e., lon, lat, alt
 geographic_points = Point3{Float64}.(JSON3.read(read(Downloads.download("https://gist.githubusercontent.com/Fil/6bc12c535edc3602813a6ef2d1c73891/raw/3ae88bf307e740ddc020303ea95d7d2ecdec0d19/points.json"), String)))
+cartesian_points = UnitCartesianFromGeographic().(geographic_points)
 z_values = last.(geographic_points)
+
 faces = spherical_triangulation(geographic_points)
 # correct the faces, since the order seems to be off
 faces = reverse.(faces)
 
-faces = spherical_triangulation(SphericalConvexHull(), geographic_points)
+faces2 = spherical_triangulation(SphericalConvexHull(), geographic_points)
 
 unique!(sort!(reduce(vcat, faces))) # so how am I getting this index?
 
-cartesian_points = UnitCartesianFromGeographic().(geographic_points)
 
 caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces))
 
 lons = -180.0:0.5:180.0
 lats = -90.0:0.5:90.0
 
-eval_laplace_coordinates(cartesian_points, faces, z_values, Point3(1.0, 0.0, 0.0))
+eval_laplace_coordinates(cartesian_points, faces2, z_values, LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0)))
 
 values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do point
     eval_laplace_coordinates(cartesian_points, faces, z_values, point)

From 6b33e28bc8479b53cf168a02016611631083120f Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Tue, 30 Jul 2024 14:23:27 +0530
Subject: [PATCH 25/29] Get the file working again

---
 docs/src/experiments/spherical_delaunay.jl    | 10 +++---
 ...herical_natural_neighbour_interpolation.jl | 34 ++++++++++++++-----
 2 files changed, 32 insertions(+), 12 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay.jl b/docs/src/experiments/spherical_delaunay.jl
index 10e9ca964..5a35197fd 100644
--- a/docs/src/experiments/spherical_delaunay.jl
+++ b/docs/src/experiments/spherical_delaunay.jl
@@ -90,14 +90,16 @@ function spherical_triangulation(::StereographicDelaunayTriangulation, input_poi
     # First, get all the "solid" faces, ie, faces not attached to boundary nodes
     original_triangles = collect(DelTri.each_solid_triangle(triangulation))
     boundary_face_inds = findall(Base.Fix1(DelTri.is_boundary_triangle, triangulation), original_triangles)
-    faces = map(facetype, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
-    
+    # Now we construct the faces.  However, we need to reverse the order of the points 
+    # because the Delaunay triangulation is constructed in the stereographic plane,
+    # and orientation there is the opposite of orientation on the sphere.
+    faces = map(facetype ∘ reverse, view(original_triangles, setdiff(axes(original_triangles, 1), boundary_face_inds)))
     for boundary_face in view(original_triangles, boundary_face_inds)
-        push!(faces, facetype(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end))
+        push!(faces, reverse(facetype(map(boundary_face) do i; first(DelTri.is_boundary_node(triangulation, i)) ? pivot_ind : i end)))
     end
     
     for ghost_face in DelTri.each_ghost_triangle(triangulation)
-        push!(faces, facetype(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end))
+        push!(faces, reverse(facetype(map(ghost_face) do i; DelTri.is_ghost_vertex(i) ? pivot_ind : i end)))
     end
     # Remove degenerate triangles
     filter!(faces) do face
diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 49788b033..179e6ae17 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -16,7 +16,6 @@ include(joinpath(@__DIR__, "spherical_delaunay.jl"))
 
 using LinearAlgebra
 using GeometryBasics
-
 struct SphericalCap{T}
     point::Point3{T}
     radius::T
@@ -104,7 +103,7 @@ function bowyer_watson_envelope!(applicable_points, query_point, points, faces,
     angles = [angle_between(three_d_points[1], query_point, point) for point in three_d_points]
     pt_inds = sortperm(angles)
     permute!(applicable_points, pt_inds)
-    push!(applicable_points, applicable_points[begin])
+    # push!(applicable_points, applicable_points[begin])
     return applicable_points
 end
 
@@ -113,8 +112,7 @@ function laplace_ratio(points, envelope, i #= current vertex index =#, interpola
     u = envelope[i]
     prev_u = envelope[previndex_circular(envelope, i)]
     next_u = envelope[nextindex_circular(envelope, i)]
-    p = points[u]
-    q, r = points[prev_u], points[next_u]
+    p, q, r = points[u], points[prev_u], points[next_u]
     g1 = circumcenter_on_unit_sphere(q, p, interpolation_point)
     g2 = circumcenter_on_unit_sphere(p, r, interpolation_point)
     ℓ = spherical_distance(g1, g2)
@@ -159,8 +157,6 @@ cartesian_points = UnitCartesianFromGeographic().(geographic_points)
 z_values = last.(geographic_points)
 
 faces = spherical_triangulation(geographic_points)
-# correct the faces, since the order seems to be off
-faces = reverse.(faces)
 
 faces2 = spherical_triangulation(SphericalConvexHull(), geographic_points)
 
@@ -169,9 +165,10 @@ unique!(sort!(reduce(vcat, faces))) # so how am I getting this index?
 
 caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces))
 
-lons = -180.0:0.5:180.0
-lats = -90.0:0.5:90.0
+lons = -180.0:1.0:180.0
+lats = -90.0:1.0:90.0
 
+eval_laplace_coordinates(cartesian_points, faces, z_values, LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0)))
 eval_laplace_coordinates(cartesian_points, faces2, z_values, LinearAlgebra.normalize(Point3(1.0, 1.0, 0.0)))
 
 values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do point
@@ -212,3 +209,24 @@ function Makie.convert_arguments(::Type{Makie.Mesh}, cap::SphericalCap)
     faces = [GeometryBasics.TriangleFace(i, i+1, 21) for i in 1:19]
     return (GeometryBasics.normal_mesh(points, faces),)
 end
+
+
+# Animate the natural neighbours of a loxodromic curve
+# To create a loxodromic curve we use a Mercator projection
+import GeometryOps as GO, GeoInterface as GI, GeoFormatTypes as GFT
+using CairoMakie, GeoMakie
+start_point = (-95.7129, 90.0)
+end_point = (-130.0, 0.0)
+
+geographic_path = GI.LineString([start_point, end_point]; crs = GFT.EPSG(4326))
+mercator_path = GO.reproject(geographic_path; target_crs = GFT.EPSG(3857))
+sampled_loxodromic_path = GO.segmentize(GO.LinearSegments(; max_distance = 1e5), mercator_path) # max distance in mercator coordinates!
+sampled_loxodromic_cartesian_path = GO.transform(UnitCartesianFromGeographic(), GO.reproject(sampled_loxodromic_path; target_crs = GFT.EPSG(4326))) |> GI.getpoint
+
+lines(sampled_loxodromic_cartesian_path)
+
+f, a, p = lines(sampled_loxodromic_path; axis = (; type = GeoAxis))
+lines!(a, GeoMakie.coastlines(); color = (:black, 0.5))
+f
+
+start_point = 
\ No newline at end of file

From 46341c1d99e536fcc4e7edb18ff7ca91b06ba4c4 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 31 Jul 2024 07:12:12 +0530
Subject: [PATCH 26/29] Fix some bugs

- Caps not propagated
- randsphere fixed
---
 docs/src/experiments/spherical_delaunay.jl                | 2 +-
 .../spherical_natural_neighbour_interpolation.jl          | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/experiments/spherical_delaunay.jl b/docs/src/experiments/spherical_delaunay.jl
index 5a35197fd..d8251c4c3 100644
--- a/docs/src/experiments/spherical_delaunay.jl
+++ b/docs/src/experiments/spherical_delaunay.jl
@@ -182,7 +182,7 @@ function randsphericalangles(n)
 end
 
 function randsphere(n)
-    θϕ = randsphericalangles(n)
+    θϕs = randsphericalangles(n)
     return Point3.(
         sin.(last.(θϕs)) .* cos.(first.(θϕs)),
         sin.(last.(θϕs)) .* sin.(first.(θϕs)),
diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 179e6ae17..0d87a2ae9 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -56,7 +56,7 @@ function angle_between(a, b, c)
     end
 end
 
-function bowyer_watson_envelope!(applicable_points, query_point, points, faces, caps = map(splat(SphericalCap), (view(cartesian_points, face) for face in faces)); applicable_cap_indices = Int64[])
+function bowyer_watson_envelope!(applicable_points, query_point, points, faces, caps = map(splat(SphericalCap), (view(points, face) for face in faces)); applicable_cap_indices = Int64[])
     # brute force for now, but try the jump and search algorithm later
     # can use e.g GeometryBasics.decompose(Point3{Float64}, GeometryBasics.Sphere(Point3(0.0), 1.0), 5) 
     # to get starting points, or similar
@@ -127,7 +127,7 @@ struct NaturalCoordinates{F, I}
     interpolation_point::Point3{Float64}
 end
 
-function laplace_nearest_neighbour_coords(points, faces, interpolation_point; envelope = Int64[], cap_idxs = Int64[])
+function laplace_nearest_neighbour_coords(points, faces, interpolation_point, caps = map(splat(SphericalCap), (view(points, face) for face in faces)); envelope = Int64[], cap_idxs = Int64[])
     envelope = bowyer_watson_envelope!(envelope, interpolation_point, points, faces, caps; applicable_cap_indices = cap_idxs)
     coords = NaturalCoordinates(Float64[], Int64[], interpolation_point)
     for i in eachindex(envelope)
@@ -139,8 +139,8 @@ function laplace_nearest_neighbour_coords(points, faces, interpolation_point; en
     return coords
 end
 
-function eval_laplace_coordinates(points, faces, zs, interpolation_point; envelope = Int64[], cap_idxs = Int64[])
-    coords = laplace_nearest_neighbour_coords(points, faces, interpolation_point; envelope, cap_idxs)
+function eval_laplace_coordinates(points, faces, zs, interpolation_point, caps = map(splat(SphericalCap), (view(points, face) for face in faces)); envelope = Int64[], cap_idxs = Int64[])
+    coords = laplace_nearest_neighbour_coords(points, faces, interpolation_point, caps; envelope, cap_idxs)
     if isempty(coords.coordinates)
         return Inf
     end

From f8fd202c7203b9e196d1cdcd3647e5c251c76acb Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 31 Jul 2024 07:12:20 +0530
Subject: [PATCH 27/29] Add a minimal test

---
 ...herical_natural_neighbour_interpolation.jl | 26 ++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 0d87a2ae9..95b3c9bfe 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -229,4 +229,28 @@ f, a, p = lines(sampled_loxodromic_path; axis = (; type = GeoAxis))
 lines!(a, GeoMakie.coastlines(); color = (:black, 0.5))
 f
 
-start_point = 
\ No newline at end of file
+
+
+# A minimal test case to debug
+randpoints = randsphere(10)
+push!(randpoints, Point3d(0.0, 0.0, 1.0))
+
+zs = [zeros(10); 1.0]
+
+faces = spherical_triangulation(randpoints)
+caps = map(splat(SphericalCap), (view(randpoints, face) for face in faces))
+values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do query_point
+    eval_laplace_coordinates(randpoints, faces, zs, query_point, caps)
+end
+
+scatter(randpoints; color = zs)
+
+heatmap(lons, lats, values)
+
+
+f = Figure();
+a = LScene(f[1, 1])
+p = meshimage!(a, lons, lats, rotl90(values); npoints = (720, 360))
+p.transformation.transform_func[] = Makie.PointTrans{3}(UnitCartesianFromGeographic())
+# scatter!(a, cartesian_points)
+f # not entirely sure what's going on here
\ No newline at end of file

From a63c3423108023974a1c3c07cb17f11d70f1520c Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 3 Aug 2024 08:17:10 -0400
Subject: [PATCH 28/29] add spherical segmentizer

---
 ...herical_natural_neighbour_interpolation.jl | 63 +++++++++++++++++--
 1 file changed, 57 insertions(+), 6 deletions(-)

diff --git a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
index 95b3c9bfe..1a7e340a9 100644
--- a/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
+++ b/docs/src/experiments/spherical_natural_neighbour_interpolation.jl
@@ -147,7 +147,10 @@ function eval_laplace_coordinates(points, faces, zs, interpolation_point, caps =
     return sum((coord * zs[idx] for (coord, idx) in zip(coords.coordinates, coords.indices)))
 end
 
-
+function get_voronoi_polygon(point_index, points, faces, caps)
+    # For a Voronoi polygon, we can simply retrieve all triangles that have
+    # point_index as a vertex, and then filter them together.  
+end
 
 
 
@@ -232,12 +235,12 @@ f
 
 
 # A minimal test case to debug
-randpoints = randsphere(10)
-push!(randpoints, Point3d(0.0, 0.0, 1.0))
+randpoints = randsphere(50)
+push!(randpoints, LinearAlgebra.normalize(Point3d(0.0, 1.0, 1.0)))
 
-zs = [zeros(10); 1.0]
+zs = [zeros(length(randpoints)-1); 1.0]
 
-faces = spherical_triangulation(randpoints)
+faces = spherical_triangulation(SphericalConvexHull(), randpoints)
 caps = map(splat(SphericalCap), (view(randpoints, face) for face in faces))
 values = map(UnitCartesianFromGeographic().(Point2.(lons, lats'))) do query_point
     eval_laplace_coordinates(randpoints, faces, zs, query_point, caps)
@@ -253,4 +256,52 @@ a = LScene(f[1, 1])
 p = meshimage!(a, lons, lats, rotl90(values); npoints = (720, 360))
 p.transformation.transform_func[] = Makie.PointTrans{3}(UnitCartesianFromGeographic())
 # scatter!(a, cartesian_points)
-f # not entirely sure what's going on here
\ No newline at end of file
+f # not entirely sure what's going on here
+
+f, a, p = wireframe(GeometryBasics.Mesh(randpoints, faces))
+
+struct UnitSphereSegments <: GO.SegmentizeMethod
+    "The maximum distance between segments, in radians (solid angle units).  A nice value here for the globe is ~0.01." 
+    max_distance::Float64
+end
+
+function GO._fill_linear_kernel!(p, q, npoints = 10)
+    # perform slerp interpolation between p and q
+    # return npoints points on the great circle between p and q
+    # using spherical linear interpolation
+    # return the points as a vector of Point3{Float64}
+    ω = acos(clamp(dot(p, q), -1.0, 1.0))
+    if isapprox(ω, 0, atol=1e-8)
+        return fill(p, (npoints,))
+    end
+    return [
+        (sin((1 - t) * ω) * p + sin(t * ω) * q) / sin(ω)
+        for t in range(0, 1, length=npoints)
+    ]
+end
+
+function _segmentize(method::UnitSphereSegments, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
+    first_coord = GI.getpoint(geom, 1)
+    x1, y1, z1 = GI.x(first_coord), GI.y(first_coord), GI.z(first_coord)
+    new_coords = NTuple{3, Float64}[]
+    sizehint!(new_coords, GI.npoint(geom))
+    push!(new_coords, (x1, y1, z1   ))
+    for coord in Iterators.drop(GI.getpoint(geom), 1)
+        x2, y2, z2 = GI.x(coord), GI.y(coord), GI.z(coord)
+        _fill_linear_kernel!(method, new_coords, x1, y1, z1, x2, y2, z2)
+        x1, y1, z1 = x2, y2, z2
+    end 
+    return rebuild(geom, new_coords)
+end
+
+function get_face_lines(points, faces)
+    return GeometryBasics.MultiLineString(map(faces) do face
+        return GeometryBasics.LineString(vcat(
+            interpolate_line_unit_sphere(points[face[1]], points[face[2]], 10),
+            interpolate_line_unit_sphere(points[face[2]], points[face[3]], 10),
+            interpolate_line_unit_sphere(points[face[3]], points[face[1]], 10),
+        ))
+    end)
+end
+
+lines(get_face_lines(cartesian_points, faces))
\ No newline at end of file

From e23c41964cb87b146a780c0ee898aa7c5ea01441 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 3 Aug 2024 08:17:28 -0400
Subject: [PATCH 29/29] Make the `_segmentize` function more generic

---
 src/transformations/segmentize.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl
index 623637bc5..8f199e4d5 100644
--- a/src/transformations/segmentize.jl
+++ b/src/transformations/segmentize.jl
@@ -182,7 +182,7 @@ _segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom))
 This is a method which performs the common functionality for both linear and geodesic algorithms, 
 and calls out to the "kernel" function which we've defined per linesegment.
 =#
-function _segmentize(method::Union{LinearSegments, GeodesicSegments}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
+function _segmentize(method::SegmentizeMethod, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait})
     first_coord = GI.getpoint(geom, 1)
     x1, y1 = GI.x(first_coord), GI.y(first_coord)
     new_coords = NTuple{2, Float64}[]