Point in polygon algorithm creates spurious lines

[danielwboyce]

Discovered at the end of fixing NanoComp/meep#864. The version of node_in_or_on_polygon (in geom.c) that was added with commit da42dcf creates small lines in polygons that shouldn't be there. For example, in these images show the error (with the correct form of the polygon on the left and the incorrect on the right):

More examples of this error can be seen in this imgur gallery. The lines only seem to form to the left of a vertex in the polygon. The even-odd algorithm that was used (from this paper) accounts for this sort of problem, and running through the algorithm on paper with problematic polygons doesn't show the lines appearing. Perhaps the occurrence of the lines is caused by a floating-point problem?

This problem could potentially be addressed by using a winding number algorithm instead of an even-odd algorithm.

[stevengj]

I would create a little .c program that creates this 12-side polygon and calls point_in_object for one of the points on those white lines (connecting two vertices), where presumably it is returning false when it should be returning true. If you can reproduce the problem this way, then it will be easy to track down — just step through the computation line-by-line with a debugger, since the computation isn't that long.

Once it is fixed, you can also include this minimal test case in make check.

[danielwboyce]

With two small changes found with above-mentioned .c program all of my test polygons are rendering perfectly, as seen in this gallery which compares renders from 4 different commits, with the last one being the commit with perfect results.

Two lines (2065 and 2088) were changed in geom.c, as it stands in commit 4dc135e:

int savedX = nodes[savedIndex].x;

else if (savedX > THRESH) {

became

double savedX = nodes[savedIndex].x;

else if (savedX > q0.x + THRESH) {

I still need to add the test case to make check, so a pull request with that addition and the updates to geom.c will get started in the next couple of days.

[danielwboyce]

I've added a unit test to test-prism.c to check ten points in this geometry

which always failed to render correctly up through the most recent version of node_in_or_on_polygon. The test reads

/************************************************************************/
/* fourth unit test: check of point in polygon test with slanted H     */
/************************************************************************/
int test_point_in_polygon(int write_log) {
  // make array of test points that should always pass
  vector3 pass[5];
  pass[0] = make_vector3(0.3, 0.5, 0.0);
  pass[1] = make_vector3(0.4, 0.4, 0.0);
  pass[2] = make_vector3(0.5, 0.7, 0.0);
  pass[3] = make_vector3(0.5, 0.5, 0.0);
  pass[4] = make_vector3(0.5, 0.3, 0.0);

  // make array of test points that should always pass
  vector3 fail[5];
  fail[0] = make_vector3(0.2, 0.2, 0.0);
  fail[1] = make_vector3(0.3, 0.3, 0.0);
  fail[2] = make_vector3(0.4, 0.6, 0.0);
  fail[3] = make_vector3(0.6, 0.4, 0.0);
  fail[4] = make_vector3(0.7, 0.7, 0.0);

  // make array of nodes for the test polygon (an H slanted by 45 degrees)
  int num_nodes = 12;
  vector3 nodes[num_nodes];
  nodes[0] = make_vector3(0.5, 0.2, 0.0);
  nodes[1] = make_vector3(0.6, 0.3, 0.0);
  nodes[2] = make_vector3(0.5, 0.4, 0.0);
  nodes[3] = make_vector3(0.6, 0.5, 0.0);
  nodes[4] = make_vector3(0.7, 0.4, 0.0);
  nodes[5] = make_vector3(0.8, 0.5, 0.0);
  nodes[6] = make_vector3(0.5, 0.8, 0.0);
  nodes[7] = make_vector3(0.4, 0.7, 0.0);
  nodes[8] = make_vector3(0.5, 0.6, 0.0);
  nodes[9] = make_vector3(0.4, 0.5, 0.0);
  nodes[10] = make_vector3(0.3, 0.6, 0.0);
  nodes[11] = make_vector3(0.2, 0.5, 0.0);

  FILE *f = write_log ? fopen("/tmp/test-prism.point-in-polygon", "w") : 0;

  boolean all_points_success = 1;
  boolean include_boundaries = 1;
  int i;
  for (i = 0; i < 5; i++) {
    boolean local_success = node_in_or_on_polygon(pass[i], nodes, num_nodes, include_boundaries);
    if (!local_success) {
        all_points_success = 0;
    }
    if (f) {
      fprintf(f, "%f %f %i\n", pass[i].x, pass[i].y, local_success);
    }
  }
  for (i = 0; i < 5; i++) {
    boolean local_success = !node_in_or_on_polygon(fail[i], nodes, num_nodes, include_boundaries);
    if (!local_success) {
      all_points_success = 0;
    }
    if (f) {
      fprintf(f, "%f %f %i\n", pass[i].x, pass[i].y, local_success);
    }
  }

  if (f) {
    if (all_points_success) {
      printf("all test points for slanted H pass and fail as expected\n");
    }
    else {
      printf("one or more test points for slanted H do not pass and fail as expected\n");
    }
    fclose(f);
  }

  int num_failed;
  if (all_points_success) {
    num_failed = 0;
    printf("all test points for slanted H pass and fail as expected\n");
  }
  else {
    num_failed = 1;
    printf("one or more test points for slanted H do not pass and fail as expected\n");
  }

  return num_failed;
}

and has also been added to the run_unit_tests function in test-prism.c. I've tested that make check shows passing and failing unit tests when some points that should be in the polygon are changed to be outside it and vice versa and make check behaves as expected. I will start a pull request for this momentarily.