0
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.