Find centroid of curve

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.

https://in.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections

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?

Answer 1

[Bjorn Gustavsson]

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

[sparsh garg]

hey thanks it works on this one,but if I want to generalize this,what ammendments should I make.

For example:

for the following set

istart is 368 and iend is 357 so R{:,istart:iend] will give NULL

i tried reversing but then the centroid vector comes out to be NULL

Same goes for the following figures

Side note tried changing the start and end values to min and max respectively but in that case the centroid vector comes out to be NULL

[sparsh garg]

i tried running the debugger,as of now when start and end are changed it doesn't enter the while loop.inpolygon returns 0,any idea why that is happening and what steps i can take from my end to generalize it.The result is in the direction which I want,so that's why I would like to discuss this further as per your conveninece.

[Bjorn Gustavsson]

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]

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.

Answer 2

[Image Analyst]

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]

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]

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.