Prism geometry bug

[smartalecH]

It looks like there is a bug within the prism geometry implementation. I'm assuming it's a bug within the node_in_or_on_polygon() method of libctl (similar to #400).

I drew a simple arc and saved the geometry as a GDS file that could easily be imported using the gdspy library (I bypassed the libGDSii library to rule out any problems there).

Using the following code:

import meep as mp
from matplotlib import pyplot as plt
import numpy as np
import gdspy

cell_size = mp.Vector3(2,2,0)

Si = mp.Medium(index=3.5)

h = 100.0


# Read in the GDS file
lib = gdspy.GdsLibrary(infile='debug.GDS',units='import',unit=1e-6)
main_cell = lib.top_level()[0]
pol_dict = main_cell.get_polygons(by_spec=True)

# Load the ring geometry
polygons = pol_dict[(0, 0)]
geometry = []
for shape in polygons:
    verts = [mp.Vector3(idx[0],idx[1],-50) for idx in shape]
    geometry.append(mp.Prism(
        vertices=verts,height=h,
        axis=mp.Vector3(0,0,1),material=Si))


sim = mp.Simulation(resolution=200,
                    geometry=geometry,
                    cell_size=cell_size,
                    eps_averaging=False)

sim.init_sim()

eps_data = sim.get_array(size=cell_size,center=mp.Vector3(), component=mp.Dielectric)
plt.figure()
plt.imshow(np.rot90(eps_data), interpolation='spline36', cmap='binary')
plt.axis('off')
plt.tight_layout()
plt.savefig('meep.png')
plt.show()

meep renders the following geometry:

If I had to guess, I would say this is one of those weird boundary cases that apply to most even-odd point-in-polygon algorithms. I've done a little digging to see if anyone has developed an algorithm resilient to these kinds of errors.

Here's a paper that compares several different even-odd methods and finds errors with all of them.

Here's another paper referencing the previous paper that implements a new even-odd algorithm that claims to be "perfect". It's not terribly different from @HomerReid's current implementation.

Here's the GDS.

Here's a tarball with all of the test cases used by the first paper (stored as text files with all of the vertices). Perhaps this could be used to generate a more comprehensive test engine.

It's probably best if @HomerReid takes a look at everything first, but I'm happy to help!

[danielwboyce]

I'm taking a look at this bug. The plan is to start by updating node_in_or_on_polygon() in the libctl library to use the "perfect" algorithm from the aforementioned paper. If rebuilding libctl and meep with that update doesn't fix it, I'll have to take a more in-depth look.

[danielwboyce]

The code from the paper seems to be having its own problems with our implementation; now almost none of the prism shows up correctly:

This is with node_in_or_on_polygon() changed from

boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes,
                              boolean include_boundaries)
{
  vector3 u = {0.0, -1.0, 0.0};
  int nn, edges_crossed=0;
  for(nn=0; nn<num_nodes; nn++)
   { int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn+1)%num_nodes], u, 0);
     if (status==IN_SEGMENT)
      return include_boundaries;
     else if (status==INTERSECTING)
      edges_crossed++;
     else if (status==ON_RAY)
      { vector3 nm1 = nodes[ (nn==0 ? num_nodes-1 : nn-1) ];
        vector3 n0  = nodes[ nn ];
        vector3 np1 = nodes[ (nn+1) % num_nodes ];
        vector3 np2 = nodes[ (nn+2) % num_nodes ];
        int last_status = intersect_ray_with_segment(q0, nm1, n0, u, 0);
        if (last_status==INTERSECTING) edges_crossed--;
        int next_status = intersect_ray_with_segment(q0, np1, np2, u, 0);
        if (next_status==INTERSECTING) edges_crossed--;
      }
   }
  return (edges_crossed % 2);
}

to

boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes,
                              boolean include_boundaries)
{
	// Create axes
	vector3 xAxisPos = {1.0, 0.0. 0.0};
	vector3 xAxisNeg = {-1.0, 0.0, 0.0};
	
	// Initial start point
	vector3 startPoint;
	vector3 endPoint;
	
	int startNodePosition = -1;
	int nn, edges_crossed=0;
	
	// Is q0 on a vertex or edge?
	// Transform coordinate system of nodes such that q0 is at 0|0
	for(nn=0; nn<num_nodes; nn++) {
		int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn+1)%%num_nodes], xAxisPos, 0);
		if (status==IN_SEGMENT) {
			return include_boundaries;
		}
		
		nodes[nn].x -= q0.x;
		nodes[nn].y -= q0.y;
		
		// Find start point which is not on the x axis (from q0)
		if (nodes[nn].y != 0) {
			startPoint.x = nodes[nn].x;
			startPoint.y = nodes[nn].y;
			startNodePosition = nn;
		}
	}
	
	// Move q0 to 0|0
	q0.x = 0;
	q0.y = 0;
	
	// No start point found and point is not on an edge or node
	// --> point is outside
	if (startNodePosition == -1) {
		return 0;
	}
	
	int checkedPoints = 0;
	nn = startNodePosition;
	
	// Consider all edges
	while (checkedPoints < num_points) {
		int savedX = nodes[ (nn+1)%num_nodes ].x;
		int savedIndex = (nn+1)%num_nodes;
		
		// Move to next point which is not on the x-axis
		do {
			nn = (nn+1)%num_nodes;
			checkedPoints++;
		} while (nodes[nn].y == 0);
		// Found end point
		endPoint.x = nodes[nn].x;
		endPoint.y = nodes[nn].y;
		
		// Only intersect lines that cross the x-axis
		if (startPoint.y * endPoint.y < 0) {
			// No nodes have been skipped and the successor node
			// has been chose as the end point
			if (savedIndex == nn) {
				status = intersect_ray_with_segment(q0, startPoint, endPoint, xAxisPos, 0);
				if (status == INTERSECTING) {
					edges_crossed++;
				}
			}
			// If at least one node on the right side has been skipped,
			// the original edge would have been intersected
			// --> intersect with full x-axis
			else if (savedX > 0) {
				int statusPos = intersect_ray_with_segment(q0, startPoint, endPoint, xAxisPos, 0);
				int statusNeg = intersect_ray_with_segment(q0, startPoint, endPoint, xAxisNeg, 0);
				if (statusPos == INTERSECTING || statusNeg == INTERSECTING) {
					edges_crossed++;
				}
			}
		}
		// End point is the next start point
		startPoint = endPoint;
	}
	
	// Odd count --> in the polygon (1)
	// Even count --> outside (0)
	return count % 2;
}

Note this code is heavily based off the demo from the paper with the "perfect" algorithm.

[smartalecH]

What polygon are you testing it on? Did you try with a simple polygon that you know works with the existing method?

[stevengj]

Note that all of the computational geometry calculations are done in https://github.com/NanoComp/libctl/blob/master/utils/geom.c — everything else in Meep calls that, so if you fix the problem in libctl (and recompile/reinstall libctl and recompile Meep) then it will be fixed everywhere.

[danielwboyce]

Some progress has been made in geom.c, but there's still some sort of bug. For reference the function as it stands is

boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes,
                              boolean include_boundaries)
{
	// Create axes
	vector3 xAxis = {1.0, 0.0, 0.0};
	
	// Initial start point
	vector3 startPoint;
	vector3 endPoint;
	
	int startNodePosition = -1;
	int nn, edges_crossed=0;
	
	// Is q0 on a vertex or edge?
	// Transform coordinate system of nodes such that q0 is at 0|0
	for(nn=0; nn<num_nodes; nn++) {
		int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn+1)%num_nodes], xAxis, 0);
		if (status==IN_SEGMENT) {
			return include_boundaries;
		}
		
		// Find start point which is not on the x axis (from q0)
		if (nodes[nn].y - q0.y != 0) {
			startPoint.x = nodes[nn].x;
			startPoint.y = nodes[nn].y;
			startNodePosition = nn;
		}
	}
	
	// No start point found and point is not on an edge or node
	// --> point is outside
	if (startNodePosition == -1) {
		return 0;
	}
	
	int checkedPoints = 0;
	nn = startNodePosition;
	
	// Consider all edges
	while (checkedPoints < num_nodes) {
		int savedX = nodes[ (nn+1)%num_nodes ].x;
		int savedIndex = (nn+1)%num_nodes;
		
		// Move to next point which is not on the x-axis
		do {
			nn = (nn+1)%num_nodes;
			checkedPoints++;
		} while (nodes[nn].y - q0.y == 0);
		// Found end point
		endPoint.x = nodes[nn].x;
		endPoint.y = nodes[nn].y;
		
		// Only intersect lines that cross the x-axis
		if ((startPoint.y - q0.y) * (endPoint.y - q0.y) < 0) {
			// No nodes have been skipped and the successor node
			// has been chose as the end point
			if (savedIndex == nn) {
				int status = intersect_ray_with_segment(q0, startPoint, endPoint, xAxis, 0);
				if (status == INTERSECTING) {
					edges_crossed++;
				}
			}
			// If at least one node on the right side has been skipped,
			// the original edge would have been intersected
			// --> intersect with full x-axis
			else if (savedX > 0) {
				int status = intersect_line_with_segment(q0, startPoint, endPoint, xAxis, 0);
				if (status == INTERSECTING) {
					edges_crossed++;
				}
			}
		}
		// End point is the next start point
		startPoint = endPoint;
	}
	
	// Odd count --> in the polygon (1)
	// Even count --> outside (0)
	return edges_crossed % 2;
}

When I build libctl and then test in meep, Alec's script from above (loading in the .gds file) gives us

However, when I write a script with the same vertices without loading in the.gds file, the results are

It's great that this looks correct, of course, but perhaps there's something wrong with gdspy?

There is still something wrong with the new version of node_in_or_on_polygon. I added some different polygons to my set of test cases. When plotted with the current master branches of Meep and libctl, almost all of them plot correctly, seen here (the notable exceptions being the 6-, 9-, and 12-pointed stars adding some weird lines). When plotted with the version of libctl containing node_in_or_on_polygon as above, we can see that many polygons don't plot correctly (the 3-, 6-, 7-, 8-, and 17-sided polygons and the 3-, 4-, 6-, 9-, 10-, and 12-pointed stars, all adding weird lines).

[danielwboyce]

There must have simply been a build error yesterday or something because today the same code builds to have the curved rectangle look correct with both gdspy and the vertices copy/pasted in. Still working on the weird lines that got added to my other test polygons.

[danielwboyce]

With some more updates to node_in_or_on_polygon I've greatly reduced the occurrence of render errors, but have not yet completely eliminated them.

For reference, this is the code for the new version of node_in_or_on_polygon:

boolean node_in_or_on_polygon(vector3 q0, vector3 *nodes, int num_nodes,
                              boolean include_boundaries)
{
	// Create axes
	vector3 xAxis = {1.0, 0.0, 0.0};
	
	// Initial start point
	vector3 startPoint;
	vector3 endPoint;
	
	int startNodePosition = -1;
	int nn, edges_crossed = 0;
	
	// Is q0 on a vertex or edge?
	// Transform coordinate system of nodes such that q0 is at 0|0
	for(nn = 0; nn < num_nodes; nn++) {
		int status = intersect_ray_with_segment(q0, nodes[nn], nodes[(nn+1)%num_nodes], xAxis, 0);
		if (status == IN_SEGMENT) {
			return include_boundaries;
		}
		
		// Find start point which is not on the x axis (from q0)
		if (fabs(nodes[nn].y - q0.y) > THRESH) {
			startNodePosition = nn;
			startPoint = nodes[startNodePosition];
		}
	}
	
	// No start point found and point is not on an edge or node
	// --> point is outside
	if (startNodePosition == -1) {
		return 0;
	}
	
	int checkedPoints = 0;
	nn = startNodePosition;
	
	// Consider all edges
	while (checkedPoints < num_nodes) {
		int savedIndex = (nn+1)%num_nodes;
		int savedX = nodes[savedIndex].x;
		
		// Move to next point which is not on the x-axis
		do {
			nn = (nn+1)%num_nodes;
			checkedPoints++;
		} while (fabs(nodes[nn].y - q0.y) < THRESH);
		// Found end point
		endPoint = nodes[nn];
		
		// Only intersect lines that cross the x-axis (don't need to correct for rounding 
		// error in the if statement because startPoint and endPoint are screened to
		// never lie on the x-axis)
		if ((startPoint.y - q0.y) * (endPoint.y - q0.y) < 0) {
			// No nodes have been skipped and the successor node
			// has been chose as the end point
			if (savedIndex == nn) {
				int status = intersect_ray_with_segment(q0, startPoint, endPoint, xAxis, 0);
				if (status == INTERSECTING) {
					edges_crossed++;
				}
			}
			// If at least one node on the right side has been skipped,
			// the original edge would have been intersected
			// --> intersect with full x-axis
			else if (savedX > THRESH) {
				int status = intersect_line_with_segment(q0, startPoint, endPoint, xAxis, 0);
				if (status == INTERSECTING) {
					edges_crossed++;
				}
			}
		}
		// End point is the next start point
		startPoint = endPoint;
	}
	
	// Odd count --> in the polygon (1)
	// Even count --> outside (0)
	return edges_crossed % 2;
}

You can see in this imgur gallery the results from comparison run between the two old and new versions of node_in_or_on_polygon for over 250 polygons. My polygons mostly come from this tarball of test polygons from the paper mentioned in the first comment of this thread, but I've also added in some regular polygons and stars. This is a spreadsheet tabulating which polygons have errors and which do not. To get a sense of the sorts of errors that are common to each of the builds it's worth looking through the imgur gallery, but overall I noticed one main type of error in each build.

In the old version, I found that the node_in_or_on_polygon function causes a lot of artifacts in a render, which we noticed at the beginning of this thread, but is also demonstrated here, for example (with the master branch render on the left):

This artifact introduction is typified in polygons with orthogonal vertices. In the batch of 80 such polygons I tested, 78 introduced artifacts to the render, but it shows up elsewhere, too. Over the 263 polygons I tested, about 90 showed this error.

By comparison, the new version of node_in_or_on_polygon occasionally omits small lines, demonstrated here:

Over the 263 polygons I tested, 49 showed this error. I think it's safe to say that, on average, the new version introduces smaller errors less often than the old version and that libctl probably wouldn't be hurt with an update containing the new version of node_in_or_on_polygon.

To solve the error with lines, I would say we should start by looking at the possibility of floating-point error. I've noticed that all or nearly all of the line errors occur to the left of another vertex. The algorithm I used is supposed to account for this possibility (see this paper), but perhaps there's a floating-point leak somewhere.

[stevengj]

Can you file a new issue (in libctl, not in Meep) with a simple case exhibiting the spurious lines?