This is similar to a problem I have been working on, but in my case I want to return the distances (or not) whether inside or outside based on prior knowledge. Below is our (still in development) ray-casting algorithm julia code from FLOWFarm that also calculates differentiable distances between a point and a polygon. We are using ForwardDiff for the derivatives. There is a discontinuity in the derivative when the points lie exactly on a vertex, but otherwise this seems to be working quite well. Note that there are some minor peripheral functions I have included, such as smooth_max() and nansafesqrt() that you can either use or replace with your own versions. abs_smooth comes from FLOWMath.
I would not be surprised if there is a better way to do this, so if anyone has any feedback I would very much appreciate it.
First, here is an example of using the function directly from FLOWFarm
using FLOWFarm; const ff=FLOWFarm
vertices = [0.0 0.0; 0.0 10.0; 10.0 10.0; 10.0 0.0]
# check if point is found inside polygon
d = ff.pointinpolygon([5.0, 5.0], vertices, return_distance=false)
And now, here is the code.
"""
pointinpolygon(point, vertices, normals=nothing; s=700, method="raycasting", shift=1E-10, return_distance=true)
Given a polygon determined by a set of vertices, determine the signed distance from the point
to the polygon.
Returns the negative (-) distance if the point is inside the polygon, positive (+) otherwise.
If return_distance is set to false, then returns -1 if in polygon or on the boundary, and 1 otherwise.
# Arguments
- `point::Vector{Number}(2)`: point of interest
- `vertices::Vector{Matrix{Number}(2)`: vertices of polygon
- `normals::Vector{Matrix{Number}(2)`: if not provided, they will be calculated
- `s::Number`: smoothing factor for ksmax function (smoothmax)
- `method::String`: currently only raycasting is available
- `shift::Float`: how far to shift point if it lies on an edge or vertex
- `return_distance::Bool`: if true, return distance. if false, return -1 if in polygon or on the boundary, and 1 otherwise.
"""
function pointinpolygon(point, vertices, normals=nothing; s=700, method="raycasting", shift=1E-10, return_distance=true)
# println(point)
if return_distance && typeof(point[1]) <: Int
throw(ArgumentError("point coordinates may not be given as Ints, must use Floats of some kind. point used $(typeof(point[1]))"))
end
if normals === nothing
normals = boundary_normals_calculator(vertices)
end
# number of turbines and boundary vertices
nvertices = size(vertices)[1]
# initialize intersection counter
intersection_counter = 0
# initialize array to hold distances from each turbine to closest boundary face
turbine_to_face_distance = zeros(typeof(point[1]), nvertices)
# get vector from turbine to the first vertex in first face
turbine_to_first_facepoint = vertices[1, :] - point # dy/dp = -1
# add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
vertices = [vertices; vertices[1,1] vertices[1,2]]
# make sure that the point is not exactly on a vertex or face
for i = 1:nvertices
if isapprox(point, vertices[i,:], atol=shift/2.0)
onvertex = true
else
onvertex = false
onface = pointonline(point, vertices[i,:], vertices[i+1,:], tol=shift/2.0)
end
if onvertex || onface
# if the point is approximately on a vertex or face, move the point slightly
# this introduces some slight error, but should be well within the error
# for actual turbine placement.
# The direction moved is perpendicular to line between the previous and
# following vertices to avoid moving along an adjacent face
if i == 1
pre_direction_vector = vertices[i+1,:] - vertices[nvertices, :]
elseif i == nvertices
pre_direction_vector = vertices[1,:] - vertices[i-1, :]
else
pre_direction_vector = vertices[i+1,:] - vertices[i-1, :]
end
# get a vector perpendicular to the pre_direction_vector
perpendicular_direction = [pre_direction_vector[2], -pre_direction_vector[1]]
# normalize perpendicular vector to make it a unit vector
# perpendicular_direction ./= norm(perpendicular_direction)
perpendicular_direction ./= nansafesqrt(sum(perpendicular_direction.^2))
# move the point by shift in the direction of the perpendicular vector
point .+= shift*perpendicular_direction
end
end
# iterate through each boundary
for j = 1:nvertices
# check if x-coordinate of turbine is between the x-coordinates of the two boundary vertices
if real(vertices[j, 1]) < real(point[1]) < real(vertices[j+1, 1]) || real(vertices[j, 1]) > real(point[1]) > real(vertices[j+1, 1])
# check to see if the turbine is below the boundary
y = (vertices[j+1, 2] - vertices[j, 2]) / (vertices[j+1, 1] - vertices[j, 1]) * (point[1] - vertices[j, 1]) + vertices[j, 2]
if real(point[2]) < real(y) #(vertices[j+1, 2] - vertices[j, 2]) / (vertices[j+1, 1] - vertices[j, 1]) * (point[1] - vertices[j, 1]) + vertices[j, 2]
# the vertical ray intersects the boundary
intersection_counter += 1
end
end
if return_distance
# define the vector from the turbine to the second point of the face
turbine_to_second_facepoint = vertices[j+1, :] - point # dy/dp = -1
# find perpendicular distance from turbine to current face (vector projection)
boundary_vector = vertices[j+1, :] - vertices[j, :]
# check if perpendicular distance is the shortest
if real(sum(boundary_vector .* -turbine_to_first_facepoint)) > 0 && real(sum(boundary_vector .* turbine_to_second_facepoint)) > 0
# if boundary_vector <= turbine_to_first_facepoint && boundary_vector <= turbine_to_second_facepoint
# perpendicular distance from turbine to face
turbine_to_face_distance[j] = abs_smooth(dot(turbine_to_first_facepoint, normals[j,:]), eps())
# check if distance to first facepoint is shortest
elseif real(sum(boundary_vector .* -turbine_to_first_facepoint)) < 0
# distance from turbine to first facepoint
# turbine_to_face_distance[j] = norm(turbine_to_first_facepoint)
turbine_to_face_distance[j] = nansafesqrt(sum(turbine_to_first_facepoint.^2))
# distance to second facepoint is shortest
else
# distance from turbine to second facepoint
# turbine_to_face_distance[j] = norm(turbine_to_second_facepoint)
turbine_to_face_distance[j] = sqrt(sum(turbine_to_second_facepoint.^2))
end
# reset for next face iteration
turbine_to_first_facepoint = turbine_to_second_facepoint # dy/dx = 1 # (for efficiency, so we don't have to recalculate for the same vertex twice)
end
end
if return_distance
# magnitude of the constraint value
c = -ff.smooth_max(-turbine_to_face_distance, s=s)
# c = -ksmax(-turbine_to_face_distance, s)
# sign of the constraint value (- is inside, + is outside)
if mod(intersection_counter, 2) == 1 #|| onvertex || onface
c = -c
end
else
if mod(intersection_counter, 2) == 1
c = -1
else
c = 1
end
end
return c
end
"""
smooth_max_ndim(x; s=100.0)
Calculate the smoothmax (a.k.a. softmax or LogSumExponential) of the elements in x.
Based on John D. Cook's writings at
(1) https://www.johndcook.com/blog/2010/01/13/soft-maximum/
and
(2) https://www.johndcook.com/blog/2010/01/20/how-to-compute-the-soft-maximum/
# Arguments
- `x::Float`: first value for comparison
- `y::Float`: second value for comparison
- `s::Float` : controls the level of smoothing used in the smooth max
"""
function smooth_max(x, y; s=10.0)
# LogSumExponential Method - used this in the past
# g = (x*exp(s*x)+y*exp(s*y))/(exp(s*x)+exp(s*y))
# non-overflowing version of Smooth Max function (see ref 2 above)
max_val = max(x, y)
min_val = min(x, y)
r = (log(1.0 + exp(s*(min_val - max_val))) + s*max_val)/s
return r
end
"""
nansafesqrt(a)
Calculate the square root of a number, but if the number is less than the given tolerance
then use the line y = a(sqrt(eps())/eps()) so that the derivative is well defined.
# Arguments
- `a::Number`: takes the square root of this value, or approximates it with a line for a < eps()
"""
function nansafesqrt(a::Number)
tol = eps()
if real(a) < tol
return a*sqrt(tol)/tol
else
return sqrt(a)
end
end
"""
single_boundary_normals_calculator(boundary_vertices)
Outputs the unit vectors perpendicular to each edge of a polygon, given the Cartesian
coordinates for the polygon's vertices.
# Arguments
- `boundary_vertices::Array{Float,1}` : m-by-2 array containing all the boundary vertices, counterclockwise
"""
function single_boundary_normals_calculator(boundary_vertices)
# get number of vertices in shape
nvertices = length(boundary_vertices[:, 1])
# add the first vertex to the end of the array to form a closed-loop
boundary_vertices = [boundary_vertices; boundary_vertices[1,1] boundary_vertices[1,2]]
# initialize array to hold boundary normals
boundary_normals = zeros(nvertices, 2)
# iterate over each boundary
for i = 1:nvertices
# create a vector normal to the boundary
boundary_normals[i, :] = [ -(boundary_vertices[i+1, 2] - boundary_vertices[i, 2]) ; boundary_vertices[i+1, 1] - boundary_vertices[i, 1] ]
# normalize the vector
boundary_normals[i, :] = boundary_normals[i, :] / norm(boundary_normals[i, :])
end
return boundary_normals
end
"""
boundary_normals_calculator(boundary_vertices; nboundaries=1)
Outputs the unit vectors perpendicular to each edge of each polygon in a set of polygons,
given the Cartesian coordinates for each polygon's vertices.
# Arguments
- `boundary_vertices::Array{Float,1}` : ragged array of arrays containing all the boundary vertices of each polygon, counterclockwise
- `nboundaries::Int` : the number of boundaries in the set
"""
function boundary_normals_calculator(boundary_vertices; nboundaries=1)
if nboundaries == 1
boundary_normals = single_boundary_normals_calculator(boundary_vertices)
else
boundary_normals = deepcopy(boundary_vertices)
for i = 1:nboundaries
boundary_normals[i] = single_boundary_normals_calculator(boundary_vertices[i])
end
end
return boundary_normals
end
