Pole of Inaccessibility and Polygons

We provide a function for computing the pole of inaccessibility of a given polygon, namely the point inside the polygon that is furthest from the boundary. Our method is primarily based on this blogpost, recursively subdividing the polygon using quadtree partitioning. The function for this is pole_of_inaccessibility:

DelaunayTriangulation.pole_of_inaccessibilityFunction
pole_of_inaccessibility(pts, boundary_nodes; precision = one(number_type(pts)))

Given a collection of points pts and a set of boundary_nodes defining the polygon connections, finds the pole of inaccessibility. This works for multiply-connected polygons, provided boundary_nodes matches the specification given in the documentation. You can control the tolerance of the returned pole using precision.

This function is also commonly called polylabel.

Note

The pole of inaccessibility is a point within a polygon that is furthest from an edge. It is useful for our purposes since it is a representative point that is guaranteed to be inside the polygon, in contrast to for example a centroid which is not always inside the polygon.

For more information about this, see e.g. this blog post or the original repo. This implementation was partially based on the python implementation and this other Julia implementation.

We needed this method since the point we need to associate ghost vertices with must be inside the domain, and so other representative points like centroids or arithmetic averages would not be sufficient if the domain is non-convex.

Below we also list some other relevant docstrings.

DelaunayTriangulation.CellType
Cell{T}

A cell in a grid. The cell is a square with side length 2half_width. The cell is centered at (x, y). The cell is assumed to live in a polygon.

Fields

• x::T

The x-coordinate of the center of the cell.

• y::T

The y-coordinate of the center of the cell.

• half_width::T

The half-width of the cell.

• dist::T

The distance from the center of the cell to the polygon.

• max_dist::T

The maximum distance from the center of the cell to the polygon. This is dist + half_width * sqrt(2).

Constructors

Cell(x::T, y::T, half_width::T, pts, boundary_nodes)

Constructs a cell with center (x, y) and half-width half_width. The cell is assumed to live in the polygon defined by pts and boundary_nodes.

DelaunayTriangulation.CellQueueType
CellQueue{T}

A struct representing the priority queue of Cells, used for sorting the cells in a grid according to their maximum distance.

Fields

• queue::PriorityQueue{Cell{T},T,typeof(Base.Order.Reverse)}

The priority queue of cells.

Constructors

CellQueue{T}()

Constructs a new CellQueue with elements of type Cell{T}.

DelaunayTriangulation.polygon_featuresMethod
polygon_features(pts, boundary_nodes)

Returns features of the polygon represented by the points pts with boundary_nodes defining the polygon connections. The features returned are (a, c), where a is the area of the polygon and c = (cx, cy) is the centroid.

Note
• The polygon is assumed to be simple, i.e. no self-intersections.
• The function works with holes, provided boundary_nodes represents these as described in the documentation.
• The polygon is assumed to have a consistent orientation for each boundary. If the orientation is positive, a > 0, and a < 0 otherwise.
DelaunayTriangulation.squared_distance_to_segmentFunction
squared_distance_to_segment(x₁, y₁, x₂, y₂, x, y)

Given a line segment (x₁, y₁) → (x₂, y₂) and a query point (x, y), returns the squared distance from (x, y) to the line segment.

DelaunayTriangulation.distance_to_polygonFunction
distance_to_polygon(q, pts, boundary_nodes)

Given a polygon represented by the points pts with boundary_nodes defining the polygon connections, and a query point q, returns the signed distance from q to the polygon. If q is outside of the polygon, then the returned distance is negative, and if it is inside then the distance is positive. Works with holes, provided boundary_nodes matches the specification of a boundary given in the documentation.

DelaunayTriangulation.polygon_boundsMethod
polygon_bounds(pts, boundary_nodes, check_all_curves = Val(false))

Given a polygon represented by the points pts with boundary_nodes defining the polygon connections, returns a bounding box of the polygon. The bounding box is returned in the order (xmin, xmax, ymin, ymax). If your polygon is not a multiple polygon, check_all_curves = Val(false) is sufficient, otherwise you might want to use Val(true).

DelaunayTriangulation.sort_convex_polygon!Function
sort_convex_polygon!(vertices, points)

Sorts the vertices of a convex polygon in counter-clockwise order. The polygon is defined by the points points and the vertices vertices. The vertices are sorted in place. It is assumed that the vertices are not circular, i.e. vertices[begin] ≠ vertices[end].

distance_to_polygon is also useful for point location, as shown in the examples below.

If you need to compute this for multiple boundaries, meaning multiple poles, use compute_representative_points!.

DelaunayTriangulation.compute_representative_points!Function
compute_representative_points!(tri::Triangulation; use_convex_hull=!has_multiple_segments(tri) && num_boundary_edges(get_boundary_nodes(tri)) == 0)

Updates get_representative_point_list(tri) to match the current position of the boundaries. If there are no boundary nodes, use_convex_hull instead represents them using the indices of the convex hull.

Example: Pole of inaccessibility

Below is a simple example of computing this pole of inaccessibility.

pts = [0.0 8.0
2.0 5.0
3.0 7.0
1.81907 8.13422
3.22963 8.865
4.24931 7.74335
4.50423 5.87393
3.67149 4.3784
2.73678 2.62795
5.50691 1.38734
8.43 2.74691
9.7046 5.53404
8.56595 7.79433
6.71353 9.03494
4.13034 9.66375
2.75378 10.3775
1.0883 10.4965
-1.138 9.83369
-2.25965 8.45712
-2.78649 5.94191
-1.39292 3.64763
0.323538 4.97322
-0.900078 6.6217
0.98633 9.68074
0.153591 9.54478
0.272554 8.66106
2.90673 8.18521
2.12497 9.42582
7.27436 2.7979
3.0 4.0
5.33697 1.88019]'
boundary_nodes = [
[[1, 4, 3, 2], [2, 9, 10, 11, 8, 7, 12], [12, 6, 13, 5, 14, 15, 16, 17, 16], [16, 17, 18, 19, 20, 21, 22, 23, 1]],
[[26, 25, 24], [24, 28, 27, 26]],
[[29, 30, 31, 29]]
]
x, y = DT.pole_of_inaccessibility(pts, boundary_nodes)

fig = Figure()
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y")
bn1 = pts[:, unique(reduce(vcat, boundary_nodes[1]))] |> x -> hcat(x, x[:, begin])
bn2 = pts[:, unique(reduce(vcat, boundary_nodes[2]))] |> x -> hcat(x, x[:, begin])
bn3 = pts[:, unique(reduce(vcat, boundary_nodes[3]))] |> x -> hcat(x, x[:, begin])
lines!(ax, bn1, color=:red, linewidth=4)
lines!(ax, bn2, color=:red, linewidth=4)
lines!(ax, bn3, color=:red, linewidth=4)
scatter!(ax, [x], [y], color=:blue, markersize=23)

Example: Querying if points are in a polygon

If we need to, distance_to_polygon is a nice way for querying whether a point is inside or outside of a given polygon. Let's use the Julia example. Let us start by defining the polygon (this is a multiple polygon, it works for even simpler polygons obviously) and placing random points inside its bounding box:

D = (14.2705719699703, 32.8530791545746)
E = (14.3, 27.2)
F = (14.1, 27.0)
G = (13.7, 27.2)
H = (13.4, 27.5)
I = (13.1, 27.6)
J = (12.7, 27.4)
K = (12.5, 27.1)
L = (12.7, 26.7)
M = (13.1, 26.5)
N = (13.6, 26.4)
O = (14.0, 26.4)
P = (14.6, 26.5)
Q = (15.1983491346581, 26.8128534095401)
R = (15.6, 27.6)
S = (15.6952958264624, 28.2344688505621)
T = (17.8088971520274, 33.1192363585346)
U = (16.3058917649589, 33.0722674401887)
V = (16.3215480710742, 29.7374742376305)
W = (16.3841732955354, 29.393035503094)
Z = (16.6190178872649, 28.9233463196351)
A1 = (17.0417381523779, 28.5319386667527)
B1 = (17.5114273358368, 28.3753756055997)
C1 = (18.1376795804487, 28.3597192994844)
D1 = (18.7169629067146, 28.5632512789833)
E1 = (19.2805899268653, 28.8920337074045)
F1 = (19.26493362075, 28.4536571361762)
G1 = (20.6426885588962, 28.4223445239456)
H1 = (20.689657477242, 33.1035800524193)
I1 = (19.2805899268653, 33.0722674401887)
J1 = (19.2962462329806, 29.7531305437458)
K1 = (19.0614016412512, 29.393035503094)
L1 = (18.7482755189452, 29.236472441941)
M1 = (18.4508057027546, 29.1425346052493)
N1 = (18.1689921926793, 29.3147539725175)
O1 = (17.7932408459121, 29.6278800948235)
P1 = (22.6466957416542, 35.4207133574833)
Q1 = (21.2219718851621, 34.9979930923702)
R1 = (21.2376281912774, 28.4693134422915)
S1 = (22.6780083538847, 28.4380008300609)
T1 = (24.5724213938357, 33.1975178891111)
U1 = (23.3512295168425, 32.8530791545746)
V1 = (23.3199169046119, 28.4380008300609)
W1 = (24.6663592305274, 28.3753756055997)
Z1 = (15.1942940307729, 35.4363696635986)
A2 = (14.7246048473139, 35.3737444391374)
B2 = (14.3645098066621, 35.1858687657538)
C2 = (14.1766341332786, 34.8570863373326)
D2 = (14.1140089088174, 34.3247719294125)
E2 = (14.2705719699703, 33.8394264398383)
F2 = (14.7246048473139, 33.6202381542241)
G2 = (15.4604512347329, 33.6045818481088)
H2 = (16.0, 34.0)
I2 = (15.9771093365377, 34.6848669700643)
J2 = (15.6170142958859, 35.2328376840997)
K2 = (24.1653574348379, 35.4520259697138)
L2 = (23.7739497819555, 35.4363696635986)
M2 = (23.4608236596496, 35.2641502963303)
N2 = (23.272947986266, 34.9040552556785)
O2 = (23.1320412312284, 34.5909291333725)
P2 = (23.1163849251131, 34.2151777866054)
Q2 = (23.2886042923813, 33.8081138276077)
R2 = (23.8209187003014, 33.6045818481088)
S2 = (24.3062641898756, 33.5576129297629)
T2 = (24.7602970672192, 33.8550827459536)
U2 = (25.010797965064, 34.4656786844502)
V2 = (24.8385785977957, 34.9666804801397)
W2 = (24.5254524754898, 35.2641502963303)
Z2 = (25.3708930057158, 37.4716894585871)
A3 = (24.7916096794498, 37.3464390096648)
B3 = (24.4471709449133, 36.9550313567823)
C3 = (24.3062641898756, 36.5636237038999)
D3 = (24.4941398632592, 35.9999966837492)
E3 = (25.0264542711793, 35.5929327247515)
F3 = (25.5587686790994, 35.5929327247515)
F3 = (25.5587686790994, 35.5929327247515)
G3 = (26.0, 36.0)
H3 = (26.1380520053653, 36.5792800100152)
I3 = (26.0, 37.0)
J3 = (25.7466443524829, 37.2838137852036)
K3 = (26.3885529032101, 35.4676822758291)
L3 = (25.9814889442124, 35.3580881330221)
M3 = (25.6840191280217, 35.1858687657538)
N3 = (25.5274560668688, 34.9040552556785)
O3 = (25.4961434546382, 34.5596165211419)
P3 = (25.5274560668688, 34.246490398836)
Q3 = (25.6683628219064, 33.8394264398383)
R3 = (26.0284578625583, 33.6358944603394)
S3 = (26.5451159643631, 33.6202381542241)
T3 = (27.0, 34.0)
U3 = (27.280962351782, 34.5596165211419)
V3 = (27.0304614539373, 35.2171813779844)
W3 = (26.1693646175959, 33.087923746304)
Z3 = (26.0, 33.0)
A4 = (25.5274560668688, 32.7278287056522)
B4 = (25.2612988629087, 32.4147025833463)
C4 = (25.1830173323322, 32.0702638488098)
D4 = (25.2299862506781, 31.7727940326191)
E4 = (25.6527065157911, 31.5222931347744)
F4 = (26.2946150665183, 31.7258251142732)
G4 = (26.5607722704784, 32.5086404200381)
H4 = (27.1557119028596, 32.7434850117675)
I4 = (27.6097447802033, 32.4929841139228)
J4 = (27.6410573924338, 32.1015764610403)
K4 = (27.7193389230103, 31.6005746653509)
L4 = (27.437525412935, 31.4283552980826)
M4 = (26.9834925355914, 31.2561359308143)
N4 = (26.5764285765937, 31.0995728696614)
O4 = (26.0441141686736, 30.7864467473554)
P4 = (25.6527065157911, 30.5672584617413)
Q4 = (25.3239240873699, 30.1915071149741)
R4 = (25.1673610262169, 29.8783809926682)
S4 = (25.1047358017558, 29.6122237887082)
T4 = (25.0890794956405, 29.1895035235952)
U4 = (25.2926114751393, 28.8294084829433)
V4 = (25.6840191280217, 28.5632512789833)
W4 = (26.1537083114806, 28.3753756055997)
Z4 = (26.8269294744384, 28.391031911715)
A5 = (27.4844943312809, 28.6102201973292)
B5 = (27.7342002330051, 28.7239579596219)
C5 = (27.7264126450755, 28.4202565942047)
D5 = (29.1825559185446, 28.3922538389457)
E5 = (29.1545531632856, 32.2146299318021)
F5 = (29.000538009361, 32.5786657501693)
G5 = (28.6785063238822, 32.9006974356481)
H5 = (28.3144705055149, 33.0827153448317)
I5 = (27.9084305542591, 33.2367304987563)
J5 = (27.3343740714492, 33.3207387645334)
K5 = (26.8303244767868, 33.2367304987563)
L5 = (27.6564057569279, 30.786489413592)
M5 = (27.6984098898165, 30.3944508399657)
N5 = (27.6984098898165, 29.7363860913787)
O5 = (27.5863988687804, 29.4143544059)
P5 = (27.2643671833016, 29.2043337414573)
Q5 = (26.9843396307114, 29.1763309861983)
R5 = (26.6903107004917, 29.3163447624934)
S5 = (26.5782996794556, 29.7503874690082)
T5 = (26.7603175886393, 30.3384453294476)
U5 = (27.3203726938197, 30.7024811478149)
J_curve = [[C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, C]]
U_curve = [[T, U, V, W, Z, A1, B1, C1, D1, E1, F1, G1, H1, I1, J1, K1, L1, M1, N1, O1, T]]
L_curve = [[P1, Q1, R1, S1, P1]]
I_curve = [[T1, U1, V1, W1, T1]]
A_curve_outline = [[
K5, W3, Z3, A4, B4, C4, D4, E4, F4, G4, H4, I4, J4, K4, L4, M4, N4,
O4, P4, Q4, R4, S4, T4, U4, V4, W4, Z4, A5, B5, C5, D5, E5, F5, G5,
H5, I5, J5, K5]]
A_curve_hole = [[L5, M5, N5, O5, P5, Q5, R5, S5, T5, U5, L5]]
dot_1 = [[Z1, A2, B2, C2, D2, E2, F2, G2, H2, I2, J2, Z1]]
dot_2 = [[Z2, A3, B3, C3, D3, E3, F3, G3, H3, I3, J3, Z2]]
dot_3 = [[K2, L2, M2, N2, O2, P2, Q2, R2, S2, T2, U2, V2, W2, K2]]
dot_4 = [[K3, L3, M3, N3, O3, P3, Q3, R3, S3, T3, U3, V3, K3]]
curves = [J_curve, U_curve, L_curve, I_curve, A_curve_outline, A_curve_hole, dot_1, dot_2, dot_3, dot_4]
nodes, points = convert_boundary_points_to_indices(curves)
xmin, xmax, ymin, ymax = DelaunayTriangulation.polygon_bounds(points, nodes, Val(true)) # Val(true) => check all parts of the polygon
query_points = [((xmax - xmin) * rand() + xmin, (ymax - ymin) * rand() + ymin) for _ in 1:1000]

fig = Figure()
ax = Axis(fig[1, 1])
scatter!(ax, query_points)
for nodes in nodes
lines!(ax, points[reduce(vcat, nodes)], color=:magenta, linewidth=3)
end

Now let's use distance_to_polygon to test if points are inside or outside of the logo. We colour points inside in blue, and points outside in red.

is_inside = [DelaunayTriangulation.distance_to_polygon(q, points, nodes) > 0 for q in query_points]
scatter!(ax, query_points[is_inside], color=:blue)
scatter!(ax, query_points[.!is_inside], color=:red)

Works perfectly!