How can I calculate the thickness of a set of strokes in an image?

Question

How can I calculate the thickness of each stroke in the attached image?

I did a segmentation for all the nerves in the image and I labelled each object in the image. What I want is to know how to calculate the thickness of each nerve? Also, how can I highlight a specific part of the nerve if there is an inflation in any nerve?

Your help in this regard is highly appreciated.

Answer 1

My approach solves this problem using the following

Idea:

For each single stroke in the picture:

1. Compute area of stroke
2. Compute length of stroke
3. Compute average thickness as area divided by stroke length

Area:

Computing the area is easy. One could simply count the pixels.

Stroke length:

The stroke length is more difficult. One may be able to use the skeleton of the image bwmorph(BW,'skel',Inf) to do this, I however use a different approach. I find the perimeter of the stroke and use some math to find a rectangle of equivalent area and perimeter, the sides of which should correspond to stroke length and average thickness.

Computing the perimeter of such a stroke, one has to pay close attention to the coastline paradox. Measuring the actual perimeter of the stroke defined by the outline of the pixels will give inaccurate results: The perimeter of a diagonal stroke would be measured as 1 pixelwidth to the right, 1 pixelwidth up, 1 pixelwidth to the right, 1 pixelwidth up, ... instead of: sqrt(2) pixelwidths diagonally, ...

We can solve this problem by using alpha shapes to triangulate our stroke. I use the implementation by Jonas Lundgren on the File Exchange: alphavol.

Final code:

image = im2bw(imread('https://i.sstatic.net/DggWW.jpg'));    
alphaRadius = 2; %// Could be chosen larger than 2 to get better results, but must not affect the overall area of the triangulation.
%%
[Labeled,N] = bwlabel(image);
thicknesses = zeros(1,N);
lengths = @(Points) sqrt(sum(Points.^2,2));                           

%// Formula for finding an equivalent rectangle:
%// Assume rectangle: Area = a*b, Perimeter = 2*(a+b),
%// solve for a and b, return minimum as thickness
solveForThickness = @(P,A) min(P/4 + (P^2/4 - 4*A)^(1/2)/2, ...
                               P/4 - (P^2/4 - 4*A)^(1/2)/2);
for idx = 1:N
    %// Get current stroke
    stroke = (Labeled == idx);
    [Y,X] = find(stroke);
    %// Get corners of pixels and generate alpha shape
    P = [X-1, Y-1; ...
         X,   Y-1;
         X-1, Y;
         X,   Y];
    [area,S] = alphavol(P,alphaRadius);
    area = nnz(stroke); %// Seems to give better results than triangulated area above
    %// Compute perimeter and thickness
    perimeter = sum(lengths(P(S.bnd(:,1),:)-P(S.bnd(:,2),:)));
    thicknesses(idx) = solveForThickness(perimeter, area);
end

Aftermath:

Keep in mind that this solution will largely depend on how well the triangulation fits the stroke data. The problems of this algorithm due to the coastline paradox will be visible by choosing strokes of uniform varying steepness. There might be a better method to compute the length of the stroke to get more accurate results. It seems to work quite well though.

An alternative to using alpha shapes could be:

perimeter = getfield(regionprops(stroke,'perimeter'),'Perimeter');

The results of which are also quite accurate.

Answer 2

I would like to point you to the following post: Measuring the average thickness of traces in an image . Essentially, there is a single stroke in an image and the objective is to determine the average thickness of this stroke. We simply have to reapply this algorithm for each stroke that exists in your image.

The basic algorithm is the following (and very ingenious):

1. For a stroke in your image, apply a distance transform to your image where the image is transformed so that all background pixels become white and the object pixels become black. Hence, we are applying the distance transform on the inverse of the image.

   The distance transform finds the distance between a point in your image to the closest non-zero pixel. As such, by applying the distance transform to the inverse, for each stroke point, this will determine the distance from that particular point to the closest boundary point on the stroke. The reason why we need the distance transform in this aspect will be made clear in the next step.

2. The points on the stroke that share the greatest distance in the distance transform define the middle of the stroke. You can think of this as the half-way point in between the stroke. Collect all of those distance values within some tolerance of this greatest distance due to floating point arithmetic.

3. Find the average or median of all of these said distances and this determines the width of half the stroke.... so if you want the full thickness, simply multiply the result by 2.

Before I continue, I'm going to assume you have access to the Image Processing Toolbox, as the distance transform is implemented in the bwdist function that is part of the toolbox. If you don't have this toolbox, then unfortunately what I have written will not work. This is the only dependency that the algorithm needs (other than using imread and im2bw of course).

I had to convert your image to binary as you uploaded the image in RGB format for some reason, so I made it into a single channel binary image. As such, your code could look something like this:

%// Read the image from StackOverflow and convert to binary
im = im2bw(imread('https://i.sstatic.net/DggWW.jpg'));

%// Find all unique contours and assign them a unique ID.  Also determine how many there are
[lbl,N] = bwlabel(im);

%// Initialize an array that determines the thickness of each stroke
thickness = zeros(N,1);

%// Tolerance to collect distances
tol = 1;

%// For each stroke...
for idx = 1 : N

    %// Step #1
    cntr = lbl == idx; %// Get a mask that segments out that stroke only
    d = bwdist(~cntr); %// Find distance transform on inverse

    %// Step #2
    max_dist = max(d(:)); %// Get the maximum distance from stroke to edge
    dists = d(abs(d - max_dist) <= tol); %// Collect those distances within a tolerance

    %// Step #3
    thickness(idx) = 2*mean(dists);  %// Find the average width and multiply by 2 for full width
end

thickness is an array of elements where it tells you the width or thickness of each stroke in your image. Each element i tells you the thickness of the contour specified by label i from the bwlabel result. Therefore, for each element i in thickness, this records the thickness for the stroke delineated at:

cntr = lbl == i;

With your image, I get the following result:

thickness =

    3.3555
    3.2494
    3.1545
    3.1773
    3.3282
    3.2607
    3.2710
    3.2772
    3.2948
    3.1607

As for determining whether or not nerves are inflamed, you would need to know what the ground truth thickness is for each stroke and determine if there is an increase. I don't have such information, so I'm going to leave that to you as an exercise. We've calculated the thickness of each stroke so from there, write some code that detects an increase and act accordingly.

---

As an additional bonus, let's make a new output image where we find the centroid of each stroke, place that stroke in this output image and print out its thickness at the centroid of the stroke. Something like this:

imshow(im); hold on;
for idx = 1 : N
    [y,x]= find(lbl == idx);
    cen = mean([x y]);
    text(cen(1), cen(2), ['T = ' num2str(thickness(idx))], 'color', 'red');
end

We get: