My goal is to fit a line through a point cloud. The point cloud is approximately cylindrical but can be curved, which is why the fitted line should not be straight. I’ve tried several things, but this is my current approach :
I use PyTorch to optimise on a surface equation and compute the loss for each point. However, this does not seem to lead to “good” results, since the plane/surface does not cut the point cloud vertically, which I would expect to lead to the least error. Since I’m new to PyTorch I don’t know whether it’s a mistake in my code or whether there’s a mathematical problem with the idea. This approach is in the code below.
Additionally, once I’ve fitted this surface, I would like to obtain a line on it, but I’m unsure on the best way to do this. My questions are : Why is the fit of the plane/surface not more centred ? How could I then obtain a curved line of best fit on the resulting surface ?
Other approaches I've tried :
- SVD, but like I said, im looking for a line that is not straight
- Computing the mean of each z-slice (i.e. each horizontal slice) in the y and x direction and using the resulting points to fit a spline (using python’s splrep). The problem with this approach is that a couple of top or bottom slices can consist of very few points that are not necessarily in the “center” of the point cloud (e.g. are on the left extreme of the top slice), which forces the spline in that direction. It’s really important for the project that the line is “central”.
The data in this code example is just toy data, the real data would have a more complex shape. To check the idea I'm fitting a plane, but in the end I'd like a more complex surface. The code results in this figure.
import numpy as np
import torch
import matplotlib.pyplot as plt
from torch import nn
from torch.functional import F
ix = np.random.uniform(0,2, 1000)
iy = np.random.uniform(0,2, 1000)
iz = np.random.uniform(0,100, 1000)
x = torch.tensor(np.array([ ix, iy ]).T).float()
y = torch.tensor(iz).float()
degree =1
class Model(nn.Module):
def __init__(self):
super().__init__()
if degree == 1 :
weights = torch.distributions.Uniform(-1, 1).sample((4,))
if degree == 3 :
weights = torch.distributions.Uniform(-1, 1).sample((8,))
self.weights = nn.Parameter(weights)
def forward(self, X):
if degree ==1 :
a_1, a_2, a_3, a_4 = self.weights
return (a_1 *X[:,0] + a_2*X[:,1] +a_3)
if degree == 3 :
a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8 = self.weights
return (a_1 * X[:,0]**3 + a_2*X[:,1]**3 + a_3*X[:,0]**2 + a_4 *X[:,1] **2 + a_5 *X[:,0] + a_6 *X[:,1] + a_7)
def training_loop(model, optimizer):
losses = []
loss = 10000
it = 0
if degree == 3 :
lim = 0.1
if degree == 1 :
lim = 0.1
while loss > lim:
it +=1
if it > 5000:
break
preds1= model(x).float()
l1 = torch.nn.L1Loss()
loss = l1(preds1, y)
loss.backward()
optimizer.step()
optimizer.zero_grad()
losses.append(loss.detach().numpy())
print(loss)
return losses
m = Model()
if degree == 1 :
opt = torch.optim.Adam(m.parameters(), lr=0.01)
losses = np.array(training_loop(m, opt))
if degree == 3 :
opt= torch.optim.Adam(m.parameters(), lr=0.001)
losses = np.array(training_loop(m, opt))
params=list(m.parameters())[0].detach().numpy()
X = np.arange(0, 2, 0.1)
Y = np.arange(0, 2, 0.1)
X, Y = np.meshgrid(X, Y)
if degree == 1 :
Z = (params[0] * X + params[1]*Y + params[2])
if degree == 3:
Z = (params[0] * X**3 + params[1]*Y**3 + params[2]*X**2 + params[3]*Y**2 + params[4]*X + params[5]*Y + params[6])
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
surf = ax.plot_surface(X, Y, Z, color='tab:orange', alpha = 0.5,linewidth=0, antialiased=False)
ax.scatter3D(ix,iy,iz, alpha = 0.3, s=2)
plt.show()
Edit : I tried the approach in this post : Fit Curve-Spline to 3D Point Cloud However this forces me to specify a source and target for the shortest path. Simply choosing the centers of the lowest and highest slice results in this :
