Created
January 2, 2023 03:15
-
-
Save tatarize/c709e7f487ec26285fbcd45087a1688a to your computer and use it in GitHub Desktop.
Fast scanbeam point in polygon code
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| def build_edge_list(polygon): | |
| edge_list = [] | |
| for i in range(0, len(polygon) - 1): | |
| if (polygon[i].imag, polygon[i].real) < (polygon[i + 1].imag, polygon[i + 1].real): | |
| edge_list.append((polygon[i], i)) | |
| edge_list.append((polygon[i + 1], ~i)) | |
| else: | |
| edge_list.append((polygon[i], ~i)) | |
| edge_list.append((polygon[i + 1], i)) | |
| def sort_key(e): | |
| return e[0].imag, e[0].real, ~e[1] | |
| edge_list.sort(key=sort_key) | |
| return edge_list | |
| def build_scanbeam(edge_list): | |
| actives = [] | |
| actives_table = [] | |
| events = [] | |
| y = -float("inf") | |
| for pt, index in edge_list: | |
| if y != pt.imag: | |
| actives_table.append(list(actives)) | |
| events.append(pt.imag) | |
| if index >= 0: | |
| actives.append(index) | |
| else: | |
| actives.remove(~index) | |
| y = pt.imag | |
| actives_table.append(list(actives)) | |
| largest_actives = max([len(a) for a in actives_table]) | |
| scan = np.zeros((len(actives_table), largest_actives), dtype=int) | |
| scan -= 1 | |
| for i, active in enumerate(actives_table): | |
| scan[i, 0: len(active)] = active | |
| return scan, events | |
| def point_in_polygon(polygon, point): | |
| edge_list = build_edge_list(polygon) | |
| scan, events = build_scanbeam(edge_list) | |
| pts_y = np.imag(point) | |
| idx = np.searchsorted(events, pts_y) | |
| actives = scan[idx] | |
| a = polygon[actives] | |
| b = polygon[actives + 1] | |
| a = np.where(actives == -1, np.nan + np.nan * 1j, a) | |
| b = np.where(actives == -1, np.nan + np.nan * 1j, b) | |
| old_np_seterr = np.seterr(invalid="ignore", divide="ignore") | |
| try: | |
| # If horizontal slope is undefined. But, all x-ints are at x since x0=x1 | |
| m = (b.imag - a.imag) / (b.real - a.real) | |
| y0 = a.imag - (m * a.real) | |
| ys = np.reshape(np.repeat(np.imag(point), y0.shape[1]), y0.shape) | |
| x_intercepts = np.where(~np.isinf(m), (ys - y0) / m, a.real) | |
| finally: | |
| np.seterr(**old_np_seterr) | |
| xs = np.reshape(np.repeat(np.real(point), y0.shape[1]), y0.shape) | |
| results = np.sum(x_intercepts <= xs, axis=1) | |
| results %= 2 | |
| return results |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment