linuxlewis / tripy

Simple polygon triangulation algorithms in pure python
MIT License
64 stars 10 forks source link

Tripy-1.0.0 unable to process polygons supported in version 0.0.3 #3

Open lelabo-m opened 5 years ago

lelabo-m commented 5 years ago

Hi,

We discovered with some colleagues by using different version of tripy that it was not able anymore to process some polygons that were supported in version 0.0.3.

Here is an example to reproduce the case.

import shapely.wkt
import numpy as np
import tripy
import matplotlib.pyplot as plt

if __name__ == '__main__':
    polygon = shapely.wkt.loads("POLYGON ((2.3374803 48.87445059972346, 2.3374866 48.87444049972347, 2.3375329 48.87436569972348, 2.3375717 48.8743029997235, 2.3376836 48.87433219972348, 2.3378011 48.87436279972348, 2.3377983 48.87436709972348, 2.3377661 48.87441459972347, 2.337753 48.87441159972347, 2.3377382 48.87443419972345, 2.337751 48.87443809972346, 2.3377694 48.87443419972345, 2.3379195 48.87447389972345, 2.3380046 48.87471509972338, 2.3379528 48.87472209972339, 2.3379554 48.87472929972338, 2.337902 48.87473809972339, 2.3379056 48.87474939972338, 2.3377288 48.87477829972337, 2.3377241 48.87476549972337, 2.3376082 48.87478359972337, 2.3376066 48.87477919972339, 2.3375998 48.87478029972337, 2.337577099999999 48.87472049972339, 2.3374907 48.87473379972338, 2.3374961 48.87459949972342, 2.3374982 48.87454969972343, 2.3375389 48.87454889972344, 2.3375331 48.87447659972346, 2.3374803 48.87445059972346))")

    poly_coords = np.array(polygon.exterior.coords)
    tris = tripy.earclip(poly_coords)
    for t in np.array(tris):
        plt.plot(t[[0, 1, 2, 0], 0], t[[0, 1, 2, 0], 1])

    plt.plot(poly_coords[:, 0], poly_coords[:, 1], 'k-', lw=3)
    plt.show()

It produces the following images:

tripy-0.0.3 tripy-1.0.0
tripy-0.0.3 tripy-1.0.0

As you can see in the second images, tripy was not able to produce any triangles.

I was reading the PR #2 and saw that you introduced some changes to handle some edge cases.

Is this a side effect ? Should the library be able to process this polygon (or not) ? Should it at least produce some of the triangles ?

AlexLVuc commented 5 years ago

I seem to be having the same problem. tripy.earclip(poly_coords) returns an empty array on release 1.0.0. The recent PR introduced the use of EPSILON into the code which is incorporated into the earclip function.

if _is_ear(prev_point, point, next_point, polygon):
  ear_vertex.append(point)

....

def _is_ear(p1, p2, p3, polygon):
    ear = _contains_no_points(p1, p2, p3, polygon) and \
            _is_convex(p1, p2, p3) and \
            _triangle_area(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y) > 0
        _is_convex(p1, p2, p3) and \
        _triangle_area(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y) > 0
    return ear

def _contains_no_points(p1, p2, p3, polygon):
    for pn in polygon:
        if pn in (p1, p2, p3):
            continue
        elif _is_point_inside(pn, p1, p2, p3):
            return False
    return True
def _is_point_inside(p, a, b, c):
    area = _triangle_area(a.x, a.y, b.x, b.y, c.x, c.y)
    area1 = _triangle_area(p.x, p.y, b.x, b.y, c.x, c.y)
    area2 = _triangle_area(p.x, p.y, a.x, a.y, c.x, c.y)
    area3 = _triangle_area(p.x, p.y, a.x, a.y, b.x, b.y)
    # OLD CODE
    return area == sum([area1, area2, area3])
    # NEW CODE
    areadiff = abs(area - sum([area1, area2, area3])) < EPSILON
    return areadiff

I ran it over this set of points:

points = [[2808.892921157646, 3427.334099047921], [2173.127752759028, 3427.334099047922], [2173.127752759025, 2784.693474047921], [2616.518377759022, 2784.693474047921], [2616.518377759022, 3063.917404293587], [2359.893377759028, 3063.917404293591], [2359.893377759028, 3144.542404293586], [2616.518377759037, 3144.542404293587], [2616.518377759022, 3063.917404293587], [2616.518377759023, 2784.693474047921], [2808.892921157641, 2784.693474047921]]

And received multiple instances where the new _is_point_inside would return True and the old _is_point_inside would return False.

Because of this change, the _is_ear check returned True 0 times on the new code and 6 times on the old code, which gave me an actual output that I was looking for. ( Or a close enough output for what I am working on )

triangulation-v0.0.3 triangulation-v1.0.0

The top image is the image produced by version 0.0.3 of code, while the bottom image is the version produced by version 1.0.0.

Was the intended solution here to create a more robust comparison between the area of the triangle and the sum of the area of the smaller triangles? This check seems to break when provided numbers such as the points provided above.

MartinFalatic commented 5 years ago

Pardon my belated reply - I only saw this issue report recently.

The goal was to reduce the potential for error when it comes to certain narrow triangles that can result from earclipping. In certain cases these can lead to very inaccurate results in v0.0.3.

An example of a working case is:

poly = np.array([
    [4732, 2255],
    [3102, 1762],
    [   0,    0],
    [ 182, 1055],
    [4732, 2255],
])

Which leads to the expected outcome:

Barycenters vs v0 0 3 larger values

But scale it all down tenfold and you get:

poly = np.array([
    [473.2, 225.5],
    [310.2, 176.2],
    [  0.0,   0.0],
    [ 18.2, 105.5],
    [473.2, 225.5],
])

Which leads to erroneous output:

Barycenters vs v0 0 3 smaller values

It was found that the strict comparison in v0.0.3 would lead to this problem in cases like these, leading to triangles being excluded that shouldn't be.

I thought that comparing within a small range epsilon would solve the problem, and indeed in all the tests I ran with the v1.0.0 code this worked. However, clearly it doesn't work in all real world cases (I learned that generating interesting, non-trivial polygons is itself quite non-trivial).

You'll notice both of the examples above are compared to a barycenter-based method of triangulation not currently used in Tripy but which is based on https://github.com/lionfish0/earclip. That method appears to be immune to the edge cases that earclipping can run into - however, it is also 3x-6x slower. (I suspect one could improve that with numpy, though in the end an optimized version may still be slower than earclipping.)

I hoped that we could keep the speed of the present method instead of adding the complexity of barycenters, but if the edge cases in earclipping cannot be contained systematically it may be good to move to that instead.

MartinFalatic commented 5 years ago

To add, I did try the barycenter method on the examples you both provided above, which gave satisfactory results for both. Modifying v1.0.0 of Tripy by increasing epsilon had a similar effect, but systematically adapting epsilon may not be tractable.

ex1

ex2

nofaceinbook commented 4 years ago

One big difference is that https://github.com/lionfish0/earclip. defines epsilon as 0.0000001 which is far higher than the value in tripy. Setting epsilon to higher value also in tripy seems to work much better.