Interpolating and resampling around closed curve to calculate centroid

I am working with bone STL data and slicing the bone to get an outline of the bone at that height. I need to find the centroid of the outline, but the points are not evenly spaced out around the curve. The points are also not plotted in order around the curve, which has made using a spline fit difficult. This is a picture of the points on the outline that I am working with. I have had a hard time finding a way to interpolate the entire closed loop. Let me know if you have any suggestions!

[Alan Stevens]

Post the data.

[Sofia Guarnieri]

This is the data for this outline

Answer 1

[David Goodmanson]

Hi Sofia,

Since your points are pretty closely spaced, there is not much need to interpolate the points. The method below picks a point in the center, sorts the points by angle about that point and does sums involving the triangles formed by two adjacent points on the perimeter and the point in the center. I used complex notation since it makes for shorter code.

I assume you want the centroid of all the enclosed mass (assumed to have constant density), not the centroid of the perimeter, which is different.

zdata = slice_proximal(:,1)+i*slice_proximal(:,2);        % x + i*y
fig(1)
plot(zdata,'o-')                       % unsorted points
z0 = mean(zdata);
zdata_centered = zdata - z0;
theta_ns = angle(zdata_centered);
r_ns = abs(zdata_centered);          
[theta ind] = sort(theta_ns);          % sort the points by angle
r = r_ns(ind);
theta = [theta; (theta(1)+2*pi)];      % close the loop
r = [r; r(1)];
z = r.*exp(i*theta);                   
fig(2)                                 % sorted points
plot(z,'o-')
dtheta = diff(theta);
rmid = (r(1:end-1) + r(2:end))/2;      % mid-segment values
thetamid = (theta(1:end-1) + theta(2:end))/2;  
A = sum((rmid.^2/2).*dtheta)                                   % the area
zcm_centered = sum(exp(i*thetamid).*rmid.^3.*dtheta)/(3*A);
% zcm_centered is relative to z0; add in z0
zcm = z0 + zcm_centered;
xcm = real(zcm)
ycm = imag(zcm) 

[Sofia Guarnieri]

Thank you so much for your help! This worked perfectly!