diff --git a/docs/src/interface/counting.md b/docs/src/interface/counting.md
index 13cf7fe29..93f9657d2 100644
--- a/docs/src/interface/counting.md
+++ b/docs/src/interface/counting.md
@@ -75,27 +75,27 @@ function DT.delete_from_triangles!(tri::Set{CustomTriangle}, triangle::CustomTri
     return nothing
 end
 function DT.orient_predicate(p::CustomPoint, q::CustomPoint, r::CustomPoint)
-    o = DT.orient_predicate(getxy(p), getxy(q), getxy(r))
+    o = DT.orient_predicate(DT._getxy(p), DT._getxy(q), DT._getxy(r))
     opstats[Base.Threads.threadid()].orient_calls += 1
     return o
 end
 function DT.incircle_predicate(p::CustomPoint, q::CustomPoint, r::CustomPoint, s::CustomPoint)
-    o = DT.incircle_predicate(getxy(p), getxy(q), getxy(r), getxy(s))
+    o = DT.incircle_predicate(DT._getxy(p), DT._getxy(q), DT._getxy(r), DT._getxy(s))
     opstats[Base.Threads.threadid()].incircle_calls += 1
     return o
 end
 function DT.parallelorder_predicate(p::CustomPoint, q::CustomPoint, r::CustomPoint, s::CustomPoint)
-    o = DT.parallelorder_predicate(getxy(p), getxy(q), getxy(r), getxy(s))
+    o = DT.parallelorder_predicate(DT._getxy(p), DT._getxy(q), DT._getxy(r), DT._getxy(s))
     opstats[Base.Threads.threadid()].parallelorder_calls += 1
     return o
 end
 function DT.sameside_predicate(p::CustomPoint, q::CustomPoint, r::CustomPoint)
-    o = DT.sameside_predicate(getxy(p), getxy(q), getxy(r))
+    o = DT.sameside_predicate(DT._getxy(p), DT._getxy(q), DT._getxy(r))
     opstats[Base.Threads.threadid()].sameside_calls += 1
     return o
 end
 function DT.meet_predicate(p::CustomPoint, q::CustomPoint, r::CustomPoint, s::CustomPoint)
-    o = DT.meet_predicate(getxy(p), getxy(q), getxy(r), getxy(s))
+    o = DT.meet_predicate(DT._getxy(p), DT._getxy(q), DT._getxy(r), DT._getxy(s))
     opstats[Base.Threads.threadid()].meet_calls += 1
     return o
 end
diff --git a/docs/src/utils.md b/docs/src/utils.md
index fcd2ce0cf..7cda275d5 100644
--- a/docs/src/utils.md
+++ b/docs/src/utils.md
@@ -45,5 +45,4 @@ convert_to_boundary_edge
 get_shared_vertex
 replace_boundary_triangle_with_ghost_triangle
 iterated_neighbourhood
-f64_getxy
 ```
\ No newline at end of file
diff --git a/src/data_structures/refinement/refinement_queue.jl b/src/data_structures/refinement/refinement_queue.jl
index 57c850088..83493f774 100644
--- a/src/data_structures/refinement/refinement_queue.jl
+++ b/src/data_structures/refinement/refinement_queue.jl
@@ -101,8 +101,8 @@ function initialise_refinement_queue(tri::Triangulation, targets::RefinementTarg
         if encroachment_flag
             u, v = edge_indices(e)
             p, q = get_point(tri, u, v)
-            px, py = getxy(p)
-            qx, qy = getxy(q)
+            px, py = _getxy(p)
+            qx, qy = _getxy(q)
             e_length² = (px - qx)^2 + (py - qy)^2
             encroachment_enqueue!(queue, e, e_length²)
         end
diff --git a/src/data_structures/statistics.jl b/src/data_structures/statistics.jl
index 773058eea..b51df2736 100644
--- a/src/data_structures/statistics.jl
+++ b/src/data_structures/statistics.jl
@@ -158,7 +158,7 @@ function statistics(tri::Triangulation)
     nconstrained_edges = num_edges(constrained_edges)
     convex_hull_indices = get_convex_hull_indices(tri)
     nconvex_hull_points = max(0, length(convex_hull_indices) - 1) # -1 because the last index is the same as the first 
-    individual_statistics = Dict{V,IndividualTriangleStatistics{F}}()
+    individual_statistics = Dict{V,IndividualTriangleStatistics{Float64}}()
     sizehint!(individual_statistics, nsolid_tris)
     smallest_angle = typemax(F)
     largest_angle = typemin(F)
@@ -297,14 +297,14 @@ triangle_perimeter(ℓmin::Number, ℓmed::Number, ℓmax::Number) = ℓmin + 
 triangle_inradius(A, perimeter) = 2A / perimeter
 triangle_aspect_ratio(inradius::Number, circumradius::Number) = inradius / circumradius
 triangle_radius_edge_ratio(circumradius::Number, ℓmin::Number) = circumradius / ℓmin
-triangle_centroid(p, q, r) = ((getx(p) + getx(q) + getx(r)) / 3, (gety(p) + gety(q) + gety(r)) / 3)
+triangle_centroid(p, q, r) = ((_getx(p) + _getx(q) + _getx(r)) / 3, (_gety(p) + _gety(q) + _gety(r)) / 3)
 
 function triangle_angles(p, q, r)
     ℓ₁², ℓ₂², ℓ₃² = squared_triangle_lengths(p, q, r)
     A = triangle_area(ℓ₁², ℓ₂², ℓ₃²)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     ax, by = px - qx, py - qy
     bx, ay = px - rx, py - ry
     dotab = ax * bx + ay * by
@@ -343,16 +343,12 @@ function triangle_angles(p, q, r)
 end
 
 function squared_triangle_area(p, q, r)
-    F = number_type(p)
-    p = f64_getxy(p)
-    q = f64_getxy(q)
-    r = f64_getxy(r) # Issue 72: Float32 is just a terrible choice for computing tessellations ...
     ℓ₁², ℓ₂², ℓ₃² = squared_triangle_lengths(p, q, r)
     A² = squared_triangle_area(ℓ₁², ℓ₂², ℓ₃²)
     if A² ≤ zero(A²)
         A² = squared_triangle_area_v2(ℓ₁², ℓ₂², ℓ₃²)
     end
-    return F(A²)
+    return number_type(p)(A²)
 end
 function triangle_area(p, q, r)
     A² = squared_triangle_area(p, q, r)
@@ -365,9 +361,9 @@ function squared_triangle_lengths(p, q, r)
 end
 
 function squared_triangle_lengths_and_smallest_index(p, q, r)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     ℓ₁² = (qx - px)^2 + (qy - py)^2
     ℓ₂² = (rx - qx)^2 + (ry - qy)^2
     ℓ₃² = (px - rx)^2 + (py - ry)^2
@@ -383,9 +379,9 @@ function triangle_lengths(p, q, r)
 end
 
 function triangle_circumcenter(p, q, r, A=triangle_area(p, q, r))
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     d11 = (px - rx)^2 + (py - ry)^2
     d12 = py - ry
     d21 = (qx - rx)^2 + (qy - ry)^2
@@ -442,9 +438,9 @@ function triangle_radius_edge_ratio(p, q, r)
 end
 
 function triangle_edge_midpoints(p, q, r)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     mx = (px + qx) / 2
     my = (py + qy) / 2
     nx = (qx + rx) / 2
diff --git a/src/data_structures/voronoi/voronoi.jl b/src/data_structures/voronoi/voronoi.jl
index 9b7687438..c9e6a2f1f 100644
--- a/src/data_structures/voronoi/voronoi.jl
+++ b/src/data_structures/voronoi/voronoi.jl
@@ -170,8 +170,8 @@ function get_polygon_coordinates(vorn::VoronoiTessellation, j, bbox=nothing)
             ghost_tri = get_circumcenter_to_triangle(vorn, C[i])
             u, v, _ = indices(ghost_tri) # w is the ghost vertex
             p, q = get_generator(vorn, u, v)
-            px, py = getxy(p)
-            qx, qy = getxy(q)
+            px, py = _getxy(p)
+            qx, qy = _getxy(q)
             @assert bbox[1] ≤ px ≤ bbox[2] && bbox[1] ≤ qx ≤ bbox[2] && bbox[3] ≤ py ≤ bbox[4] && bbox[3] ≤ qy ≤ bbox[4] "The bounding box is not large enough to contain the circumcenter."
             m = (px + qx) / 2, (py + qy) / 2
             is_first = is_first_boundary_index(C, i)
@@ -183,7 +183,7 @@ function get_polygon_coordinates(vorn::VoronoiTessellation, j, bbox=nothing)
                 r = get_polygon_point(vorn, C[next_index])
             end
             if r == m # It's possible for the circumcenter to lie on the edge and exactly at the midpoint (e.g. [(0.0,1.0),(-1.0,2.0),(-2.0,-1.0)]). In this case, just rotate 
-                mx, my = getxy(m)
+                mx, my = _getxy(m)
                 dx, dy = qx - mx, qy - my
                 rotated_dx, rotated_dy = -dy, dx
                 r = mx + rotated_dx, my + rotated_dy
@@ -192,7 +192,7 @@ function get_polygon_coordinates(vorn::VoronoiTessellation, j, bbox=nothing)
                     r = mx + rotated_dx, my + rotated_dy
                 end
             end
-            r = getxy(r)
+            r = _getxy(r)
             if is_left(point_position_relative_to_line(p, q, r))
                 intersection = intersection_of_ray_with_bounding_box(m, r, bbox[1], bbox[2], bbox[3], bbox[4])
             else
@@ -487,14 +487,14 @@ function polygon_bounds(vorn::VoronoiTessellation, unbounded_extension_factor=0.
     ymin = typemax(F)
     ymax = typemin(F)
     for i in each_polygon_vertex(vorn)
-        x, y = getxy(get_polygon_point(vorn, i))
+        x, y = _getxy(get_polygon_point(vorn, i))
         xmin = min(xmin, x)
         xmax = max(xmax, x)
         ymin = min(ymin, y)
         ymax = max(ymax, y)
     end
     for i in each_generator(vorn)
-        x, y = getxy(get_generator(vorn, i))
+        x, y = _getxy(get_generator(vorn, i))
         xmin = min(xmin, x)
         xmax = max(xmax, x)
         ymin = min(ymin, y)
@@ -518,17 +518,17 @@ end
 
 function jump_to_voronoi_polygon(tri::Triangulation, q; kwargs...)
     V = jump_and_march(tri, q; kwargs...)
-    qx, qy = getxy(q)
+    qx, qy = _getxy(q)
     V = rotate_triangle_to_standard_form(V)
     i, j, k = V
     a, b = get_point(tri, i, j)
-    ax, ay = getxy(a)
-    bx, by = getxy(b)
+    ax, ay = _getxy(a)
+    bx, by = _getxy(b)
     daq = (qx - ax)^2 + (qy - ay)^2
     dbq = (qx - bx)^2 + (qy - by)^2
     if !is_boundary_index(k)
         c = get_point(tri, k)
-        cx, cy = getxy(c)
+        cx, cy = _getxy(c)
         dcq = (qx - cx)^2 + (qy - cy)^2
     else
         dcq = typemax(number_type(tri))
diff --git a/src/geometry_utils/intersections.jl b/src/geometry_utils/intersections.jl
index 006a91256..11849605e 100644
--- a/src/geometry_utils/intersections.jl
+++ b/src/geometry_utils/intersections.jl
@@ -12,8 +12,8 @@ function intersection_of_ray_with_boundary(points, boundary_nodes, p, q, tol=1e-
     # TODO: Write this in terms of an angle θ rather than q, computing θ from pq. 
     # Probably don't need bisection if we do that.
     p == q && throw(ArgumentError("p and q must be distinct."))
-    px, py = getxy(p)
-    qx, qy = getxy(q)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
 
     ## Start by finding a bracketing for the intersection point 
     t1 = zero(px)
@@ -67,7 +67,7 @@ of the rectangle that the point is on:
 - `:top`: `r` is on the top side of the rectangle.
 """
 function identify_side(r, a, b, c, d)
-    rx, ry = getxy(r)
+    rx, ry = _getxy(r)
     if rx == a
         return :left
     elseif rx == b
@@ -89,8 +89,8 @@ of the ray with the bounding box `[a, b] × [c, d]`. It is assumed that `p` is i
 the bounding box, but `q` can be inside or outside.
 """
 function intersection_of_ray_with_bounding_box(p, q, a, b, c, d)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
     pℓbx, pℓby = a, c
     pℓtx, pℓty = a, d
     prtx, prty = b, d
@@ -134,10 +134,10 @@ Finds the coordinates of the intersection of the line segment from `a` to `b`
 with the line segment from `c` to `d`.
 """
 function segment_intersection_coordinates(a, b, c, d)
-    ax, ay = getxy(a)
-    bx, by = getxy(b)
-    cx, cy = getxy(c)
-    dx, dy = getxy(d)
+    ax, ay = _getxy(a)
+    bx, by = _getxy(b)
+    cx, cy = _getxy(c)
+    dx, dy = _getxy(d)
     num = (cx - ax) * (dy - ay) - (cy - ay) * (dx - ax)
     den = (bx - ax) * (dy - cy) - (by - ay) * (dx - cx)
     α = num / den
@@ -157,8 +157,8 @@ to the left of the edge.
 function intersection_of_edge_and_bisector_ray(a, b, c)
     cert = point_position_relative_to_line(a, b, c)
     if !is_left(cert)
-        ax, ay = getxy(a)
-        bx, by = getxy(b)
+        ax, ay = _getxy(a)
+        bx, by = _getxy(b)
         m = (ax + bx) / 2, (ay + by) / 2
         return cert, m
     else
diff --git a/src/geometry_utils/polygons.jl b/src/geometry_utils/polygons.jl
index 01086eba8..8b93bc9a6 100644
--- a/src/geometry_utils/polygons.jl
+++ b/src/geometry_utils/polygons.jl
@@ -28,11 +28,11 @@ function polygon_features_single_segment(pts, boundary_nodes; scale=Val(true))
     n_edge = num_boundary_edges(boundary_nodes)
     vᵢ = get_boundary_nodes(boundary_nodes, 1)
     pᵢ = get_point(pts, vᵢ)
-    xᵢ, yᵢ = getxy(pᵢ)
+    xᵢ, yᵢ = _getxy(pᵢ)
     for j in 2:(n_edge+1)
         vᵢ₊₁ = get_boundary_nodes(boundary_nodes, j)
         pᵢ₊₁ = get_point(pts, vᵢ₊₁)
-        xᵢ₊₁, yᵢ₊₁ = getxy(pᵢ₊₁)
+        xᵢ₊₁, yᵢ₊₁ = _getxy(pᵢ₊₁)
         area_contribution = xᵢ * yᵢ₊₁ - xᵢ₊₁ * yᵢ
         cx += (xᵢ + xᵢ₊₁) * area_contribution
         cy += (yᵢ + yᵢ₊₁) * area_contribution
@@ -116,17 +116,17 @@ function distance_to_polygon(q, pts, boundary_nodes)
     end
 end
 function distance_to_polygon_single_segment(q, pts, boundary_nodes, is_in_outer=false; return_sqrt=Val(true))
-    x, y = getxy(q)
+    x, y = _getxy(q)
     F = number_type(pts)
     dist = typemax(F)
     n_edge = num_boundary_edges(boundary_nodes)
     vᵢ = get_boundary_nodes(boundary_nodes, 1)
     pᵢ = get_point(pts, vᵢ)
-    xᵢ, yᵢ = getxy(pᵢ)
+    xᵢ, yᵢ = _getxy(pᵢ)
     for j in 2:(n_edge+1)
         vᵢ₊₁ = get_boundary_nodes(boundary_nodes, j)
         pᵢ₊₁ = get_point(pts, vᵢ₊₁)
-        xᵢ₊₁, yᵢ₊₁ = getxy(pᵢ₊₁)
+        xᵢ₊₁, yᵢ₊₁ = _getxy(pᵢ₊₁)
         if (yᵢ₊₁ > y) ≠ (yᵢ > y)
             intersect_x = (xᵢ - xᵢ₊₁) * (y - yᵢ₊₁) / (yᵢ - yᵢ₊₁) + xᵢ₊₁
             if x < intersect_x
@@ -209,7 +209,7 @@ function polygon_bounds_single_segment(pts, boundary_nodes)
     for i in 1:n_edge
         vᵢ = get_boundary_nodes(boundary_nodes, i)
         pᵢ = get_point(pts, vᵢ)
-        xᵢ, yᵢ = getxy(pᵢ)
+        xᵢ, yᵢ = _getxy(pᵢ)
         xmin = min(xᵢ, xmin)
         xmax = max(xᵢ, xmax)
         ymin = min(yᵢ, ymin)
@@ -241,7 +241,7 @@ assumed that the vertices are not circular, i.e. `vertices[begin] ≠ vertices[e
 """
 function sort_convex_polygon!(vertices, points)
     cx, cy = mean_points(points, vertices)
-    to_angle = p -> atan(gety(p) - cy, getx(p) - cx)
+    to_angle = p -> atan(gety(p) - cy, _getx(p) - cx)
     vert_to_angle = v -> (to_angle ∘ get_point)(points, v)
     sort!(vertices, by=vert_to_angle)
     return vertices
diff --git a/src/geometry_utils/polylabel.jl b/src/geometry_utils/polylabel.jl
index 9e71b56d0..862c64a4b 100644
--- a/src/geometry_utils/polylabel.jl
+++ b/src/geometry_utils/polylabel.jl
@@ -37,10 +37,10 @@ function pole_of_inaccessibility(pts, boundary_nodes; precision=one(number_type(
 
     ## Initialise the priority queue and decide if the polygon centroid of bounding box centroid is the best initial guess
     _, centroid = polygon_features(pts, boundary_nodes)
-    centroid_cell = Cell(getx(centroid), gety(centroid), zero(half_width), pts, boundary_nodes)
-    bounding_box_cell = Cell(xmin + width / 2, ymin + height / 2, zero(half_width), pts, boundary_nodes)
+    centroid_cell = Cell(_getx(centroid), _gety(centroid), zero(half_width), pts, boundary_nodes)
+    bounding_box_cell = Cell(Float64(xmin + width / 2), Float64(ymin + height / 2), zero(half_width), pts, boundary_nodes)
     best_cell = centroid_cell.dist > bounding_box_cell.dist ? centroid_cell : bounding_box_cell
-    queue = CellQueue{F}()
+    queue = CellQueue{Float64}()
     insert_cell!(queue, best_cell)
 
     ## Now fill the bounding box with more cells
@@ -74,8 +74,8 @@ function process_cell!(queue::CellQueue{F}, best_cell::Cell{F}, pts, boundary_no
     end
     # We now break the cell into four, since we have not yet found the largest radius. This is the "quadtree" part of the approach.
     h = next_cell.half_width / 2
-    x = getx(next_cell)
-    y = gety(next_cell)
+    x = _getx(next_cell)
+    y = _gety(next_cell)
     cell_1 = Cell(x - h, y - h, h, pts, boundary_nodes)
     cell_2 = Cell(x + h, y - h, h, pts, boundary_nodes)
     cell_3 = Cell(x - h, y + h, h, pts, boundary_nodes)
diff --git a/src/interfaces/points.jl b/src/interfaces/points.jl
index 041e9783f..10f03ffaa 100644
--- a/src/interfaces/points.jl
+++ b/src/interfaces/points.jl
@@ -41,6 +41,10 @@ Given a point `p`, returns `(getx(p), gety(p))`.
 """
 getxy(p) = (getx(p), gety(p))
 
+@inline _getx(p) = Float64(getx(p))
+@inline _gety(p) = Float64(gety(p))
+@inline _getxy(p) = (_getx(p), _gety(p))
+
 """
     getpoint(pts::P, i)
 
diff --git a/src/point_location/initialisers.jl b/src/point_location/initialisers.jl
index 4a561453c..fec7e5091 100644
--- a/src/point_location/initialisers.jl
+++ b/src/point_location/initialisers.jl
@@ -8,7 +8,7 @@ default_num_samples(num_points::I) where {I} = ceil(I, cbrt(num_points))
 
 function compare_distance(current_dist, current_idx::I, pts, i, qx, qy) where {I}
     p = get_point(pts, i)
-    px, py = getxy(p)
+    px, py = _getxy(p)
     sq_dist = (px - qx)^2 + (py - qy)^2
     if sq_dist < current_dist
         current_dist = sq_dist
@@ -48,7 +48,7 @@ function select_initial_point(tri::Triangulation{P,Ts,I}, q;
     F = number_type(tri)
     current_dist = typemax(F)
     current_idx = I(first(point_indices))
-    qx, qy = getxy(q)
+    qx, qy = _getxy(q)
     for _ in 1:m # Not using replacement, but probability of duplicates is approximately 0.5num_solid_vertices(tri)^(-1/3)
         i = I(rand(rng, point_indices))
         is_boundary_index(i) && continue
diff --git a/src/predicates/general.jl b/src/predicates/general.jl
index 30e461464..1bc6b5acb 100644
--- a/src/predicates/general.jl
+++ b/src/predicates/general.jl
@@ -15,7 +15,7 @@ Returns `ExactPredicates.orient(p, q, r)`, in particular we return:
     \\text{orient}(p, q, r) = \\text{sgn} \\det \\begin{vmatrix} p_x & p_y & 1 \\\\ q_x & q_y & 1 \\\\ r_x & r_y & 1 \\end{vmatrix} = \\text{sgn} \\det \\begin{vmatrix} p_x-r_x & p_y-r_y \\\\ q_x-r_x & q_y-r_y \\end{vmatrix}.
     ```
 """
-orient_predicate(p, q, r) = orient(f64_getxy(p), f64_getxy(q), f64_getxy(r))
+orient_predicate(p, q, r) = orient(_getxy(p), _getxy(q), _getxy(r))
 
 """
     incircle_predicate(a, b, c, p)
@@ -34,7 +34,7 @@ Returns `ExactPredicates.incircle(a, b, c, p)`, in particular we return:
     \\text{incircle}(a, b, c, d) = \\text{sgn} \\det \\begin{vmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\\\ b_x & b_y & b_x62 + b_y^2 & 1 \\\\ c_x & c_y & c_x^2 + c_y^2 & 1 \\\\ d_x & d_y & d_x^2 + d_y^2 & 1 \\end{vmatrix} = \\text{sgn} \\det \\begin{vmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\\\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\\\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \\end{vmatrix}.
     ```
 """
-incircle_predicate(a, b, c, p) = incircle(f64_getxy(a), f64_getxy(b), f64_getxy(c), f64_getxy(p))
+incircle_predicate(a, b, c, p) = incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p))
 
 """
     parallelorder_predicate(a, b, p, q)
@@ -49,7 +49,7 @@ Returns `ExactPredicates.parallelorder(a, b, p, q)`, in particular we return:
 
     The parallelorder predicate is the same as `orient_predicate(b-a, q-p, 0)`.
 """
-parallelorder_predicate(a, b, p, q) = parallelorder(f64_getxy(a), f64_getxy(b), f64_getxy(p), f64_getxy(q))
+parallelorder_predicate(a, b, p, q) = parallelorder(_getxy(a), _getxy(b), _getxy(p), _getxy(q))
 
 """
     sameside_predicate(a, b, p)
@@ -62,9 +62,9 @@ Returns `ExactPredicates.sameside(p, a, b)` (but we redefine it here).
     main point being tested is the last argument.
 """
 function sameside_predicate(a, b, p)
-    _p = f64_getxy(p)
-    _a = f64_getxy(a)
-    _b = f64_getxy(b)
+    _p = _getxy(p)
+    _a = _getxy(a)
+    _b = _getxy(b)
     if _a < _p && _b < _p || _a > _p && _b > _p
         return 1
     elseif _a < _p && _b > _p || _a > _p && _b < _p
diff --git a/src/refinement/encroachment.jl b/src/refinement/encroachment.jl
index 6ed24791f..b6421faf4 100644
--- a/src/refinement/encroachment.jl
+++ b/src/refinement/encroachment.jl
@@ -13,9 +13,9 @@ The test is done by checking if `dot(p - r, q - r) ≤ 0`, where `< 0` means obt
     we do not employ this here.
 """
 function opposite_angle_is_obtuse(p, q, r) # https://math.stackexchange.com/a/701682/861404 
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     d = (px - rx) * (qx - rx) + (py - ry) * (qy - ry)
     if d < zero(d)
         return Certificate.Obtuse
@@ -115,8 +115,8 @@ Ruppert's original paper [here](https://www.nas.nasa.gov/assets/pdf/techreports/
 and the book by Cheng, Dey, and Shewchuk.
 """
 function compute_concentric_shell_ternary_split_position(p, q)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
     ℓ² = (px - qx)^2 + (py - qy)^2
     ℓ = sqrt(ℓ²)
     #=
@@ -145,8 +145,8 @@ function compute_concentric_shell_ternary_split_position(p, q)
     return t
 end
 function compute_concentric_shell_quarternary_split_position(p, q)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
     ℓ² = (px - qx)^2 + (py - qy)^2
     ℓ = sqrt(ℓ²)
     s = balanced_power_of_two_ternary_split(ℓ)
diff --git a/src/refinement/quality_assessment.jl b/src/refinement/quality_assessment.jl
index b2e6aa312..943ae0812 100644
--- a/src/refinement/quality_assessment.jl
+++ b/src/refinement/quality_assessment.jl
@@ -53,9 +53,9 @@ function edge_is_seditious(tri::Triangulation, u, v, w, smallest_idx, ratio_flag
     !shared_flag && return false
     # ρ ≥ 5.0 && return true
     p, q, r = get_point(tri, i, j, common_vertex)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     ℓrp² = (px - rx)^2 + (py - ry)^2
     ℓrq² = (qx - rx)^2 + (qy - ry)^2
     # Kept running into precision errors when only checking if ℓrp² == ℓrq² or e.g. (1-ε) < ℓrp² / ℓrq² < (1+ε)
@@ -113,8 +113,8 @@ function assess_added_triangles!(tri::Triangulation, queue, events, targets)
                 if contains_constrained_edge(tri, e) && is_encroached(tri, e)
                     u, v = edge_indices(e)
                     p, q = get_point(tri, u, v)
-                    px, py = getxy(p)
-                    qx, qy = getxy(q)
+                    px, py = _getxy(p)
+                    qx, qy = _getxy(q)
                     ℓ² = (qx - px)^2 + (qy - py)^2
                     encroachment_enqueue!(queue, e, ℓ²)
                 end
diff --git a/src/refinement/refinement_operations.jl b/src/refinement/refinement_operations.jl
index b39b80441..db1a7787f 100644
--- a/src/refinement/refinement_operations.jl
+++ b/src/refinement/refinement_operations.jl
@@ -25,9 +25,9 @@ function try_circumcenter_insertion!(tri::Triangulation, T, events::InsertionEve
     A² ≤ ε^2 && return Cert.PrecisionFailure
 
     cx, cy = triangle_circumcenter(p, q, r)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
-    rx, ry = getxy(r)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
+    rx, ry = _getxy(r)
     # Just as we need when checking for precision issues when splitting a segment, we need to check that the circumcenter is not one of the vertices.
     if ((px == cx) && (py == cy)) || ((qx == cx) && (qy == cy)) || ((rx == cx) && (ry == cy)) || # equality is actually valid to check here, sometimes the precison really does push them to be equal. We do it separately as it's faster.
        (!iszero(cx) && !iszero(cy) && ((abs((px - cx) / cx) ≤ ε) && (abs((py - cy) / cy) ≤ ε)) || ((abs((qx - cx) / cx) ≤ ε) && (abs((qy - cy) / cy) ≤ ε)) || ((abs((rx - cx) / cx) ≤ ε) && (abs((ry - cy) / cy) ≤ ε))) ||
@@ -73,9 +73,9 @@ function try_circumcenter_insertion!(tri::Triangulation, T, events::InsertionEve
     end
     v1, v2, v3 = indices(V)
     p1, p2, p3 = get_point(tri, v1, v2, v3)
-    p1x, p1y = getxy(p1)
-    p2x, p2y = getxy(p2)
-    p3x, p3y = getxy(p3)
+    p1x, p1y = _getxy(p1)
+    p2x, p2y = _getxy(p2)
+    p3x, p3y = _getxy(p3)
     if ((p1x == cx) && (p1y == cy)) || ((p2x == cx) && (p2y == cy)) || ((p3x == cx) && (p3y == cy)) ||
        (!iszero(cx) && !iszero(cy) && ((abs((p1x - cx) / cx) ≤ ε) && (abs((p1y - cy) / cy) ≤ ε)) || ((abs((p2x - cx) / cx) ≤ ε) && (abs((p2y - cy) / cy) ≤ ε)) || ((abs((p3x - cx) / cx) ≤ ε) && (abs((p3y - cy) / cy) ≤ 1e-14))) ||
        ((abs(p1x - cx) ≤ ε) && (abs(p1y - cy) ≤ ε)) || ((abs(p2x - cx) ≤ ε) && (abs(p2y - cy) ≤ ε)) || ((abs(p3x - cx) ≤ ε) && (abs(p3y - cy) ≤ ε))
@@ -104,8 +104,8 @@ function try_circumcenter_insertion!(tri::Triangulation, T, events::InsertionEve
             any_encroached = true
             u, v = edge_indices(edge)
             p, q = get_point(tri, u, v)
-            px, py = getxy(p)
-            qx, qy = getxy(q)
+            px, py = _getxy(p)
+            qx, qy = _getxy(q)
             ℓ² = (qx - px)^2 + (qy - py)^2
             encroachment_enqueue!(queue, edge, ℓ²)
         end
@@ -169,8 +169,8 @@ end
 function compute_split_position(tri::Triangulation, p, q, e, segment_list)
     if e ∈ segment_list
         # Split at the midpoint when splitting for the first time
-        px, py = getxy(p)
-        qx, qy = getxy(q)
+        px, py = _getxy(p)
+        qx, qy = _getxy(q)
         mx, my = px + (qx - px) / 2, py + (qy - py) / 2
         return mx, my
     else
@@ -187,8 +187,8 @@ function compute_split_position(tri::Triangulation, p, q, e, segment_list)
             t = inv(2.0)
         end
         # Need to place the split on the side that is closer to the midpoint
-        px, py = getxy(p)
-        qx, qy = getxy(q)
+        px, py = _getxy(p)
+        qx, qy = _getxy(q)
         mx, my = px + (qx - px) / 2, py + (qy - py) / 2
         if abs(t - 1 / 2) < 1e-6 # close enough to the midpoint, so just bisect - every third split on a segment should be a bisection anyway. 
             return mx, my
@@ -213,8 +213,8 @@ function split_subsegment!(tri::Triangulation, queue, events, targets, segment_l
     mx, my = compute_split_position(tri, p, q, e, segment_list)
     # We can run into precision issues (e.g. I found this was needed when reproducing Figure 6.14 in the Delaunay 
     # Mesh Generation book). Let's ensure that the computed splitting point isn't just one of p or q.
-    px, py = getxy(p)
-    qx, qy = getxy(q)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
     if ((mx == px) && (my == py)) || ((mx == qx) && (my == qy)) ||
        (abs(mx - px) ≤ ε) && (abs(my - py) ≤ ε) || (abs(mx - qx) ≤ ε) && (abs(my - qy) ≤ ε)
         return nothing
diff --git a/src/utils.jl b/src/utils.jl
index 49edc1496..d98a8443a 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -735,13 +735,13 @@ function edge_lies_on_two_distinct_segments(tri::Triangulation, i, j)
     else
         F = number_type(tri)
         p, q = get_point(tri, i, j)
-        px, py = getxy(p)
-        qx, qy = getxy(q)
+        px, py = _getxy(p)
+        qx, qy = _getxy(q)
         min_dist = typemax(F)
         min_idx = I(0)
         for k in i_segments
             r = get_point(tri, k)
-            rx, ry = getxy(r)
+            rx, ry = _getxy(r)
             δ = squared_distance_to_segment(px, py, qx, qy, rx, ry)
             if δ < min_dist
                 min_dist = δ
@@ -888,14 +888,4 @@ function iterated_neighbourhood!(neighbours, tri, i, d)
         union!(neighbours, new_neighbours)
     end
     return neighbours
-end
-
-"""
-    f64_getxy(p)
-
-Returns the coordinates of the point `p` as `Float64`s.
-"""
-function f64_getxy(p)
-    x, y = getxy(p)
-    return Float64(x), Float64(y)
 end
\ No newline at end of file
diff --git a/src/voronoi/clipped_construction.jl b/src/voronoi/clipped_construction.jl
index 4e3a1b083..6a44c0633 100644
--- a/src/voronoi/clipped_construction.jl
+++ b/src/voronoi/clipped_construction.jl
@@ -196,8 +196,8 @@ Enqueue the edge `e` of the boundary to be processed.
 function enqueue_new_edge!(polygon_edge_queue, vorn::VoronoiTessellation, e)
     u, v = edge_indices(e)
     p, q = get_generator(vorn, u, v)
-    px, py = getxy(p)
-    qx, qy = getxy(q)
+    px, py = _getxy(p)
+    qx, qy = _getxy(q)
     m = (px + qx) / 2, (py + qy) / 2
     incident_polygon = jump_and_march(vorn, m; k=u)
     enqueue!(polygon_edge_queue, (e, incident_polygon))
diff --git a/src/voronoi/lloyd.jl b/src/voronoi/lloyd.jl
index 58a82e61b..95c378006 100644
--- a/src/voronoi/lloyd.jl
+++ b/src/voronoi/lloyd.jl
@@ -1,8 +1,8 @@
 function move_generator_to_centroid!(points, vorn::VoronoiTessellation, generator)
     c = get_centroid(vorn, generator)
-    cx, cy = getxy(c)
+    cx, cy = _getxy(c)
     p = get_generator(vorn, generator)
-    px, py = getxy(p)
+    px, py = _getxy(p)
     dist = sqrt((cx - px)^2 + (cy - py)^2)
     set_point!(points, generator, cx, cy)
     return dist
diff --git a/src/voronoi/unbounded_construction.jl b/src/voronoi/unbounded_construction.jl
index 76fd2577c..066f922f1 100644
--- a/src/voronoi/unbounded_construction.jl
+++ b/src/voronoi/unbounded_construction.jl
@@ -7,7 +7,7 @@ function initialise_voronoi_tessellation(tri::Tr) where {Tr<:Triangulation}
     I = integer_type(tri)
     T = triangle_type(tri)
     F = number_type(tri)
-    P = NTuple{2,F}
+    P = NTuple{2,Float64}
     polygon_points = Vector{P}()
     circumcenter_to_triangle = Dict{I,T}()
     triangle_to_circumcenter = Dict{T,I}()
diff --git a/test/utils.jl b/test/utils.jl
index 03e4de2a5..cfadc3303 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -909,7 +909,7 @@ end
       ))
 end
 
-@testset "f64_getxy" begin 
-      @test DT.f64_getxy((0.3,0.5)) == (0.3,0.5)
-      @test DT.f64_getxy((0.3f0, 0.7f0)) == (Float64(0.3f0), Float64(0.7f0))
+@testset "_getxy" begin 
+      @test DT._getxy((0.3,0.5)) == (0.3,0.5)
+      @test DT._getxy((0.3f0, 0.7f0)) == (Float64(0.3f0), Float64(0.7f0))
 end
\ No newline at end of file