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!
You did not specify a programming language, so if C++ works for you, PCL could be what you are looking for. PCL has an implementation of the RANSAC method to find a geometric model in noisy data. Many different models are supported, most interesting for you would be CIRCLE2D, CIRLCE3D, and CYLINDER (see https://pointclouds.org/documentation/group__sample__consensus.html)
Here is a tutorial on how to fit a cylinder: https://pcl.readthedocs.io/projects/tutorials/en/master/cylinder_segmentation.html