From ebf50094c181d585bb2e5c0464402fc317f1be00 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sat, 8 Apr 2023 10:58:41 -0700 Subject: [PATCH 01/25] mostly working --- stltovoxel/playground.py | 187 ++++++++++++++++++++++++++++++++++++ stltovoxel/visualization.py | 57 +++++++++++ 2 files changed, 244 insertions(+) create mode 100644 stltovoxel/playground.py create mode 100644 stltovoxel/visualization.py diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py new file mode 100644 index 0000000..a5f9357 --- /dev/null +++ b/stltovoxel/playground.py @@ -0,0 +1,187 @@ +import matplotlib.pyplot as plt +import numpy as np +import math +import pdb + +def normalize(num): + return ((num + math.pi) % (2*math.pi)) - math.pi + +def atansum(f1, f2): + y, x = f1 + z, w = f2 + return (y*w + x*z, x*w - y*z) + +def atandiff(f1, f2): + y1, x1 = f1 + y2, x2 = f2 + return (y1*x2 - x1*y2, x1*x2 + y1*y2) + +def negatan(f1): + y,x = f1 + return -y, x + +# calculate baseline from angle of last two line segments of all incomplete loops +## Find a more pleasing way to represent baseline. Angles are good, but the asymmetry is uncomfortable. +# Look at neighbors distance and angle, incorporate baseline, decide who to connect to. + +# +# +def is_right(line, point): + a, b = line + x1, y1 = a + x2, y2 = b + x0, y0 = point + num = ((x2 - x1)*(y1 - y0) - (x1 - x0)*(y2 - y1)) + return num > 0 + +def corners(segs): + im = np.zeros((100,100)) + for i in range(len(segs) - 1): + # pdb.set_trace() + + s1 = segs[i] + s2 = segs[i+1] + assert s1[1] == s2[0] + a = np.array(list(s1[0])) + b = np.array(list(s1[1])) + c = np.array(list(s2[1])) + l1 = a - b + l2 = b - c + # pdb.set_trace() + + ang = math.atan2(*atandiff((l2[1], l2[0]), (l1[1], l1[0]))) + + for x in range(-50, 50): + for y in range(-50, 50): + + if ang < 0 : + # Concave case + if not is_right((a, b), (x,y)) and not is_right((c, b), (x,y)): + im[y+50][x+50] += ang + 2 * math.pi + else: + im[y+50][x+50] += ang + else : + # Convex case (normal) + if is_right((a, b), (x,y)) and is_right((c, b), (x,y)): + im[y+50][x+50] += ang - 2 * math.pi + else: + im[y+50][x+50] += ang + return im + + +def edge_start(me_seg, them_pt): + # Starter unpaired point + s1, s2 = me_seg + start = s1 + angle = (s1[1] - s2[1], s1[0] - s2[0]) + p1 = (start[1] - them_pt[1], start[0] - them_pt[0]) + return math.atan2(*atansum(negatan(p1), angle)) + +def edge_end(me_seg, them_pt): + # Starter unpaired point + s1, s2 = me_seg + start = s2 + angle = (s1[1] - s2[1], s1[0] - s2[0]) + p1 = (start[1] - them_pt[1], start[0] - them_pt[0]) + return math.atan2(*atansum(p1, negatan(angle))) + + +def line(segs, inc): + im = np.zeros((100,100)) + for x in range(-50, 50): + for y in range(-50, 50): + answers = [] + for p1, p2 in segs: + x1, y1 = p1 + x2, y2 = p2 + angle = (y2 - y1, x2 - x1) + # angle = (0,0) + p1 = (y-y1, x-x1) + + p2 = (y-y2, x-x2) + answers.append(atansum(negatan(p1), angle)) + answers.append(atansum(p2, negatan(angle))) + # answers.append(atansum( + # atansum(p1, negatan(angle)), + # atansum(negatan(p2), angle), + # ) + # ) + # answers.append(-normalize(math.atan2(y-y2, x-x2) - angle)) + accum = 0 + for i in inc: + accum += math.atan2(*answers[i]) + + im[y+50][x+50] = accum + return im + +def grad(x, y): + dist2 = (x**2 + y**2) + if dist2 == 0: + return np.array([0,0]) + dx = -y/dist2 + dy = x/dist2 + return np.array([dx, dy]) + +def vecs(segs): + Xs = np.zeros((100,100)) + Ys = np.zeros((100,100)) + for x in range(-50, 50): + for y in range(-50, 50): + answers = [] + for p1, p2 in segs: + x1, y1 = p1 + x2, y2 = p2 + angle = grad(y2-y1, x2-x1) + answers.append(grad(y-y1, x-x1)-angle) + # answers.append(-grad(y-y2, x-x2)-angle) + for i in range(len(answers)): + # pdb.set_trace() + # print(answers[i]) + dx, dy = answers[i].tolist() + Xs[y+50][x+50] += dx + Ys[y+50][x+50] += dy + + return Xs, Ys + + +segs = [ + ((-30,0),(-10,-10)), + ((-10,-10), (10,-15)), + ((10,-15), (0,-30)), + # ((20,-10), (40,-10)) + # ((0,20), (30,20)), +] +fig, axs = plt.subplots(2, 4) +print(edge_start(segs[0], segs[-1][-1])) +print(edge_end(segs[-1], segs[0][0])) + +axs[0,0].imshow(line(segs, [0]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +axs[0,1].imshow(line(segs, [5]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +axs[0,2].imshow(line(segs, [1,2,3,4]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +axs[0,3].imshow(line(segs, [0,1,2,3,4,5]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +axs[1,2].imshow(corners(segs), origin='lower', vmin=-math.pi, vmax=2*math.pi) +axs[1,3].imshow(corners(segs) + line(segs, [0,-1]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +# axs[0,1].imshow(line(antiseg, [0,1]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +# axs[1,0].imshow(line(segs, [0,1,2,3,4,5]) - line(segs, [0,5]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +# axs[1,1].imshow(line(segs, [0,1]), origin='lower', vmin=-math.pi, vmax=2*math.pi) +# gx, gy = np.gradient(line(segs, [0,1,2,3])) +# plt.quiver(gx, gy) + +# axs[1, 0].imshow(line(segs, [2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[2, 0].imshow(line(segs, [4,5]), origin='lower', vmin=-math.pi, vmax=math.pi) + +# axs[0, 1].imshow(line(segs, [0, 5]), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[1, 1].imshow(line(segs, [1, 2]), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[2, 1].imshow(line(segs, [3, 4]), origin='lower', vmin=-math.pi, vmax=math.pi) + +# axs[0, 2].imshow(line(segs, [0,1,2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) + +# axs[0, 2].imshow(line(segs, [1,2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[1, 2].imshow(line(segs, [0,1,2]), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[2, 2].imshow(line(segs, [0,1,2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) + +# axs[1].imshow(line((0,0),(0,30)), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[2].imshow(line((0,0),(30,30)), origin='lower', vmin=-math.pi, vmax=math.pi) +# axs[3].imshow(line((0,0),(30,-30)), origin='lower', vmin=-math.pi, vmax=math.pi) + +plt.show() 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() From 0d5ebad8fe0e13ffb12cc98549ea0cffbb066ed1 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Wed, 27 Dec 2023 22:50:18 -0800 Subject: [PATCH 02/25] working --- stltovoxel/playground.py | 239 +++++++++++++++++++++++---------------- 1 file changed, 139 insertions(+), 100 deletions(-) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index a5f9357..d7f52ea 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -23,8 +23,6 @@ def negatan(f1): # calculate baseline from angle of last two line segments of all incomplete loops ## Find a more pleasing way to represent baseline. Angles are good, but the asymmetry is uncomfortable. # Look at neighbors distance and angle, incorporate baseline, decide who to connect to. - -# # def is_right(line, point): a, b = line @@ -34,11 +32,15 @@ def is_right(line, point): num = ((x2 - x1)*(y1 - y0) - (x1 - x0)*(y2 - y1)) return num > 0 -def corners(segs): - im = np.zeros((100,100)) +def corners(segs, pt, step_towards): + x, y = pt + # Without a small movement toward the other point the value is undefined + # because we are at the center of a sweep + dx, dy = step_towards + x += dx/1000 + y += dy/1000 + ret = 0 for i in range(len(segs) - 1): - # pdb.set_trace() - s1 = segs[i] s2 = segs[i+1] assert s1[1] == s2[0] @@ -47,26 +49,22 @@ def corners(segs): c = np.array(list(s2[1])) l1 = a - b l2 = b - c - # pdb.set_trace() ang = math.atan2(*atandiff((l2[1], l2[0]), (l1[1], l1[0]))) - for x in range(-50, 50): - for y in range(-50, 50): - - if ang < 0 : - # Concave case - if not is_right((a, b), (x,y)) and not is_right((c, b), (x,y)): - im[y+50][x+50] += ang + 2 * math.pi - else: - im[y+50][x+50] += ang - else : - # Convex case (normal) - if is_right((a, b), (x,y)) and is_right((c, b), (x,y)): - im[y+50][x+50] += ang - 2 * math.pi - else: - im[y+50][x+50] += ang - return im + if ang < 0 : + # Concave case + if not is_right((a, b), (x,y)) and not is_right((c, b), (x,y)): + ret += ang + 2 * math.pi + else: + ret += ang + else : + # Convex case (normal) + if is_right((a, b), (x,y)) and is_right((c, b), (x,y)): + ret += ang - 2 * math.pi + else: + ret += ang + return ret def edge_start(me_seg, them_pt): @@ -84,7 +82,12 @@ def edge_end(me_seg, them_pt): angle = (s1[1] - s2[1], s1[0] - s2[0]) p1 = (start[1] - them_pt[1], start[0] - them_pt[0]) return math.atan2(*atansum(p1, negatan(angle))) - + +def edge_end_baseline(me_seg): + s1, s2 = me_seg + start = s2 + angle = (s1[1] - s2[1], s1[0] - s2[0]) + return math.atan2(*angle) def line(segs, inc): im = np.zeros((100,100)) @@ -95,18 +98,10 @@ def line(segs, inc): x1, y1 = p1 x2, y2 = p2 angle = (y2 - y1, x2 - x1) - # angle = (0,0) p1 = (y-y1, x-x1) - p2 = (y-y2, x-x2) answers.append(atansum(negatan(p1), angle)) answers.append(atansum(p2, negatan(angle))) - # answers.append(atansum( - # atansum(p1, negatan(angle)), - # atansum(negatan(p2), angle), - # ) - # ) - # answers.append(-normalize(math.atan2(y-y2, x-x2) - angle)) accum = 0 for i in inc: accum += math.atan2(*answers[i]) @@ -114,74 +109,118 @@ def line(segs, inc): im[y+50][x+50] = accum return im -def grad(x, y): - dist2 = (x**2 + y**2) - if dist2 == 0: - return np.array([0,0]) - dx = -y/dist2 - dy = x/dist2 - return np.array([dx, dy]) - -def vecs(segs): - Xs = np.zeros((100,100)) - Ys = np.zeros((100,100)) - for x in range(-50, 50): - for y in range(-50, 50): - answers = [] - for p1, p2 in segs: - x1, y1 = p1 - x2, y2 = p2 - angle = grad(y2-y1, x2-x1) - answers.append(grad(y-y1, x-x1)-angle) - # answers.append(-grad(y-y2, x-x2)-angle) - for i in range(len(answers)): - # pdb.set_trace() - # print(answers[i]) - dx, dy = answers[i].tolist() - Xs[y+50][x+50] += dx - Ys[y+50][x+50] += dy - - return Xs, Ys - + +def get_direction(pos, segs, dangling_end): + x, y = pos + accum = 0 + for seg in segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + # accum += edge_start(dangling_end, pos) + + base = edge_end_baseline(dangling_end) + return -accum + base + + +def angle_to_delta(theta): + delta_x = math.cos(theta) + delta_y = math.sin(theta) + return np.array([delta_x, delta_y]) segs = [ - ((-30,0),(-10,-10)), - ((-10,-10), (10,-15)), - ((10,-15), (0,-30)), - # ((20,-10), (40,-10)) - # ((0,20), (30,20)), + ((0,0),(5,-5)), + ((10,10),(-10,10)), + ((5,-5), (10,0)), + ((-10,8), (-10,2)), ] -fig, axs = plt.subplots(2, 4) -print(edge_start(segs[0], segs[-1][-1])) -print(edge_end(segs[-1], segs[0][0])) - -axs[0,0].imshow(line(segs, [0]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -axs[0,1].imshow(line(segs, [5]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -axs[0,2].imshow(line(segs, [1,2,3,4]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -axs[0,3].imshow(line(segs, [0,1,2,3,4,5]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -axs[1,2].imshow(corners(segs), origin='lower', vmin=-math.pi, vmax=2*math.pi) -axs[1,3].imshow(corners(segs) + line(segs, [0,-1]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -# axs[0,1].imshow(line(antiseg, [0,1]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -# axs[1,0].imshow(line(segs, [0,1,2,3,4,5]) - line(segs, [0,5]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -# axs[1,1].imshow(line(segs, [0,1]), origin='lower', vmin=-math.pi, vmax=2*math.pi) -# gx, gy = np.gradient(line(segs, [0,1,2,3])) -# plt.quiver(gx, gy) - -# axs[1, 0].imshow(line(segs, [2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[2, 0].imshow(line(segs, [4,5]), origin='lower', vmin=-math.pi, vmax=math.pi) - -# axs[0, 1].imshow(line(segs, [0, 5]), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[1, 1].imshow(line(segs, [1, 2]), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[2, 1].imshow(line(segs, [3, 4]), origin='lower', vmin=-math.pi, vmax=math.pi) - -# axs[0, 2].imshow(line(segs, [0,1,2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) - -# axs[0, 2].imshow(line(segs, [1,2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[1, 2].imshow(line(segs, [0,1,2]), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[2, 2].imshow(line(segs, [0,1,2,3]), origin='lower', vmin=-math.pi, vmax=math.pi) - -# axs[1].imshow(line((0,0),(0,30)), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[2].imshow(line((0,0),(30,30)), origin='lower', vmin=-math.pi, vmax=math.pi) -# axs[3].imshow(line((0,0),(30,-30)), origin='lower', vmin=-math.pi, vmax=math.pi) - -plt.show() + +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 + +print(find_polyline_endpoints([((0,0),(0,1)),((0,1),(0,0))])) + +# while there are some ends that need repair +tries = 3 +while find_polyline_endpoints(segs) and tries > 0: + tries -= 1 + start_to_end = find_polyline_endpoints(segs) + print(start_to_end) + # Get a value in the map + new_start = start_to_end[next(iter(start_to_end))] + + lines = [] + lines.extend(segs) + + + pos = new_start + # find the segment I am a part of + my_seg = next(filter(lambda seg: seg[1] == pos, segs)) + # find all other segments + other_segs = list(filter(lambda seg: seg[1] != pos, segs)) + + accum = 0 + for seg in other_segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + # target = accum + math.pi*0 + target = math.pi + + delta = None + for i in range(50): + angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi + delta_new = angle_to_delta(angle_forward) + if delta is None: + delta = delta_new + closest_dist = 100000 + closest_pt = None + for seg in segs: + for x,y in seg: + d = math.sqrt((y-pos[1])**2 + (x-pos[0])**2) + if d < closest_dist and d != 0: + closest_dist = d + closest_pt = (x,y) + multiplier = closest_dist / 2 + if closest_dist < 0.001: + print("breaking after ", i) + break + # multiplier = 1 + delta = delta_new + lines.append((pos, pos+(delta * multiplier))) + pos = pos + (delta * multiplier) + + + for p1, p2 in lines: + x_values = [p1[0], p2[0]] + y_values = [p1[1], p2[1]] + plt.plot(x_values, y_values, 'bo', linestyle="-") + plt.show() + + if closest_dist < 0.001: + segs.append((new_start, closest_pt)) + else: + raise "could not find pt" + From a463f2285b178fa4ae0ff7042aca4897ed98d49e Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Wed, 27 Dec 2023 23:44:53 -0800 Subject: [PATCH 03/25] working better --- stltovoxel/playground.py | 39 +++++++++++++++++++++++++++------------ 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index d7f52ea..9746ca6 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -119,7 +119,21 @@ def get_direction(pos, segs, dangling_end): # accum += edge_start(dangling_end, pos) base = edge_end_baseline(dangling_end) - return -accum + base + return base - accum + +def get_direction_iter(pos, segs, dangling_end, target): + x, y = pos + accum = 0 + for seg in segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + + # accum += edge_start(dangling_end, pos) + + base = edge_end_baseline(dangling_end) + # target - accum is a compensation factor to point back toward the ideal winding number + print(accum) + return (base - accum) #+ (target - accum) def angle_to_delta(theta): @@ -132,6 +146,7 @@ def angle_to_delta(theta): ((10,10),(-10,10)), ((5,-5), (10,0)), ((-10,8), (-10,2)), + ((-5, 4), (-10, -4)), ] def find_polyline_endpoints(segs): @@ -164,7 +179,7 @@ def find_polyline_endpoints(segs): print(find_polyline_endpoints([((0,0),(0,1)),((0,1),(0,0))])) # while there are some ends that need repair -tries = 3 +tries = 4 while find_polyline_endpoints(segs) and tries > 0: tries -= 1 start_to_end = find_polyline_endpoints(segs) @@ -191,10 +206,11 @@ def find_polyline_endpoints(segs): delta = None for i in range(50): - angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi - delta_new = angle_to_delta(angle_forward) - if delta is None: - delta = delta_new + if i == 0: + angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi + else: + angle_forward = get_direction_iter(pos, segs, (new_start, pos), target) + target + math.pi + delta = angle_to_delta(angle_forward) closest_dist = 100000 closest_pt = None for seg in segs: @@ -203,12 +219,11 @@ def find_polyline_endpoints(segs): if d < closest_dist and d != 0: closest_dist = d closest_pt = (x,y) - multiplier = closest_dist / 2 + multiplier = closest_dist / 4 if closest_dist < 0.001: print("breaking after ", i) break # multiplier = 1 - delta = delta_new lines.append((pos, pos+(delta * multiplier))) pos = pos + (delta * multiplier) @@ -219,8 +234,8 @@ def find_polyline_endpoints(segs): plt.plot(x_values, y_values, 'bo', linestyle="-") plt.show() - if closest_dist < 0.001: - segs.append((new_start, closest_pt)) - else: - raise "could not find pt" + # if closest_dist < 0.001: + segs.append((new_start, closest_pt)) + # else: + # raise "could not find pt" From aba2aba63cbc50e8856b68d9d2c79ce22cb216cf Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sun, 25 Feb 2024 20:28:32 -0800 Subject: [PATCH 04/25] working one way --- stltovoxel/playground.py | 82 +++++++++++++++++----------------------- 1 file changed, 34 insertions(+), 48 deletions(-) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index 9746ca6..ad0f54a 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -176,40 +176,29 @@ def find_polyline_endpoints(segs): return start_to_end -print(find_polyline_endpoints([((0,0),(0,1)),((0,1),(0,0))])) # while there are some ends that need repair -tries = 4 -while find_polyline_endpoints(segs) and tries > 0: - tries -= 1 +while find_polyline_endpoints(segs): start_to_end = find_polyline_endpoints(segs) print(start_to_end) - # Get a value in the map - new_start = start_to_end[next(iter(start_to_end))] - - lines = [] - lines.extend(segs) - - - pos = new_start - # find the segment I am a part of - my_seg = next(filter(lambda seg: seg[1] == pos, segs)) - # find all other segments - other_segs = list(filter(lambda seg: seg[1] != pos, segs)) - - accum = 0 - for seg in other_segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - # target = accum + math.pi*0 - target = math.pi - - delta = None - for i in range(50): - if i == 0: - angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi - else: - angle_forward = get_direction_iter(pos, segs, (new_start, pos), target) + target + math.pi + for start in start_to_end.keys(): + end = start_to_end[start] + lines = [] + lines.extend(segs) + pos = end + # find the segment I am a part of + my_seg = next(filter(lambda seg: seg[1] == pos, segs)) + # find all other segments + other_segs = list(filter(lambda seg: seg[1] != pos, segs)) + + # accum = 0 + # for seg in other_segs: + # accum += edge_start(seg, pos) + # accum += edge_end(seg, pos) + # target = accum + math.pi*0 + target = math.pi + + angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi delta = angle_to_delta(angle_forward) closest_dist = 100000 closest_pt = None @@ -219,23 +208,20 @@ def find_polyline_endpoints(segs): if d < closest_dist and d != 0: closest_dist = d closest_pt = (x,y) - multiplier = closest_dist / 4 + # multiplier = closest_dist / 4 + multiplier = 1 if closest_dist < 0.001: - print("breaking after ", i) - break - # multiplier = 1 - lines.append((pos, pos+(delta * multiplier))) - pos = pos + (delta * multiplier) - - - for p1, p2 in lines: - x_values = [p1[0], p2[0]] - y_values = [p1[1], p2[1]] - plt.plot(x_values, y_values, 'bo', linestyle="-") - plt.show() - - # if closest_dist < 0.001: - segs.append((new_start, closest_pt)) - # else: - # raise "could not find pt" + segs.append((pos, closest_pt)) + continue + else: + newpos = pos+(delta * multiplier) + segs.append((pos, tuple(newpos))) + + + for p1, p2 in lines: + x_values = [p1[0], p2[0]] + y_values = [p1[1], p2[1]] + plt.plot(x_values, y_values, 'bo', linestyle="-") + plt.show() + From d6618156b5014293bcff3fa34e23a18fe2531660 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sun, 25 Feb 2024 20:57:36 -0800 Subject: [PATCH 05/25] working both ends --- stltovoxel/playground.py | 88 +++++++++++++++++++++++++++------------- 1 file changed, 60 insertions(+), 28 deletions(-) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index ad0f54a..cdc0a59 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -70,22 +70,24 @@ def corners(segs, pt, step_towards): def edge_start(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg - start = s1 angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (start[1] - them_pt[1], start[0] - them_pt[0]) + p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) return math.atan2(*atansum(negatan(p1), angle)) def edge_end(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg - start = s2 angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (start[1] - them_pt[1], start[0] - them_pt[0]) + p1 = (s2[1] - them_pt[1], s2[0] - them_pt[0]) return math.atan2(*atansum(p1, negatan(angle))) +def edge_start_baseline(me_seg): + s1, s2 = me_seg + angle = (s2[1] - s1[1], s2[0] - s1[0]) + return math.atan2(*angle) + def edge_end_baseline(me_seg): s1, s2 = me_seg - start = s2 angle = (s1[1] - s2[1], s1[0] - s2[0]) return math.atan2(*angle) @@ -111,7 +113,6 @@ def line(segs, inc): def get_direction(pos, segs, dangling_end): - x, y = pos accum = 0 for seg in segs: accum += edge_start(seg, pos) @@ -121,8 +122,17 @@ def get_direction(pos, segs, dangling_end): base = edge_end_baseline(dangling_end) return base - accum +def get_direction_backwards(pos, segs, dangling_start): + accum = 0 + for seg in segs: + accum -= edge_start(seg, pos) + accum -= edge_end(seg, pos) + # accum += edge_start(dangling_end, pos) + + base = edge_start_baseline(dangling_start) + return base - accum + def get_direction_iter(pos, segs, dangling_end, target): - x, y = pos accum = 0 for seg in segs: accum += edge_start(seg, pos) @@ -183,33 +193,26 @@ def find_polyline_endpoints(segs): print(start_to_end) for start in start_to_end.keys(): end = start_to_end[start] - lines = [] - lines.extend(segs) + + # add the end point pos = end # find the segment I am a part of my_seg = next(filter(lambda seg: seg[1] == pos, segs)) # find all other segments other_segs = list(filter(lambda seg: seg[1] != pos, segs)) - # accum = 0 - # for seg in other_segs: - # accum += edge_start(seg, pos) - # accum += edge_end(seg, pos) - # target = accum + math.pi*0 target = math.pi - angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi delta = angle_to_delta(angle_forward) closest_dist = 100000 closest_pt = None - for seg in segs: - for x,y in seg: - d = math.sqrt((y-pos[1])**2 + (x-pos[0])**2) - if d < closest_dist and d != 0: - closest_dist = d - closest_pt = (x,y) - # multiplier = closest_dist / 4 - multiplier = 1 + for x,y in start_to_end.keys(): + d = math.sqrt((y-pos[1])**2 + (x-pos[0])**2) + if d < closest_dist and d != 0: + closest_dist = d + closest_pt = (x,y) + multiplier = closest_dist / 4 + # multiplier = 1 if closest_dist < 0.001: segs.append((pos, closest_pt)) continue @@ -217,11 +220,40 @@ def find_polyline_endpoints(segs): newpos = pos+(delta * multiplier) segs.append((pos, tuple(newpos))) + # add the start point + pos = start + # find the segment I am a part of + my_seg = next(filter(lambda seg: seg[0] == pos, segs)) + # find all other segments + other_segs = list(filter(lambda seg: seg[0] != pos, segs)) + + print(get_direction_backwards(pos, other_segs, my_seg)) + target = math.pi + angle_forward = get_direction_backwards(pos, other_segs, my_seg) + target + math.pi + print(angle_forward) + delta = angle_to_delta(angle_forward) + closest_dist = 100000 + closest_pt = None + for x,y in start_to_end.values(): + d = math.sqrt((y-pos[1])**2 + (x-pos[0])**2) + if d < closest_dist and d != 0: + closest_dist = d + closest_pt = (x,y) + multiplier = closest_dist / 4 + # multiplier = 1 + if closest_dist < 0.001: + segs.append((closest_pt, pos)) + continue + else: + newpos = pos+(delta * multiplier) + segs.append((tuple(newpos), pos)) + + for p1, p2 in segs: + x_values = [p1[0], p2[0]] + y_values = [p1[1], p2[1]] + plt.plot(x_values, y_values, 'bo', linestyle="-") + - for p1, p2 in lines: - x_values = [p1[0], p2[0]] - y_values = [p1[1], p2[1]] - plt.plot(x_values, y_values, 'bo', linestyle="-") - plt.show() + plt.show() From fc4099cee5d4d6131cd9f825a4c467b31f3e3fba Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sun, 25 Feb 2024 21:10:15 -0800 Subject: [PATCH 06/25] working with background --- stltovoxel/playground.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index cdc0a59..a80d1df 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -111,6 +111,12 @@ def line(segs, inc): im[y+50][x+50] = accum return im +def get_winding(pos, segs): + accum = 0 + for seg in segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + return accum def get_direction(pos, segs, dangling_end): accum = 0 @@ -186,11 +192,21 @@ def find_polyline_endpoints(segs): return start_to_end +background = np.zeros((200,200)) +for posxa in range(200): + for posya in range(200): + posx = (posxa - 100) * 0.2 + posy = (posya - 100) * 0.2 + winding = get_winding((posx, posy), segs) + # if abs(winding - math.pi) < .1: + # winding = math.pi*2 + background[199-posya][posxa] = winding # while there are some ends that need repair while find_polyline_endpoints(segs): start_to_end = find_polyline_endpoints(segs) print(start_to_end) + for start in start_to_end.keys(): end = start_to_end[start] @@ -254,6 +270,10 @@ def find_polyline_endpoints(segs): plt.plot(x_values, y_values, 'bo', linestyle="-") + + + plt.imshow(background, aspect='auto', cmap='viridis', extent=[-20, 20, -20, 20], vmin=-math.pi, vmax=math.pi*2) + plt.show() From f3a45727a1c98732387f4c4a58a446868d7f7e71 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Wed, 19 Jun 2024 11:56:55 -0700 Subject: [PATCH 07/25] working!!! --- stltovoxel/playground.py | 157 ++++++++++++++++++++------------------- 1 file changed, 80 insertions(+), 77 deletions(-) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index a80d1df..20337db 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -74,6 +74,13 @@ def edge_start(me_seg, them_pt): p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) return math.atan2(*atansum(negatan(p1), angle)) +def edge_start_no_atan(me_seg, them_pt): + # Starter unpaired point + s1, s2 = me_seg + angle = (s1[1] - s2[1], s1[0] - s2[0]) + p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) + return atansum(negatan(p1), angle) + def edge_end(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg @@ -158,11 +165,9 @@ def angle_to_delta(theta): return np.array([delta_x, delta_y]) segs = [ - ((0,0),(5,-5)), - ((10,10),(-10,10)), - ((5,-5), (10,0)), - ((-10,8), (-10,2)), - ((-5, 4), (-10, -4)), + ((0,0),(20,0)), + ((20,10),(0,10)), + # ((0,2.5),(0,7.5)), ] def find_polyline_endpoints(segs): @@ -198,82 +203,80 @@ def find_polyline_endpoints(segs): posx = (posxa - 100) * 0.2 posy = (posya - 100) * 0.2 winding = get_winding((posx, posy), segs) - # if abs(winding - math.pi) < .1: - # winding = math.pi*2 + if abs(winding - math.pi) < .1: + winding = math.pi*2 background[199-posya][posxa] = winding -# while there are some ends that need repair -while find_polyline_endpoints(segs): - start_to_end = find_polyline_endpoints(segs) - print(start_to_end) - - for start in start_to_end.keys(): - end = start_to_end[start] - - # add the end point - pos = end - # find the segment I am a part of - my_seg = next(filter(lambda seg: seg[1] == pos, segs)) - # find all other segments - other_segs = list(filter(lambda seg: seg[1] != pos, segs)) - - target = math.pi - angle_forward = get_direction(pos, other_segs, my_seg) + target + math.pi - delta = angle_to_delta(angle_forward) - closest_dist = 100000 - closest_pt = None - for x,y in start_to_end.keys(): - d = math.sqrt((y-pos[1])**2 + (x-pos[0])**2) - if d < closest_dist and d != 0: - closest_dist = d - closest_pt = (x,y) - multiplier = closest_dist / 4 - # multiplier = 1 - if closest_dist < 0.001: - segs.append((pos, closest_pt)) - continue - else: - newpos = pos+(delta * multiplier) - segs.append((pos, tuple(newpos))) - - # add the start point - pos = start - # find the segment I am a part of - my_seg = next(filter(lambda seg: seg[0] == pos, segs)) - # find all other segments - other_segs = list(filter(lambda seg: seg[0] != pos, segs)) - - print(get_direction_backwards(pos, other_segs, my_seg)) - target = math.pi - angle_forward = get_direction_backwards(pos, other_segs, my_seg) + target + math.pi - print(angle_forward) - delta = angle_to_delta(angle_forward) - closest_dist = 100000 - closest_pt = None - for x,y in start_to_end.values(): - d = math.sqrt((y-pos[1])**2 + (x-pos[0])**2) - if d < closest_dist and d != 0: - closest_dist = d - closest_pt = (x,y) - multiplier = closest_dist / 4 - # multiplier = 1 - if closest_dist < 0.001: - segs.append((closest_pt, pos)) - continue - else: - newpos = pos+(delta * multiplier) - segs.append((tuple(newpos), pos)) - - for p1, p2 in segs: - x_values = [p1[0], p2[0]] - y_values = [p1[1], p2[1]] - plt.plot(x_values, y_values, 'bo', linestyle="-") - +def grad(ox,oy,pt): + px,py = pt + x = ox - px + y = oy - py + gx = -((y)/(x**2 + y**2)) + gy = x/(x**2 + y**2) + return (gx, gy) + +def accum_grad_90(x, y, segs): + ax = 0 + ay = 0 + for start, end in segs: + # for pt in [start, end]: + gx, gy = grad(x+0.00001,y+0.00001,start) + ax += gy + ay -= gx + gx, gy = grad(x+0.00001,y+0.00001,end) + ax -= gy + ay += gx + return (ax, ay) + +def vecnorm(pt): + x, y = pt + dist = math.sqrt(x**2 + y**2) + return (x / dist, y / dist) - +# cx, cy = -.5,0 +# for i in range(200): +# # calculate gradient for current position +# mgx, mgy = vecnorm(accum_grad_90(cx, cy, segs)) +# # move forward +# mgxnorm = mgx * 0.1 +# mgynorm = mgy * 0.1 +# plt.plot([cx, cx + mgxnorm], [cy, cy + mgynorm], 'bo', linestyle="-") +# cx += mgxnorm +# cy += mgynorm +# # add back segment +for x in range(-10,10): + for y in range(-10,10): + ax, ay = accum_grad_90(x,y,segs) + plt.quiver(x,y, ax, ay, color=['r','b','g'], scale=21) - plt.imshow(background, aspect='auto', cmap='viridis', extent=[-20, 20, -20, 20], vmin=-math.pi, vmax=math.pi*2) - plt.show() + + + +# while there are some ends that need repair +# while find_polyline_endpoints(segs): +start_to_end = find_polyline_endpoints(segs) +# print(start_to_end) + +start = list(start_to_end.keys())[0] +end = start_to_end[start] +# find the segment I am a part of +my_seg = next(filter(lambda seg: seg[0] == start, segs)) +# find all other segments +other_segs = list(filter(lambda seg: seg[0] != start, segs)) +target = math.pi +angle_forward = get_direction_backwards(start, other_segs, my_seg) + target + math.pi +delta = angle_to_delta(angle_forward) +pos = start + (delta * 0.1) +# plt.quiver(*start, *delta, color=['r','b','g'], scale=21) +for i in range(100): + delt = vecnorm(accum_grad_90(*pos, segs)) + plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") + pos += delt + + +plt.imshow(background, aspect='auto', cmap='viridis', extent=[-20, 20, -20, 20], vmin=-math.pi, vmax=math.pi*2) + +plt.show() \ No newline at end of file From 2965291de7c8895dd8da7a2303ad3337305918d0 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 12:36:27 -0700 Subject: [PATCH 08/25] working --- stltovoxel/playground.py | 24 +-- stltovoxel/winding_query.py | 384 +++++++++++++++++++++++++----------- test/test_winding_query.py | 47 ----- 3 files changed, 270 insertions(+), 185 deletions(-) diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index 20337db..9417ae9 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -1,7 +1,6 @@ import matplotlib.pyplot as plt import numpy as np import math -import pdb def normalize(num): return ((num + math.pi) % (2*math.pi)) - math.pi @@ -165,9 +164,9 @@ def angle_to_delta(theta): return np.array([delta_x, delta_y]) segs = [ - ((0,0),(20,0)), - ((20,10),(0,10)), - # ((0,2.5),(0,7.5)), + ((0,10),(0,0)), + ((0,0),(10,0)), + ((10,5),(5,10)), ] def find_polyline_endpoints(segs): @@ -233,28 +232,11 @@ def vecnorm(pt): dist = math.sqrt(x**2 + y**2) return (x / dist, y / dist) -# cx, cy = -.5,0 -# for i in range(200): -# # calculate gradient for current position -# mgx, mgy = vecnorm(accum_grad_90(cx, cy, segs)) -# # move forward -# mgxnorm = mgx * 0.1 -# mgynorm = mgy * 0.1 -# plt.plot([cx, cx + mgxnorm], [cy, cy + mgynorm], 'bo', linestyle="-") -# cx += mgxnorm -# cy += mgynorm -# # add back segment for x in range(-10,10): for y in range(-10,10): ax, ay = accum_grad_90(x,y,segs) plt.quiver(x,y, ax, ay, color=['r','b','g'], scale=21) - - - - - - # while there are some ends that need repair # while find_polyline_endpoints(segs): start_to_end = find_polyline_endpoints(segs) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index 18e7fdf..a3dbc2a 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -1,7 +1,4 @@ import numpy as np -import math -import functools -from queue import PriorityQueue def find_polylines(segments): # noqa: C901 @@ -63,39 +60,278 @@ def find_polylines(segments): # noqa: C901 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) - +import matplotlib.pyplot as plt +import numpy as np +import math def normalize(num): return ((num + math.pi) % (2*math.pi)) - math.pi - -def signed_point_line_dist(line, point): +def atansum(f1, f2): + y, x = f1 + z, w = f2 + return (y*w + x*z, x*w - y*z) + +def atandiff(f1, f2): + y1, x1 = f1 + y2, x2 = f2 + return (y1*x2 - x1*y2, x1*x2 + y1*y2) + +def negatan(f1): + y,x = f1 + return -y, x + +# calculate baseline from angle of last two line segments of all incomplete loops +## Find a more pleasing way to represent baseline. Angles are good, but the asymmetry is uncomfortable. +# Look at neighbors distance and angle, incorporate baseline, decide who to connect to. +# +def is_right(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 - + return num > 0 + +def corners(segs, pt, step_towards): + x, y = pt + # Without a small movement toward the other point the value is undefined + # because we are at the center of a sweep + dx, dy = step_towards + x += dx/1000 + y += dy/1000 + ret = 0 + for i in range(len(segs) - 1): + s1 = segs[i] + s2 = segs[i+1] + assert s1[1] == s2[0] + a = np.array(list(s1[0])) + b = np.array(list(s1[1])) + c = np.array(list(s2[1])) + l1 = a - b + l2 = b - c + + ang = math.atan2(*atandiff((l2[1], l2[0]), (l1[1], l1[0]))) + + if ang < 0 : + # Concave case + if not is_right((a, b), (x,y)) and not is_right((c, b), (x,y)): + ret += ang + 2 * math.pi + else: + ret += ang + else : + # Convex case (normal) + if is_right((a, b), (x,y)) and is_right((c, b), (x,y)): + ret += ang - 2 * math.pi + else: + ret += ang + return ret + + +def edge_start(me_seg, them_pt): + # Starter unpaired point + s1, s2 = me_seg + angle = (s1[1] - s2[1], s1[0] - s2[0]) + p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) + return math.atan2(*atansum(negatan(p1), angle)) + +def edge_start_no_atan(me_seg, them_pt): + # Starter unpaired point + s1, s2 = me_seg + angle = (s1[1] - s2[1], s1[0] - s2[0]) + p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) + return atansum(negatan(p1), angle) + +def edge_end(me_seg, them_pt): + # Starter unpaired point + s1, s2 = me_seg + angle = (s1[1] - s2[1], s1[0] - s2[0]) + p1 = (s2[1] - them_pt[1], s2[0] - them_pt[0]) + return math.atan2(*atansum(p1, negatan(angle))) + +def edge_start_baseline(me_seg): + s1, s2 = me_seg + angle = (s2[1] - s1[1], s2[0] - s1[0]) + return math.atan2(*angle) + +def edge_end_baseline(me_seg): + s1, s2 = me_seg + angle = (s1[1] - s2[1], s1[0] - s2[0]) + return math.atan2(*angle) + +def line(segs, inc): + im = np.zeros((100,100)) + for x in range(-50, 50): + for y in range(-50, 50): + answers = [] + for p1, p2 in segs: + x1, y1 = p1 + x2, y2 = p2 + angle = (y2 - y1, x2 - x1) + p1 = (y-y1, x-x1) + p2 = (y-y2, x-x2) + answers.append(atansum(negatan(p1), angle)) + answers.append(atansum(p2, negatan(angle))) + accum = 0 + for i in inc: + accum += math.atan2(*answers[i]) + + im[y+50][x+50] = accum + return im + +def get_winding(pos, segs): + accum = 0 + for seg in segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + return accum + +def get_direction(pos, segs, dangling_end): + accum = 0 + for seg in segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + # accum += edge_start(dangling_end, pos) + + base = edge_end_baseline(dangling_end) + return base - accum + +def get_direction_backwards(pos, segs, dangling_start): + accum = 0 + for seg in segs: + accum -= edge_start(seg, pos) + accum -= edge_end(seg, pos) + # accum += edge_start(dangling_end, pos) + + base = edge_start_baseline(dangling_start) + return base - accum + +def get_direction_iter(pos, segs, dangling_end, target): + accum = 0 + for seg in segs: + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) + + # accum += edge_start(dangling_end, pos) + + base = edge_end_baseline(dangling_end) + # target - accum is a compensation factor to point back toward the ideal winding number + print(accum) + return (base - accum) #+ (target - accum) + + +def angle_to_delta(theta): + delta_x = math.cos(theta) + delta_y = math.sin(theta) + return np.array([delta_x, delta_y]) + +segs = [ + ((0,10),(0,0)), + ((0,0),(10,0)), + ((10,5),(5,10)), +] + +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 + +background = np.zeros((200,200)) +for posxa in range(200): + for posya in range(200): + posx = (posxa - 100) * 0.2 + posy = (posya - 100) * 0.2 + winding = get_winding((posx, posy), segs) + if abs(winding - math.pi) < .1: + winding = math.pi*2 + background[199-posya][posxa] = winding + +def grad(ox,oy,pt): + px,py = pt + x = ox - px + y = oy - py + gx = -((y)/(x**2 + y**2)) + gy = x/(x**2 + y**2) + return (gx, gy) + +def accum_grad_90(x, y, segs): + ax = 0 + ay = 0 + for start, end in segs: + # for pt in [start, end]: + gx, gy = grad(x+0.00001,y+0.00001,start) + ax += gy + ay -= gx + gx, gy = grad(x+0.00001,y+0.00001,end) + ax -= gy + ay += gx + return (ax, ay) + +def vecnorm(pt): + x, y = pt + dist = math.sqrt(x**2 + y**2) + return (x / dist, y / dist) -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 +def dist(p1, p2): + x1, y1 = p1 + x2, y2 = p2 + return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) +def find_initial_angle(start, other_segs, my_seg): + target = math.pi + angle_forward = get_direction_backwards(start, other_segs, my_seg) + target + math.pi + delta = angle_to_delta(angle_forward) + return delta + +def grad_90_norm(pos, segs): + delt = vecnorm(accum_grad_90(*pos, segs)) + return delt + +def find_flow(start, ends, segs): + for seg in segs: + plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") + + # find the segment I am a part of + my_seg = next(filter(lambda seg: seg[0] == start, segs)) + # find all other segments + other_segs = list(filter(lambda seg: seg[0] != start, segs)) + delta = find_initial_angle(start, other_segs, my_seg) + pos = start + (delta * 0.1) + # plt.quiver(*start, *delta, color=['r','b','g'], scale=21) + for _ in range(200): + delt = grad_90_norm(pos, segs) + plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") + pos += delt + for end in ends: + print(pos, end, dist(pos, end)) + if dist(pos, end) < 1: + plt.show() + return end + + plt.show() + raise "Flow not found" class WindingQuery(): def __init__(self, segments): @@ -125,98 +361,12 @@ def repair_all(self): def repair_segment(self): # Search starts at the end of a polyline - start = self.polylines[0][-1] + start = self.polylines[0][0] # 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) + endpoints = [polyline[-1] for polyline in self.polylines] + end = find_flow(start, endpoints, self.original_segments) + print(self.original_segments) + new_segment = (end, start) + print(new_segment) 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_winding_query.py b/test/test_winding_query.py index 51f2598..a6f5246 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -5,20 +5,6 @@ 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))] @@ -46,39 +32,6 @@ def test_find_polylines_out_of_order(self): 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() From 18574fd68c2364f1d86e19237d0b2bc7ebbe84f2 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 13:13:53 -0700 Subject: [PATCH 09/25] tests --- test/test_winding_query.py | 61 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/test/test_winding_query.py b/test/test_winding_query.py index a6f5246..04c12b4 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -1,10 +1,14 @@ import unittest import math +import numpy as np from stltovoxel import winding_query class TestWindingQuery(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))] @@ -32,6 +36,63 @@ def test_find_polylines_out_of_order(self): actual_polylines = winding_query.find_polylines(line_segments) self.assertEqual(actual_polylines, expected_polylines) + def test_get_direction(self): + my_seg = ((0,0),(0,1)) + other_segs = [((0,1),(1,1)),((1,1),(1,0))] + actual = winding_query.find_initial_angle(my_seg[0], other_segs, my_seg) + actual = tuple(actual) + expected = (1.0,0) + self.tuples_almost_equal(actual, expected) + + def test_get_direction2(self): + my_seg = ((1,1),(0,1)) + other_segs = [((0,0),(1,0))] + actual = winding_query.find_initial_angle(my_seg[0], other_segs, my_seg) + actual = tuple(actual) + expected = winding_query.vecnorm((-1,-1)) + self.tuples_almost_equal(actual, expected) + + def test_get_direction3(self): + my_seg = ((0,0),(1,0)) + other_segs = [((1,0),(1,1))] + actual = winding_query.find_initial_angle(my_seg[0], other_segs, my_seg) + actual = tuple(actual) + expected = winding_query.vecnorm((1,1)) + self.tuples_almost_equal(actual, expected) + + def test_grad_90_norm(self): + segs = [((0,0),(1,0)), ((1,0),(1,1))] + pos = np.array((0,0)) + np.array((0.1,0.1)) + # Treats starts as repellers and ends as attractors + actual = winding_query.grad_90_norm(pos, segs) + expected = winding_query.vecnorm((1,1)) + self.tuples_almost_equal(actual, expected) + + + def test_grad_90_norm2(self): + segs = [((0,0),(1,0)), ((1,0),(1,1)), ((1,1),(0,1))] + pos = (0,0.5) + actual = winding_query.grad_90_norm(pos, segs) + expected = winding_query.vecnorm((0,1)) + self.tuples_almost_equal(actual, expected) + + pos = (-.5,0.5) + actual = winding_query.grad_90_norm(pos, segs) + expected = winding_query.vecnorm((0,1)) + self.tuples_almost_equal(actual, expected) + + pos = (.5,0.5) + actual = winding_query.grad_90_norm(pos, segs) + expected = winding_query.vecnorm((0,1)) + self.tuples_almost_equal(actual, expected) + + def test_grad_90_norm3(self): + segs = [((0,0),(1,0)), ((1,1),(0,1))] + pos = (0.00001, 0.00001) + actual = winding_query.grad_90_norm(pos, segs) + expected = winding_query.vecnorm((1,1)) + self.tuples_almost_equal(actual, expected) + if __name__ == '__main__': unittest.main() From d0841002b107f12bfc2a7ffa899101e3aa54828a Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 13:22:29 -0700 Subject: [PATCH 10/25] works --- stltovoxel/winding_query.py | 40 ++++++++----------------------------- 1 file changed, 8 insertions(+), 32 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index a3dbc2a..fe52141 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -197,14 +197,12 @@ def get_direction(pos, segs, dangling_end): return base - accum def get_direction_backwards(pos, segs, dangling_start): - accum = 0 + accum = edge_start_baseline(dangling_start) for seg in segs: - accum -= edge_start(seg, pos) - accum -= edge_end(seg, pos) - # accum += edge_start(dangling_end, pos) + accum += edge_start(seg, pos) + accum += edge_end(seg, pos) - base = edge_start_baseline(dangling_start) - return base - accum + return accum def get_direction_iter(pos, segs, dangling_end, target): accum = 0 @@ -219,18 +217,6 @@ def get_direction_iter(pos, segs, dangling_end, target): print(accum) return (base - accum) #+ (target - accum) - -def angle_to_delta(theta): - delta_x = math.cos(theta) - delta_y = math.sin(theta) - return np.array([delta_x, delta_y]) - -segs = [ - ((0,10),(0,0)), - ((0,0),(10,0)), - ((10,5),(5,10)), -] - def find_polyline_endpoints(segs): start_to_end = dict() end_to_start = dict() @@ -258,16 +244,6 @@ def find_polyline_endpoints(segs): return start_to_end -background = np.zeros((200,200)) -for posxa in range(200): - for posya in range(200): - posx = (posxa - 100) * 0.2 - posy = (posya - 100) * 0.2 - winding = get_winding((posx, posy), segs) - if abs(winding - math.pi) < .1: - winding = math.pi*2 - background[199-posya][posxa] = winding - def grad(ox,oy,pt): px,py = pt x = ox - px @@ -300,9 +276,10 @@ def dist(p1, p2): return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) def find_initial_angle(start, other_segs, my_seg): - target = math.pi - angle_forward = get_direction_backwards(start, other_segs, my_seg) + target + math.pi - delta = angle_to_delta(angle_forward) + angle_forward = get_direction_backwards(start, other_segs, my_seg) + delta_x = math.cos(angle_forward) + delta_y = math.sin(angle_forward) + delta = np.array([delta_x, delta_y]) return delta def grad_90_norm(pos, segs): @@ -319,7 +296,6 @@ def find_flow(start, ends, segs): other_segs = list(filter(lambda seg: seg[0] != start, segs)) delta = find_initial_angle(start, other_segs, my_seg) pos = start + (delta * 0.1) - # plt.quiver(*start, *delta, color=['r','b','g'], scale=21) for _ in range(200): delt = grad_90_norm(pos, segs) plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") From 651dd3131cfcfb688854aa35ac1944d1b66c345d Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 13:25:31 -0700 Subject: [PATCH 11/25] works --- stltovoxel/winding_query.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index fe52141..e132369 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -133,7 +133,7 @@ def edge_start(me_seg, them_pt): s1, s2 = me_seg angle = (s1[1] - s2[1], s1[0] - s2[0]) p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) - return math.atan2(*atansum(negatan(p1), angle)) + return atansum(negatan(p1), angle) def edge_start_no_atan(me_seg, them_pt): # Starter unpaired point @@ -147,12 +147,12 @@ def edge_end(me_seg, them_pt): s1, s2 = me_seg angle = (s1[1] - s2[1], s1[0] - s2[0]) p1 = (s2[1] - them_pt[1], s2[0] - them_pt[0]) - return math.atan2(*atansum(p1, negatan(angle))) + return atansum(p1, negatan(angle)) def edge_start_baseline(me_seg): s1, s2 = me_seg angle = (s2[1] - s1[1], s2[0] - s1[0]) - return math.atan2(*angle) + return angle def edge_end_baseline(me_seg): s1, s2 = me_seg @@ -199,10 +199,10 @@ def get_direction(pos, segs, dangling_end): def get_direction_backwards(pos, segs, dangling_start): accum = edge_start_baseline(dangling_start) for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) + accum = atansum(accum, edge_start(seg, pos)) + accum = atansum(accum, edge_end(seg, pos)) - return accum + return math.atan2(*accum) def get_direction_iter(pos, segs, dangling_end, target): accum = 0 From 4523c031ab3ee55702ed97f251fa67a9fa3f70ea Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 13:28:24 -0700 Subject: [PATCH 12/25] works --- stltovoxel/winding_query.py | 112 +----------------------------------- 1 file changed, 3 insertions(+), 109 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index e132369..bcd9b4c 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -81,53 +81,6 @@ def negatan(f1): y,x = f1 return -y, x -# calculate baseline from angle of last two line segments of all incomplete loops -## Find a more pleasing way to represent baseline. Angles are good, but the asymmetry is uncomfortable. -# Look at neighbors distance and angle, incorporate baseline, decide who to connect to. -# -def is_right(line, point): - a, b = line - x1, y1 = a - x2, y2 = b - x0, y0 = point - num = ((x2 - x1)*(y1 - y0) - (x1 - x0)*(y2 - y1)) - return num > 0 - -def corners(segs, pt, step_towards): - x, y = pt - # Without a small movement toward the other point the value is undefined - # because we are at the center of a sweep - dx, dy = step_towards - x += dx/1000 - y += dy/1000 - ret = 0 - for i in range(len(segs) - 1): - s1 = segs[i] - s2 = segs[i+1] - assert s1[1] == s2[0] - a = np.array(list(s1[0])) - b = np.array(list(s1[1])) - c = np.array(list(s2[1])) - l1 = a - b - l2 = b - c - - ang = math.atan2(*atandiff((l2[1], l2[0]), (l1[1], l1[0]))) - - if ang < 0 : - # Concave case - if not is_right((a, b), (x,y)) and not is_right((c, b), (x,y)): - ret += ang + 2 * math.pi - else: - ret += ang - else : - # Convex case (normal) - if is_right((a, b), (x,y)) and is_right((c, b), (x,y)): - ret += ang - 2 * math.pi - else: - ret += ang - return ret - - def edge_start(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg @@ -135,13 +88,6 @@ def edge_start(me_seg, them_pt): p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) return atansum(negatan(p1), angle) -def edge_start_no_atan(me_seg, them_pt): - # Starter unpaired point - s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) - return atansum(negatan(p1), angle) - def edge_end(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg @@ -159,63 +105,13 @@ def edge_end_baseline(me_seg): angle = (s1[1] - s2[1], s1[0] - s2[0]) return math.atan2(*angle) -def line(segs, inc): - im = np.zeros((100,100)) - for x in range(-50, 50): - for y in range(-50, 50): - answers = [] - for p1, p2 in segs: - x1, y1 = p1 - x2, y2 = p2 - angle = (y2 - y1, x2 - x1) - p1 = (y-y1, x-x1) - p2 = (y-y2, x-x2) - answers.append(atansum(negatan(p1), angle)) - answers.append(atansum(p2, negatan(angle))) - accum = 0 - for i in inc: - accum += math.atan2(*answers[i]) - - im[y+50][x+50] = accum - return im - -def get_winding(pos, segs): - accum = 0 - for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - return accum - -def get_direction(pos, segs, dangling_end): - accum = 0 - for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - # accum += edge_start(dangling_end, pos) - - base = edge_end_baseline(dangling_end) - return base - accum - def get_direction_backwards(pos, segs, dangling_start): accum = edge_start_baseline(dangling_start) for seg in segs: accum = atansum(accum, edge_start(seg, pos)) accum = atansum(accum, edge_end(seg, pos)) - - return math.atan2(*accum) - -def get_direction_iter(pos, segs, dangling_end, target): - accum = 0 - for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - - # accum += edge_start(dangling_end, pos) - - base = edge_end_baseline(dangling_end) - # target - accum is a compensation factor to point back toward the ideal winding number - print(accum) - return (base - accum) #+ (target - accum) + accum = vecnorm(accum) + return accum def find_polyline_endpoints(segs): start_to_end = dict() @@ -276,9 +172,7 @@ def dist(p1, p2): return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) def find_initial_angle(start, other_segs, my_seg): - angle_forward = get_direction_backwards(start, other_segs, my_seg) - delta_x = math.cos(angle_forward) - delta_y = math.sin(angle_forward) + delta_y, delta_x = get_direction_backwards(start, other_segs, my_seg) delta = np.array([delta_x, delta_y]) return delta From 1b44cea8c7f5e49b26993d3142334bebc2a1be3d Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 13:47:02 -0700 Subject: [PATCH 13/25] working --- stltovoxel/winding_query.py | 56 +++++++++++-------------------------- 1 file changed, 17 insertions(+), 39 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index bcd9b4c..b171205 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -1,5 +1,6 @@ import numpy as np - +import matplotlib.pyplot as plt +import math def find_polylines(segments): # noqa: C901 polylines = [] @@ -60,53 +61,31 @@ def find_polylines(segments): # noqa: C901 return polylines -import matplotlib.pyplot as plt -import numpy as np -import math - -def normalize(num): - return ((num + math.pi) % (2*math.pi)) - math.pi - def atansum(f1, f2): - y, x = f1 - z, w = f2 - return (y*w + x*z, x*w - y*z) - -def atandiff(f1, f2): - y1, x1 = f1 - y2, x2 = f2 - return (y1*x2 - x1*y2, x1*x2 + y1*y2) - -def negatan(f1): - y,x = f1 - return -y, x + x,y = f1 + w,z = f2 + return ( x*w - y*z, y*w + x*z,) def edge_start(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) - return atansum(negatan(p1), angle) + angle = ( s1[0] - s2[0], s1[1] - s2[1],) + p1 = ( s1[0] - them_pt[0],them_pt[1] - s1[1],) + return atansum(p1, angle) def edge_end(me_seg, them_pt): # Starter unpaired point s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (s2[1] - them_pt[1], s2[0] - them_pt[0]) - return atansum(p1, negatan(angle)) + angle = ( s1[0] - s2[0], s2[1] - s1[1],) + p1 = (s2[0] - them_pt[0], s2[1] - them_pt[1], ) + return atansum(p1, angle) -def edge_start_baseline(me_seg): - s1, s2 = me_seg - angle = (s2[1] - s1[1], s2[0] - s1[0]) +def subtract(s1, s2): + angle = (s1[0] - s2[0], s1[1] - s2[1]) return angle -def edge_end_baseline(me_seg): - s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - return math.atan2(*angle) - def get_direction_backwards(pos, segs, dangling_start): - accum = edge_start_baseline(dangling_start) + accum = subtract(dangling_start[1], dangling_start[0]) for seg in segs: accum = atansum(accum, edge_start(seg, pos)) accum = atansum(accum, edge_end(seg, pos)) @@ -152,11 +131,10 @@ def accum_grad_90(x, y, segs): ax = 0 ay = 0 for start, end in segs: - # for pt in [start, end]: - gx, gy = grad(x+0.00001,y+0.00001,start) + gx, gy = grad(x,y,start) ax += gy ay -= gx - gx, gy = grad(x+0.00001,y+0.00001,end) + gx, gy = grad(x,y,end) ax -= gy ay += gx return (ax, ay) @@ -172,7 +150,7 @@ def dist(p1, p2): return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) def find_initial_angle(start, other_segs, my_seg): - delta_y, delta_x = get_direction_backwards(start, other_segs, my_seg) + delta_x, delta_y = get_direction_backwards(start, other_segs, my_seg) delta = np.array([delta_x, delta_y]) return delta From 04e279287fde8e942f00e583ad2b81728556db1f Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 14:15:40 -0700 Subject: [PATCH 14/25] simplified --- stltovoxel/winding_query.py | 71 ++++++++++++++----------------------- test/test_winding_query.py | 10 +++--- 2 files changed, 32 insertions(+), 49 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index b171205..b9078b4 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -66,29 +66,24 @@ def atansum(f1, f2): w,z = f2 return ( x*w - y*z, y*w + x*z,) -def edge_start(me_seg, them_pt): - # Starter unpaired point - s1, s2 = me_seg - angle = ( s1[0] - s2[0], s1[1] - s2[1],) - p1 = ( s1[0] - them_pt[0],them_pt[1] - s1[1],) - return atansum(p1, angle) - -def edge_end(me_seg, them_pt): - # Starter unpaired point - s1, s2 = me_seg - angle = ( s1[0] - s2[0], s2[1] - s1[1],) - p1 = (s2[0] - them_pt[0], s2[1] - them_pt[1], ) - return atansum(p1, angle) +def negatan(f1): + x,y = f1 + return x,-y def subtract(s1, s2): - angle = (s1[0] - s2[0], s1[1] - s2[1]) - return angle + return (s1[0] - s2[0], s1[1] - s2[1]) + +def add(s1, s2): + return (s1[0] + s2[0], s1[1] + s2[1]) def get_direction_backwards(pos, segs, dangling_start): accum = subtract(dangling_start[1], dangling_start[0]) for seg in segs: - accum = atansum(accum, edge_start(seg, pos)) - accum = atansum(accum, edge_end(seg, pos)) + s1, s2 = seg + p1 = subtract(s1, pos) + p2 = subtract(s2, pos) + accum = atansum(accum, p2) + accum = atansum(accum, negatan(p1)) accum = vecnorm(accum) return accum @@ -119,26 +114,13 @@ def find_polyline_endpoints(segs): return start_to_end -def grad(ox,oy,pt): - px,py = pt - x = ox - px - y = oy - py - gx = -((y)/(x**2 + y**2)) - gy = x/(x**2 + y**2) +def winding_contour(pos,pt): + x,y = subtract(pos, pt) + dist = (x**2 + y**2) + gx = x/dist + gy = y/dist return (gx, gy) -def accum_grad_90(x, y, segs): - ax = 0 - ay = 0 - for start, end in segs: - gx, gy = grad(x,y,start) - ax += gy - ay -= gx - gx, gy = grad(x,y,end) - ax -= gy - ay += gx - return (ax, ay) - def vecnorm(pt): x, y = pt dist = math.sqrt(x**2 + y**2) @@ -150,13 +132,16 @@ def dist(p1, p2): return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) def find_initial_angle(start, other_segs, my_seg): - delta_x, delta_y = get_direction_backwards(start, other_segs, my_seg) - delta = np.array([delta_x, delta_y]) - return delta + return np.array(get_direction_backwards(start, other_segs, my_seg)) -def grad_90_norm(pos, segs): - delt = vecnorm(accum_grad_90(*pos, segs)) - return delt +def total_winding_contour(pos, segs): + accum = (0,0) + for start, end in segs: + start_grad = winding_contour(pos, start) + end_grad = winding_contour(pos, end) + accum = add(accum, start_grad) + accum = subtract(accum, end_grad) + return vecnorm(accum) def find_flow(start, ends, segs): for seg in segs: @@ -169,7 +154,7 @@ def find_flow(start, ends, segs): delta = find_initial_angle(start, other_segs, my_seg) pos = start + (delta * 0.1) for _ in range(200): - delt = grad_90_norm(pos, segs) + delt = total_winding_contour(pos, segs) plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") pos += delt for end in ends: @@ -214,7 +199,5 @@ def repair_segment(self): # Search will conclude when it finds the beginning of any polyline (including itself) endpoints = [polyline[-1] for polyline in self.polylines] end = find_flow(start, endpoints, self.original_segments) - print(self.original_segments) new_segment = (end, start) - print(new_segment) self.original_segments.append(new_segment) diff --git a/test/test_winding_query.py b/test/test_winding_query.py index 04c12b4..1d32d0c 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -64,7 +64,7 @@ def test_grad_90_norm(self): segs = [((0,0),(1,0)), ((1,0),(1,1))] pos = np.array((0,0)) + np.array((0.1,0.1)) # Treats starts as repellers and ends as attractors - actual = winding_query.grad_90_norm(pos, segs) + actual = winding_query.total_winding_contour(pos, segs) expected = winding_query.vecnorm((1,1)) self.tuples_almost_equal(actual, expected) @@ -72,24 +72,24 @@ def test_grad_90_norm(self): def test_grad_90_norm2(self): segs = [((0,0),(1,0)), ((1,0),(1,1)), ((1,1),(0,1))] pos = (0,0.5) - actual = winding_query.grad_90_norm(pos, segs) + actual = winding_query.total_winding_contour(pos, segs) expected = winding_query.vecnorm((0,1)) self.tuples_almost_equal(actual, expected) pos = (-.5,0.5) - actual = winding_query.grad_90_norm(pos, segs) + actual = winding_query.total_winding_contour(pos, segs) expected = winding_query.vecnorm((0,1)) self.tuples_almost_equal(actual, expected) pos = (.5,0.5) - actual = winding_query.grad_90_norm(pos, segs) + actual = winding_query.total_winding_contour(pos, segs) expected = winding_query.vecnorm((0,1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm3(self): segs = [((0,0),(1,0)), ((1,1),(0,1))] pos = (0.00001, 0.00001) - actual = winding_query.grad_90_norm(pos, segs) + actual = winding_query.total_winding_contour(pos, segs) expected = winding_query.vecnorm((1,1)) self.tuples_almost_equal(actual, expected) From 0583b699b62375fd5e7d6baed49bcc242d47d8c8 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 15:37:56 -0700 Subject: [PATCH 15/25] more polish --- stltovoxel/winding_query.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index b9078b4..a8882ae 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -144,26 +144,23 @@ def total_winding_contour(pos, segs): return vecnorm(accum) def find_flow(start, ends, segs): - for seg in segs: - plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") - # find the segment I am a part of my_seg = next(filter(lambda seg: seg[0] == start, segs)) # find all other segments other_segs = list(filter(lambda seg: seg[0] != start, segs)) delta = find_initial_angle(start, other_segs, my_seg) pos = start + (delta * 0.1) + seg_outs = [(tuple(pos), tuple(start), )] for _ in range(200): delt = total_winding_contour(pos, segs) - plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") + seg_outs.append((tuple(pos + delt), tuple(pos), )) pos += delt for end in ends: print(pos, end, dist(pos, end)) if dist(pos, end) < 1: - plt.show() - return end + seg_outs.append((tuple(end), tuple(pos))) + return seg_outs - plt.show() raise "Flow not found" class WindingQuery(): @@ -190,6 +187,9 @@ def repair_all(self): old_seg_length = len(self.polylines) self.collapse_segments() assert old_seg_length - 1 == len(self.polylines) + # for seg in self.original_segments: + # plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") + # plt.show() assert len(self.polylines) == 0 def repair_segment(self): @@ -198,6 +198,5 @@ def repair_segment(self): # Search will conclude when it finds the beginning of any polyline (including itself) endpoints = [polyline[-1] for polyline in self.polylines] - end = find_flow(start, endpoints, self.original_segments) - new_segment = (end, start) - self.original_segments.append(new_segment) + new_segs = find_flow(start, endpoints, self.original_segments) + self.original_segments.extend(new_segs) From 86e2089e5c2e68e292c37b42211de12ed6144c66 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 16:11:37 -0700 Subject: [PATCH 16/25] polish --- stltovoxel/winding_query.py | 29 +++++++++++++---------------- test/test_winding_query.py | 6 +++--- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index a8882ae..63a108e 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -76,17 +76,6 @@ def subtract(s1, s2): def add(s1, s2): return (s1[0] + s2[0], s1[1] + s2[1]) -def get_direction_backwards(pos, segs, dangling_start): - accum = subtract(dangling_start[1], dangling_start[0]) - for seg in segs: - s1, s2 = seg - p1 = subtract(s1, pos) - p2 = subtract(s2, pos) - accum = atansum(accum, p2) - accum = atansum(accum, negatan(p1)) - accum = vecnorm(accum) - return accum - def find_polyline_endpoints(segs): start_to_end = dict() end_to_start = dict() @@ -131,8 +120,17 @@ def dist(p1, p2): x2, y2 = p2 return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) -def find_initial_angle(start, other_segs, my_seg): - return np.array(get_direction_backwards(start, other_segs, my_seg)) +def find_initial_angle(my_seg, other_segs): + start, end = my_seg + accum = subtract(end, start) + for seg in other_segs: + s1, s2 = seg + p1 = subtract(s1, start) + p2 = subtract(s2, start) + accum = atansum(accum, p2) + accum = atansum(accum, negatan(p1)) + accum = vecnorm(accum) + return np.array(accum) def total_winding_contour(pos, segs): accum = (0,0) @@ -148,15 +146,14 @@ def find_flow(start, ends, segs): my_seg = next(filter(lambda seg: seg[0] == start, segs)) # find all other segments other_segs = list(filter(lambda seg: seg[0] != start, segs)) - delta = find_initial_angle(start, other_segs, my_seg) + delta = find_initial_angle(my_seg, other_segs) pos = start + (delta * 0.1) seg_outs = [(tuple(pos), tuple(start), )] for _ in range(200): delt = total_winding_contour(pos, segs) - seg_outs.append((tuple(pos + delt), tuple(pos), )) + seg_outs.append((tuple(pos + delt), tuple(pos))) pos += delt for end in ends: - print(pos, end, dist(pos, end)) if dist(pos, end) < 1: seg_outs.append((tuple(end), tuple(pos))) return seg_outs diff --git a/test/test_winding_query.py b/test/test_winding_query.py index 1d32d0c..7f90b14 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -39,7 +39,7 @@ def test_find_polylines_out_of_order(self): def test_get_direction(self): my_seg = ((0,0),(0,1)) other_segs = [((0,1),(1,1)),((1,1),(1,0))] - actual = winding_query.find_initial_angle(my_seg[0], other_segs, my_seg) + actual = winding_query.find_initial_angle(my_seg, other_segs) actual = tuple(actual) expected = (1.0,0) self.tuples_almost_equal(actual, expected) @@ -47,7 +47,7 @@ def test_get_direction(self): def test_get_direction2(self): my_seg = ((1,1),(0,1)) other_segs = [((0,0),(1,0))] - actual = winding_query.find_initial_angle(my_seg[0], other_segs, my_seg) + actual = winding_query.find_initial_angle(my_seg, other_segs) actual = tuple(actual) expected = winding_query.vecnorm((-1,-1)) self.tuples_almost_equal(actual, expected) @@ -55,7 +55,7 @@ def test_get_direction2(self): def test_get_direction3(self): my_seg = ((0,0),(1,0)) other_segs = [((1,0),(1,1))] - actual = winding_query.find_initial_angle(my_seg[0], other_segs, my_seg) + actual = winding_query.find_initial_angle(my_seg, other_segs) actual = tuple(actual) expected = winding_query.vecnorm((1,1)) self.tuples_almost_equal(actual, expected) From 00686cd4cacdffea4ec70dc90cf9b45a5423e434 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 16:22:50 -0700 Subject: [PATCH 17/25] flip start/end --- stltovoxel/winding_query.py | 34 +++++++++++++++++----------------- test/test_winding_query.py | 24 ++++++++++++------------ 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index 63a108e..2967baf 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -122,13 +122,13 @@ def dist(p1, p2): def find_initial_angle(my_seg, other_segs): start, end = my_seg - accum = subtract(end, start) + accum = subtract(start, end) for seg in other_segs: s1, s2 = seg - p1 = subtract(s1, start) - p2 = subtract(s2, start) - accum = atansum(accum, p2) - accum = atansum(accum, negatan(p1)) + p1 = subtract(s1, end) + p2 = subtract(s2, end) + accum = atansum(accum, p1) + accum = atansum(accum, negatan(p2)) accum = vecnorm(accum) return np.array(accum) @@ -137,25 +137,25 @@ def total_winding_contour(pos, segs): for start, end in segs: start_grad = winding_contour(pos, start) end_grad = winding_contour(pos, end) - accum = add(accum, start_grad) - accum = subtract(accum, end_grad) + accum = subtract(accum, start_grad) + accum = add(accum, end_grad) return vecnorm(accum) def find_flow(start, ends, segs): # find the segment I am a part of - my_seg = next(filter(lambda seg: seg[0] == start, segs)) + my_seg = next(filter(lambda seg: seg[1] == start, segs)) # find all other segments - other_segs = list(filter(lambda seg: seg[0] != start, segs)) + other_segs = list(filter(lambda seg: seg[1] != start, segs)) delta = find_initial_angle(my_seg, other_segs) pos = start + (delta * 0.1) - seg_outs = [(tuple(pos), tuple(start), )] + seg_outs = [(tuple(start), tuple(pos))] for _ in range(200): delt = total_winding_contour(pos, segs) - seg_outs.append((tuple(pos + delt), tuple(pos))) + seg_outs.append((tuple(pos), tuple(pos + delt))) pos += delt for end in ends: if dist(pos, end) < 1: - seg_outs.append((tuple(end), tuple(pos))) + seg_outs.append((tuple(pos), tuple(end))) return seg_outs raise "Flow not found" @@ -184,16 +184,16 @@ def repair_all(self): old_seg_length = len(self.polylines) self.collapse_segments() assert old_seg_length - 1 == len(self.polylines) - # for seg in self.original_segments: - # plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") - # plt.show() + for seg in self.original_segments: + plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") + plt.show() assert len(self.polylines) == 0 def repair_segment(self): # Search starts at the end of a polyline - start = self.polylines[0][0] + start = self.polylines[0][-1] # Search will conclude when it finds the beginning of any polyline (including itself) - endpoints = [polyline[-1] for polyline in self.polylines] + endpoints = [polyline[0] for polyline in self.polylines] new_segs = find_flow(start, endpoints, self.original_segments) self.original_segments.extend(new_segs) diff --git a/test/test_winding_query.py b/test/test_winding_query.py index 7f90b14..ac1278d 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -37,24 +37,24 @@ def test_find_polylines_out_of_order(self): self.assertEqual(actual_polylines, expected_polylines) def test_get_direction(self): - my_seg = ((0,0),(0,1)) - other_segs = [((0,1),(1,1)),((1,1),(1,0))] + other_segs = [((1,0),(1,1)),((1,1),(0,1)),] + my_seg = ((0,1),(0,0)) actual = winding_query.find_initial_angle(my_seg, other_segs) actual = tuple(actual) expected = (1.0,0) self.tuples_almost_equal(actual, expected) def test_get_direction2(self): - my_seg = ((1,1),(0,1)) - other_segs = [((0,0),(1,0))] + other_segs = [((1,0),(0,0))] + my_seg = ((0,1),(1,1)) actual = winding_query.find_initial_angle(my_seg, other_segs) actual = tuple(actual) expected = winding_query.vecnorm((-1,-1)) self.tuples_almost_equal(actual, expected) def test_get_direction3(self): - my_seg = ((0,0),(1,0)) - other_segs = [((1,0),(1,1))] + other_segs = [((1,1),(1,0))] + my_seg = ((1,0),(0,0)) actual = winding_query.find_initial_angle(my_seg, other_segs) actual = tuple(actual) expected = winding_query.vecnorm((1,1)) @@ -62,10 +62,10 @@ def test_get_direction3(self): def test_grad_90_norm(self): segs = [((0,0),(1,0)), ((1,0),(1,1))] - pos = np.array((0,0)) + np.array((0.1,0.1)) + pos = np.array((1,1)) - np.array((0.1,0.1)) # Treats starts as repellers and ends as attractors actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((1,1)) + expected = winding_query.vecnorm((-1,-1)) self.tuples_almost_equal(actual, expected) @@ -73,24 +73,24 @@ def test_grad_90_norm2(self): segs = [((0,0),(1,0)), ((1,0),(1,1)), ((1,1),(0,1))] pos = (0,0.5) actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((0,1)) + expected = winding_query.vecnorm((0,-1)) self.tuples_almost_equal(actual, expected) pos = (-.5,0.5) actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((0,1)) + expected = winding_query.vecnorm((0,-1)) self.tuples_almost_equal(actual, expected) pos = (.5,0.5) actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((0,1)) + expected = winding_query.vecnorm((0,-1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm3(self): segs = [((0,0),(1,0)), ((1,1),(0,1))] pos = (0.00001, 0.00001) actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((1,1)) + expected = winding_query.vecnorm((-1,-1)) self.tuples_almost_equal(actual, expected) From b17444853e3ceeb4292862deca28b74d6dc82368 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Thu, 20 Jun 2024 20:37:34 -0700 Subject: [PATCH 18/25] cleanup --- data/cube_hole.stl | Bin 0 -> 584 bytes stltovoxel/playground.py | 225 ++++++++---------------------------- stltovoxel/winding_query.py | 123 ++++++++++++-------- test/test_winding_query.py | 30 ++--- 4 files changed, 141 insertions(+), 237 deletions(-) create mode 100644 data/cube_hole.stl diff --git a/data/cube_hole.stl b/data/cube_hole.stl new file mode 100644 index 0000000000000000000000000000000000000000..597c630b1c407f86c64538dae14f1e3dfcc4b97d GIT binary patch literal 584 zcmb7=Q3?Vv3`G0wDfS{FMf?y}78dW*+-$FACh@UZ#Sa6?w3C_8^t#N;dLCm;kN&)^ z_ddr(g!1n=wFrB|C2{JlkzzQ5sZup}h*ebNT=@_++Y+gDW_5YrQX@xI;RKJanD0(l x@P}F(yIO%#NW$-TZA|R!w;N4jxa>Gr9VS%t2Yvie&d7=qcd~4avZ~-f^ 0 - -def corners(segs, pt, step_towards): - x, y = pt - # Without a small movement toward the other point the value is undefined - # because we are at the center of a sweep - dx, dy = step_towards - x += dx/1000 - y += dy/1000 - ret = 0 - for i in range(len(segs) - 1): - s1 = segs[i] - s2 = segs[i+1] - assert s1[1] == s2[0] - a = np.array(list(s1[0])) - b = np.array(list(s1[1])) - c = np.array(list(s2[1])) - l1 = a - b - l2 = b - c - - ang = math.atan2(*atandiff((l2[1], l2[0]), (l1[1], l1[0]))) - - if ang < 0 : - # Concave case - if not is_right((a, b), (x,y)) and not is_right((c, b), (x,y)): - ret += ang + 2 * math.pi - else: - ret += ang - else : - # Convex case (normal) - if is_right((a, b), (x,y)) and is_right((c, b), (x,y)): - ret += ang - 2 * math.pi - else: - ret += ang - return ret - - def edge_start(me_seg, them_pt): - # Starter unpaired point s1, s2 = me_seg angle = (s1[1] - s2[1], s1[0] - s2[0]) p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) - return math.atan2(*atansum(negatan(p1), angle)) - -def edge_start_no_atan(me_seg, them_pt): - # Starter unpaired point - s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) - return atansum(negatan(p1), angle) + return math.atan2(*atan_sum(atan_neg(p1), angle)) def edge_end(me_seg, them_pt): - # Starter unpaired point s1, s2 = me_seg angle = (s1[1] - s2[1], s1[0] - s2[0]) p1 = (s2[1] - them_pt[1], s2[0] - them_pt[0]) - return math.atan2(*atansum(p1, negatan(angle))) + return math.atan2(*atan_sum(p1, atan_neg(angle))) def edge_start_baseline(me_seg): s1, s2 = me_seg @@ -97,26 +33,6 @@ def edge_end_baseline(me_seg): angle = (s1[1] - s2[1], s1[0] - s2[0]) return math.atan2(*angle) -def line(segs, inc): - im = np.zeros((100,100)) - for x in range(-50, 50): - for y in range(-50, 50): - answers = [] - for p1, p2 in segs: - x1, y1 = p1 - x2, y2 = p2 - angle = (y2 - y1, x2 - x1) - p1 = (y-y1, x-x1) - p2 = (y-y2, x-x2) - answers.append(atansum(negatan(p1), angle)) - answers.append(atansum(p2, negatan(angle))) - accum = 0 - for i in inc: - accum += math.atan2(*answers[i]) - - im[y+50][x+50] = accum - return im - def get_winding(pos, segs): accum = 0 for seg in segs: @@ -124,107 +40,56 @@ def get_winding(pos, segs): accum += edge_end(seg, pos) return accum -def get_direction(pos, segs, dangling_end): +def get_angle(pos, segs, dangling_end, target): accum = 0 for seg in segs: accum += edge_start(seg, pos) accum += edge_end(seg, pos) - # accum += edge_start(dangling_end, pos) base = edge_end_baseline(dangling_end) - return base - accum + return (base - accum) + target + math.pi -def get_direction_backwards(pos, segs, dangling_start): +def get_angle_backwards(pos, segs, dangling_start, target): accum = 0 for seg in segs: accum -= edge_start(seg, pos) accum -= edge_end(seg, pos) - # accum += edge_start(dangling_end, pos) base = edge_start_baseline(dangling_start) - return base - accum - -def get_direction_iter(pos, segs, dangling_end, target): - accum = 0 - for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - - # accum += edge_start(dangling_end, pos) - - base = edge_end_baseline(dangling_end) - # target - accum is a compensation factor to point back toward the ideal winding number - print(accum) - return (base - accum) #+ (target - accum) - + return (base - accum) + target + math.pi def angle_to_delta(theta): delta_x = math.cos(theta) delta_y = math.sin(theta) return np.array([delta_x, delta_y]) -segs = [ - ((0,10),(0,0)), - ((0,0),(10,0)), - ((10,5),(5,10)), -] - -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 - -background = np.zeros((200,200)) -for posxa in range(200): - for posya in range(200): - posx = (posxa - 100) * 0.2 - posy = (posya - 100) * 0.2 - winding = get_winding((posx, posy), segs) - if abs(winding - math.pi) < .1: - winding = math.pi*2 - background[199-posya][posxa] = winding - def grad(ox,oy,pt): px,py = pt x = ox - px y = oy - py - gx = -((y)/(x**2 + y**2)) + gx = -y/(x**2 + y**2) gy = x/(x**2 + y**2) return (gx, gy) -def accum_grad_90(x, y, segs): +def contour_left(ox, oy, pt): + gx, gy = grad(ox, oy, pt) + return (-gy, gx) + +def contour_right(ox, oy, pt): + gx, gy = grad(ox, oy, pt) + return (gy, -gx) + +def sum_contour(x, y, segs): ax = 0 ay = 0 for start, end in segs: - # for pt in [start, end]: - gx, gy = grad(x+0.00001,y+0.00001,start) - ax += gy - ay -= gx - gx, gy = grad(x+0.00001,y+0.00001,end) - ax -= gy - ay += gx + # Adding 0.00001 to avoid undefined errors due to finding contour eactly on a segment point + cx, cy = contour_right(x+0.00001,y+0.00001,start) + ax += cx + ay += cy + cx, cy = contour_left(x+0.00001,y+0.00001,end) + ax += cx + ay += cy return (ax, ay) def vecnorm(pt): @@ -232,29 +97,41 @@ def vecnorm(pt): dist = math.sqrt(x**2 + y**2) return (x / dist, y / dist) -for x in range(-10,10): - for y in range(-10,10): - ax, ay = accum_grad_90(x,y,segs) - plt.quiver(x,y, ax, ay, color=['r','b','g'], scale=21) +segs = [ + ((0,10),(0,0)), + ((0,0),(10,0)), + ((10,5),(5,10)), +] + +background = np.zeros((200,200)) +for posxa in range(200): + for posya in range(200): + posx = (posxa - 100) * 0.2 + posy = (posya - 100) * 0.2 + winding = get_winding((posx, posy), segs) + if abs(winding - math.pi) < .1: + winding = math.pi*2 + background[199-posya][posxa] = winding -# while there are some ends that need repair -# while find_polyline_endpoints(segs): -start_to_end = find_polyline_endpoints(segs) -# print(start_to_end) +for x in range(-20,20, 1): + for y in range(-20,20, 1): + ax, ay = sum_contour(x,y,segs) + if abs(ax) > 1000 or abs(ay) > 1000: + continue + plt.quiver(x,y, ax, ay, color=['r'], scale=21) -start = list(start_to_end.keys())[0] -end = start_to_end[start] +start = segs[0][0] +end = segs[0][1] # find the segment I am a part of my_seg = next(filter(lambda seg: seg[0] == start, segs)) # find all other segments other_segs = list(filter(lambda seg: seg[0] != start, segs)) -target = math.pi -angle_forward = get_direction_backwards(start, other_segs, my_seg) + target + math.pi +angle_forward = get_angle_backwards(start, other_segs, my_seg, math.pi) delta = angle_to_delta(angle_forward) pos = start + (delta * 0.1) -# plt.quiver(*start, *delta, color=['r','b','g'], scale=21) +plt.quiver(*start, *delta, color=['g'], scale=21) for i in range(100): - delt = vecnorm(accum_grad_90(*pos, segs)) + delt = vecnorm(sum_contour(*pos, segs)) plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") pos += delt diff --git a/stltovoxel/winding_query.py b/stltovoxel/winding_query.py index 2967baf..89673c4 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/winding_query.py @@ -61,12 +61,15 @@ def find_polylines(segments): # noqa: C901 return polylines -def atansum(f1, f2): - x,y = f1 - w,z = f2 - return ( x*w - y*z, y*w + x*z,) - -def negatan(f1): +def atan_sum(f1, f2): + # Angle sum and difference 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_neg(f1): + # atan2(atan_neg(f1)) == -atan2(f1) x,y = f1 return x,-y @@ -103,62 +106,86 @@ def find_polyline_endpoints(segs): return start_to_end -def winding_contour(pos,pt): - x,y = subtract(pos, pt) - dist = (x**2 + y**2) - gx = x/dist - gy = y/dist - return (gx, gy) - -def vecnorm(pt): +def winding_contour_monopole(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 dist(p1, p2): +def distance(p1, p2): x1, y1 = p1 x2, y2 = p2 return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) -def find_initial_angle(my_seg, other_segs): - start, end = my_seg - accum = subtract(start, end) - for seg in other_segs: - s1, s2 = seg - p1 = subtract(s1, end) - p2 = subtract(s2, end) - accum = atansum(accum, p1) - accum = atansum(accum, negatan(p2)) - accum = vecnorm(accum) +def initial_direction(my_seg, other_segs): + # Computing the winding number of a given point requires 2 atan2 calls per line segment (one per point). + # This method takes a line segment (my_seg) and makes 1 atan2 call for that and 2 atan2 calls for all other line segments (other_segs). + # This produces an angle which if followed out of the end of my_seg would produce a winding number of the target value. + # In theory this target winding number can be any value, but here pi radians / 180 degrees is used for simplicity. + pt = my_seg[1] + accum = subtract(my_seg[0], my_seg[1]) + for seg_start, seg_end in other_segs: + accum = atan_sum(accum, subtract(seg_start, pt)) + accum = atan_sum(accum, atan_neg(subtract(seg_end, pt))) + # Without this accum can get arbitrarily large which causes problems with lots of segments + accum = normalize(accum) return np.array(accum) -def total_winding_contour(pos, segs): +def winding_contour(pos, segs): accum = (0,0) for start, end in segs: - start_grad = winding_contour(pos, start) - end_grad = winding_contour(pos, end) - accum = subtract(accum, start_grad) - accum = add(accum, end_grad) - return vecnorm(accum) - -def find_flow(start, ends, segs): + # Starting points attract the vector field + start_vec = winding_contour_monopole(pos, start, repel=False) + accum = add(accum, start_vec) + # Ending points repel the vector field + end_vec = winding_contour_monopole(pos, end, repel=True) + accum = add(accum, end_vec) + return normalize(accum) + +def find_flow(start, ends, segs, polyline_endpoints): # find the segment I am a part of my_seg = next(filter(lambda seg: seg[1] == start, segs)) # find all other segments other_segs = list(filter(lambda seg: seg[1] != start, segs)) - delta = find_initial_angle(my_seg, other_segs) - pos = start + (delta * 0.1) + # Find the initial direction to start marching towards. + direction = initial_direction(my_seg, other_segs) + # Move slightly toward that direction to pick the contour that we will use below + pos = start + (direction * 0.001) seg_outs = [(tuple(start), tuple(pos))] for _ in range(200): - delt = total_winding_contour(pos, segs) - seg_outs.append((tuple(pos), tuple(pos + delt))) - pos += delt + # 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 each line segment pair + # which has a start and end at the same point will cancel out. + direction = winding_contour(pos, polyline_endpoints) + # March along the flowline using euler's method 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 dist(pos, end) < 1: + if distance(pos, end) < 1: seg_outs.append((tuple(pos), tuple(end))) return seg_outs - raise "Flow not found" + raise Exception("Max iteration number exceeded to repair mesh") class WindingQuery(): def __init__(self, segments): @@ -184,16 +211,16 @@ def repair_all(self): old_seg_length = len(self.polylines) self.collapse_segments() assert old_seg_length - 1 == len(self.polylines) - for seg in self.original_segments: - plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") - plt.show() + # for seg in self.original_segments: + # plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") + # plt.show() assert len(self.polylines) == 0 def repair_segment(self): - # Search starts at the end of a polyline - start = self.polylines[0][-1] - + # 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) - endpoints = [polyline[0] for polyline in self.polylines] - new_segs = find_flow(start, endpoints, self.original_segments) + search_ends = [polyline[0] for polyline in self.polylines] + polyline_endpoints = [(polyline[0], polyline[-1]) for polyline in self.polylines] + new_segs = find_flow(search_start, search_ends, self.original_segments, polyline_endpoints) self.original_segments.extend(new_segs) diff --git a/test/test_winding_query.py b/test/test_winding_query.py index ac1278d..4831df7 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -39,7 +39,7 @@ def test_find_polylines_out_of_order(self): def test_get_direction(self): other_segs = [((1,0),(1,1)),((1,1),(0,1)),] my_seg = ((0,1),(0,0)) - actual = winding_query.find_initial_angle(my_seg, other_segs) + actual = winding_query.initial_direction(my_seg, other_segs) actual = tuple(actual) expected = (1.0,0) self.tuples_almost_equal(actual, expected) @@ -47,50 +47,50 @@ def test_get_direction(self): def test_get_direction2(self): other_segs = [((1,0),(0,0))] my_seg = ((0,1),(1,1)) - actual = winding_query.find_initial_angle(my_seg, other_segs) + actual = winding_query.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = winding_query.vecnorm((-1,-1)) + expected = winding_query.normalize((-1,-1)) self.tuples_almost_equal(actual, expected) def test_get_direction3(self): other_segs = [((1,1),(1,0))] my_seg = ((1,0),(0,0)) - actual = winding_query.find_initial_angle(my_seg, other_segs) + actual = winding_query.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = winding_query.vecnorm((1,1)) + expected = winding_query.normalize((1,1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm(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 = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((-1,-1)) + actual = winding_query.winding_contour(pos, segs) + expected = winding_query.normalize((-1,-1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm2(self): segs = [((0,0),(1,0)), ((1,0),(1,1)), ((1,1),(0,1))] pos = (0,0.5) - actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((0,-1)) + actual = winding_query.winding_contour(pos, segs) + expected = winding_query.normalize((0,-1)) self.tuples_almost_equal(actual, expected) pos = (-.5,0.5) - actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((0,-1)) + actual = winding_query.winding_contour(pos, segs) + expected = winding_query.normalize((0,-1)) self.tuples_almost_equal(actual, expected) pos = (.5,0.5) - actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((0,-1)) + actual = winding_query.winding_contour(pos, segs) + expected = winding_query.normalize((0,-1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm3(self): segs = [((0,0),(1,0)), ((1,1),(0,1))] pos = (0.00001, 0.00001) - actual = winding_query.total_winding_contour(pos, segs) - expected = winding_query.vecnorm((-1,-1)) + actual = winding_query.winding_contour(pos, segs) + expected = winding_query.normalize((-1,-1)) self.tuples_almost_equal(actual, expected) From 9432aa61721eb921cce9984f48f8650cea4f04e1 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sat, 22 Jun 2024 16:38:16 -0700 Subject: [PATCH 19/25] stash --- .gitignore | 5 +- stltovoxel/perimeter.py | 5 +- stltovoxel/playground.py | 33 +++----- .../{winding_query.py => polygon_repair.py} | 78 +++++++++---------- test/test_winding_query.py | 38 ++++----- 5 files changed, 75 insertions(+), 84 deletions(-) rename stltovoxel/{winding_query.py => polygon_repair.py} (83%) 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/stltovoxel/perimeter.py b/stltovoxel/perimeter.py index 6fd97c3..de5e9cb 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): 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) wq.repair_all() new_line_list = [] for polyline in wq.loops: diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py index 0a40d56..2c6fcfc 100644 --- a/stltovoxel/playground.py +++ b/stltovoxel/playground.py @@ -28,11 +28,6 @@ def edge_start_baseline(me_seg): angle = (s2[1] - s1[1], s2[0] - s1[0]) return math.atan2(*angle) -def edge_end_baseline(me_seg): - s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - return math.atan2(*angle) - def get_winding(pos, segs): accum = 0 for seg in segs: @@ -40,23 +35,14 @@ def get_winding(pos, segs): accum += edge_end(seg, pos) return accum -def get_angle(pos, segs, dangling_end, target): - accum = 0 - for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - - base = edge_end_baseline(dangling_end) - return (base - accum) + target + math.pi - -def get_angle_backwards(pos, segs, dangling_start, target): +def get_angle(pos, segs, dangling_start, target): accum = 0 for seg in segs: accum -= edge_start(seg, pos) accum -= edge_end(seg, pos) base = edge_start_baseline(dangling_start) - return (base - accum) + target + math.pi + return base - accum + math.pi - target def angle_to_delta(theta): delta_x = math.cos(theta) @@ -98,8 +84,8 @@ def vecnorm(pt): return (x / dist, y / dist) segs = [ - ((0,10),(0,0)), - ((0,0),(10,0)), + ((-10,15),(-10,-15)), + ((-10,-15),(10,-15)), ((10,5),(5,10)), ] @@ -113,12 +99,13 @@ def vecnorm(pt): winding = math.pi*2 background[199-posya][posxa] = winding -for x in range(-20,20, 1): - for y in range(-20,20, 1): +for x in range(-20,20): + for y in range(-20,20): ax, ay = sum_contour(x,y,segs) if abs(ax) > 1000 or abs(ay) > 1000: + # Too big, likely near a segment point continue - plt.quiver(x,y, ax, ay, color=['r'], scale=21) + plt.quiver(x,y, ax, ay, color=['r'], scale=20) start = segs[0][0] end = segs[0][1] @@ -126,10 +113,10 @@ def vecnorm(pt): my_seg = next(filter(lambda seg: seg[0] == start, segs)) # find all other segments other_segs = list(filter(lambda seg: seg[0] != start, segs)) -angle_forward = get_angle_backwards(start, other_segs, my_seg, math.pi) +angle_forward = get_angle(start, other_segs, my_seg, math.pi) delta = angle_to_delta(angle_forward) pos = start + (delta * 0.1) -plt.quiver(*start, *delta, color=['g'], scale=21) +plt.quiver(*start, *(delta * 10), color=['g'], scale=40) for i in range(100): delt = vecnorm(sum_contour(*pos, segs)) plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") diff --git a/stltovoxel/winding_query.py b/stltovoxel/polygon_repair.py similarity index 83% rename from stltovoxel/winding_query.py rename to stltovoxel/polygon_repair.py index 89673c4..e3019eb 100644 --- a/stltovoxel/winding_query.py +++ b/stltovoxel/polygon_repair.py @@ -1,5 +1,4 @@ import numpy as np -import matplotlib.pyplot as plt import math def find_polylines(segments): # noqa: C901 @@ -61,24 +60,6 @@ def find_polylines(segments): # noqa: C901 return polylines -def atan_sum(f1, f2): - # Angle sum and difference 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_neg(f1): - # atan2(atan_neg(f1)) == -atan2(f1) - x,y = f1 - return x,-y - -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 find_polyline_endpoints(segs): start_to_end = dict() end_to_start = dict() @@ -106,14 +87,32 @@ def find_polyline_endpoints(segs): return start_to_end -def winding_contour_monopole(pos, pt, repel): +def atan_sum(f1, f2): + # Angle sum and difference 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_neg(f1): + # atan2(atan_neg(f1)) == -atan2(f1) + x,y = f1 + return x,-y + +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) 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) @@ -136,16 +135,20 @@ def distance(p1, p2): return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) def initial_direction(my_seg, other_segs): - # Computing the winding number of a given point requires 2 atan2 calls per line segment (one per point). - # This method takes a line segment (my_seg) and makes 1 atan2 call for that and 2 atan2 calls for all other line segments (other_segs). - # This produces an angle which if followed out of the end of my_seg would produce a winding number of the target value. - # In theory this target winding number can be any value, but here pi radians / 180 degrees is used for simplicity. + # 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. pt = my_seg[1] accum = subtract(my_seg[0], my_seg[1]) for seg_start, seg_end in other_segs: accum = atan_sum(accum, subtract(seg_start, pt)) accum = atan_sum(accum, atan_neg(subtract(seg_end, pt))) - # Without this accum can get arbitrarily large which causes problems with lots of segments + # Without this accum can get arbitrarily large and lead to floating point problems accum = normalize(accum) return np.array(accum) @@ -153,14 +156,14 @@ def winding_contour(pos, segs): accum = (0,0) for start, end in segs: # Starting points attract the vector field - start_vec = winding_contour_monopole(pos, start, repel=False) + start_vec = winding_contour_pole(pos, start, repel=False) accum = add(accum, start_vec) # Ending points repel the vector field - end_vec = winding_contour_monopole(pos, end, repel=True) + end_vec = winding_contour_pole(pos, end, repel=True) accum = add(accum, end_vec) return normalize(accum) -def find_flow(start, ends, segs, polyline_endpoints): +def winding_number_search(start, ends, segs, polyline_endpoints): # find the segment I am a part of my_seg = next(filter(lambda seg: seg[1] == start, segs)) # find all other segments @@ -172,10 +175,10 @@ def find_flow(start, ends, segs, polyline_endpoints): seg_outs = [(tuple(start), tuple(pos))] for _ in range(200): # 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 each line segment pair - # which has a start and end at the same point will cancel out. + # 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 using euler's method to find where it terminates. + # March along the flowline to find where it terminates. new_pos = pos + direction seg_outs.append((tuple(pos), tuple(new_pos))) pos = new_pos @@ -187,7 +190,7 @@ def find_flow(start, ends, segs, polyline_endpoints): raise Exception("Max iteration number exceeded to repair mesh") -class WindingQuery(): +class PolygonRepair(): def __init__(self, segments): # Maps endpoints to the polygon they form self.loops = [] @@ -207,20 +210,17 @@ def collapse_segments(self): def repair_all(self): while self.polylines: - self.repair_segment() + self.repair_polyline() old_seg_length = len(self.polylines) self.collapse_segments() assert old_seg_length - 1 == len(self.polylines) - # for seg in self.original_segments: - # plt.plot([seg[0][0], seg[1][0]], [seg[0][1], seg[1][1]], 'bo', linestyle="-") - # plt.show() assert len(self.polylines) == 0 - def repair_segment(self): + 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] - new_segs = find_flow(search_start, search_ends, self.original_segments, polyline_endpoints) + new_segs = winding_number_search(search_start, search_ends, self.original_segments, polyline_endpoints) self.original_segments.extend(new_segs) diff --git a/test/test_winding_query.py b/test/test_winding_query.py index 4831df7..131ce31 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -2,7 +2,7 @@ import math import numpy as np -from stltovoxel import winding_query +from stltovoxel import polygon_repair class TestWindingQuery(unittest.TestCase): @@ -13,7 +13,7 @@ def tuples_almost_equal(self, 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 = winding_query.find_polylines(line_segments) + actual_polylines = polygon_repair.find_polylines(line_segments) self.assertEqual(actual_polylines, expected_polylines) def test_find_polylines_no_cycle(self): @@ -27,19 +27,19 @@ def test_find_polylines_no_cycle(self): ((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) + 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 = winding_query.find_polylines(line_segments) + actual_polylines = polygon_repair.find_polylines(line_segments) self.assertEqual(actual_polylines, expected_polylines) def test_get_direction(self): other_segs = [((1,0),(1,1)),((1,1),(0,1)),] my_seg = ((0,1),(0,0)) - actual = winding_query.initial_direction(my_seg, other_segs) + actual = polygon_repair.initial_direction(my_seg, other_segs) actual = tuple(actual) expected = (1.0,0) self.tuples_almost_equal(actual, expected) @@ -47,50 +47,50 @@ def test_get_direction(self): def test_get_direction2(self): other_segs = [((1,0),(0,0))] my_seg = ((0,1),(1,1)) - actual = winding_query.initial_direction(my_seg, other_segs) + actual = polygon_repair.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = winding_query.normalize((-1,-1)) + expected = polygon_repair.normalize((-1,-1)) self.tuples_almost_equal(actual, expected) def test_get_direction3(self): other_segs = [((1,1),(1,0))] my_seg = ((1,0),(0,0)) - actual = winding_query.initial_direction(my_seg, other_segs) + actual = polygon_repair.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = winding_query.normalize((1,1)) + expected = polygon_repair.normalize((1,1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm(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 = winding_query.winding_contour(pos, segs) - expected = winding_query.normalize((-1,-1)) + actual = polygon_repair.winding_contour(pos, segs) + expected = polygon_repair.normalize((-1,-1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm2(self): segs = [((0,0),(1,0)), ((1,0),(1,1)), ((1,1),(0,1))] pos = (0,0.5) - actual = winding_query.winding_contour(pos, segs) - expected = winding_query.normalize((0,-1)) + actual = polygon_repair.winding_contour(pos, segs) + expected = polygon_repair.normalize((0,-1)) self.tuples_almost_equal(actual, expected) pos = (-.5,0.5) - actual = winding_query.winding_contour(pos, segs) - expected = winding_query.normalize((0,-1)) + actual = polygon_repair.winding_contour(pos, segs) + expected = polygon_repair.normalize((0,-1)) self.tuples_almost_equal(actual, expected) pos = (.5,0.5) - actual = winding_query.winding_contour(pos, segs) - expected = winding_query.normalize((0,-1)) + actual = polygon_repair.winding_contour(pos, segs) + expected = polygon_repair.normalize((0,-1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm3(self): segs = [((0,0),(1,0)), ((1,1),(0,1))] pos = (0.00001, 0.00001) - actual = winding_query.winding_contour(pos, segs) - expected = winding_query.normalize((-1,-1)) + actual = polygon_repair.winding_contour(pos, segs) + expected = polygon_repair.normalize((-1,-1)) self.tuples_almost_equal(actual, expected) From 2999e3131d8f0ba627541c1533f8e47cf36e6afb Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sat, 22 Jun 2024 16:38:24 -0700 Subject: [PATCH 20/25] remove playground --- stltovoxel/playground.py | 128 --------------------------------------- 1 file changed, 128 deletions(-) delete mode 100644 stltovoxel/playground.py diff --git a/stltovoxel/playground.py b/stltovoxel/playground.py deleted file mode 100644 index 2c6fcfc..0000000 --- a/stltovoxel/playground.py +++ /dev/null @@ -1,128 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -import math - -def atan_sum(f1, f2): - y, x = f1 - z, w = f2 - return (y*w + x*z, x*w - y*z) - -def atan_neg(f1): - y,x = f1 - return -y, x - -def edge_start(me_seg, them_pt): - s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (s1[1] - them_pt[1], s1[0] - them_pt[0]) - return math.atan2(*atan_sum(atan_neg(p1), angle)) - -def edge_end(me_seg, them_pt): - s1, s2 = me_seg - angle = (s1[1] - s2[1], s1[0] - s2[0]) - p1 = (s2[1] - them_pt[1], s2[0] - them_pt[0]) - return math.atan2(*atan_sum(p1, atan_neg(angle))) - -def edge_start_baseline(me_seg): - s1, s2 = me_seg - angle = (s2[1] - s1[1], s2[0] - s1[0]) - return math.atan2(*angle) - -def get_winding(pos, segs): - accum = 0 - for seg in segs: - accum += edge_start(seg, pos) - accum += edge_end(seg, pos) - return accum - -def get_angle(pos, segs, dangling_start, target): - accum = 0 - for seg in segs: - accum -= edge_start(seg, pos) - accum -= edge_end(seg, pos) - - base = edge_start_baseline(dangling_start) - return base - accum + math.pi - target - -def angle_to_delta(theta): - delta_x = math.cos(theta) - delta_y = math.sin(theta) - return np.array([delta_x, delta_y]) - -def grad(ox,oy,pt): - px,py = pt - x = ox - px - y = oy - py - gx = -y/(x**2 + y**2) - gy = x/(x**2 + y**2) - return (gx, gy) - -def contour_left(ox, oy, pt): - gx, gy = grad(ox, oy, pt) - return (-gy, gx) - -def contour_right(ox, oy, pt): - gx, gy = grad(ox, oy, pt) - return (gy, -gx) - -def sum_contour(x, y, segs): - ax = 0 - ay = 0 - for start, end in segs: - # Adding 0.00001 to avoid undefined errors due to finding contour eactly on a segment point - cx, cy = contour_right(x+0.00001,y+0.00001,start) - ax += cx - ay += cy - cx, cy = contour_left(x+0.00001,y+0.00001,end) - ax += cx - ay += cy - return (ax, ay) - -def vecnorm(pt): - x, y = pt - dist = math.sqrt(x**2 + y**2) - return (x / dist, y / dist) - -segs = [ - ((-10,15),(-10,-15)), - ((-10,-15),(10,-15)), - ((10,5),(5,10)), -] - -background = np.zeros((200,200)) -for posxa in range(200): - for posya in range(200): - posx = (posxa - 100) * 0.2 - posy = (posya - 100) * 0.2 - winding = get_winding((posx, posy), segs) - if abs(winding - math.pi) < .1: - winding = math.pi*2 - background[199-posya][posxa] = winding - -for x in range(-20,20): - for y in range(-20,20): - ax, ay = sum_contour(x,y,segs) - if abs(ax) > 1000 or abs(ay) > 1000: - # Too big, likely near a segment point - continue - plt.quiver(x,y, ax, ay, color=['r'], scale=20) - -start = segs[0][0] -end = segs[0][1] -# find the segment I am a part of -my_seg = next(filter(lambda seg: seg[0] == start, segs)) -# find all other segments -other_segs = list(filter(lambda seg: seg[0] != start, segs)) -angle_forward = get_angle(start, other_segs, my_seg, math.pi) -delta = angle_to_delta(angle_forward) -pos = start + (delta * 0.1) -plt.quiver(*start, *(delta * 10), color=['g'], scale=40) -for i in range(100): - delt = vecnorm(sum_contour(*pos, segs)) - plt.plot([pos[0], pos[0] + delt[0]], [pos[1], pos[1] + delt[1]], 'bo', linestyle="-") - pos += delt - - -plt.imshow(background, aspect='auto', cmap='viridis', extent=[-20, 20, -20, 20], vmin=-math.pi, vmax=math.pi*2) - -plt.show() \ No newline at end of file From bb7a70a97652156cd256b314cad8f92f5fae7f13 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sat, 22 Jun 2024 16:54:51 -0700 Subject: [PATCH 21/25] fixes --- stltovoxel/perimeter.py | 4 ++-- stltovoxel/polygon_repair.py | 18 ++++++++++-------- stltovoxel/slice.py | 2 +- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/stltovoxel/perimeter.py b/stltovoxel/perimeter.py index de5e9cb..497034d 100644 --- a/stltovoxel/perimeter.py +++ b/stltovoxel/perimeter.py @@ -1,11 +1,11 @@ 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 segments = [[tuple(pt.tolist())[:2] for pt in seg] for seg in line_list] - wq = polygon_repair.PolygonRepair(segments) + 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 index e3019eb..25636be 100644 --- a/stltovoxel/polygon_repair.py +++ b/stltovoxel/polygon_repair.py @@ -163,7 +163,7 @@ def winding_contour(pos, segs): accum = add(accum, end_vec) return normalize(accum) -def winding_number_search(start, ends, segs, polyline_endpoints): +def winding_number_search(start, ends, segs, polyline_endpoints, max_iterations): # find the segment I am a part of my_seg = next(filter(lambda seg: seg[1] == start, segs)) # find all other segments @@ -171,9 +171,9 @@ def winding_number_search(start, ends, segs, polyline_endpoints): # Find the initial direction to start marching towards. direction = initial_direction(my_seg, other_segs) # Move slightly toward that direction to pick the contour that we will use below - pos = start + (direction * 0.001) + pos = start + (direction * 0.1) seg_outs = [(tuple(start), tuple(pos))] - for _ in range(200): + 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. @@ -186,17 +186,18 @@ def winding_number_search(start, ends, segs, polyline_endpoints): for end in ends: if distance(pos, end) < 1: seg_outs.append((tuple(pos), tuple(end))) - return seg_outs + return end - raise Exception("Max iteration number exceeded to repair mesh") + raise Exception("Failed to repair mesh") class PolygonRepair(): - def __init__(self, segments): + 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): @@ -222,5 +223,6 @@ def repair_polyline(self): # 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] - new_segs = winding_number_search(search_start, search_ends, self.original_segments, polyline_endpoints) - self.original_segments.extend(new_segs) + 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 From 507d2c6c750b0c0b6784c066cad93143b88926ba Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Sat, 22 Jun 2024 16:57:32 -0700 Subject: [PATCH 22/25] lint --- stltovoxel/polygon_repair.py | 30 ++++++++++++++++++-------- test/test_winding_query.py | 42 ++++++++++++++++++------------------ 2 files changed, 42 insertions(+), 30 deletions(-) diff --git a/stltovoxel/polygon_repair.py b/stltovoxel/polygon_repair.py index 25636be..8eb9cf1 100644 --- a/stltovoxel/polygon_repair.py +++ b/stltovoxel/polygon_repair.py @@ -60,6 +60,7 @@ def find_polylines(segments): # noqa: C901 return polylines + def find_polyline_endpoints(segs): start_to_end = dict() end_to_start = dict() @@ -87,30 +88,35 @@ def find_polyline_endpoints(segs): return start_to_end + def atan_sum(f1, f2): # Angle sum and difference identity # atan2(atan_sum(f1, f2)) == atan2(f1) + atan2(f2) - x1,y1 = f1 - x2,y2 = f2 + x1, y1 = f1 + x2, y2 = f2 return ( x1*x2 - y1*y2, y1*x2 + x1*y2) + def atan_neg(f1): # atan2(atan_neg(f1)) == -atan2(f1) - x,y = f1 - return x,-y + x, y = f1 + return x, -y + 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 + # 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) + # (-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) @@ -124,24 +130,27 @@ def winding_contour_pole(pos, pt, repel): 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(my_seg, other_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. + # 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. + # can be avoided through use of atan2 identities. # Since the final angle would be converted back into a vector, no atan2 call is required. pt = my_seg[1] accum = subtract(my_seg[0], my_seg[1]) @@ -152,8 +161,9 @@ def initial_direction(my_seg, other_segs): accum = normalize(accum) return np.array(accum) + def winding_contour(pos, segs): - accum = (0,0) + accum = (0, 0) for start, end in segs: # Starting points attract the vector field start_vec = winding_contour_pole(pos, start, repel=False) @@ -163,6 +173,7 @@ def winding_contour(pos, segs): accum = add(accum, end_vec) return normalize(accum) + def winding_number_search(start, ends, segs, polyline_endpoints, max_iterations): # find the segment I am a part of my_seg = next(filter(lambda seg: seg[1] == start, segs)) @@ -190,6 +201,7 @@ def winding_number_search(start, ends, segs, polyline_endpoints, max_iterations) raise Exception("Failed to repair mesh") + class PolygonRepair(): def __init__(self, segments, dimensions): # Maps endpoints to the polygon they form diff --git a/test/test_winding_query.py b/test/test_winding_query.py index 131ce31..9f03956 100644 --- a/test/test_winding_query.py +++ b/test/test_winding_query.py @@ -37,60 +37,60 @@ def test_find_polylines_out_of_order(self): self.assertEqual(actual_polylines, expected_polylines) def test_get_direction(self): - other_segs = [((1,0),(1,1)),((1,1),(0,1)),] - my_seg = ((0,1),(0,0)) + other_segs = [((1, 0), (1, 1)), ((1, 1), (0, 1)), ] + my_seg = ((0, 1), (0, 0)) actual = polygon_repair.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = (1.0,0) + expected = (1.0, 0) self.tuples_almost_equal(actual, expected) def test_get_direction2(self): - other_segs = [((1,0),(0,0))] - my_seg = ((0,1),(1,1)) + other_segs = [((1, 0), (0, 0))] + my_seg = ((0, 1), (1, 1)) actual = polygon_repair.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = polygon_repair.normalize((-1,-1)) + expected = polygon_repair.normalize((-1, -1)) self.tuples_almost_equal(actual, expected) def test_get_direction3(self): - other_segs = [((1,1),(1,0))] - my_seg = ((1,0),(0,0)) + other_segs = [((1, 1), (1, 0))] + my_seg = ((1, 0), (0, 0)) actual = polygon_repair.initial_direction(my_seg, other_segs) actual = tuple(actual) - expected = polygon_repair.normalize((1,1)) + expected = polygon_repair.normalize((1, 1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm(self): - segs = [((0,0),(1,0)), ((1,0),(1,1))] - pos = np.array((1,1)) - np.array((0.1,0.1)) + 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)) + expected = polygon_repair.normalize((-1, -1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm2(self): - segs = [((0,0),(1,0)), ((1,0),(1,1)), ((1,1),(0,1))] - pos = (0,0.5) + 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)) + expected = polygon_repair.normalize((0, -1)) self.tuples_almost_equal(actual, expected) - pos = (-.5,0.5) + pos = (-.5, 0.5) actual = polygon_repair.winding_contour(pos, segs) - expected = polygon_repair.normalize((0,-1)) + expected = polygon_repair.normalize((0, -1)) self.tuples_almost_equal(actual, expected) - pos = (.5,0.5) + pos = (.5, 0.5) actual = polygon_repair.winding_contour(pos, segs) - expected = polygon_repair.normalize((0,-1)) + expected = polygon_repair.normalize((0, -1)) self.tuples_almost_equal(actual, expected) def test_grad_90_norm3(self): - segs = [((0,0),(1,0)), ((1,1),(0,1))] + 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)) + expected = polygon_repair.normalize((-1, -1)) self.tuples_almost_equal(actual, expected) From 0f387fd5ced533c7a8a33e8b9c1ec1b04fd5a923 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Tue, 25 Jun 2024 18:17:42 -0700 Subject: [PATCH 23/25] lint --- stltovoxel/polygon_repair.py | 3 ++- test/{test_winding_query.py => test_polygon_repair.py} | 8 +++----- 2 files changed, 5 insertions(+), 6 deletions(-) rename test/{test_winding_query.py => test_polygon_repair.py} (98%) diff --git a/stltovoxel/polygon_repair.py b/stltovoxel/polygon_repair.py index 8eb9cf1..92c3784 100644 --- a/stltovoxel/polygon_repair.py +++ b/stltovoxel/polygon_repair.py @@ -1,6 +1,7 @@ import numpy as np import math + def find_polylines(segments): # noqa: C901 polylines = [] segment_forward_dict = {} @@ -94,7 +95,7 @@ def atan_sum(f1, f2): # atan2(atan_sum(f1, f2)) == atan2(f1) + atan2(f2) x1, y1 = f1 x2, y2 = f2 - return ( x1*x2 - y1*y2, y1*x2 + x1*y2) + return (x1*x2 - y1*y2, y1*x2 + x1*y2) def atan_neg(f1): diff --git a/test/test_winding_query.py b/test/test_polygon_repair.py similarity index 98% rename from test/test_winding_query.py rename to test/test_polygon_repair.py index 9f03956..f52f76f 100644 --- a/test/test_winding_query.py +++ b/test/test_polygon_repair.py @@ -1,11 +1,10 @@ -import unittest -import math import numpy as np +import unittest from stltovoxel import polygon_repair -class TestWindingQuery(unittest.TestCase): +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}") @@ -68,7 +67,6 @@ def test_grad_90_norm(self): expected = polygon_repair.normalize((-1, -1)) self.tuples_almost_equal(actual, expected) - def test_grad_90_norm2(self): segs = [((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (0, 1))] pos = (0, 0.5) @@ -92,7 +90,7 @@ def test_grad_90_norm3(self): actual = polygon_repair.winding_contour(pos, segs) expected = polygon_repair.normalize((-1, -1)) self.tuples_almost_equal(actual, expected) - + if __name__ == '__main__': unittest.main() From 19364ac8abba3b70449f8124e2ab1e23e3f43389 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Tue, 25 Jun 2024 20:33:24 -0700 Subject: [PATCH 24/25] fixes --- stltovoxel/polygon_repair.py | 35 +++++++++++++++------------- test/test_polygon_repair.py | 44 ++++++++++++++++++++++++------------ 2 files changed, 48 insertions(+), 31 deletions(-) diff --git a/stltovoxel/polygon_repair.py b/stltovoxel/polygon_repair.py index 92c3784..a6fab8f 100644 --- a/stltovoxel/polygon_repair.py +++ b/stltovoxel/polygon_repair.py @@ -91,17 +91,20 @@ def find_polyline_endpoints(segs): def atan_sum(f1, f2): - # Angle sum and difference identity + # 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_neg(f1): - # atan2(atan_neg(f1)) == -atan2(f1) - x, y = f1 - return x, -y +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): @@ -144,7 +147,7 @@ def distance(p1, p2): return math.sqrt((y2 - y1)**2 + (x2 - x1)**2) -def initial_direction(my_seg, other_segs): +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]. @@ -153,11 +156,15 @@ def initial_direction(my_seg, other_segs): # 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. - pt = my_seg[1] - accum = subtract(my_seg[0], my_seg[1]) - for seg_start, seg_end in other_segs: - accum = atan_sum(accum, subtract(seg_start, pt)) - accum = atan_sum(accum, atan_neg(subtract(seg_end, pt))) + 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) @@ -176,12 +183,8 @@ def winding_contour(pos, segs): def winding_number_search(start, ends, segs, polyline_endpoints, max_iterations): - # find the segment I am a part of - my_seg = next(filter(lambda seg: seg[1] == start, segs)) - # find all other segments - other_segs = list(filter(lambda seg: seg[1] != start, segs)) # Find the initial direction to start marching towards. - direction = initial_direction(my_seg, other_segs) + 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))] diff --git a/test/test_polygon_repair.py b/test/test_polygon_repair.py index f52f76f..5c149bc 100644 --- a/test/test_polygon_repair.py +++ b/test/test_polygon_repair.py @@ -35,31 +35,31 @@ def test_find_polylines_out_of_order(self): actual_polylines = polygon_repair.find_polylines(line_segments) self.assertEqual(actual_polylines, expected_polylines) - def test_get_direction(self): - other_segs = [((1, 0), (1, 1)), ((1, 1), (0, 1)), ] - my_seg = ((0, 1), (0, 0)) - actual = polygon_repair.initial_direction(my_seg, other_segs) + 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_get_direction2(self): - other_segs = [((1, 0), (0, 0))] - my_seg = ((0, 1), (1, 1)) - actual = polygon_repair.initial_direction(my_seg, other_segs) + 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_get_direction3(self): - other_segs = [((1, 1), (1, 0))] - my_seg = ((1, 0), (0, 0)) - actual = polygon_repair.initial_direction(my_seg, other_segs) + 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_grad_90_norm(self): + 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 @@ -67,7 +67,7 @@ def test_grad_90_norm(self): expected = polygon_repair.normalize((-1, -1)) self.tuples_almost_equal(actual, expected) - def test_grad_90_norm2(self): + 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) @@ -84,13 +84,27 @@ def test_grad_90_norm2(self): expected = polygon_repair.normalize((0, -1)) self.tuples_almost_equal(actual, expected) - def test_grad_90_norm3(self): + 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() From 94d2b4407ceb6e352d6ceec6ca1f3010e0bd9df3 Mon Sep 17 00:00:00 2001 From: Christian Pederkoff <6548797+cpederkoff@users.noreply.github.com> Date: Fri, 28 Jun 2024 23:32:16 -0700 Subject: [PATCH 25/25] stash --- stltovoxel/polygon_repair.py | 1 - 1 file changed, 1 deletion(-) diff --git a/stltovoxel/polygon_repair.py b/stltovoxel/polygon_repair.py index a6fab8f..05598bc 100644 --- a/stltovoxel/polygon_repair.py +++ b/stltovoxel/polygon_repair.py @@ -106,7 +106,6 @@ def atan_diff(f1, f2): return (x1*x2 + y1*y2, y1*x2 - x1*y2) - def subtract(s1, s2): return (s1[0] - s2[0], s1[1] - s2[1])