Cylinder fit to a point cloud data

40 次查看（过去 30 天）
Feri
Feri 2020-2-13
评论： Stefan Süssmaier 2022-2-21
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 个评论
显示 4更早的评论
Feri
Feri 2020-2-13
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2020-2-14
编辑：John D'Errico 2020-2-14
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 个评论
显示 4更早的评论
Feri
Feri 2020-2-15
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 2022-2-21
this helped a lot. thanks for sharing your functions!

请先登录，再进行评论。

更多回答（2 个）

darova
darova 2020-2-13
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 个评论
显示 10更早的评论
darova
darova 2020-2-14
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 2020-2-14
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...

请先登录，再进行评论。

Jacob Wood
Jacob Wood 2020-2-14
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 个评论
Feri
Feri 2020-2-14
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Multirate Signal Processing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by