[SOLVED] A way to fit circles on 3D data

What I am trying to do is to fit a circle across all cross-sections of 3D point cloud data that resembles a pipe.

The data resembles a cylindrical shape with disconnected parts in between, which can be ignored at this point. These data are composed of X, Y, and Z component.

At first I thought it would be better to take a cross section of the data, and find the best circle to fit on it, then apply it across the entire data to get many circles.

I first was able to obtain the curve-fitting after swapping X and Z axis to make the data explicit(thanks to jlandercy for answering my previous question(Describing the data points with arc-shaped graph). Then extracted a circle from the point where slope is 0. This is the idea coming from the fact that a circle should be able to fit where it is 90 degrees between a horizontal and a vertical lines, as depicted with the following picture: .

To explain it in more detail, what I did was finding the point along the curve of the cross-section data where the slope is 0 and finding another point where slope is the infinite. But as there is no point along the curve where the slope is infinite, I have assumed that a perfect circle could not fit along the entire curve, but the best-fitting one could be drawn on the point where slope is 0. So I found the coordinates of the point whose slope is 0 and compared it with the lowest point of the cross-section data. I used their coordinates to find the center of the circle, and I used the difference of X component between the point of slope 0 and the center to find the radius. As a result, I was able to draw a circle that seems to fit best with the curve on the cross section as displayed in yellow in the picture below.

I ignored Y components at that time, because I found that their variance and standard deviation were very low compared to those of the X and Z components, meaning that they don't really affect my analysis.

Now my plan is trying to apply it across the entire data so that I can create circles for every cross section. Somewhat like below.

The ultimate goal is to connect the centers of these circles on every cross section of the data, to obtain a non-straight line, which I believe wouldn't be too difficult once I am able to obtain the circles on all the cross sections.

I hope that what I am trying to achieve is clear. I just have no idea how to apply my cross-section method across the original data at this point. I was thinking about taking a for loop to apply it across the entire data but I couldn't quite imagine how to specify the loop for it. I was trying to get some help of Chat-GPT, but it didn't really help me.

I also have tried other methods of finding certerline of the point cloud such as combination of Delanay Triangulation and Voronoi Tessellation as suggested by another post Fit Curve-Spline to 3D Point Cloud. But I was not able to succeed.

If anyone can give me some hint or advice, I would appreciate it.

The original data is too large so I am linking on my OneDrive with a link: Original Data

Here is the piece of my code for how I approached with my idea of cross-section circle fitting.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Load the point cloud data
data = pd.read_csv('data/selected_2.csv')
points = data.to_numpy()
x = points[:, 0]
y = points[:, 1]
z = points[:, 2]

# Calculate polynomial
a = np.polyfit(z, x, 7)
f_poly = np.poly1d(a)

# Find new x and y data along the curve
z_new = np.linspace(z[0], z[-1], 50)
x_new = f_poly(z_new)
data_new = np.vstack((z_new, x_new)).T

# Define the derivative function
def FindingDerivatives(data_new):
    slopes = []
    for i in range(len(data_new) - 1):
        z1, x1 = data_new[i]
        z2, x2 = data_new[i + 1]
        slope = (x2 - x1) / (z2 - z1)
        slopes.append((z1, slope))
    return slopes

# Define the function to find perpendicular line slopes
def FindingPerpendicular(slopes):
    perpendicular = []
    for z1, slope in slopes:
        if slope != 0:
            inverse = -1 / slope # Perpendicular line is the inverse of the original slope
        else:
            inverse = float('inf')  # Perpendicular to a horizontal line is a vertical line
        perpendicular.append((z1, inverse))
    return perpendicular

# Find the maximum slope element
def FindMinMaxSlopeElement(slopes):
    target_slope_elements = []
    max_slope_elements = 0
    min_slope_elements = 0
    max_slope = max(slopes)[1]
    min_slope = min(slopes)[1]
    print(max_slope, min_slope)
    for i in range(len(slopes)):
        if slopes[i][1] == min_slope:
            max_slope_elements = slopes[i]
        if slopes[i][1] == max_slope:
            min_slope_elements = slopes[i]
    target_slope_elements.append(min_slope_elements)
    target_slope_elements.append(max_slope_elements)
    return target_slope_elements

# Find the corresponding coordinates to the maximum slope element
def FindCorrespondingCoordinates(target_slope_elements, data):
    min_coordinates = []
    max_coordinates = []
    target_coordinates = []
    for i in range(len(data)):
        if data[:, 0][i] == target_slope_elements[0][0]:
            min_coordinates.append(data[:, 0][i])
            min_coordinates.append(data[:, 1][i])
        if data[:, 0][i] == target_slope_elements[1][0]:
            max_coordinates.append(data[:, 0][i])
            max_coordinates.append(data[:, 1][i])
    print(min_coordinates)
    print(max_coordinates)
    target_coordinates.append(min_coordinates)
    target_coordinates.append(max_coordinates)
    return target_coordinates

# Find the circle based on the coordinate values corresponding to the maximum and minimum slopes
def FindCircle(target_coordinates):
    circle_info = []
    center_z = target_coordinates[0][1] 
    center_x = target_coordinates[1][0]
    center = [center_z, center_x]
    radius = target_coordinates[0][0] - target_coordinates[1][0]
    print(center, radius)
    circle_info.append(radius)
    circle_info.append(center)
    return circle_info

# Draw the circle
def CircleDraw(circle, test_data, z_new, x_new):
    radius = circle[0]
    fig, ax = plt.subplots()
    ax.scatter(test_data[:, 1], test_data[:, 0])
    ax.plot(z_new, x_new, 'y')
    cir = plt.Circle((circle[1][0], circle[1][1]), radius, color='r',fill=False)
    center_label = f"center: \n {round(circle[1][0], 2), round(circle[1][1], 2)}"
    plt.annotate(center_label, xy=(circle[1][0], circle[1][1]), ha='center') 
    ax.set_aspect('equal', adjustable='datalim')
    ax.add_patch(cir)
    ax.set_xlabel('Z Data')  # X axis data label
    ax.set_ylabel('X Data')  # Y axis data label
    plt.show()


# Organize the test data (only taking 'X' and 'Z' coordinates of the original data)
test_data = points[:, [0, 2]]
# Find the slope of the original data
slopes = FindingDerivatives(test_data)
# Find perpendicular lines of the original data
perpendicular = FindingPerpendicular(slopes)
# Find minimum and maximum slope elements 
target_slope_elements = FindMinMaxSlopeElement(slopes) 
print("Target Slopes", target_slope_elements)
# Find the target coordinates
target_coordinates = FindCorrespondingCoordinates(target_slope_elements, test_data) 
print("Target Coordinates: ", target_coordinates)
# Find the circle information
circle = FindCircle(target_coordinates)
print("radius: %s center: %s" % (circle[0], circle[1]))
# Call the circle function
CircleDraw(circle, test_data, z_new, x_new)

# Here is the cross-section of the data
# ('Y' components are not really used as I decided that they don't really affect my analysis)
X,Y,Z
1702113.3399999738,2434223.5100000203,9080.16
1702132.699999988,2434227.2700000107,9079.86
1702172.7199999988,2434233.629999995,9106.22
1702153.0400000215,2434230.4699999997,9087.98
1702143.6899999976,2434227.5699999933,9083.15
1702167.079999983,2434224.6899999976,9100.11
1702159.8799999952,2434236.3600000143,9093.52
1702177.9900000098,2434225.169999987,9113.97
1702192.0399999917,2434220.0899999742,9155.37
1702186.7199999988,2434228.7599999905,9128.54
1702191.330000013,2434221.2299999893,9174.77
1702182.9099999962,2434234.8899999857,9120.3
1702182.580000013,2434217.7899999917,9121.94
1702192.1200000048,2434219.9599999785,9164.910000000002
1702189.5399999917,2434224.2199999997,9137.12
1702191.229999989,2434221.3899999857,9146.08
1702188.9300000072,2434225.100000024,9185.21
1702183.5,2434216.300000012,9198.22
1702175.5,2434229.150000006,9220.28
1702178.099999994,2434224.9799999893,9209.75
1702188.080000013,2434226.4600000083,9223.87
1702184.669999987,2434231.9600000083,9196.33
1702214.9900000098,2434218.3599999845,9220.34

Please comment if my question is not clear enough. I will update and revise the question as much as necessary.

Thanks a lot!