find radius of xyz data

@Image Analyst

I have a 3d point cloud that contains xyz data that I would like to find the radius of. I would like to use inclircle function to find the diameter that contains all points but it is not giving me my desired output. I dont know if it has to do with the fact that my point cloud is tilted.

Alternatively, I was thinking to find a way to find the length of the axis that would give me the diameter in this example. Thank you

https://www.mathworks.com/matlabcentral/answers/413107-i-have-3d-points-and-i-should-determine-a-method-to-produce-and-visualize-a-center-line-down-the-c

I included zip containing my data

Code:

% Read data.
cloud = pcread('data.ply');
xyz = double(cloud.Location);
xyz(:,3) = [];
%k = convhull(xyz);
x = xyz(:,1);
y = xyz(:,2);
[center,radius] = minboundcircle(x,y)
theta = linspace(0,2*pi,100);
xc = center(1) + radius*cos(theta);
yc = center(2) + radius*sin(theta);
% Sometimes the minimal radius circle may pass
% through only two points.
plot(x,y,'ro',xc,yc,'b-')
set(gcf,'color','white')
% Create ylabel
ylabel({'Y'});
% Create xlabel
xlabel({'X'});
% Create title
title({'Minimum radius enclosing circle'},'Editing','on');
figure 
pcshow(cloud);

[jordan Ashton]

@Image Analyst

Answer 1

[Image Analyst]

Try this:

% Create a set of points in 3-D
xyz = rand(9, 3);
% Find distances between all of them with pdist2() in the stats toolbox.
distances = pdist2(xyz, xyz)
maxDistance = max(distances(:))
[index1, index2] = find(distances == maxDistance)
% Just take one singe they're symmetric.
index1 = index1(1);
index2 = index2(1);
fprintf('Furthest apart points:\n');
fprintf('(x1,y1,z1) = (%f, %f, %f)\n', xyz(index1, 1), xyz(index1, 2), xyz(index1, 3))
fprintf('(x2,y2,z2) = (%f, %f, %f)\n', xyz(index2, 1), xyz(index2, 2), xyz(index2, 3))
% Plot all of them.
plot3(xyz(:, 1), xyz(:, 2), xyz(:, 3), 'b.', 'MarkerSize', 20);
grid on;
xlabel('x');
ylabel('y');
zlabel('z');
% Plot the two that are farthest away and define the diameter.
hold on;
plot3([xyz(index1, 1), xyz(index2, 1)],...
    [xyz(index1, 2), xyz(index2, 2)],...
    [xyz(index1, 3), xyz(index2, 3)],...
    'r-', 'LineWidth', 3);

[Image Analyst]

On second thought, the two farthest away points in 3-D won't necessarily give you the diameter - imagine 3 equidistant points every 60 degrees along the perimeter.

I suppose if there are not too many points, you cold always brute force it with a triple for loop. John D'Errico (pronounced DEER - ih - coe I think, or maybe der-REE-coe) may know of a better way.

Answer 2

[Torsten]

It's an optimization problem:

min: r^2

s.c.

(x_i-x_s)^2 + (y_i-y_s)^2 + (z_i-z_s)^2 <= r^2 (i=1,...,number of points in 3d cloud)

where the unknowns are x_s, y_s, z_s and r and the triples (x_i,y_i,z_i) are your point cloud data.

You can try "fmincon" to solve.

Here is a test case:

xyz = randn(100,3);
X = xyz(:,1);
Y = xyz(:,2);
Z = xyz(:,3);
xyz0(1) = sum(X)/numel(X);
xyz0(2) = sum(Y)/numel(Y);
xyz0(3) = sum(Z)/numel(Z);
xyz0(4) = max(sqrt(((X-xyz(1)).^2+(Y-xyz(2)).^2+(Z-xyz(3)).^2)));
xyz0
%xyz0 = [1 1 1 100];
sol = fmincon(@fun,xyz0,[],[],[],[],[],[],@(p)nonlcon(p,X,Y,Z))
end

function obj = fun(p)
obj = p(4).^2
end

function [c,ceq] = nonlcon(p,X,Y,Z)
ceq = [];
c = (X-p(1)).^2 + (Y-p(2)).^2 + (Z-p(3)).^2 - p(4)^2;
end

Maybe the center of gravity x_s = 1/n *sum_x_i, y_s = 1/n *sum_y_i , z_s = 1/n * sum_z_i is the solution.

I'm not sure.

After the test runs I can say: it's not the center of gravity.

[Matt J]

That would be trivial to solve. It's just,

max_{i,j} ( norm([xi,yi,zi]-[xj,yj,zj]))

The radius that is really being sought is that of the tightest cylinder that will fit around the data, or in other words the tightest circle radius after projecting the (x,y,z) data onto a 2D plane. So, what you really need is to define the projection P_{theta,phi} onto a 2D plane that is angled with unknown parameters theta and phi. The constraint is then,

[ xp,yp]=P_{theta,phi} (x,y,z)

(xp_i-xp_s)^2 + (yp_i-yp_s)^2 <= r^2

[Torsten]

That would be trivial to solve. It's just,

max_{i,j} ( norm([xi,yi,zi]-[xj,yj,zj]))

Maybe you mean max_{i,j} ( norm([xi,yi,zi]-[xj,yj,zj]))/2, but it's not correct.

Consider an equal-sided triangle in 2d with sidelength 1.

The radius of the circumcircle is 1/sqrt(3) > 1/2.

The radius that is really being sought is that of the tightest cylinder that will fit around the data, or in other words the tightest circle radius after projecting the (x,y,z) data onto a 2D plane.

Is it so ? If this is the usual definition of the radius of a point cloud, I didn't know.

[Matt J]

Consider an equal-sided triangle in 2d with sidelength 1.

Yep. Never mind.

Answer 3

[Matt J]

I would like to use inclircle [sic.] function to find the diameter

Are you trying to enclose the points in a sphere or a cylinder? If a sphere, how would a function having to do with circles be applicable to you?

And where did you get the 3rd party code routine "incircle"? If you got it here, note that the package also contains a routine called minboundsphere() which might be more applicable.

[jordan Ashton]

Thank you for all your help @Image Analyst @Torsten @Matt J. I am trying to find the diameter of the point cloud. My though process was to set z to zero and find the minimum enclosing circle to get the diameter using this function linked, however I am not getting the result I am looking for. I think becuase my point cloud is at an angle so its giving a diameter that is too large.

https://www.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects?s_tid=srchtitle

I really liked the example by @Image Analyst here: https://www.mathworks.com/matlabcentral/answers/413107-i-have-3d-points-and-i-should-determine-a-method-to-produce-and-visualize-a-center-line-down-the-c

It seems like those two lines in the center (graph) could provide me my diameter. However, I am just not sure how to find the length of those lines

[jordan Ashton]

What I am asking specifically is how I can find the length of

plot3(xyz2Limit(:, 1) + xMean, xyz2Limit(:, 2) + yMean, xyz2Limit(:, 3) + zMean, 'k-', 'LineWidth', 4);

plot3(xyz3Limit(:, 1) + xMean, xyz3Limit(:, 2) + yMean, xyz3Limit(:, 3) + zMean, 'k-', 'LineWidth', 4);

from the example I linked

[Torsten]

But the formulae for all the parameters in the plot command are given in the code - you just have to copy the lines.