# Find centroid of curve

10 views (last 30 days)
sparsh garg on 8 Sep 2021
Commented: sparsh garg on 9 Sep 2021
Given this curve,I would like to find the centroid of the curve.
Currently I am selecting the piece wise splines and using the Inter X module to find out where the projection of one spline will intersect the other spline.
However the issue with this method is that some points are not matched properly.
Therefore is there a way to find the centroid of the following set of spline points.
Result figure
As you can see at some points interX returns Null I tried doing it from right to left,but that gives indice exceeds error.
Also given the set of points can I skeletonize them?

Bjorn Gustavsson on 8 Sep 2021
Something like this ought to work:
[~,iEnd] = max(spline_pts(1,:));
[~,iStart] = max(spline_pts(2,:));
R1 = spline_pts(:,iStart:iEnd);
R2 = fliplr([spline_pts(1,[iEnd:end,1:iStart]);spline_pts(2,[iEnd:end,1:iStart])]);
xC = spline_pts(1,:);
yC = spline_pts(2,:);
r0 = R1(:,1)';
phi0 = -pi/2;
h = 0.5;
r_next = @(phi,r0,h) r0 + h*[cos(phi),sin(phi)]; % function for dr in x-y plane
r_n = r_next(phi0,r0,h);
r_N = [];
clf
plot(R1(1,:),R1(2,:),'r','linewidth',2)
hold on
plot(R2(1,:),R2(2,:),'g','linewidth',2)
plot(R1(1,1),R1(2,1),'r.','markersize',15)
plot(R2(1,1),R2(2,1),'g.','markersize',12)
plot(r0(1),r0(2),'k.')
OPS = optimoptions('fsolve');
OPS.Display = 'off';
while inpolygon(r_n(1),r_n(2),xC,yC) % while the next step is inside
% look for the h-long step with the same shortest distance to R1 and R2
phi_n = fsolve(@(phi) min(distance2curveBG(R1',r_next(phi,r0,h))) - ...
min(distance2curveBG(R2',r_next(phi,r0,h))),...
phi0,OPS);
% calculate
r_n = r_next(phi_n,r0,h);
r_N = [r_N;r_n]; % yeah, yeah, dynamically appending etc...
r0 = r_n; % move r0 one step forward
phi0 = phi_n;
plot(r_n(1),r_n(2),'k.')
drawnow
end
This worked on a very simple test-case I tried, and for your curve. Here I used distance2curve from the file exchange, but in order to use it here I had to swap the order of the first and second output argument of that function, hence the name-change.
HTH
Bjorn Gustavsson on 8 Sep 2021
The choise of the starting-point and initial step-direction was simply to find "the natural end-point" of the 'h' character at the "red +". That just conveniently happened to be the max y-value of the perimeter, and the initial step ought to be perpendicular to the curve at that point (I think it mathematically has to, but are too lazy to prove this) - therefore it should be in the downward direction at -pi/2. So for your "M" character you will have to pick some point on either left or right foot and set a phi0 corresponding to a stepping-direction perpendicular to the perimeter there. However, I'm not entirely sure this algorithm will be general enough to handle the branching of "M" and "R" seems to be an even tricker problem. For that case you might need a completely different approach, I recall Steve Eddins making blog-posts about this, what I find in a hurry are his series on detecting the US continental divide (continental-divide-1-intro) and river-drainage basins (upslope-area). Exactly how complicated curves you will be determined by how complicated characters you need to simplify...
HTH
Bjorn Gustavsson on 8 Sep 2021
If you check some combinations of where to put the end-points of the cuts to split the curve in two and build up the two sides, and then make sure that the first initial step points into the closed curve by chosing a proper phi0 then the method works as well as can be expected for all characters - however the "R" will for obvious reasons be "challenging" for this method. For that some ridge-detection or skelletoning seems more likely to give you the solution you want.

Image Analyst on 8 Sep 2021
Maybe you can convert the points to a digital image and then skeletonize it. Something like
binaryImage = poly2mask(x, y, rows, columns);
skelImage = bwskel(binaryImage);
% Get rows and columns of the skeleton
[rows, columns] = find(skelImage);
Requires the Image Processing Toolbox.
sparsh garg on 9 Sep 2021
Edited: sparsh garg on 9 Sep 2021
Ok so let's say that my current set of points is a 2XN matrix,when I do polymask(M(1,:),M(2,:),2,N);
the result of this is a 2XN logical where every entry is set to 0.
Secondly when I do [rows,columns]=find(skelImage) for the below img
so first of all as expected it drew all the points here except the skeletal points ,secondly why is the figure coming out like this
let's say skeletal image is BW_Thinned
and I do [r,c]=find(~BW_Thinned)
then plot(r(1,:),c(1,:),'or');
BW_Thinned is 88*65 double r and c is 5524*1
Ok one part got taken care I was able to get the points corresponding to R by doing [r,c]=find(BW_Thinned)
but the plot is still coming out like this
plus the line plot is unordered.
sparsh garg on 9 Sep 2021
I was able to order them,but the ordered plot is flipped by 360 degree
R_skeleton is BW_Thinned
R_exact_skeleton is the points shown in the previous comment
and ordred_R_skeleton is the plot in this comment.
I would appreciate it if we can discuss this as per your conveninece.