JCash / voronoi

A C implementation for creating 2D voronoi diagrams
MIT License
622 stars 92 forks source link

Incorrect number of edges in specific case #22

Closed archimedes4000 closed 5 years ago

archimedes4000 commented 6 years ago

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

david94133 commented 6 years ago

I'm pretty sure this is the same problem I ran into. When a site is on the boundary, the edges are not created properly. The open pull request I have fixes this issue; try applying it to your branch.

archimedes4000 commented 6 years ago

Sorry for the delay in responding, @david94133

This did not fix the issue I have run across.

JCash commented 5 years ago

Fixed in #27 (v0.5.0)