JCash / voronoi

A C implementation for creating 2D voronoi diagrams
MIT License
632 stars 94 forks source link

Duplicate points for site #63

Closed williamleong closed 2 years ago

williamleong commented 2 years ago

Given 2 points, 144.15753183116175, 148.76233737400082 and -140.05365988451959, 143.19845406198522 with the polygon:

{  250,  200 },
{  250, -200 },
{ -250, -200 },
{ -250,  200 }

The pos[0].x and pos[0].y of the edges in site 0 is:

{ 0.994417, 200 }
{ 8.82505, -200 }
{ 250, -200 }
{ 250, 200 }
{ 0.994417, 200 }
{ 8.82505, -200 }
{ 250, -200 }
{ 250, 200 }

The first 4 points are correct, but they are duplicated, resulting in a total of 8 points.

williamleong commented 2 years ago

Manage to replicate the issue, it is something to do with double precision and it seems to work if I remove double precision defines below:

#define JCV_REAL_TYPE double
#define JCV_FLT_MAX std::numeric_limits<JCV_REAL_TYPE>::max()
#define JCV_ATAN2 atan2
#define JCV_SQRT sqrt
#define JCV_PI 3.141592653589793115997963468544185161590576171875

then the behaviour is correct and there will be 4 points. With the double precision defines, there are 8 points.

Minimal example below:

#include "doctest/doctest.h"

#include <iostream>
#include <memory>
#include <vector>
#include <cstring>
#include <limits>

#define JC_VORONOI_CLIP_IMPLEMENTATION
#define JC_VORONOI_IMPLEMENTATION

//Use double precision
#define JCV_REAL_TYPE double
#define JCV_FLT_MAX std::numeric_limits<JCV_REAL_TYPE>::max();
#define JCV_ATAN2 atan2
#define JCV_SQRT sqrt
#define JCV_PI 3.141592653589793115997963468544185161590576171875

#include "voronoi/src/jc_voronoi_clip.h"

TEST_CASE("TestDuplicateVertices")
{
    std::vector<jcv_point> polygonVector =
    {
        jcv_point{  250,  200 },
        jcv_point{ -250,  200 },
        jcv_point{ -250, -200 },
        jcv_point{  250, -200 }
    };

    std::vector<jcv_point> pointsVector =
    {
        jcv_point{ 144.15753183116175, 148.76233737400082},
        jcv_point{ -140.05365988451959, 143.19845406198522 }
    };

    std::vector<jcv_point> partition;

    //polygon
    jcv_clipping_polygon polygon;
    polygon.num_points = (int)polygonVector.size();
    polygon.points = polygonVector.data();

    //clipper
    jcv_clipper polygonclipper;
    polygonclipper.test_fn = jcv_clip_polygon_test_point;
    polygonclipper.clip_fn = jcv_clip_polygon_clip_edge;
    polygonclipper.fill_fn = jcv_clip_polygon_fill_gaps;
    polygonclipper.ctx = &polygon;

    jcv_rect rect;
    rect.min = { -250, -200 };
    rect.max = { 250, 200 };

    //Compute diagram
    jcv_diagram diagram;
    memset(&diagram, 0, sizeof(jcv_diagram));
    jcv_diagram_generate((int)pointsVector.size(), (const jcv_point*)pointsVector.data(), &rect, &polygonclipper, &diagram);

    //Get sites
    const jcv_site* sites = jcv_diagram_get_sites(&diagram);
    jcv_graphedge* currentEdge = nullptr;

    for(int i = 0; i < diagram.numsites; i++)
    {
        const jcv_site* currentSite = &sites[i];
        if (currentSite->index == 0)
        {
            currentEdge = currentSite->edges;
            break;
        }
    }

    size_t i = 0;

    while (currentEdge)
    {
        i++;
        std::cout << currentEdge->pos[0].x << "\t" << currentEdge->pos[0].y << std::endl;
        currentEdge = currentEdge->next;
    }

    REQUIRE(i == 4);

    jcv_diagram_free(&diagram);
} //TEST_CASE
williamleong commented 2 years ago
// check that we didn't accidentally add a duplicate (rare), then remove it
if( ge->next && ge->angle == ge->next->angle )
{
    if( jcv_point_eq( &ge->pos[0], &ge->next->pos[0] ) && jcv_point_eq( &ge->pos[1], &ge->next->pos[1] ) )
    {
        ge->next = ge->next->next; // Throw it away, they're so few anyways
    }
}

jcv_point_eq( &ge->pos[1], &ge->next->pos[1] ) is returning false, causing the edge not to be pruned:

ge->pos[1]
0.9944168342190074
200

ge->next->pos[1]
0.99441683421900784
200

I'm guessing it's a side effect of using double precision.

Possible solution:

static inline int jcv_point_eq( const jcv_point* pt1, const jcv_point* pt2 )
{
    return (fabs(pt1->y - pt2->y) < JCV_EPSILON) && (fabs(pt1->x - pt2->x) < JCV_EPSILON);
}

@JCash, what are your thoughts?