Cylinder fit to a point cloud data

35 views (last 30 days)
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.).
  6 Comments
Feri
Feri on 13 Feb 2020
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.

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 14 Feb 2020
Edited: John D'Errico on 14 Feb 2020
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.
  6 Comments
Feri
Feri on 15 Feb 2020
Well, dumb question, I down sampled the data and got the results in 2 seconds.
I really appreciate your help.
Stefan Süssmaier
Stefan Süssmaier on 21 Feb 2022
this helped a lot. thanks for sharing your functions!

Sign in to comment.

More Answers (2)

darova
darova on 13 Feb 2020
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
  12 Comments
darova
darova on 14 Feb 2020
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
Feri on 14 Feb 2020
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...

Sign in to comment.


Jacob Wood
Jacob Wood on 14 Feb 2020
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
  1 Comment
Feri
Feri on 14 Feb 2020
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!