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. 232
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(point1) <: Int throw(ArgumentError("point coordinates may not be given as Ints, must use Floats of some kind. point used $(typeof(point1))")) 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(point1), 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. 232