Splitting a skeleton in a binary object to equal lengths!

Hello everyone,

I have a skeleton image obtained from running bwskel on a binary image. Here's couple of images as a sample:

I want find to specific number of points (for example 10 points) on these skeletons with equal distances from each other. Basically like splitting the skeleton to 10 equal parts.

I tried using bwconncomp to find the connected components in the the image and their pixel locations. But the locations are not in an order that I could use to split these to equal distances.

I would appreciate any kinds of suggestions!

Thank You

Answer 1

[Walter Roberson]

See https://www.mathworks.com/matlabcentral/fileexchange/34874-interparc "A common request is to interpolate a set of points at fixed distances along some curve in space (2 or more dimensions.) The user typically has a set of points along a curve, some of which are closely spaced, others not so close, and they wish to create a new set which is uniformly spaced along the same curve."

[Farshad Bolouri]

Hello Walter!

Thank You for your answer, I have actually tried this function before. I was wondering if this is possible without interpolation.

But even with interparc() I think I am not understanding something correctly. When the points are sorted with respect to 'x' this is what I get:

skel = bwskel(worm); %worm is the original image
[y, x] = find(skel);
pt = interparc(100,x,y,'spline');
plot(x,y,'r*',pt(:,1),pt(:,2),'b-o')

And when I sort with respect to 'y' values:

skel = bwskel(worm); %worm is the original image
[y, x] = find(skel);
[sortedY, sortIndex] = sort(y)
sortedX = x(sortIndex)
pt = interparc(100,sortedX,sortedY,'spline');
plot(x,y,'r*',pt(:,1),pt(:,2),'b-o')

It seems like that the order of 'x' and 'y' values matter. Could you help me on how to fix this problem?

Thank You

[Image Analyst]

You need to get the worm sorted. Using find() does not do that. It has no concept of what pixel is next to any other pixel. It just goes down rows and moves over column-by-column from left to right. You need to use regionprops(), which I think will sort it right, or at least better, for you. Try

props = regionprops(mask, 'PixelList', 'PixelIdxList');

props.PixelList is the list of (x,y) coordinates, hopefully going from one end to the other and not going column by column which is what find() does.

props.PixelIdxList is the linear Index of the pixel in the original mask image.

Attach 'skel' in a .mat file if you need more help.

save('answers.mat', 'skel');

[Walter Roberson]

https://www.mathworks.com/help/images/ref/bwtraceboundary.html

[Farshad Bolouri]

Thank you @Image Analyst and @Walter Roberson for quick responses!

The bwtraceboundary() function worked perfectly in this case! Here are the results for interpolation:

Again, thank you!

[Walter Roberson]

Did you end up using interparc() or did you just take every N'th point from the boundary trace?

[Farshad Bolouri]

Well the results that I was showing were from interparc()!

But here are the results without using interparc() and just by using N'th point.

As you can see they basically look identical which is what I was interested in. I was trying to eliminate the step of interpolation for efficiency.

Now I am trying to see if I want to show the normals with respect to these points (blue circles), can I still work with the N'th point approach or do I need to interpolate.

This is what I currently have for the normal approximation:

skel_L = length(find(skel));

[endpointsX, endpointsY] = find( bwmorph(skel == 1, 'endpoints') );
trace = bwtraceboundary(skel, [endpointsX(1), endpointsY(1)],'N');
y = trace(1:skel_L,1);
x = trace(1:skel_L,2);

imshow(skel)

hold on;
norms = zeros(length(4: floor(skel_L/15) :skel_L-3),2);
i = 1;

for k=4: floor(skel_L/15) :skel_L-3
%plot(x(k),y(k),'bo')
% step 1 dv/dt
dx = (mean(x(k+1:k+3))-mean(x(k-3:k-1)));
dy =(mean(y(k+1:k+3))-mean(y(k-3:k-1)));

dvdt = [dx;dy];

% step2 rotate 90%
dvdtRot = [-dy ;dx];

% Step 3: Scale it to magnitude 1
% unit vector
dvdtRotNorm = dvdtRot/norm(dvdtRot);

scale=50;
pNorm = [x(k);y(k)]+scale*dvdtRotNorm;
%plot(pNorm(1),pNorm(2),'gs')
%line([x(k);pNorm(1)],[y(k);pNorm(2)])
norms(i,:) = pNorm;
i = i + 1;
end
k=4: floor(skel_L/15) :length(x)-3;
quiver(x(k),y(k), abs(x(k) - norms(:,1))*3, abs(y(k) - norms(:,2))*3, 'r' ,...
    'LineWidth', 3, 'MaxHeadSize', 3)

I am not sure if this is a good approach to this task or not. In my opinion the normals don't look great but don't look bad either.

Do you have any suggestion on how I should approach this to possibly get better results?

PS: I attached the 2 mat files to this comment which are the two skeletons shown above.

Thank You

[Jhe Mag]

Hello @Farshad Bolouri can I please ask for the code you used here?

I tried to run your code in the recent comment but it gives me error such as:

Index in position 1 exceeds array bounds (must not exceed 3).

y = trace(1:skel_L,1);

Hoping for help