# Loss function for points inside polygon

I am trying to optimize some parameters that used to transform 2d points from a place to another (you may think of that as rotation & translation parameter for simplicity)

The parameters are considered optimal if the transformed points lay inside a pre-defined convex polygon. Otherwise, the parameters should be adjusted till all points lay inside that polygon. I do not care how the points are arranged inside that polygon, my only concern that they are inside..

How can I formulate my loss function in smooth way that the farer the point outside the polygon the more the loss function is and in the same time if they are inside the polygon, the distance does not matter.

P.S. the loss has to be differentiable since I am using derivative-based optimizer.

The first thing came in mind is to test if the point is inside the polygon (traditional point to polygon test), if it is inside, I will return 0, otherwise I will return the distance from the point to the nearest segment of the polygon. Is it valid/good? if yes, is there better?

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