Fix GetOrientation

[SamuelReithmeir]

There is a Problem in GetOrientation when the Normal is Suppsed to have one Coordinate 0, but due to floating point Input Coordinates the points are slightly non coplanar. and this the Normal has a axis value of eg 0.0000001... instead of 0. This leads to the sign beeing unequal In cases where Crossproduct is exactly 0 but the normal not. I switched logic to test the sign of the crossproduct and thus allow for errors in the normal

[SamuelReithmeir]

This also fixes Issue #4

[NMO13]

Ah sorry I totally forgot to investigate your issue. Thanks for the PR. I will check it the next couple of days!

[NMO13]

Floating point inaccuracies should not be present because of the rational datatype. Considering that the input vertices are perfectly coplanar, the normal will be correct. Can you show a concrete example where this error happens?

[SamuelReithmeir]

Okay, I think i wrapped my head around the problem finally. There are two conflicting problems. First one is that while the numbers one in your libary have perfect precision, they need to be defined/calculated before. This normally happens with floating point logic (in my case reading from a XML File) or like polygons from a 3d model.

This is why it makes it really difficult to expect the user to provide perfect coplanar points, it simply is not possible without first fixing coplanarity by looking for the "intended plane" projecting points that are not on there on there and then passing them to the triangulation. This imho should be done by this libarary itself and not be passed as responisbility to the user.

The problem in issue #4 arises exactly from this precision error. At one point it checks for direction of these three points:

While they are theoretically perfectly in one line, practically there is like a 1e-12 degree angle between them. Big Arrow is normal and center point is slightly moved backwards. => Squaredlength of CrossProduct is not zero. =>but since the crossProduct it is perpendicular to the normalVector, taking the dot from both results in (a perfect) 0 again

your version can not return 0 once the length check fails. the version using the dot product can.

So i see 2 option to solve this:

• change the logic and GetDirection => allows for non coplanar points
• If you do not want to change the logic, i would suggest adding a function inbetween SetPoints() and Triangulate() called EnforceCoplanirty() that would work like this: get the plane that contains the most points. project the other points on that plane. then after triangulate one could map the vectors back to their original value so there is no change from the input data

[NMO13]

Thanks for this comment. I checked #4, to me it looks like this polygon is self intersecting and that is the reason why its not working. But I will check that again in more detail. OK so you want to allow non-perfect coplanar input. I don't see a reason why we should not allow that. Can you please provide a simple polygon that has your mentioned problem above? I will then create a testcase.

[SamuelReithmeir]

The polygon from #4 is the H seen in my screenshots, it has exactly this issue. And it is not self intersecting. It is just that when the three marked points are being tested for direction, that they are theoretically in one line (when projected onto the plane) but actually have a slight angle perpendiclat to the normal (because the middle point seems to be slightly behind the plane)

[NMO13]

I just checked #4, it is self-intersecting.

Check that:

 var points = new List<Vector3m>()
 {
new Vector3m(7197, -131, -6003),
 new Vector3m(7197, 131, -6003),
 new Vector3m(7103, 131, -6115),
 new Vector3m(7103, 145, -6115),
 //new Vector3m(7296, 145, -5884),
 //new Vector3m(7296, 131, -5884),
 //new Vector3m(7202, 131, -5996),
 //new Vector3m(7202, -131, -5996),
 //new Vector3m(7296, -131, -5884),
 //new Vector3m(7296, -145, -5884),
 //new Vector3m(7103, -145, -6115),
 //new Vector3m(7103, -131, -6115)
    };

 var A = points[0];
 var B = points[1];
 var C = points[2];
 var D = points[3];

 var isCoplanar = (D - A).Dot((B - A).Cross(C - A));
 Console.WriteLine(isCoplanar); // outputs 0

Your suggested fix won't help in this case either, because the sign of the dot product stays -1.

[SamuelReithmeir]

I just tested it again and you are right that my change does not fix the problem concerning the problem on hand in #4. First of all i beliefe that the given polygon is correct, accoring to my understanding of "self intersection" at least. When we look at the entire polygon: We see that there are no intersecting lines, but when we look at the first 4 vertices alone, of course they intersect:

Then i tried to run earclipper with my change to GetDirection and you were right, it is missing one triangle:

Which i find interesting, because this problem shoudl occur on all 4 quadrants, and not just bottom right.

[SamuelReithmeir]

Btw here is some python code i wrote to visualize this stuff:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot_3d_polygon(polygons):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    for vertices in polygons:
        # Extract x, y, and z coordinates from vertices
        x_coords = [vertex[0] for vertex in vertices]
        y_coords = [vertex[1] for vertex in vertices]
        z_coords = [vertex[2] for vertex in vertices]

        # Connect the vertices to form the polygon
        ax.plot(x_coords + [x_coords[0]], y_coords + [y_coords[0]], z_coords + [z_coords[0]], color='b')
        ax.scatter(x_coords, y_coords, z_coords, color='r')

    ax.scatter(x_coords, y_coords, z_coords, color='r')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('3D Polygon')
    # ax.set_xlim([0, 0.94248])
    # ax.set_ylim([0, 0.12])
    # ax.set_zlim([0, 1.8753])

    plt.show()

# Example vertices of a polygon
polygons = [
    # np.array([[0.942475438117981, 0, 2],
    #          [0.942475438117981, 0.1199876144528389, 2],
    #          [0.04999000206589699, 0.1199876144528389, 2],
    #          [0.04999000206589699, 0.06931385397911072, 2],
    #          [0.043292585760354996, 0.05998999997973442, 2],
    #          [0, 0.05998999997973442, 2],
    #          [0, 0, 2]]),
    # np.array([
    #    [7197.2792162972564, -131, -6003.2649404376552],
    #    [7197.2792162972564, 131, -6003.2649404376552],
    #    [7103.9723400461444, 131, -6115.2331919389362],
    #    [7103.9723400461444, 145, -6115.2331919389362],
    #    #[7296.0276599455174, 145, -5884.7668080597978],
    #    #[7296.0276599455174, 131, -5884.7668080597978],
    #    #[7202.7207836944053, 131, -5996.7350595610787],
    #    #[7202.7207836944053, -131, -5996.7350595610787],
    #    #[7296.0276599455174, -131, -5884.7668080597978],
    #    #[7296.0276599455174, -145, -5884.7668080597978],
    #    #[7103.9723400461444, -145, -6115.2331919389362],
    #    #[7103.9723400461444, -131, -6115.2331919389362]
    #          ]),
    # np.array([
    #    [7197.2792162972564, 131, -6003.2649404376552],
    #    [72#02.7207836944053, 131, -5996.7350595610787],
    #    [7103.9723400461444, 131, -6115.2331919389362],
    #         ]),
    np.array([[7197.279216297256, 131, -6003.264940437655], [7103.972340046144, 145, -6115.233191938936],
              [7296.027659945517, 145, -5884.766808059798], ]),
    np.array([[7296.027659945517, 145, -5884.766808059798], [7296.027659945517, 131, -5884.766808059798],
              [7202.720783694405, 131, -5996.735059561079], ]),
    np.array([[7197.279216297256, 131, -6003.264940437655], [7296.027659945517, 145, -5884.766808059798],
              [7202.720783694405, 131, -5996.735059561079], ]),
    np.array([[7197.279216297256, -131, -6003.264940437655], [7197.279216297256, 131, -6003.264940437655],
              [7202.720783694405, 131, -5996.735059561079], ]),
    np.array([[7197.279216297256, -131, -6003.264940437655], [7202.720783694405, 131, -5996.735059561079],
              [7202.720783694405, -131, -5996.735059561079], ]),
    np.array([[7202.720783694405, -131, -5996.735059561079], [7296.027659945517, -131, -5884.766808059798],
              [7296.027659945517, -145, -5884.766808059798], ]),
    np.array([[7197.279216297256, -131, -6003.264940437655], [7202.720783694405, -131, -5996.735059561079],
              [7296.027659945517, -145, -5884.766808059798], ]),
    np.array([[7103.972340046144, -131, -6115.233191938936], [7197.279216297256, -131, -6003.264940437655],
              [7296.027659945517, -145, -5884.766808059798], ]),
    np.array([[7103.972340046144, -145, -6115.233191938936], [7103.972340046144, -131, -6115.233191938936],
              [7296.027659945517, -145, -5884.766808059798], ])

]

plot_3d_polygon(polygons)

[NMO13]

Yes you are right, its not self intersecting if you take all points. Thanks for pointing that out.

Regarding the core problem, inaccurate input brings multiple problems in my opinion:

1. If the input points are not exact, self intersections might occur. This brakes this algorithm for sure as I have already shown.
2. If the input points are not self intersecting but are not on the same plane, then non-convex corners might be considered as convex and vice versa. This might break the algorithm, depending which direction the points are shifted.

Your fix changes the result of the IsConvex() function in the case that the polygon normal and triangle normal are normal to each other. So when your dot product return 0. Before, isConvex() returned false, now it returns true. This case is very odd anyway. If the points are "a bit" off then the current implementation might work. (I added a testcase for that). But thats also not guaranteed because that depends on which direction the error shifts the vertices.

So long story short: I highly suppose that your fix might work in this case but will break in other cases.

[SamuelReithmeir]

You are right, this change is not the real fix. I created these methods:

public static List<Vector3m> GetCoplanarMapping(List<Vector3m> points,
    out Dictionary<Vector3m, Vector3m> mappingDict)
{
  var normals = points.Select((x, i) => (points[(i - 1 + points.Count) % points.Count] - x)
      .Cross(points[(i + 1) % points.Count] - x)).ToList();
  var directionGroups = normals.GroupBy(x => x, new DirectionEqualComparer()).ToList();
  var dominantDirection = directionGroups.OrderByDescending(x => x.Count()).First().Key;
  var pointOnPlane = normals.Zip(points, (x, y) => new { x, y }).First(x => x.x.Equals(dominantDirection)).y;
  var mappedPoints = points
      .Select(x => x.ProjectOntoPlane(dominantDirection, pointOnPlane)).ToList();
  mappingDict = points
      .Zip(mappedPoints, (x, y) => new { x, y })
      .Where(arg => !arg.x.Equals(arg.y))
      .ToDictionary(arg => arg.x, arg => arg.y);
  return mappedPoints;
}

public static List<Vector3m> RevertCoplanarityMapping(List<Vector3m> mappedPoints,
    Dictionary<Vector3m, Vector3m> mapping)
{
    var reverseMapping = mapping.ToDictionary(kvp => kvp.Value, kvp => kvp.Key);
    return mappedPoints.Select(p => reverseMapping.TryGetValue(p, out var mapped) ? mapped : p).ToList();
}

it checks for the "dominant normal" by grouping the normals based on direction and taking the most occuring it should not matter if it finds the cw or ccwnormal, as they both describe the same plane. it then takes one of the points that had this normal and uses these to describe a plane. it then projects all points on this plane and returns a dictionary to revert the mapping

it can be used like this:

// example 5
points = new List<Vector3m>()
{
    new Vector3m(7197, -131, -6003),
    new Vector3m(7197, 131, -6003),
    new Vector3m(7103, 131, -6115),
    new Vector3m(7103, 145, -6115),
    new Vector3m(7296, 145, -5884),
    new Vector3m(7296, 131, -5884),
    new Vector3m(7202, 131, -5996),
    new Vector3m(7202, -131, -5996),
    new Vector3m(7296, -131, -5884),
    new Vector3m(7296, -145, -5884),
    new Vector3m(7103, -145, -6115),
    new Vector3m(7103, -131, -6115)
};
points= EarClipping.GetCoplanarMapping(points,out var reverseMapping);
earClipping = new EarClipping();
earClipping.SetPoints(points);
earClipping.Triangulate();
res = earClipping.Result;
res=EarClipping.RevertCoplanarityMapping(res,reverseMapping);
PrintTriangles(res);

i decided against automatically calling it inside setpoints, a it should not be used in cases that are already coplaner because it will do unneccessary calculation

[SamuelReithmeir]

you can see it working correctly here: