Поиск всех возможных линий внутри набора точек

Я хочу обнаружить все возможные комбинации «прямых» линий (при условии соблюдения критериев (условие внутри функции extend_line(): if x_sign == new_x_sign and y_sign == new_y_sign and np.sign(slope) == np.sign(new_slope) and abs(slope - new_slope) <= 2)) из моего выборочного набора точек. Поэтому я попробовал подход с вложением функции расширения линии extend_line(), чтобы функция запускалась столько раз, сколько необходимо, чтобы объединить все сегменты линии, соответствующие критериям. Чтобы добиться этого, я вычисляю ближайшие точки для каждой отдельной точки, а затем проверяю комбинации линий.

Основная идея работает для первой «итерации», и обнаруживается первая «прямая» линия, отвечающая критериям:

Но как только последнее выполнение вложенного extend_line() завершено, у меня возникает проблема: я не могу избавиться от всех сегментов линии, которые принадлежат предыдущей строке, но не принадлежат текущей строке (коричневый и фиолетовый сегменты линий должны быть удалены:

Я пробовал манипулировать разными переменными (например, helper_list), но это не сработало, и как только вложенная функция была выполнена, переменные были «сброшены» до переменных из исходной функции extend_line(). Я также попробовал усечь helper_list, где я сохраняю сегмент линии, чтобы он соответствовал количеству посещенных точек. Это казалось довольно многообещающим, но мой способ/идея не сработала должным образом.

import geopandas as gpd
from shapely.geometry import Point, LineString, MultiLineString
from shapely.ops import linemerge
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist

def calc_slope(line):
    x1,y1 = line.coords[0]
    x2,y2 = line.coords[1]
    try:
        s = (y2-y1)/(x2-x1)
    except:
        s = 0
    return s

def extend_line(start_line, current_line, current_point, x_sign, y_sign, slope, points, visited_points, nearby_points_dict, helper_list=[]):
    # Add current line to helper_list if first "iteration"
    if len(visited_points) == 1:
        helper_list.append(current_line)
    # Add the current point to visited points
    visited_points.append(tuple(current_point))
    
    # Find nearby points that have not been visited yet
    key = np.where(points == current_point)
    nearby_points = [p for p in nearby_points_dict[key[0][0]] if tuple(p) not in visited_points]
    
    # Iterate over nearby points to extend the line
    for next_point in nearby_points:
        # # Truncate list to match amount of visited points
        # if len(visited_points) < len(helper_list):
        #     helper_list = helper_list[:len(visited_points)]
        # Create a new line from current point to next point
        new_line = LineString([current_point, next_point])
        
        # Check for signs in order to extent only points "in front"
        new_dx = current_point[0] - next_point[0]
        new_x_sign = np.sign(new_dx)
        new_dy = current_point[1] - next_point[1]
        new_y_sign = np.sign(new_dy)
        
        # Calculate slope of new line segment
        new_slope = calc_slope(new_line)
        
        # Check if conditions are met (same sign -> point in front, same lope sign -> same tilted line segment, slope witin threshold -> keeps "straight" line)
        if x_sign == new_x_sign and y_sign == new_y_sign and np.sign(slope) == np.sign(new_slope) and abs(slope - new_slope) <= 2:
            # Append valid line to the list of valid lines
            helper_list.append(new_line)
            
            # Recursively extend the line from the next point
            extend_line(start_line, new_line, next_point, x_sign, y_sign, new_slope, points, visited_points.copy(), nearby_points_dict, helper_list)
            
            if len(helper_list) > 1:
                mls = MultiLineString(helper_list)
                ls = linemerge(mls)
                valid_lines.append(ls)
                fig, ax = plt.subplots(figsize=(8, 8))
                gdf.plot(ax=ax, marker='o', color='gray', markersize=50, label='Original Points')
                try:
                    ax.plot(*ls.xy)
                except:
                    for geom in mls.geoms:
                        ax.plot(*geom.xy)
                plt.show()

# Sample data for GeoDataFrame
data = {
    'id': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
           11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
           21, 22, 23, 24, 25, 26, 27, 28, 29, 30],
    'x': [0.0, 1.0, 3.0, 2.0, 5.0, 4.0, 7.0, 6.0, 9.0, 8.0,
          1.5, 2.5, 4.5, 3.5, 6.5, 5.5, 8.5, 7.5, 10.0, 9.0,
          3.0, 4.0, 6.0, 5.0, 8.0, 7.0, 10.0, 9.0, 11.0, 10.0],
    'y': [0.0, 2.0, 1.0, 3.0, 0.5, 4.0, 2.5, 5.0, 1.5, 6.0,
          3.5, 5.0, 2.0, 6.0, 4.5, 8.0, 3.0, 7.0, 5.5, 9.0,
          7.0, 8.0, 6.0, 9.0, 8.5, 10.0, 7.5, 11.0, 9.0, 12.0]
}

# Create GeoDataFrame
geometry = [Point(xy) for xy in zip(data['x'], data['y'])]
gdf = gpd.GeoDataFrame(data, geometry=geometry)
gdf["id"] = gdf["id"] - 1

points = gdf[['x', 'y']].values

# Calculate distance matrix
pairwise_distances = cdist(points, points)

# Find indices of nearby points based on the threshold distance
nearby_indices = np.where((pairwise_distances > 1) & (pairwise_distances < 3))

# Exclude self-distances (where distance is zero)
nearby_indices = np.array([i for i in zip(*nearby_indices) if i[0] != i[1]])

# Create a list to store nearby points for each point
nearby_points_dict = {}

# Loop over each point to collect nearby points
for i in range(len(points)):
    nearby_points = points[nearby_indices[nearby_indices[:,0]==i,1]]
    nearby_points_dict[i] = nearby_points
all_lines = []
for i, nearby_points in nearby_points_dict.items():
    start_point = points[i]
    
    # Plot sample points and area to be considered nearby
    fig, ax = plt.subplots(figsize=(8, 8))
    gdf.plot(ax=ax, marker='o', color='gray', markersize=50, label='Original Points')
    for k, v in gdf.iterrows():
        ax.annotate(f'{v["id"]}\n({v["x"]},{v["y"]})',(v["x"],v["y"]), ha='center')
    ax.plot(points[i,0],points[i,1], marker='o', color='blue', label='Current Point')
    circle = plt.Circle(points[i], 3, color = "green", alpha=0.25)
    ax.add_patch(circle)
    
    # Iterate over all nearby_points of start_point
    for point in nearby_points:
        end_point = point
        
        # Check signs of points in order to only extend points "in front" of start_point
        dx = start_point[0] - end_point[0]
        x_sign = np.sign(dx)
        dy = start_point[1] - end_point[1]
        y_sign = np.sign(dy)
        
        line = LineString([start_point, end_point])
        ax.plot(point[0],point[1], marker='o', color='red', linestyle='', label='Nearby Points')
        
        # Calculate slope of line in order to only extend with line segments within slope threshold
        slope = calc_slope(line)
        
        visited_points = [tuple(start_point)]
        valid_lines = []
        extend_line(line, line, end_point, x_sign, y_sign, slope, points, visited_points, nearby_points_dict)
        
    all_lines.append(valid_lines)
    break

Обновлено: Также интересно, почему мне нужно получить доступ к функции keyinside extend_line() с помощью key[0][0], если каждая точка в points существует только один раз? Почему key распечатанный возврат (array([ 1, 1, 12], dtype=int64), array([0, 1, 1], dtype=int64)) при первом выполнении extent_line() и не только, например, (array([1,], dtype=int64)?