Skip to content

Commit

Permalink
Merge pull request #75 from DanielVandH/f64-conv
Browse files Browse the repository at this point in the history
Convert all numbers to Float64 internally
  • Loading branch information
DanielVandH authored Aug 6, 2023
2 parents 4167b51 + c498c41 commit dc18f8b
Show file tree
Hide file tree
Showing 19 changed files with 104 additions and 115 deletions.
10 changes: 5 additions & 5 deletions docs/src/interface/counting.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 0 additions & 1 deletion docs/src/utils.md
Original file line number Diff line number Diff line change
Expand Up @@ -45,5 +45,4 @@ convert_to_boundary_edge
get_shared_vertex
replace_boundary_triangle_with_ghost_triangle
iterated_neighbourhood
f64_getxy
```
4 changes: 2 additions & 2 deletions src/data_structures/refinement/refinement_queue.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
34 changes: 15 additions & 19 deletions src/data_structures/statistics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
= squared_triangle_area(ℓ₁², ℓ₂², ℓ₃²)
if zero(A²)
= squared_triangle_area_v2(ℓ₁², ℓ₂², ℓ₃²)
end
return F(A²)
return number_type(p)(A²)
end
function triangle_area(p, q, r)
= squared_triangle_area(p, q, r)
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
20 changes: 10 additions & 10 deletions src/data_structures/voronoi/voronoi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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))
Expand Down
22 changes: 11 additions & 11 deletions src/geometry_utils/intersections.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
14 changes: 7 additions & 7 deletions src/geometry_utils/polygons.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions src/geometry_utils/polylabel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
4 changes: 4 additions & 0 deletions src/interfaces/points.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions src/point_location/initialisers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

2 comments on commit dc18f8b

@DanielVandH
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/89159

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.8.3 -m "<description of version>" dc18f8b89cecdbbcb4dbbe5d62a66fb5bfb83652
git push origin v0.8.3

Please sign in to comment.