unexpected/wrong output of Polygon::offsetInward function

[Arjan-TNO]

When using the Polygon::offsetInward function I get an unexpected result for the following input:

Polygon polygon(
	{Position(0.0, 0.0),
	Position(1.0, 1.0),
	Position(2.0, 1.0),
	Position(3.0, 0.0),
	Position(3.0, 1.0),
	Position(3.0, 2.0),
	Position(2.0, 3.0),
	Position(1.0, 3.0),
	Position(0.0, 2.0),
	Position(0.0, 1.0)});
polygon.offsetInward(0.1);

The problem here is two-fold, the function does not seem to work well :

1. for non-convex polygons (see sequence of points: 0, 1, 2, 3 in picture below),
2. with points in the polygon that are positioned exactly (on a straight line) between their neighbors (see points 4 & 9, in picture below).

I have noticed the statement on non-convex polygons here, but that same statement is in a TODO, so this might be the moment to address convex part of the TODO... :-)

A picture of the result is below (red dashed line is output of the offsetInward function):

A small analysis (in Matlab which for me is a tool I'm most familiar with) shows some intermediate steps in the process:

Here, the purple and green vectors at each point, are oriented towards their respective previous and next neighbor. Theses correspond with the v1 (purple) and the v2 (green) vectors of the function.
The black vectors is the summation of the green and purple vectors and is the direction over which the offset is applied.

Especially for the case of the 2nd problem I mention (points 4 & 9) it is obvious: the summation of these (v1+v2) equals zero, and therefore any offset value will be placed exactly on the point itself.

[Arjan-TNO]

A possible solution for this issue is to make use of the cross products of the v1 and v2 vectors with an up-vector. Basically the vectors between a point and its next neighbor, form the x-axis.

Together with some up-vector (which I take as z=[0,0,1]) the cross products result in two y-vectors: $y_1 = z \times v_1$ and $y_2 = z \times -v_2$.

The vector summation of both y-vectors is the vector along which the offset should be applied.
$y_0 = y_1 + y_2 $

The unit-vector for this is:
$\hat{n} = \frac{y_0}{|y_0|}$

Note that for the assumption that all points in the polygon are in counterclockwise order (which is the assumption that the function makes), the unit vectors point inwards, otherwise they will point outward.

The offset itself should be calculated, using:

$\alpha = cos^{-1}(y_1 \cdot y_2)$

and

$s = \frac{margin}{cos(\frac{1}{2}\alpha)}$

Now, the point can be offset correctly using:
$p_{new} = p_{orig} + s . \hat{n}$

A result of this is presented below (red dashed is original output of function, blue dashed line is the output of the suggested approach of above):

To generate the picture of above, I have changed the original function to:

bool Polygon::offsetInward(const double margin) {
    // Assuming an ordered set of vertices (clockwise or counterclockwise)

    const unsigned int n = nVertices();

    // vector with offset-vectors
    std::vector<Eigen::Vector2d> offset_vectors;
    offset_vectors.resize(n);

    // integrate angles to determine if polygon is clockwise or not
    double sum_angles = 0;

    // loop over all points
    for (unsigned int i = 0; i < n; ++i) {
        int i_prev = (i > 0 ? (i-1)%n : n-1);
        int i_next = (i + 1) % n;

        // define incoming/outgoing edges as x-axis
        Eigen::Vector2d x1 = vertices_[i] - vertices_[i_prev];
        Eigen::Vector2d x2 = vertices_[i_next] - vertices_[i];

        // normalize the vectors
        x1.normalize();
        x2.normalize();
    
        // create 'y-axis' using ISO definition (cross-product of z=[0;0;1] with x1 and x2)
        Eigen::Vector2d y1(-x1[1],x1[0]);
        Eigen::Vector2d y2(-x2[1],x2[0]);

        // determine vector over which to apply the offset
        Eigen::Vector2d v = y1 + y2;
        v.normalize();

        // determine angle between lines and add them to sum
        // using the atan2 approach here, to avoid y1.dot(y2)>1 which sometimes occurs!
        const double angle1 = atan2(y1[1],y1[0]);
        const double angle2 = atan2(y2[1],y2[0]);
        const double angle = atan2(sin(angle1-angle2),cos(angle1-angle2));
        sum_angles += angle;

        // store offset vector
        offset_vectors[i] = (margin / cos(angle/2)) * v;
    }

    if (sum_angles!=sum_angles){
        // Somehow, a NaN has entered the scene
        return false;
    } else {
        // check result of sum_angles to see if polygon is (counter-)clockwise
        bool counter_clockwise = (sum_angles<0);
        for (unsigned int i = 0; i < n; ++i) {
            if (counter_clockwise){
                vertices_[i] = vertices_[i] + offset_vectors[i];
            } else {
                vertices_[i] = vertices_[i] - offset_vectors[i];
            }
        }
    }
        
    return true;
}