diff --git a/.gitignore b/.gitignore
index c62d86f..a7bd822 100644
--- a/.gitignore
+++ b/.gitignore
@@ -143,4 +143,7 @@ dmypy.json
 cython_debug/
 
 # Mac Desktop Services Store file
-.DS_STORE
\ No newline at end of file
+.DS_STORE
+
+# vscode
+.vscode/
\ No newline at end of file
diff --git a/data/cube_hole.stl b/data/cube_hole.stl
new file mode 100644
index 0000000..597c630
Binary files /dev/null and b/data/cube_hole.stl differ
diff --git a/stltovoxel/perimeter.py b/stltovoxel/perimeter.py
index 6fd97c3..497034d 100644
--- a/stltovoxel/perimeter.py
+++ b/stltovoxel/perimeter.py
@@ -1,10 +1,11 @@
-from . import winding_query
+from . import polygon_repair
 
 
-def repaired_lines_to_voxels(line_list, pixels):
+def repaired_lines_to_voxels(line_list, pixels, plane_shape):
     if not line_list:
         return
-    wq = winding_query.WindingQuery([[tuple(pt.tolist())[:2] for pt in seg] for seg in line_list])
+    segments = [[tuple(pt.tolist())[:2] for pt in seg] for seg in line_list]
+    wq = polygon_repair.PolygonRepair(segments, plane_shape)
     wq.repair_all()
     new_line_list = []
     for polyline in wq.loops:
diff --git a/stltovoxel/polygon_repair.py b/stltovoxel/polygon_repair.py
new file mode 100644
index 0000000..05598bc
--- /dev/null
+++ b/stltovoxel/polygon_repair.py
@@ -0,0 +1,243 @@
+import numpy as np
+import math
+
+
+def find_polylines(segments):  # noqa: C901
+    polylines = []
+    segment_forward_dict = {}
+    segment_backward_dict = {}
+
+    # create a dictionary where each endpoint is a key, and the value is a list of segments
+    for segment in segments:
+        start, end = segment
+        if start not in segment_forward_dict:
+            segment_forward_dict[start] = []
+        segment_forward_dict[start].append(end)
+        if end not in segment_backward_dict:
+            segment_backward_dict[end] = []
+        segment_backward_dict[end].append(start)
+    # loop through each segment, and recursively add connected segments to form polylines
+    while len(segment_forward_dict) > 0:
+        start = list(segment_forward_dict.keys())[0]
+        polyline = [start]
+
+        # Iterate through forward dict
+        while True:
+            if start not in segment_forward_dict:
+                break
+            next_points = segment_forward_dict[start]
+            end = next_points[0]
+
+            segment_forward_dict[start].remove(end)
+            if len(segment_forward_dict[start]) == 0:
+                del segment_forward_dict[start]
+            segment_backward_dict[end].remove(start)
+            if len(segment_backward_dict[end]) == 0:
+                del segment_backward_dict[end]
+            # Append points to the end of the list
+            polyline.append(end)
+            start = end
+
+        # Reset the start point
+        start = polyline[0]
+        # Iterate through backward dict
+        while True:
+            if start not in segment_backward_dict:
+                break
+            next_points = segment_backward_dict[start]
+            end = next_points[0]
+
+            segment_backward_dict[start].remove(end)
+            if len(segment_backward_dict[start]) == 0:
+                del segment_backward_dict[start]
+            segment_forward_dict[end].remove(start)
+            if len(segment_forward_dict[end]) == 0:
+                del segment_forward_dict[end]
+
+            # Insert points at the front of the list
+            polyline.insert(0, end)
+            start = end
+        polylines.append(polyline)
+
+    return polylines
+
+
+def find_polyline_endpoints(segs):
+    start_to_end = dict()
+    end_to_start = dict()
+
+    for start, end in segs:
+        # Update connections for the new segment
+        actual_start = end_to_start.get(start, start)
+        actual_end = start_to_end.get(end, end)
+
+        # Check for loops or redundant segments
+        if actual_start == actual_end:
+            del end_to_start[actual_start]
+            del start_to_end[actual_start]
+            continue  # This segment forms a loop or is redundant, so skip it
+
+        # Merge polylines or add new segment
+        start_to_end[actual_start] = actual_end
+        end_to_start[actual_end] = actual_start
+
+        # Remove old references if they are now internal points of a polyline
+        if start in end_to_start:
+            del end_to_start[start]
+        if end in start_to_end:
+            del start_to_end[end]
+
+    return start_to_end
+
+
+def atan_sum(f1, f2):
+    # Atan sum identity
+    # atan2(atan_sum(f1, f2)) == atan2(f1) + atan2(f2)
+    x1, y1 = f1
+    x2, y2 = f2
+    return (x1*x2 - y1*y2, y1*x2 + x1*y2)
+
+
+def atan_diff(f1, f2):
+    # Atan difference identity
+    # atan2(atan_diff(f1, f2)) == atan2(f1) - atan2(f2)
+    x1, y1 = f1
+    x2, y2 = f2
+    return (x1*x2 + y1*y2, y1*x2 - x1*y2)
+
+
+def subtract(s1, s2):
+    return (s1[0] - s2[0], s1[1] - s2[1])
+
+
+def add(s1, s2):
+    return (s1[0] + s2[0], s1[1] + s2[1])
+
+
+def winding_contour_pole(pos, pt, repel):
+    # The total winding number of a system is composed of a sum of atan2 functions (one at each point of each line segment)
+    # The gradient of atan2(y, x) is
+    # Gx = -y/(x^2+y^2); Gy = x/(x^2+y^2)
+    # The contour (direction of no increase) is orthogonal to the gradient, either
+    # (-Gy, Gx) or (Gy, -Gx)
+    # This is represented by:
+    # Cx =  x/(x^2+y^2); Cy =  y/(x^2+y^2) or
+    # Cx = -x/(x^2+y^2); Cy = -y/(x^2+y^2)
+    # In practice, one of each is used per line segment which repels & attracts the vector field.
+    x, y = subtract(pos, pt)
+    dist2 = (x**2 + y**2)
+    cx = x/dist2
+    cy = y/dist2
+    if repel:
+        return (cx, cy)
+    else:
+        return (-cx, -cy)
+
+
+def normalize(pt):
+    x, y = pt
+    dist = math.sqrt(x**2 + y**2)
+    return (x / dist, y / dist)
+
+
+def distance(p1, p2):
+    x1, y1 = p1
+    x2, y2 = p2
+    return math.sqrt((y2 - y1)**2 + (x2 - x1)**2)
+
+
+def initial_direction(pt, segs):
+    # Computing the winding number of a given point requires 2 angle computations per line segment (one per point).
+    # This method makes 2 angle computations for every segment in other_segs to compute the winding number at my_seg[1].
+    # It also makes one angle computation from my_seg[0] to my_seg[1].
+    # The result is an angle out of my_seg which leads to a winding number of a target value.
+    # In theory this target winding number can be any value, but here pi radians (180 degrees) is used for simplicity.
+    # Also, generally angle computation would take the form of [atan2(start-pt)-atan2(end-pt)]+... , but multiple atan2 calls
+    # can be avoided through use of atan2 identities.
+    # Since the final angle would be converted back into a vector, no atan2 call is required.
+    accum = (1, 0)
+    for seg_start, seg_end in segs:
+        # Low quality meshes may have multiple segments with the same start or end point, so we should not attempt
+        # to compute the angle from pt because the result is undefined.
+        if seg_start != pt:
+            accum = atan_sum(accum, subtract(seg_start, pt))
+        if seg_end != pt:
+            # This will trigger for at least one segment because pt comes from this same list of segments.
+            accum = atan_diff(accum, subtract(seg_end, pt))
+        # Without this accum can get arbitrarily large and lead to floating point problems
+        accum = normalize(accum)
+    return np.array(accum)
+
+
+def winding_contour(pos, segs):
+    accum = (0, 0)
+    for start, end in segs:
+        # Starting points attract the vector field
+        start_vec = winding_contour_pole(pos, start, repel=False)
+        accum = add(accum, start_vec)
+        # Ending points repel the vector field
+        end_vec = winding_contour_pole(pos, end, repel=True)
+        accum = add(accum, end_vec)
+    return normalize(accum)
+
+
+def winding_number_search(start, ends, segs, polyline_endpoints, max_iterations):
+    # Find the initial direction to start marching towards.
+    direction = initial_direction(start, segs)
+    # Move slightly toward that direction to pick the contour that we will use below
+    pos = start + (direction * 0.1)
+    seg_outs = [(tuple(start), tuple(pos))]
+    for _ in range(max_iterations):
+        # Flow lines in this vector field have equivalent winding numbers
+        # As an optimization, polyline endpoints are used instead of original_segments for winding_contour because
+        # start and end points in the same place cancel out.
+        direction = winding_contour(pos, polyline_endpoints)
+        # March along the flowline to find where it terminates.
+        new_pos = pos + direction
+        seg_outs.append((tuple(pos), tuple(new_pos)))
+        pos = new_pos
+        # It should terminate at an endpoint
+        for end in ends:
+            if distance(pos, end) < 1:
+                seg_outs.append((tuple(pos), tuple(end)))
+                return end
+
+    raise Exception("Failed to repair mesh")
+
+
+class PolygonRepair():
+    def __init__(self, segments, dimensions):
+        # Maps endpoints to the polygon they form
+        self.loops = []
+        # Populate initially
+        self.polylines = []
+        self.original_segments = segments
+        self.dimensions = dimensions
+        self.collapse_segments()
+
+    def collapse_segments(self):
+        self.loops = []
+        self.polylines = []
+        for polyline in find_polylines(self.original_segments):
+            if polyline[0] == polyline[-1]:
+                self.loops.append(polyline)
+            else:
+                self.polylines.append(polyline)
+
+    def repair_all(self):
+        while self.polylines:
+            self.repair_polyline()
+            old_seg_length = len(self.polylines)
+            self.collapse_segments()
+            assert old_seg_length - 1 == len(self.polylines)
+        assert len(self.polylines) == 0
+
+    def repair_polyline(self):
+        # Search starts at the end of an arbitrary polyline
+        search_start = self.polylines[0][-1]
+        # Search will conclude when it finds the beginning of any polyline (including itself)
+        search_ends = [polyline[0] for polyline in self.polylines]
+        polyline_endpoints = [(polyline[0], polyline[-1]) for polyline in self.polylines]
+        max_iterations = self.dimensions[0] + self.dimensions[1]
+        end = winding_number_search(search_start, search_ends, self.original_segments, polyline_endpoints, max_iterations)
+        self.original_segments.append((search_start, end))
diff --git a/stltovoxel/slice.py b/stltovoxel/slice.py
index 127d2b0..9f481a6 100644
--- a/stltovoxel/slice.py
+++ b/stltovoxel/slice.py
@@ -66,7 +66,7 @@ def paint_z_plane(mesh, height, plane_shape):
                 pt2 = points[(i + 1) % 3]
                 lines.append((pt, pt2))
 
-    perimeter.repaired_lines_to_voxels(lines, pixels)
+    perimeter.repaired_lines_to_voxels(lines, pixels, plane_shape)
 
     return height, pixels
 
diff --git a/stltovoxel/visualization.py b/stltovoxel/visualization.py
new file mode 100644
index 0000000..8090417
--- /dev/null
+++ b/stltovoxel/visualization.py
@@ -0,0 +1,57 @@
+import matplotlib.pyplot as plt
+
+
+def plot_line_segments_pixels(line_segments, pixels):
+    # Loop over each line segment and plot it using the plot function
+    for p1, p2 in line_segments:
+        x1, y1 = p1[:2]
+        x2, y2 = p2[:2]
+        dx = x2 - x1
+        dy = y2 - y1
+        plt.plot([x1, x2], [y1, y2])
+        plt.arrow(x1, y1, dx/2, dy/2, head_width=0.1, head_length=0.1, fc='k', ec='k')
+    for y in range(pixels.shape[0]):
+        for x in range(pixels.shape[1]):
+            xoff = 0
+            yoff = 0
+            plt.gca().add_patch(plt.Rectangle((x + xoff, y + yoff), 1, 1, fill=False))
+            if pixels[y][x]:
+                plt.plot(x + .5 + xoff, y + .5 + yoff, 'ro')
+            else:
+                plt.plot(x + .5 + xoff, y + .5 + yoff, 'bo')
+
+    # Show the plot
+    plt.show()
+
+
+def plot_polylines(polylines):
+    """
+    Plots a polyline using pyplot.
+
+    Parameters:
+    polyline (list): A list of (x,y) tuples representing the vertices of the polyline.
+
+    Returns:
+    None
+    """
+    # Create a new figure
+    _, ax = plt.subplots()
+    for polyline in polylines:
+        # Extract the x and y coordinates from the polyline
+        x_coords, y_coords = zip(*polyline)
+
+        # Plot the polyline using the plot function
+        ax.plot(x_coords, y_coords)
+
+        middle_ind = int(len(polyline)/2)
+
+        if len(polyline) <= middle_ind+1:
+            middle_ind = 0
+        x1, y1 = polyline[middle_ind]
+        x2, y2 = polyline[middle_ind+1]
+        dx = x2 - x1
+        dy = y2 - y1
+        ax.arrow(x1, y1, dx/2, dy/2, head_width=0.1, head_length=0.1, fc='k', ec='k')
+
+    # Show the plot
+    plt.show()
diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py
deleted file mode 100644
index 18e7fdf..0000000
--- a/stltovoxel/winding_query.py
+++ /dev/null
@@ -1,222 +0,0 @@
-import numpy as np
-import math
-import functools
-from queue import PriorityQueue
-
-
-def find_polylines(segments):  # noqa: C901
-    polylines = []
-    segment_forward_dict = {}
-    segment_backward_dict = {}
-
-    # create a dictionary where each endpoint is a key, and the value is a list of segments
-    for segment in segments:
-        start, end = segment
-        if start not in segment_forward_dict:
-            segment_forward_dict[start] = []
-        segment_forward_dict[start].append(end)
-        if end not in segment_backward_dict:
-            segment_backward_dict[end] = []
-        segment_backward_dict[end].append(start)
-    # loop through each segment, and recursively add connected segments to form polylines
-    while len(segment_forward_dict) > 0:
-        start = list(segment_forward_dict.keys())[0]
-        polyline = [start]
-
-        # Iterate through forward dict
-        while True:
-            if start not in segment_forward_dict:
-                break
-            next_points = segment_forward_dict[start]
-            end = next_points[0]
-
-            segment_forward_dict[start].remove(end)
-            if len(segment_forward_dict[start]) == 0:
-                del segment_forward_dict[start]
-            segment_backward_dict[end].remove(start)
-            if len(segment_backward_dict[end]) == 0:
-                del segment_backward_dict[end]
-            # Append points to the end of the list
-            polyline.append(end)
-            start = end
-
-        # Reset the start point
-        start = polyline[0]
-        # Iterate through backward dict
-        while True:
-            if start not in segment_backward_dict:
-                break
-            next_points = segment_backward_dict[start]
-            end = next_points[0]
-
-            segment_backward_dict[start].remove(end)
-            if len(segment_backward_dict[start]) == 0:
-                del segment_backward_dict[start]
-            segment_forward_dict[end].remove(start)
-            if len(segment_forward_dict[end]) == 0:
-                del segment_forward_dict[end]
-
-            # Insert points at the front of the list
-            polyline.insert(0, end)
-            start = end
-        polylines.append(polyline)
-
-    return polylines
-
-
-def closest_distance(point, goals):
-    return min([dist(point, goal) for goal in goals])
-
-
-def dist(p1, p2):
-    x1, y1 = p1
-    x2, y2 = p2
-    return math.sqrt((y2 - y1)**2 + (x2 - x1)**2)
-
-
-def normalize(num):
-    return ((num + math.pi) % (2*math.pi)) - math.pi
-
-
-def signed_point_line_dist(line, point):
-    a, b = line
-    x1, y1 = a
-    x2, y2 = b
-    x0, y0 = point
-    num = ((x2 - x1)*(y1 - y0) - (x1 - x0)*(y2 - y1))
-    denom = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
-    return num / denom
-
-
-def close_to_goal(start, goals):
-    sx, sy = start
-    for goal in goals:
-        gx, gy = goal
-        if abs(sx-gx) <= 1 and abs(sy-gy) <= 1:
-            return goal
-    return False
-
-
-class WindingQuery():
-    def __init__(self, segments):
-        # Maps endpoints to the polygon they form
-        self.loops = []
-        # Populate initially
-        self.polylines = []
-        self.original_segments = segments
-        self.collapse_segments()
-
-    def collapse_segments(self):
-        self.loops = []
-        self.polylines = []
-        for polyline in find_polylines(self.original_segments):
-            if polyline[0] == polyline[-1]:
-                self.loops.append(polyline)
-            else:
-                self.polylines.append(polyline)
-
-    def repair_all(self):
-        while self.polylines:
-            self.repair_segment()
-            old_seg_length = len(self.polylines)
-            self.collapse_segments()
-            assert old_seg_length - 1 == len(self.polylines)
-        assert len(self.polylines) == 0
-
-    def repair_segment(self):
-        # Search starts at the end of a polyline
-        start = self.polylines[0][-1]
-
-        # Search will conclude when it finds the beginning of any polyline (including itself)
-        endpoints = [polyline[0] for polyline in self.polylines]
-        end = self.a_star(start, endpoints)
-
-        new_segment = (start, end)
-        self.original_segments.append(new_segment)
-
-    def a_star(self, start, goals):
-        frontier = PriorityQueue()
-        frontier.put((0, start))
-        cost_so_far = {start: 0}
-
-        current = None
-        while not frontier.empty():
-            score, current = frontier.get()
-            close_goal = close_to_goal(current, goals)
-            if close_goal:
-                current = close_goal
-                break
-
-            for dx in range(-1, 2):
-                for dy in range(-1, 2):
-                    next_point = (current[0] + dx, current[1] + dy)
-                    heuristic_cost = closest_distance(next_point, goals)
-                    new_cost = cost_so_far[current] + abs(self.query_winding(next_point) - math.pi) * 200
-
-                    if next_point not in cost_so_far or new_cost < cost_so_far[next_point]:
-                        cost_so_far[next_point] = new_cost
-                        priority = new_cost + heuristic_cost
-                        frontier.put((priority, next_point))
-        assert current is not None
-        return current
-
-    def query_winding(self, point):
-        total = 0
-        for polyline in self.polylines:
-            total += self.winding_segment(polyline, point)
-        return total
-
-    def winding_segment(self, polyline, point):
-        collapsed = (polyline[0], polyline[-1])
-        inner_line, outer_line = self.get_lines(tuple(polyline))
-
-        if len(polyline) == 2:
-            # This is the actual segment so okay to be behind it.
-            return self.winding(collapsed, point)
-        elif signed_point_line_dist(inner_line, point) < 0:
-            # We are inside and beyond any concavity
-            return self.winding(collapsed, point)
-        elif signed_point_line_dist(outer_line, point) > 0:
-            # We are outside beyond any convexity
-            return self.winding(collapsed, point)
-        else:
-            split = int(len(polyline) / 2)
-            # Otherwise, split the segment
-            return self.winding_segment(polyline[:split+1], point) + self.winding_segment(polyline[split:], point)
-
-    def winding(self, segment, point):
-        start, end = segment
-        start = np.array(start)
-        end = np.array(end)
-        point = np.array(point)
-        # Side 1
-        offset = point - start
-        ang1 = math.atan2(offset[1], offset[0])
-        # Side 2
-        offset = point - end
-        ang2 = math.atan2(offset[1], offset[0])
-        return normalize(ang2 - ang1)
-
-    @functools.lru_cache(maxsize=None)
-    def get_lines(self, polyline):
-        start = np.array(polyline[0])
-        end = np.array(polyline[-1])
-        slope = end - start
-
-        furthest_in = 0
-        innermost = polyline[0]
-        furthest_out = 0
-        outermost = polyline[0]
-        for pt in polyline:
-            dist = signed_point_line_dist((polyline[0], polyline[-1]), pt)
-            if dist < furthest_in:
-                innermost = pt
-                furthest_in = dist
-            elif dist > furthest_out:
-                outermost = pt
-                furthest_out = dist
-        innermost = np.array(innermost)
-        outermost = np.array(outermost)
-        inner_line = (innermost, innermost + slope)
-        outer_line = (outermost, outermost + slope)
-        return inner_line, outer_line
diff --git a/test/test_polygon_repair.py b/test/test_polygon_repair.py
new file mode 100644
index 0000000..5c149bc
--- /dev/null
+++ b/test/test_polygon_repair.py
@@ -0,0 +1,110 @@
+import numpy as np
+import unittest
+
+from stltovoxel import polygon_repair
+
+
+class TestPolygonRepair(unittest.TestCase):
+    def tuples_almost_equal(self, actual, expected):
+        for ac_val, exp_val in zip(actual, expected):
+            self.assertAlmostEqual(ac_val, exp_val, 3, f"{actual} != {expected}")
+
+    def test_find_polylines_cycle(self):
+        line_segments = [((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (0, 1)), ((0, 1), (0, 0))]
+        expected_polylines = [[(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]]
+        actual_polylines = polygon_repair.find_polylines(line_segments)
+        self.assertEqual(actual_polylines, expected_polylines)
+
+    def test_find_polylines_no_cycle(self):
+        line_segments = [
+            ((0, 0), (1, 0)),
+            ((1, 0), (2, 1)),
+            ((3, 1), (4, 0)),
+            ((4, 0), (5, 0)),
+            ((5, 0), (6, 1)),
+            ((7, 1), (8, 0)),
+            ((8, 0), (9, 0))
+        ]
+        expected_polylines = [[(0, 0), (1, 0), (2, 1)], [(3, 1), (4, 0), (5, 0), (6, 1)], [(7, 1), (8, 0), (9, 0)]]
+        actual_polylines = polygon_repair.find_polylines(line_segments)
+        self.assertEqual(actual_polylines, expected_polylines)
+
+    def test_find_polylines_out_of_order(self):
+        line_segments = [((0, 0), (1, 0)), ((-1, 0), (0, 0)), ((1, 0), (1, 1))]
+        expected_polylines = [[(-1, 0), (0, 0), (1, 0), (1, 1)]]
+        actual_polylines = polygon_repair.find_polylines(line_segments)
+        self.assertEqual(actual_polylines, expected_polylines)
+
+    def test_initial_direction(self):
+        segs = [((0, 1), (0, 0)), ((1, 0), (1, 1)), ((1, 1), (0, 1)), ]
+        pt = (0, 0)
+        actual = polygon_repair.initial_direction(pt, segs)
+        actual = tuple(actual)
+        expected = (1.0, 0)
+        self.tuples_almost_equal(actual, expected)
+
+    def test_initial_direction2(self):
+        segs = [((0, 1), (1, 1)), ((1, 0), (0, 0))]
+        pt = (1, 1)
+        actual = polygon_repair.initial_direction(pt, segs)
+        actual = tuple(actual)
+        expected = polygon_repair.normalize((-1, -1))
+        self.tuples_almost_equal(actual, expected)
+
+    def test_initial_direction3(self):
+        segs = [((1, 0), (0, 0)), ((1, 1), (1, 0))]
+        pt = (0, 0)
+        actual = polygon_repair.initial_direction(pt, segs)
+        actual = tuple(actual)
+        expected = polygon_repair.normalize((1, 1))
+        self.tuples_almost_equal(actual, expected)
+
+    def test_winding_contour(self):
+        segs = [((0, 0), (1, 0)), ((1, 0), (1, 1))]
+        pos = np.array((1, 1)) - np.array((0.1, 0.1))
+        # Treats starts as repellers and ends as attractors
+        actual = polygon_repair.winding_contour(pos, segs)
+        expected = polygon_repair.normalize((-1, -1))
+        self.tuples_almost_equal(actual, expected)
+
+    def test_winding_contour2(self):
+        segs = [((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (0, 1))]
+        pos = (0, 0.5)
+        actual = polygon_repair.winding_contour(pos, segs)
+        expected = polygon_repair.normalize((0, -1))
+        self.tuples_almost_equal(actual, expected)
+
+        pos = (-.5, 0.5)
+        actual = polygon_repair.winding_contour(pos, segs)
+        expected = polygon_repair.normalize((0, -1))
+        self.tuples_almost_equal(actual, expected)
+
+        pos = (.5, 0.5)
+        actual = polygon_repair.winding_contour(pos, segs)
+        expected = polygon_repair.normalize((0, -1))
+        self.tuples_almost_equal(actual, expected)
+
+    def test_winding_contour3(self):
+        segs = [((0, 0), (1, 0)), ((1, 1), (0, 1))]
+        pos = (0.00001, 0.00001)
+        actual = polygon_repair.winding_contour(pos, segs)
+        expected = polygon_repair.normalize((-1, -1))
+        self.tuples_almost_equal(actual, expected)
+
+    def test_repair_all_close_case(self):
+        segs = [((0, 0), (1, 0)), ((1, .5), (0, .5))]
+        repair = polygon_repair.PolygonRepair(segs, (10, 10))
+        repair.repair_all()
+        expected = [[(0, 0), (1, 0), (1, 0.5), (0, 0.5), (0, 0)]]
+        self.assertEqual(repair.loops, expected)
+
+    def test_repair_all_far_case(self):
+        segs = [((0, 0), (1, 0)), ((1, 1.5), (0, 1.5))]
+        repair = polygon_repair.PolygonRepair(segs, (10, 10))
+        repair.repair_all()
+        expected = [[(0, 0), (1, 0), (0, 0)], [(1, 1.5), (0, 1.5), (1, 1.5)]]
+        self.assertEqual(repair.loops, expected)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/test/test_winding_query.py b/test/test_winding_query.py
deleted file mode 100644
index 51f2598..0000000
--- a/test/test_winding_query.py
+++ /dev/null
@@ -1,84 +0,0 @@
-import unittest
-import math
-
-from stltovoxel import winding_query
-
-
-class TestWindingQuery(unittest.TestCase):
-    def test_winding_query(self):
-        segments = [
-            [(10, 10), (20, 10)],
-            [(20, 10), (30, 10)],
-            [(30, 10), (40, 10)],
-            [(40, 10), (30, 20)],
-            [(30, 20), (40, 40)],
-            [(40, 40), (10, 40)],
-            [(10, 40), (10, 9.9)],
-        ]
-        wq = winding_query.WindingQuery(segments)
-        self.assertAlmostEqual(wq.query_winding((35, 35)), math.pi*2, places=2)
-        self.assertAlmostEqual(wq.query_winding((35, 24)), 0, places=2)
-        self.assertAlmostEqual(wq.query_winding((35, 12)), math.pi*2, places=2)
-
-    def test_find_polylines_cycle(self):
-        line_segments = [((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (0, 1)), ((0, 1), (0, 0))]
-        expected_polylines = [[(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]]
-        actual_polylines = winding_query.find_polylines(line_segments)
-        self.assertEqual(actual_polylines, expected_polylines)
-
-    def test_find_polylines_no_cycle(self):
-        line_segments = [
-            ((0, 0), (1, 0)),
-            ((1, 0), (2, 1)),
-            ((3, 1), (4, 0)),
-            ((4, 0), (5, 0)),
-            ((5, 0), (6, 1)),
-            ((7, 1), (8, 0)),
-            ((8, 0), (9, 0))
-        ]
-        expected_polylines = [[(0, 0), (1, 0), (2, 1)], [(3, 1), (4, 0), (5, 0), (6, 1)], [(7, 1), (8, 0), (9, 0)]]
-        actual_polylines = winding_query.find_polylines(line_segments)
-        self.assertEqual(actual_polylines, expected_polylines)
-
-    def test_find_polylines_out_of_order(self):
-        line_segments = [((0, 0), (1, 0)), ((-1, 0), (0, 0)), ((1, 0), (1, 1))]
-        expected_polylines = [[(-1, 0), (0, 0), (1, 0), (1, 1)]]
-        actual_polylines = winding_query.find_polylines(line_segments)
-        self.assertEqual(actual_polylines, expected_polylines)
-
-    def test_regression(self):
-        segments3 = [
-            ((55, 0), (84, 16)),
-            ((84, 16), (95, 48)),
-            ((95, 48), (84, 79)),
-            ((84, 79), (84, 80)),
-            ((84, 80), (83, 80)),
-            ((83, 80), (55, 97)),
-            ((0, 65), (0, 32)),
-            ((0, 32), (0, 31)),
-            ((0, 31), (21, 6)),
-            ((21, 6), (21, 5)),
-            ((22, 5), (55, 0)),
-        ]
-        wq = winding_query.WindingQuery(segments3)
-        wq.repair_all()
-        self.assertEqual(wq.loops, [[
-                (55, 0),
-                (84, 16),
-                (95, 48),
-                (84, 79),
-                (84, 80),
-                (83, 80),
-                (55, 97),
-                (0, 65),
-                (0, 32),
-                (0, 31),
-                (21, 6),
-                (21, 5),
-                (22, 5),
-                (55, 0)
-            ]])
-
-
-if __name__ == '__main__':
-    unittest.main()