Cylinder fit to a point cloud data

35 visualizaciones (últimos 30 días)
Feri
Feri el 13 de Feb. de 2020
Comentada: Stefan Süssmaier el 21 de Feb. de 2022
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 comentarios
darova
darova el 13 de Feb. de 2020
Attach the data
Feri
Feri el 13 de Feb. de 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.

Iniciar sesión para comentar.

Respuesta aceptada

John D'Errico
John D'Errico el 14 de Feb. de 2020
Editada: John D'Errico el 14 de Feb. de 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 comentarios
Feri
Feri el 15 de Feb. de 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 el 21 de Feb. de 2022
this helped a lot. thanks for sharing your functions!

Iniciar sesión para comentar.

Más respuestas (2)

darova
darova el 13 de Feb. de 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 comentarios
darova
darova el 14 de Feb. de 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 el 14 de Feb. de 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...

Iniciar sesión para comentar.


Jacob Wood
Jacob Wood el 14 de Feb. de 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 comentario
Feri
Feri el 14 de Feb. de 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.

Iniciar sesión para comentar.

Community Treasure Hunt

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

Start Hunting!

Translated by