Cylinder fit to a point cloud data
24 views (last 30 days)
Show older comments
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
Accepted Answer
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/271588/image.jpeg)
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
More Answers (2)
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
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')
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
See Also
Categories
Find more on Multirate Signal Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!