Fit a cylinder to scattered 3D XYZ point data

Question

As in the title, I want to fit a cylinder to a group of 3D points with Python. This is a nice solution with MATLAB. How can we do it with Python?

Wildon

You could convert the matlab function? – Bonny

The Matlab function is quite complicated. I hope we have a library that supports this function (or with minimal code) using Python. – Wildon

Answer 1

Using scipy.optimize.leastsq, we can create an error function in which the difference between the observed cylinder radius and the modelled radius is minimized. The following is an example of fitting a vertical cylinder

import numpy as np
from scipy.optimize import leastsq


def cylinderFitting(xyz,p,th):

    """
    This is a fitting for a vertical cylinder fitting
    Reference:
    http://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XXXIX-B5/169/2012/isprsarchives-XXXIX-B5-169-2012.pdf

    xyz is a matrix contain at least 5 rows, and each row stores x y z of a cylindrical surface
    p is initial values of the parameter;
    p[0] = Xc, x coordinate of the cylinder centre
    P[1] = Yc, y coordinate of the cylinder centre
    P[2] = alpha, rotation angle (radian) about the x-axis
    P[3] = beta, rotation angle (radian) about the y-axis
    P[4] = r, radius of the cylinder

    th, threshold for the convergence of the least squares

    """   
    x = xyz[:,0]
    y = xyz[:,1]
    z = xyz[:,2]

    fitfunc = lambda p, x, y, z: (- np.cos(p[3])*(p[0] - x) - z*np.cos(p[2])*np.sin(p[3]) - np.sin(p[2])*np.sin(p[3])*(p[1] - y))**2 + (z*np.sin(p[2]) - np.cos(p[2])*(p[1] - y))**2 #fit function
    errfunc = lambda p, x, y, z: fitfunc(p, x, y, z) - p[4]**2 #error function 

    est_p , success = leastsq(errfunc, p, args=(x, y, z), maxfev=1000)

    return est_p

if __name__=="__main__":

    np.set_printoptions(suppress=True)    
    xyz = np.loadtxt('cylinder11.xyz')
    #print xyz
    print "Initial Parameters: "
    p = np.array([-13.79,-8.45,0,0,0.3])
    print p
    print " "

    print "Performing Cylinder Fitting ... "
    est_p =  cylinderFitting(xyz,p,0.00001)
    print "Fitting Done!"
    print " "


    print "Estimated Parameters: "
    print est_p

Blab

Just a remark but the th parameter is ignored inside the function. – Gerik

Can you elaborate on how exactly did you come up with this function? – Veedis

Answer 2

There is paper at David Eberly site "Fitting 3D Data with a Cylinder" that describes math basics and shows pseudocode.

You can also refer to C++ code in Geometric Tools Engine at the same site. I think that some auxiliary math functions like matrix inverse etc could be implemented in NymPy.

Nobelium

Answer 3

Cylinder Image

I had a similar situation of fitting a cylinder through several points. we measure the gap between 2 cylinders using a gap sensor at several points and I had to visualize how gap varies in comparison with the cylinder.

I used ax.plot_surface(x, y, z, alpha=0.5) where x , y, z are numpy arrays of 3D location of all the points. Take a look at image below.

code snippet,

# Extract X,Y,Z values from the sensor data 
for i in range(num_of_sensors):
  ax.scatter(x[:,i], y[:,i], z[:,i], color = "k", marker=".", s=5, cmap='hot')
  ax.text(x[0,i]+ 10, y[0,i] + 10, z[0,i]+100, '%s' % (idx[i]), size=5, zorder=1, color='b') 

# plot the surface
ax.plot_surface(x, y, z, alpha=0.5)

Aletheaalethia

Answer 4

This repo by xingjiepan allows you to compute the best fit cylinder using Python.

The algorithm is by David Eberly.

As stated by David Eberly, the main assumption is that the underlying data is modelled by a cylinder and that errors have caused the points not to be exactly on the cylinder.

It's very easy to use:

from cylinder_fitting import fit
w_fit, C_fit, r_fit, fit_err = fit(data)

I have also developed an object-oriented implementation of xingjiepan's repo with input validation and type hints. You can find it here and can be used as follows:

>>> from py_cylinder_fitting import BestFitCylinder
>>> from skspatial.objects import Points
>>> points = [[2, 0, 0], [0, 2, 0], [0, -2, 0], [2, 0, 4], [0, 2, 4], [0, -2, 4]]
>>> best_fit_cylinder = BestFitCylinder(Points(points))
>>> best_fit_cylinder.point
[0., 0., 0.]
>>> best_fit_cylinder.vector
[0.  0.  4.]
>>> best_fit_cylinder.radius
2.0

We are also working to add this feature to the scikit-spatial library.

Sherrard

Now available also in the scikit-spatial library. – Sherrard