abellgithub/delaunator-cpp

A really fast C++ library for Delaunay triangulation of 2D points

delaunator-cpp's Issues

Incorrect bounding box height calculation

While porting this code to Rust I noticed that when the bounding box height is being calculated, the code does max_x - min_y instead of max_y - min_y:

delaunator-cpp/include/delaunator.cpp

Line 223 in 2e12c1f

double height = max_x - min_y;

I assume this is incorrect, right? I would submit a pull request, but since all tests are passing currently, I am unsure if I should.

Support using float instead of double

Hi,

It looks like the library currently only supports double. But for platforms with only a single-precision FPU it would be nice to allow using float for performance reasons. From usability point of view, it would be great to make the classes a template so we could pass either float or double as template argument. Alternatively it would also be OK to replace all occurrences of double with an own typedef and e.g. allow switching to float with a define passed to the compiler, or something like that.

Are there any plans in this direction, or would this be totally out of scope of this library?

Allow use as header-only

triangulation fails with large grids of points (solved?)

Hello,

I'm considering using delaunator for speeding up GDL's TRIANGULATE function, that badly needs it!

I've made the changes in a branch of my fork of GDL
It happens that GDL's TRIANGULATE function is frequently used on regular grids of points (ex: pixels positions on a camera). Certainly not the primary intended use of a triangulation, but "it is the use that makes the function" 😄 . A beta-tester reported failures, so I investigated.

I found that weird triangles were produced in a manner similar to issue #8 when the number of points was large, example: 3017x2018 points in a normalized box of size ~1 (that is, each x and y values are divided by the absolute max of the x and y range. This trick is intended to make the (mathematical) triangulation independent of the magnitude of the x and y values, be it kilometers or microns -- a trick necessary when using legacy triangulation programs of the 60's...)

The following image is the triangulation obtained on a regular grid 3017x2018 , each of the 7340942 triangles painted with a random color.

The movie of the triangles being painted is something to be seen. Everything is fine until the sweeping hull is near the border and then it goes awry.

I have tracked back the culprit to the size of the hash_key. I changed the size of the haskey by adding 3000 entries, and it worked.

I am afraid delaunator is unfortunately not robust to non-square high number of grid-aligned points:

• The size m_hash_size = static_cast<std::size_t>(std::ceil(std::sqrt(n))); is too small for large "non-square" problems. At least with the way the hashkey is computed. Of course, speed is related to the hash size so it is important to maintain it small...
• Possibly, there is an adverse effect of the dependency of pseudo_angle() to the actual values of dx and dy. It will make the angle comparison uncertain for close points (very small values of dx-dy)

This drove me to try smaller non-square arrays and I found that wrong triangles were possible for rather small values, the smallest I found was 1602x408.

By multiplying the size of the hastable by the golden ratio (after all this seems related to the ratio of the problem's rectangle to the implicit squareness of sqrt()), all the tests I made worked beautifully, with almost no speed loss.

Logic error in bounding rectangle initialization.

std::numeric_limits::min() is used to initialize maximum value variables that should grow to the highest encountered value. It represents the smallest positive normal double, so this value can never grow negative (e. g. for only negative inputs). Use std::numeric_limits::lowest() istead.

Lines 199 and 200 of delaunator.cpp:
double max_x = (std::numeric_limits<double>::min)();
double max_y = (std::numeric_limits<double>::min)();

should be:
double max_x = std::numeric_limits<double>::lowest();
double max_y = std::numeric_limits<double>::lowest();

sort is slow on GCC

GCC 7.5+ copies the comparator in the process of sorting. Since we cache distances in a vector in the comparator, we have to copy those distances and it's very slow if the numbers are large.

Iterating hull in wrong direction leads to infinite loop

See: mapbox/delaunator#49

I don't think there's any way a race condition is involved, but whatever...

Has anyone been able to write the delaunator-js demo in C++?

The demo here lists a bunch of code that works with the JS version of this lib, but because of what I'm assuming is type differences across the implementations, porting it isn't straight forward. Has anyone been able to render the voronoi edges like the js demo does?

Cache distance from center when sorting IDs

INVALID_INDEX problem

Hello,
I don't have any other means to contact you, so I resorted to creating new issue. I am trying to figure out what am I doing wrong - I have ~1000 points defined and I am trying to run triangulation on them. Then i am iterating halfedges for some further processing but some of them are pointing to "INVALID_INDEX" (std::numeric_limitsstd::size_t::max()). It occours close the end of halfedges vector. What can cause this behavior? Positions of points that I am feeding into Delaunator are in magnitude of 100-s.

Thanks in advance and sorry for this form of contact.

Triangulation error on tilted grids.

I am processing point clouds around the size of 100 k to 1.000 k Points. The points' outer hull is a near perfect rectangle and their distribution is grid like. In pre-processing the point cloud gets slightly tilted towards a plane for error compensation.

While I have not yet observed the issue below in our raw data, occasionally there was a number of big triangles (hundreds in a cloud of a million points) spanning across the whole pre-processed (tilted) point cloud. To reproduce this, I wrote a test. It generates simple square grid, rotates it by a few degrees and then calls delaunator to triangulate.
Validation is done knowing the grid consists of 1-by-1 unit squares, thus the triangulation should consist of right triangles with two unit-length edges and one sqrt(2) unit-length edge. Degenerate triangles on the outer hull (due to point rotation rounding errors) are discarded by checking for very small triangle areas.

The smallest point set I could find to let the test fail is the one supplied below.
Compiler is Visual C++ 14.0, tested platforms are Win32 and x64. Testing framework is googletest, but you should have no difficulties to modify it for Catch (I did not because I do not fully understand the test structure in your repo).

A plot of the triangulated grid is attached here:

#include "gtest/gtest.h"
#include "delaunator/delaunator.hpp"

struct Point2D {
    double x, y;

    Point2D()
        : x(std::numeric_limits<double>::quiet_NaN())
        , y(std::numeric_limits<double>::quiet_NaN())
    {}

    Point2D(double x, double y)
        : x(x)
        , y(y)
    {}
};

Point2D rotate(Point2D const& p_in, Point2D const& center, double const& angle_deg) {
    const double pi        = 3.14159265358979323846;
    const double angle_rad = angle_deg * pi / 180.0;
    const double s         = sin(angle_rad);
    const double c         = cos(angle_rad);
    Point2D p;

    // translate point back to origin:
    p.x = p_in.x - center.x;
    p.y = p_in.y - center.y;

    // rotate point
    double tmp = p.x * c - p.y * s;
    p.y = p.x * s + p.y * c;
    p.x = tmp;

    // translate point back:
    p.x += center.x;
    p.y += center.y;

    return p;
}

double distance(Point2D const& p1, Point2D const& p2) {
    return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}

bool approx(double v1, double v2) {
    const double eps = 0.01;
    return fabs(v1 - v2) < eps;
}

TEST(External_Lib_Delaunator_Test, rotate_grid_6_8_8_207_0) {
    const size_t  grid_size = 6;
    const Point2D center(8, 8);
    const double  angle = 207.0;

    std::vector<double> points;

    // Generate rotated point grid of 1-by-1 squares.
    for (size_t x = 0; x < grid_size; ++x) {
        for (size_t y = 0; y < grid_size; ++y) {
            Point2D p = rotate(Point2D(x, y), center, angle);
            points.push_back(p.x);
            points.push_back(p.y);
        }
    }

    delaunator::Delaunator delaunator(points);

    // Validate.
    // Assume:
    // (1) since the grid is 1-by-1 square, that each triangle has to have two edges with length 1 and one edge wit length sqrt(2) - give or take epsilon.
    // (2) due to rounding errors during rotation the hull points of the grid are not in a numerically perfect line, thus degenerate triangles will appear on the outer edges. Their area mitght be rather small (< 0.01) compared to a regular grid triangle (== 0.5).
    for (size_t i = 0; i < delaunator.triangles.size(); i += 3) {
        Point2D p1(points[2 * delaunator.triangles[i]],     points[1 + 2 * delaunator.triangles[i]]);
        Point2D p2(points[2 * delaunator.triangles[i + 1]], points[1 + 2 * delaunator.triangles[i + 1]]);
        Point2D p3(points[2 * delaunator.triangles[i + 2]], points[1 + 2 * delaunator.triangles[i + 2]]);

        // Determine edge lengths and triangle area.
        double len_1 = distance(p1, p2);
        double len_2 = distance(p2, p3);
        double len_3 = distance(p3, p1);
        double area  = fabs((p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y)) / 2.0);

        // Check edge lengths according to (1).
        bool cond_1 = false;
        if (approx(len_1, sqrt(2.0))) {
            cond_1 = approx(len_2, 1.0) && approx(len_3, 1.0);
        }
        else if (approx(len_2, sqrt(2.0))) {
            cond_1 = approx(len_1, 1.0) && approx(len_3, 1.0);
        }
        else if (approx(len_3, sqrt(2.0))) {
            cond_1 = approx(len_1, 1.0) && approx(len_2, 1.0);
        }

        // Check triangle area according to (2).
        bool cond_2 = false;
        if (cond_1) {
            // Trianlge seems regular half 1-b-1 square. Just to be sure, check the area.
            cond_2 = approx(area, 0.5);
        }
        else {
            // Irregular triangle. Accept degenerate ones. (Assume they are on the outer egde, but do not check.)
            cond_2 = area < 0.01;
        }

        EXPECT_TRUE(cond_2) << "triangle = [" << p1.x << ", " << p1.y << "]; [" << p2.x << ", " << p2.y << "]; " << "[" << p3.x << ", " << p3.y << "]";
    }
}

Windows max-guarding missing

delaunator.cpp line 267 should be:

if (!(min_radius < (std::numeric_limits::max)())) {