Convert all numbers to Float64 internally

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

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
 ```

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

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

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))

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

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

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)

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)
 

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