If I want to create a mesh from the sites can I use edges and assume they are connected, i.e. iterate through edges and add each starting point to a list of vertices ?

Not sure where to put these comments / requests for info ...

(Thanks for creating this library, btw ! Has saved me a lot of time and brain power)

Is there an easy way to determine which site within the rectangle belongs to an x,y co-ordinate ? All I can think to do is loop through each one and do a polygonal collision check using each edge but I figured Voronoi is sort-of designed to know which points belong inside (that's basically the definition of a Voronoi diagram).

Hi, I'm using this code to generate voronoi with 1036 points in my project. but the result is it't not a voronoi map and many polygons intersects. I'm using this function: jcv_diagram_generate(1036, points, &bounding_box, nullptr, &diagram) bounding_box is jcv_rect bounding_box = {{-45.8605, -653.969}, {746.861, 142.3}}; points is in the points.csv file. points.csv

diagram is definited as this: jcv_diagram diagram; memset(&diagram, 0, sizeof(jcv_diagram)); as a result, I'm using OGRGeometry library to render that diagram as above and I'm pretty sure it's not a problem of ogrGeometry library. I also confirmed that all points are within the bounding_box and there is no abnormal point.

can you help me to handle this problem? thank you very much @JCash @dgavedissian @williamleong

There is a bug in this library that occasionally causes the assertion at line 1134 (internal->numsites == 1) to fail. I do not know what the issue is, but I have attached a test program and data file that faithfully reproduce the error. Apologies for the size of the data file—this is an extremely infrequent error, and I've only been able to trigger it with sets this large.

Test file: mem.bin.zip

Test code:

#include <stdio.h>

#define JC_VORONOI_IMPLEMENTATION
#define JCV_REAL_TYPE double
#define JCV_ATAN2 atan2
#define JCV_FLT_MAX 1.7976931348623157E+308
#include "jc_voronoi.h"

int main() {
  jcv_point *points = (jcv_point *)malloc(9216 * sizeof(jcv_point));

  FILE *in = fopen("mem.bin", "rb");
  fread(points, sizeof(jcv_point), 9216, in);
  fclose(in);

  /* should start with
     (x = 40.232121213226684, y = 13.714460519854523)
     (x = 168.23212121322669, y = 13.714460519854523)
     (x = -87.767878786773309, y = 13.714460519854523)
     (x = 40.232121213226684, y = 29.714460519854523)
     (x = 40.232121213226684, y = -2.2855394801454771)
     (x = 168.23212121322669, y = 29.714460519854523)
     (x = -87.767878786773309, y = 29.714460519854523)
     (x = 168.23212121322669, y = -2.2855394801454771)
     (x = -87.767878786773309, y = -2.2855394801454771)
     (x = 123.81366674520085, y = 1.1403291016984274)
     */
  for (unsigned int i = 0; i < 10; i++) {
    printf("(x = %.14f, y = %.14f)\n", points[i].x, points[i].y);
  }

  jcv_diagram diagram;
  memset(&diagram, 0, sizeof(jcv_diagram));

  jcv_rect rect = {{-128.0, -16.0}, {256.0, 32.0}};

  jcv_diagram_generate(9216, points, &rect, &diagram);

  return 0;
}

Example use:

> unzip mem.bin.zip
Archive:  mem.bin.zip
  inflating: mem.bin 
> clang -lm test.c
> ./a.out
(x = 40.23212121322668, y = 13.71446051985452)
(x = 168.23212121322669, y = 13.71446051985452)
(x = -87.76787878677331, y = 13.71446051985452)
(x = 40.23212121322668, y = 29.71446051985452)
(x = 40.23212121322668, y = -2.28553948014548)
(x = 168.23212121322669, y = 29.71446051985452)
(x = -87.76787878677331, y = 29.71446051985452)
(x = 168.23212121322669, y = -2.28553948014548)
(x = -87.76787878677331, y = -2.28553948014548)
(x = 123.81366674520085, y = 1.14032910169843)
a.out: ./jc_voronoi.h:1134: void jcv_fillgaps(jcv_diagram *): Assertion `internal->numsites == 1' failed.
[1]    21138 abort (core dumped)  ./a.out

One more issue :/

Because I want to create a mesh, I need a list of unique vertices. These are stored as x,y points. Any idea how to do this ?

My naive algorithm would go as follows:

• Loop through all sites and their edges and add them all together (i.e. 5 sites each with, say, 5 edges gives a maximum of 25 vertices)
• Now we have a maximum. Allocate an integer for each (to be used as boolean). Set all to 1.
• Loop through again, this time with an inner loop starting from the outer loop +1
• In the inner loop, check if the outer loop position as the same as inner. If so, set the value of the integer to 0 (not unique).

Now we have a list of unique vertices ...

Not very elegant. Any better ideas ?

Also - can I just use an equals check on the positions (x1==x2) even though they are floats ?

It is very rare, and saw it only once, but I was able to hit this assert:

jc_voronoi.h:887: int jcv_pq_push(jcv_priorityqueue *, void *): Assertion `pq->numitems < pq->maxnumitems' failed.

Unfortunately, I didn't get a core dump to investigate further.

I feed it between 2 and 16 points, with no duplicates.

From now on, I'll run my game in a debugger, so I can catch it and report back.

First of all, thank you very much for your great work!

We're are using your implementation in a RoboCup soccer league and believe to have encountered a bug.

The points are the player's positions {-4050,0}, {-3300,500}, {0,1000}, {0,-1000}, {2250,0}. The bounding box to clip against is the field's corners {-4500,-3000}, {4500,3000}.

When iterating over the graph edges of the last point, the top right ({4500, 3000}) and bottom right ({4500,-3000}) corners of the field are missing entirely. Interestingly, changing the point {2250,0} to {2250,1} will fix the issue and the Voronoi diagram is constructed correctly.

Please find my screenshot attached.

Any help would be much appreciated. I'd be happy to submit a PR fixing the bug when found.

Hey, congrats on the great library.

I am trying to generate a voronoi diagram inside one of the cells of another voronoi diagram. I am converting the bounding polygon to a jcv_polygon and then making a clipper and setting it as a ctx.

The generation asserts at line 218 of jc_voronoi_clip.h at this assertion assert(min_edge >= 0);

Can you please help me understand what is happening? Thanks

when i iterate through a site's halfedges and look at their neighbors, they do not correspond to actual neighbors. and even if i look at the half edge's corresponding edge, its two sites are not neighboring in the diagram.

When using jc_voronoi for a natural neighbor interpolation problem, I was testing jc_voronoi to make sure it was getting linked correctly by giving it a square of points at {0,0}, {val, 0}, {0, val}, {-val, 0}, and {0, -val}.

There appears to be a bug when val = 2 in this case. I went through and was checking the number of edges that each site said it had, and the {0,0} site is always supposed to say 4. When val = 2, it says it has two. This doesn't seem to happen when I rotate or shift the square. This code runs through these cases and some other vals other than 2

#define JC_VORONOI_IMPLEMENTATION
#include "jc_voronoi.h"
#include <cmath>
#include <stdlib.h>
#include <iostream>

void printSquare( float val, int mode )
{
  std::cout << "\nStart Of Test\n";
  int pointCount = 5;
  jcv_point list[pointCount];
  // list[0] = {  0,   0 };

  if( mode == 1)
  {
    //45 degree rotation on unit circle
    list[0].x = 0;
    list[0].y = 0;
    float valUpdated = val/(std::sqrt(2));
    // list[1] = { valUpdated, valUpdated };
    list[1].x = valUpdated;
    list[1].y = valUpdated;
    // list[2] = { -valUpdated, valUpdated };
    list[2].x = -valUpdated;
    list[2].y = valUpdated;
    // list[3] = { valUpdated, -valUpdated };
    list[3].x = valUpdated;
    list[3].y = -valUpdated;
    // list[4] = { -valUpdated, -valUpdated };
    list[4].x = -valUpdated;
    list[4].y = -valUpdated;
  }else if ( mode == 2)
  {
    //offset from origin
    float xOffset = std::rand() % 100 + 1;
    float yOffset = std::rand() % 100 + 1;

    list[0].x = 0 + xOffset;
    list[0].y = 0 + yOffset;

    // list[1] = { val + xOffset , 0 + yOffset };
    list[1].x = val + xOffset;
    list[1].y = 0 + yOffset;
    // list[2] = { 0 + xOffset, val + yOffset  };
    list[2].x = 0 + xOffset;
    list[2].y = val + yOffset;
    // list[3] = { -val + xOffset, 0 + yOffset  };
    list[3].x = -val + xOffset;
    list[3].y = 0 + yOffset;
    // list[4] = { 0 + xOffset, -val + yOffset };
    list[4].x = 0 + xOffset;
    list[4].y = -val + yOffset;
  }else
  {
    list[0].x = 0;
    list[0].y = 0;
    // list[1] = { val, 0 };
    list[1].x = val;
    list[1].y = 0;
    // list[2] = { 0, val };
    list[2].x = 0;
    list[2].y = val;
    // list[3] = { -val, 0 };
    list[3].x = -val;
    list[3].y = 0;
    // list[4] = { 0, -val };
    list[4].x = 0;
    list[4].y = -val;
  }
  jcv_diagram diagram;
  memset(&diagram, 0, sizeof(jcv_diagram));
  jcv_diagram_generate(pointCount, list, nullptr, &diagram);

  const jcv_site *sites = jcv_diagram_get_sites(&diagram);
  for( int i = 0; i < diagram.numsites; ++i )
  {
      const jcv_site* site = &sites[i];
      std::cout << "\nAt site index " << site->index << " with position (" << site->p.x << ", " << site->p.y << ")\n";
      const jcv_graphedge* e = site->edges;
      int edgeCount = 0;
			while( e )
			{
			// 	std::cout << "\nSite pos  = " << site->p.x << ", " << site->p.y << "\n";
      //   std::cout << "Edge 1 pos = " << e->pos[0].x << ", " << e->pos[0].y << "\n";
      //   std::cout << "Edge 2 pos = " << e->pos[1].x << ", " << e->pos[1].y << "\n";
        edgeCount++;
				e = e->next;
			}
      // std::cout << "\nSite pos  = " << site->p.x << ", " << site->p.y << "\n";
      std::cout << "Number of edges is " << edgeCount << "\n";
      if(site->p.x == 0 && site->p.y == 0)
      {
        if(edgeCount != 4)
        {
          std::cout << "At (0,0) the cell does not have 4 edges!!!!!!!!!!!\n";
        }else{
          std::cout << "As expected, the cell at (0,0) has 4 edges\n";
        }
      }
  }
  std::cout << "Done printing sites and edge counts for this test\n\n\n";
  jcv_diagram_free( &diagram );

}

int main( int argc, const char *argv[])
{

  float eps = std::numeric_limits<float>::epsilon();
  std::cout << "\nDemonstration of bug at 2\n\n";

  printSquare( 1, 0 );
  printSquare( 3, 0 );
  std::cout << "\nThis is the buggy one\n";
  printSquare( 2, 0 ); //issue is here
  std::cout << "\nEnd of the buggy one\n";
  printSquare( 2, 1 );
  printSquare( 2, 2 );
  printSquare( 2+2*eps, 0 );
  printSquare( 2-2*eps, 0 );

  std::cout << "\n\n\nEnd of Tests\n";

  return 0;

}

The result from this is

Demonstration of bug at 2


Start Of Test

At site index 4 with position (0, -1)
Number of edges is 4

At site index 3 with position (-1, 0)
Number of edges is 4

At site index 0 with position (0, 0)
Number of edges is 4
As expected, the cell at (0,0) has 4 edges

At site index 1 with position (1, 0)
Number of edges is 4

At site index 2 with position (0, 1)
Number of edges is 4
Done printing sites and edge counts for this test



Start Of Test

At site index 4 with position (0, -3)
Number of edges is 4

At site index 3 with position (-3, 0)
Number of edges is 4

At site index 0 with position (0, 0)
Number of edges is 4
As expected, the cell at (0,0) has 4 edges

At site index 1 with position (3, 0)
Number of edges is 4

At site index 2 with position (0, 3)
Number of edges is 4
Done printing sites and edge counts for this test



This is the buggy one

Start Of Test

At site index 4 with position (0, -2)
Number of edges is 3

At site index 3 with position (-2, 0)
Number of edges is 2

At site index 0 with position (0, 0)
Number of edges is 2
At (0,0) the cell does not have 4 edges!!!!!!!!!!!

At site index 1 with position (2, 0)
Number of edges is 4

At site index 2 with position (0, 2)
Number of edges is 4
Done printing sites and edge counts for this test



End of the buggy one

Start Of Test

At site index 4 with position (-1.41421, -1.41421)
Number of edges is 5

At site index 3 with position (1.41421, -1.41421)
Number of edges is 5

At site index 0 with position (0, 0)
Number of edges is 4
As expected, the cell at (0,0) has 4 edges

At site index 2 with position (-1.41421, 1.41421)
Number of edges is 5

At site index 1 with position (1.41421, 1.41421)
Number of edges is 5
Done printing sites and edge counts for this test



Start Of Test

At site index 4 with position (8, 48)
Number of edges is 4

At site index 3 with position (6, 50)
Number of edges is 4

At site index 0 with position (8, 50)
Number of edges is 4

At site index 1 with position (10, 50)
Number of edges is 4

At site index 2 with position (8, 52)
Number of edges is 4
Done printing sites and edge counts for this test



Start Of Test

At site index 4 with position (0, -2)
Number of edges is 4

At site index 3 with position (-2, 0)
Number of edges is 4

At site index 0 with position (0, 0)
Number of edges is 4
As expected, the cell at (0,0) has 4 edges

At site index 1 with position (2, 0)
Number of edges is 4

At site index 2 with position (0, 2)
Number of edges is 4
Done printing sites and edge counts for this test



Start Of Test

At site index 4 with position (0, -2)
Number of edges is 4

At site index 3 with position (-2, 0)
Number of edges is 4

At site index 0 with position (0, 0)
Number of edges is 4
As expected, the cell at (0,0) has 4 edges

At site index 1 with position (2, 0)
Number of edges is 4

At site index 2 with position (0, 2)
Number of edges is 4
Done printing sites and edge counts for this test





End of Tests

I am not sure why this is happening, especially since it works when val = 2 + 2 *numeric_limits::epsilon.

I am hoping someone more familiar with the code can find it the reason

Thanks for the great library!

It would be helpful to provide a simple example program of usage. The main.c is a little verbose and has a lot of cruft dealing with coloring and saving to an image file.

Here is an example program:

// to compile:
//
// gcc jc_voronoi_example.c -I../src -o jc_voronoi_example  -lm
//

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define JC_VORONOI_IMPLEMENTATION
// If you wish to use doubles
//#define JCV_REAL_TYPE double
//#define JCV_FABS fabs
//#define JCV_ATAN2 atan2
#include "jc_voronoi.h"

#define NPOINT 10

int main(int argc, char **argv) {
  int i;
  jcv_rect bounding_box = { { 0.0, 0.0 }, { 1.0, 1.0 } };
  jcv_diagram diagram;
  jcv_point points[NPOINT];
  const jcv_site *sites;
  jcv_graphedge *graph_edge;

  memset(&diagram, 0, sizeof(jcv_diagram));

  srand(0);
  for (i=0; i<NPOINT; i++) {
    points[i].x = (float)rand()/(1.0 + RAND_MAX);
    points[i].y = (float)rand()/(1.0 + RAND_MAX);
  }

  // seed sites
  //
  for (i=0; i<NPOINT; i++) {
    printf("%f %f\n\n", points[i].x, points[i].y);
  }

  jcv_diagram_generate(NPOINT, (const jcv_point *)points, &bounding_box, &diagram);

  sites = jcv_diagram_get_sites(&diagram);
  for (i=0; i<diagram.numsites; i++) {

    graph_edge = sites[i].edges;
    while (graph_edge) {
      printf("%f %f\n", graph_edge->pos[0].x, graph_edge->pos[0].y);
      printf("%f %f\n", graph_edge->pos[1].x, graph_edge->pos[1].y);
      printf("\n");

      graph_edge = graph_edge->next;
    }

  }

  jcv_diagram_free(&diagram);
}

The generated seed sites:

0.840188 0.394383
0.783099 0.798440
0.911647 0.197551
0.335223 0.768230
0.277775 0.553970
0.477397 0.628871
0.364784 0.513401
0.952230 0.916195
0.635712 0.717297
0.141603 0.606969

Visualizing the seed sites and edges with gnuplot:

This assumes you're in an example subdirectory, say, to compile and run. The number of points is hard coded and it creates the points randomly in a 1x1 box, but the above example gets across clearly how to set up, use and get useful information out of the library. Presumably this 'double covers' the half edges but for illustration purposes I don't think that's a problem.

I'd be happy to submit a pull request if that's helpful.

I was running simulations using your library, and I found some errors. Here's an example:

jcv_diagram d{};

// example that fails using double
//jcv_point points[]
//{
//	{888.19238281250000, 377.82843017578125},
//	{914.00000000000000, 341.00000000000000},
//};

// example that fails using the standard float version of the library
jcv_point points[]
	{
		{883.382263f, 340.749908f},
		{850.622253f, 378.323486f},
	};

jcv_rect rect;
rect.min = { 600, 250 };
rect.max = { 1000, 650 };
const auto count = sizeof(points) / sizeof(*points);
jcv_diagram_generate(count, points, &rect, 0, &d);
const jcv_site* sites = jcv_diagram_get_sites(&d);
for (int i = 0; i != d.numsites; i++)
{
	const jcv_site* site = &sites[i];
	const jcv_graphedge* e = site->edges;
	int cnt = 0;
	while (e)
	{
		cnt++;
		e = e->next;
	}
	std::cout << cnt << " sides\n";
}

/* output:
4 sides
2 sides
Obviously wrong. One can clearly see the voronoi should have 5 sides and 3 sides in each cell
*/

I believe it is associated with this issue, but this one is easier to reproduce because of only 2 vertices.

Now, to get access to a separate site by its index, you need to either make a separate array sorted by site indexes, or use a chart search. Maybe there are other ways that I can't see?

Hi, is there any simple method to enable the library to support clip by nonconvex boundary? I tested when use concave polygon as clip polygon, the result is blank.

Hi. first of all, thank you so much for your great OSS. This library is very useful and easy to use. Thanks!!

I might find an issue of this OSS. When I input a certain data into jcv_diagram_generate function, an assertion was reported:

Assertion failed: (internal->numsites == 1), function jcv_fillgaps, file src/jc_voronoi.h, line 1143.

This repo is a fork of your voronoi repo and I added a test code for reproducing the issue. https://github.com/AtsushiSakai/voronoi If you have time, please take a look the file test/assert_text.c file. https://github.com/AtsushiSakai/voronoi/blob/master/test/assert_test.c test/invalid_data.h includes input x-y point data for the issue.

#38 might be a same issue.

In my knowledge, the time complexity of Fortune's Sweepline algorithm is O(n log n). This algorithm uses a balanced binary search tree(BBST) to insert/delete parabola and to do a binary search in O(log n).

I found that this code uses a linked list, instead of BBST. The linked list makes this code O(n^2), and it means this code will take lots of time to calculate Voronoi Diagram in specific inputs.

Generator of test input is here.

#!/usr/bin/ruby
n = 1000000
n.times do |i|
	puts "%d %d" % [(i+1), -(i+1)]
	puts "%d %d" % [-(i+1), -(i+1)]
end

You can check that your program is almost stopped at this part or this part.

Actually, an implementation used linked list will work well in the average case.