Cylinder fit to a point cloud data

I'm trying to fit a cylinder with a known radius of curvature (R = 25.86) to a set of data points in space. The data point has orientations with respect to X, Y, and Z axes. I am looking for the minimum residual after subtracting the form. How can I do that (without using Computer Vision Toolbox etc.).

[Rik]

It is probably the easiest to try to fit the rotation and scaling matrix that converts your Cartesian coordinates to cylinder coordinates.

[Feri]

Thank you for your response. Would you be able to provide more details about the rotation and scaling matrix? How can I obtain such matices?

[Feri]

What I have includes X, Y, Z coordinates of the cylinder.

[Star Strider]

I am not certain what approach you want to take. Consider the cart2pol function if you want to transform your Cartesian coordinates to cylindrical coordinates to fit them.

[darova]

Attach the data

[Feri]

I have attached three .mat files here, [X, Y, Z], the whole data file is very huge and cannot be uploaded here, my problem is that how can I subtract the best cylinder fit from this data, and obtain the residual.

Thank you in advance.

Answer 1

[John D'Errico]

You seem to be way overthinking this.

xyz = REAL_data;
xyz0 = mean(xyz)
xyz0 =
    -0.024621       -18.89      0.75025
xyzhat = xyz - xyz0;
[V,D] = eig(xyzhat'*xyzhat)
V =
   1.3309e-05   0.00074361           -1
      0.17374      0.98479   0.00073461
      0.98479     -0.17374  -0.00011609
D =
       1387.3            0            0
            0   6.4145e+06            0
            0            0   1.0888e+08
trans = V(:,1:2)
trans =
   1.3309e-05   0.00074361
      0.17374      0.98479
      0.98479     -0.17374
uv = xyzhat*trans;
plot(uv(:,1),uv(:,2),'.')

I've used what is essentially principal components to compute the transformation here.

So projecting the data into a plane that is perpendicular to the axis of the cylinder, we see:

Now, can we find the circle parameters in the subspace uv? The simplest way is...

N = 100000;
ind1 = randperm(size(uv,1),N);
ind2 = randperm(size(uv,1),N);

uvcenter = (2*uv(ind1,:) - uv(ind2,:)) \ sum(uv(ind1,:).^2 - uv(ind2,:).^2,2)
uvcenter =
       130.65
    -0.014346

The center was found by another trick. One way to visualize it is ... Pick any two points at random on the "circle". If we look at the line segment through the points, the perpendicular to that line through the midpoint goes through the center of the circle. What I actually did was to take the equations of a circle with unknown center and radius through each point, then subtract. The difference yields a set of LINEAR equations in the unknown center of the circle.

uvcenter is in the uv plane. Now, compute the estimated radius, as the average distance from the center to every one of the points in the circle.

mean(sqrt(sum((uvcenter.' - uv).^2,2)))
ans =
       130.69

We can convert this back into the original coordinate system, thus xyz as...

xyzcenter = xyz0 + uvcenter'*trans'
xyzcenter =
    -0.022893       3.7941       129.41

That is a point on the centerline of the "cylinder".

V(:,3)
ans =
           -1
   0.00073461
  -0.00011609

The centerline of the cylinder moves along the vector V(:,3). And the cylinder is seen to have radius 130.69.

Actually, pretty easy, as I said.

[darova]

V = [      -1
   0.00073461
  -0.00011609];
acosd(-V(1)/norm(V))
ans =
    0.0426

angle is 1.7 degree, not 0.04

Am im missing something?

[Feri]

Thank you John,

I tried the eigen vectores that the covarriance matrix gives, in order to find the principal axies. However, I don't know how to get to the minimum residual (below right) from the raw data (below left). I have attached a reduced data file of a section of a cylinder with R=25.86 mm.

How can I get the best fit cylinder to obtain the minimum residual?

[John D'Errico]

I'll make it even easier. I wrote some code a few years ago. It takes a set of data, then uses a scheme not unlike like what I showed above, but then pushes it into an optimizer to refine the fit. It uses either fminunc (if you have the optimization toolbox) or fminsearch if not.

[axispoint,axisvec,radius,rmse] = cylinderfit(xyz)
axispoint =
   -0.0012239    -0.085028      -43.939
axisvec =
      0.99956     0.029657  -0.00014082
radius =
       25.856
rmse =
   0.00029309

Thus cylinderfit does the hard work, internally calling circlefit along the way, another tool I wrote for fitting a circle in some number of dimensions.

As you can see here, on Reduced_data, it tells you a point on the axis of the cylinder, a vector that points along the axis, the radius of the cylinder (which matches what you said) and the RMSE, which here is pretty small, since the data is very nearly perfectly cylindrical.

[Feri]

Thank you so much John, this works perfect, and the arccos of axisvec(1) is equal to -1.7, great.

My last question I promis, can I get this major axis vector using the eigen vectors of the covarriance matrix, just to speed things up? (or how can I make this process faster if possible?).

Running your previous script won't give the same major axis, is this because of the lack of optimization?

Thanks a million.

[Feri]

Well, dumb question, I down sampled the data and got the results in 2 seconds.

I really appreciate your help.

[Stefan Süssmaier]

this helped a lot. thanks for sharing your functions!

Answer 2

[darova]

Here is an attempt

    % load data
load('Data.mat');
x = REAL_data(:,1);
y = REAL_data(:,2);
z = REAL_data(:,3);
ix = 1:1000:length(x);          % extract only every 100th point
plot3(x(ix),y(ix),z(ix),'.r')

    % fit data in YZ plane
f1 = @(a,b,x) -sqrt(109^2 - (x-a).^2) + b;
f = fit(y(ix),z(ix),f1)
z1 = f(y(ix));
hold on
plot3(z1*0,y(ix),z1,'-b')
hold off
axis equal
view(90,0)

figure
plot(z(ix)-z1)                  % residuals

[Feri]

Thank you for the code.

Running the whole data points the residual looks like this, there is still a cylinderical form present. How can I get rid of this?

[darova]

I dont understand. What data do you use?

[Feri]

I used entire data point (too big to be able to upload here).

When running your code on it and plotting the residual in 3D I get the above figure.

The residual of this data set (at each point) should be in the order of ~10e-4.

[darova]

Can you upload reduced data?

ix = 1:1000:length(x);          % extract only every 100th point
x1 = x(ix);
y1 = y(ix);
z1 = z(ix);
save rdata

[Feri]

my apologies, the radius of curvature is 25.86, but still doesn't give the minimum residual.

Is there a way to fit the data to both X and Y in the fit command?

[Feri]

Attached is the reduced data.

Thank you so much.

[darova]

I rotated data around Z axis and sorted along Y

The result

Check the attached script

[Feri]

Thank you.

How did you figure out the rotation angle? How can this be done automatically?

[darova]

I think it will be too complicated to calculate the angle. I just tried every angle to get pure circle (manually)

Then experimented with residuals (to get minimum)

[Feri]

This angle that you got is exactly what I get when putting the whole precess in a loop to change that angle, calculate the residual and continue until the minimum residual is achieved.

Do you have any idea how to get this angle (I know you said it's complicated), can this be done in the fitting command, as C constant something like "f1 = @(a,b,c, x) sqrt(R^2 - (C*x-a).^2) + b;" or somethinng like it?

Thank you.

[darova]

Maybe for loop? Just loop through every angle and see max residual

a = -10:10;
max_residual = a*0;
for i = 1:length(a)
    % a = 1.7;                      % rotate data around Z axis 2 degree
    v = [cosd(a(i)) 1;-sind(a(i)) 1]*[x y]';
    x1 = v(1,:)';
    y1 = v(2,:)';
    [~,ix1] = sort(y1);            % sorted Y data
    ix2 = ix1(1:10:end);          % extract only every 1000th point
    % plot3(x1(ix2),y1(ix2),z(ix2),'.r')
    
        % fit data in YZ plane
    R = 25.86;
    f1 = @(a,b,x) sqrt(R^2 - (x-a).^2) + b;
    f = fit(y1(ix2),z(ix2),f1);
    z1 = f(y1(ix2));
    
    max_residual(i) = max(abs(z1-z(ix2)));
end
plot(a,max_residual)
title('max residual')

[Feri]

This appears to be the only option left. and time consuming.

I am comparing my MATLAB results with a software that takes less than a second to remove the best fit cylinder and output the minimum residual. ANd it's killing me that what is it that I don't get...

Answer 3

[Jacob Wood]

Not the fastest running solution, but I think rather straightforward. This uses fminsearch() to find the 6 parameters that define the cylinder's centerline and 1 parameter that defines the cylinder radius. If you wish to use a known radius, rather than fit one, you can tailor the solution.

%% Define inputs
d = Reduced_data;
r0 = 25.86;

%X(1,2,3) - point on the cylinder centerline
%X(4,5,6) - vector that points along centerline
%X(7) - radius of cylinder
X0 = [0,0,-r0,1,0,0,r0]; %solution starting point


%% Do minimization
%Minimize the average residual of all the points
Xr = fminsearch(@(x) get_avg_res(x,d),X0);


%% Plot results
dis = dist_point_line(Xr(1:3),Xr(4:6),d);
signed_res = dis-Xr(7);

figure
scatter3(d(:,1),d(:,2),d(:,3),[],d(:,3),'filled');
colormap('jet')
colorbar

figure
scatter3(d(:,1),d(:,2),res,[],res,'filled');
colormap('jet')
colorbar
view(0,90)


%% Helpers
function avg_res = get_avg_res(X,d)
    dis = dist_point_line(X(1:3),X(4:6),d);
    avg_res = sum(abs(dis-X(7)))/size(d,1);
end

function d = dist_point_line(line_p,line_s,points) 
%https://onlinemschool.com/math/library/analytic_geometry/p_line/
    M0M1 = line_p - points;
    cp = cross(M0M1,repmat(line_s,size(points,1),1));
    d = sqrt(cp(:,1).^2+cp(:,2).^2+cp(:,3).^2) / norm(line_s);
end

[Feri]

Thank you Jacob,

Running you script the residual still appearas to have some sort of form left in it, due to its initial orientation I guess.