Walk along path of discrete line segments evenly distributing points

Question

I am trying to write a program that given a list of points indicating a path and given a desired number of marks, it should distribute these marks exactly evenly along the path. As it happens the path is cyclical but given an arbitrary point to both start and end at I don't think it affects the algorithm at all.

The first step is to sum up the length of the line segments to determine the total length of the path, and then dividing that by the number of marks to get the desired distance between marks. Easy enough.

The next step is to walk along the path, storing the coordinates of each mark each time you traverse another even multiple's worth of the distance between marks.

In my code, the traversal seems correct but the distribution of marks is not even and does not exactly follow the path. I have created a visualization in matplotlib to plot where the marks are landing showing this (see last section).

path data

point_data = [
 (53.8024, 50.4762), (49.5272, 51.8727), (45.0118, 52.3863), (40.5399, 53.0184), (36.3951, 54.7708),
 (28.7127, 58.6807), (25.5306, 61.4955), (23.3828, 65.2082), (22.6764, 68.3316), (22.6945, 71.535),
 (24.6674, 77.6427), (28.8279, 82.4529), (31.5805, 84.0346), (34.7024, 84.8875), (45.9183, 84.5739),
 (57.0529, 82.9846), (64.2141, 79.1657), (71.089, 74.802), (76.7944, 69.8429), (82.1092, 64.4783),
 (83.974, 63.3605), (85.2997, 61.5455), (85.7719, 59.4206), (85.0764, 57.3729), (82.0979, 56.0247),
 (78.878, 55.1062), (73.891, 53.0987), (68.7101, 51.7283), (63.6943, 51.2997), (58.6791, 51.7438),
 (56.1255, 51.5243), (53.8024, 50.4762), (53.8024, 50.4762)]

traversal

import math

number_of_points = 20

def euclid_dist(x1, y1, x2, y2):
  return ((x1-x2)**2 + (y1-y2)**2)**0.5

def move_point(x0, y0, d, theta_rad):
  return x0 + d*math.cos(theta_rad), y0 + d*math.sin(theta_rad)

total_dist = 0
for i in range(1, len(point_data), 1):
  x1, y1 = point_data[i - 1]
  x2, y2 = point_data[i]
  total_dist += euclid_dist(x1, y1, x2, y2)

dist_per_point = total_dist / number_of_points

length_left_over = 0  # distance left over from the last segment
# led_id = 0

results = []

for i in range(1, len(point_data), 1):
  x1, y1 = point_data[i - 1]
  x2, y2 = point_data[i]

  angle_rads  = math.atan2(y1-y2, x1-x2)
  extra_rotation = math.pi / 2  # 90deg
  angle_output = math.degrees((angle_rads + extra_rotation + math.pi) % (2*math.pi) - math.pi)
  length_of_segment = euclid_dist(x1, y1, x2, y2)
    
  distance_to_work_with = length_left_over + length_of_segment
  
  current_dist = dist_per_point - length_left_over

  while distance_to_work_with > dist_per_point:

    new_point = move_point(x1, y1, current_dist, angle_rads)
    results.append((new_point[0], new_point[1], angle_output))

    current_dist += dist_per_point
    
    distance_to_work_with -= dist_per_point

  length_left_over = distance_to_work_with

visualization code

import matplotlib.pyplot as plt
from matplotlib import collections  as mc
import numpy as np

X = np.array([x for x, _, _ in results])
Y = np.array([y for _, y, _ in results])

plt.scatter(X, Y)

for i, (x, y) in enumerate(zip(X, Y)):
    plt.text(x, y, str(i), color="red", fontsize=12)

possible_colors = [(1, 0, 0, 1), (0, 1, 0, 1), (0, 0, 1, 1)]

lines = []
colors = []
for i in range(len(point_data) -1 , 0, -1):
  x1, y1 = point_data[i - 1]
  x2, y2 = point_data[i]
  lines.append(((x1, y1), (x2, y2)))
  colors.append(possible_colors[i % 3])

lc = mc.LineCollection(lines, colors = colors, linewidths=2)
fig, ax = plt.subplots()
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)

plt.show()

visualization result

Answer 1

The key here is to find the segment on the path for each of the points we want to distribute along the path based on the cumulative distance (across segments) from the starting point on the path. Then, interpolate for the point based on the distance between the two end points of the segment in which the point is on the path. The following code does this using a mixture of numpy array processing and list comprehension:

    point_data = [
        (53.8024, 50.4762), (49.5272, 51.8727), (45.0118, 52.3863), (40.5399, 53.0184), (36.3951, 54.7708),
        (28.7127, 58.6807), (25.5306, 61.4955), (23.3828, 65.2082), (22.6764, 68.3316), (22.6945, 71.535),
        (24.6674, 77.6427), (28.8279, 82.4529), (31.5805, 84.0346), (34.7024, 84.8875), (45.9183, 84.5739),
        (57.0529, 82.9846), (64.2141, 79.1657), (71.089, 74.802), (76.7944, 69.8429), (82.1092, 64.4783),
        (83.974, 63.3605), (85.2997, 61.5455), (85.7719, 59.4206), (85.0764, 57.3729), (82.0979, 56.0247),
        (78.878, 55.1062), (73.891, 53.0987), (68.7101, 51.7283), (63.6943, 51.2997), (58.6791, 51.7438),
        (56.1255, 51.5243), (53.8024, 50.4762), (53.8024, 50.4762)
    ]

    number_of_points = 20

    def euclid_dist(x1, y1, x2, y2):
        return ((x1-x2)**2 + (y1-y2)**2)**0.5

    # compute distances between segment end-points (padded with 0. at the start)
    # I am using the OP supplied function and list comprehension, but this
    # can also be done using numpy
    dist_between_points = [0.] + [euclid_dist(p0[0], p0[1], p1[0], p1[1])
                                  for p0, p1 in zip(point_data[:-1], point_data[1:])]
    cum_dist_to_point = np.cumsum(dist_between_points)
    total_dist = sum(dist_between_points)

    cum_dist_per_point = np.linspace(0., total_dist, number_of_points, endpoint=False)
    # find the segment that the points will be in
    point_line_segment_indices = np.searchsorted(cum_dist_to_point, cum_dist_per_point, side='right').astype(int)
    # then do linear interpolation for the point based on distance between the two end points of the segment
    #   d0s: left end-point cumulative distances from start for segment containing point
    #   d1s: right end-point cumulative distances from start for segment containing point
    #   alphas: the interpolation distance in the segment
    #   p0s: left end-point for segment containing point
    #   p1s: right end-point for segment containing point    
    d0s = cum_dist_to_point[point_line_segment_indices - 1]
    d1s = cum_dist_to_point[point_line_segment_indices]
    alphas = (cum_dist_per_point - d0s) / (d1s - d0s)
    p0s = [point_data[segment_index - 1] for segment_index in point_line_segment_indices]
    p1s = [point_data[segment_index] for segment_index in point_line_segment_indices]
    results = [(p0[0] + alpha * (p1[0] - p0[0]), p0[1] + alpha * (p1[1] - p0[1]))
               for p0, p1, alpha in zip(p0s, p1s, alphas)]

The array cum_dist_to_point is the cumulative (across segments) distance along the path from the start to each point in point_data, and the array cum_dist_per_point is the cumulative distance along the path for the number of points we want to evenly distribute along the path. Note that we use np.searchsorted to identify the segment on the path (by cumulative distance from start) that the point, with a given distance from the start, lies in. According to the documentation, searchsorted:

Find the indices into a sorted array (first argument) such that, if the corresponding elements in the second argument were inserted before the indices, the order would be preserved.

Then, using the OP's plot function (slightly modified because results no longer has an angle component):

def plot_me(results):
    X = np.array([x for x, _ in results])
    Y = np.array([y for _, y in results])

    plt.scatter(X, Y)

    for i, (x, y) in enumerate(zip(X, Y)):
        plt.text(x, y, str(i), color="red", fontsize=12)

    possible_colors = [(1, 0, 0, 1), (0, 1, 0, 1), (0, 0, 1, 1)]

    lines = []
    colors = []
    for i in range(len(point_data) -1, 0, -1):
        x1, y1 = point_data[i - 1]
        x2, y2 = point_data[i]
        lines.append(((x1, y1), (x2, y2)))
        colors.append(possible_colors[i % 3])

    lc = mc.LineCollection(lines, colors=colors, linewidths=2)
    fig, ax = plt.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(0.1)

    plt.show()

We have:

    plot_me(results)