Calculating paths around complicated shape

I have a several datasets for irregular 2D shapes. These datasets contain unordered points that make up the outline of shapes, but occasionally contain loops or islands of points. An example is shown below in fig 1. Coordinates for the shown shape are available at https://pastebin.com/0cHxPnua.

Figure 1

Purpose: I am trying to calculate a continuous line around these outlines in order to measure their perimeter, and eventually the distance between 2 points in that perimeter.

Problem: In the instances of complex shapes like the one shown in figure 1, the looped regions and disconnected regions can cause issues in generating a reasonable outline. I handle the islands by disregarding any connections that are much longer than the average distance between nearest point, but I haven't figured out what to do with the loops yet.

If you look below to fig 2, my pathfinding so far is able to handle most of the shape and successfully ignores island points, but fails after picking a bad initial path around the loop. It has to go to a new unlinked point, so it doubles back to a previously missed point near the intersection of the loop, but is then stuck as the next nearest unlinked points are beyond the distance threshold for ignoring island points.

Figure 2 Blue start, Red end. Note: The last point ALWAYS connects back to the first (red)

Goal: My goal here is either to chage my pathfinding method so that it can consistently handle complicated loops like these, or to find a wholly alternative method for pathfinding that is known to be applicable to this issue.

A test version of my code is available below to trial. You need to copy the shape's points from pastebin and paste them on the points = line

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.spatial
import sys
import os
import numpy as np
import pandas as pd
import re
from scipy.interpolate import CubicSpline
from scipy.signal import savgol_filter
from matplotlib.collections import LineCollection


def calculate_path_length(path):
    """
    Calculate the path length.

    Parameters:
    path (array): An array of 2D points in the path.

    Returns:
    float: The total length of the path.
    """
    return np.sum(np.sqrt(np.sum(np.diff(path, axis=0)**2, axis=1)))


def get_shortest_path(full_path, start_point, end_point):
    """
    Get the shortest path between start and end points in a circular path.

    Parameters:
    full_path (array): An array of 2D points in the circular path.
    start_point (array): The start point.
    end_point (array): The end point.

    Returns:
    array: The shortest path.
    """
    # Get the indices of the start and end points in the full path
    start_index = np.where((full_path == start_point).all(axis=1))[0][0]
    end_index = np.where((full_path == end_point).all(axis=1))[0][0]

    # Calculate the paths
    if start_index < end_index:
        direct_path = full_path[start_index:end_index + 1]
        indirect_path = np.concatenate((full_path[end_index:], full_path[:start_index + 1]))
    else:
        direct_path = full_path[end_index:start_index + 1]
        indirect_path = np.concatenate((full_path[:end_index + 1], full_path[start_index:]))

    # Calculate the lengths of the paths
    direct_path_length = calculate_path_length(direct_path)
    indirect_path_length = calculate_path_length(indirect_path)

    # Return the shortest path
    print(direct_path_length, indirect_path_length)
    if direct_path_length < indirect_path_length:
        return direct_path
    else:
        return indirect_path


def GetPointDistance(p1, p2):
    return np.sqrt( ((p2[0] - p1[0])**2) + ((p2[1] - p1[1])**2) )


def AverageDistanceWithBuffer(points, buffer_size=20):
    """
    Calculate the average distance of a group of points and add buffer_size standard deviations to 
    add a large ceiling for distance variances. This should still be able to exclude extremely 
    large deviations that occur when outliers are present.

    Parameters:
    points (list): 2D list of 2xN points.
    buffer_size (int): Padding multiplier to use for outlier rejection.

    Returns:
    float: The average distance with buffer.
    """
    distances = []
    for i, point in enumerate(points):
        if i != len(points) - 1:
            distances.append(GetPointDistance(point, points[i+1]))
    return np.mean(distances) + buffer_size * np.std(distances)
    

def random_downsample(arr, percentage):
    """
    Randomly downsample an array by the given percentage.

    Parameters:
    arr (array): The input array.
    percentage (float): The percentage of data to retain.

    Returns:
    array: The downsampled array.
    """
    if percentage <= 0 or percentage > 100:
        raise ValueError("Percentage should be between 0 and 100 (inclusive).")

    sample_size = int(len(arr) * percentage / 100)
    random_indices = np.random.choice(len(arr), sample_size, replace=False)
    sampled_arr = arr[random_indices]

    return sampled_arr



def shortest_path(points, start_point, end_point):
    """
    Find the shortest path between start and end points.

    Parameters:
    points (array): An array of 2D points.
    start_point (array): The start point.
    end_point (array): The end point.

    Returns:
    array: The full path.
    array: The shortest path between the start and end points.
    """
    points = points.tolist()
    start_point = start_point.tolist()
    end_point = end_point.tolist()
    points = points + [start_point, end_point]
    current_point = start_point
    full_path = [current_point]
    points.remove(current_point)

    use_cutoff = False
    while points:
        if len(full_path) == 1000:
            buffered_avg_dist = AverageDistanceWithBuffer(full_path)
            use_cutoff = True

        if len(points) == 1 and points[0] == end_point:
            full_path.append(end_point)
            points.remove(end_point)
        else:
            closest_point = min(points, key=lambda x: np.linalg.norm(np.array(current_point) - np.array(x)))
            closest_point_dist = GetPointDistance(current_point, closest_point)
            if use_cutoff:
                if closest_point_dist < buffered_avg_dist:
                    full_path.append(closest_point)
                    current_point = closest_point
            else:
                full_path.append(closest_point)
                current_point = closest_point
            points.remove(closest_point)
        

    # To form a closed loop for the full perimeter path
    full_path.append(start_point)

    # Shortest path from start to end
    shortest_path = get_shortest_path(np.array(full_path), start_point, end_point)

    plt.scatter(shortest_path[:, 0], shortest_path[:, 1], c=range(len(shortest_path)), cmap='inferno')
    plt.close()

    # Convert lists into numpy arrays
    full_path = np.array(full_path)
    shortest_path = np.array(shortest_path)

    return full_path, shortest_path


points = # Copy points from https://pastebin.com/0cHxPnua
points = np.array(points)

# Returns a full path and the shortest path between the start and end points
full, shortest = shortest_path(points, points[48], points[int(len(points)/2)])  # Arbitrary Start and End points

# Plot the outline
plt.scatter(points[0, 0], points[0, 1], s=200, facecolor='red', zorder=1000)  # Plot start point
plt.scatter(points[int(len(points)/2), 0], points[int(len(points)/2), 1], s=200, facecolor='green', zorder=1000)  # Plot end point
plt.scatter(points[:, 0], points[:, 1])  # Plot all points
plt.plot(full[:, 0], full[:, 1], color='k')  # Plot path line
plt.axis('equal')
plt.show()