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

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

.

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

— Type`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.CellQueue`

— Type`CellQueue{T}`

A struct representing the priority queue of `Cell`

s, 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_features`

— Method`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.

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

— Function`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_polygon`

— Function`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_bounds`

— Method`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!