How to calculate the curvature of a boundaries in binary images?

Dear all,

I have a binary image that contains one shape as I show below:

I want to represent the curvature using color coding of the boundariessimilar to image I got from the video in this link http://en.wikipedia.org/wiki/Curvature

Does anyone knows how to do that?

Any help or suggestion will be appreciated.

Meshoo

Answer 1

[Image Analyst]

Use bwboundaries to get a list of the (x,y) boundary coordinates. Run around that border taking a section of the curve, say 9 elements, and fit a polynomial to it. Knowing the derivative of the curve model you used, and the coefficients polyfit() gave you, you can get the curvature at each point. Then run around the list of (x,y) edge coordinates using a marker color that corresponds to the curvature. Try to program that up. If you can't do it all , then post what you have and I'll try to help with the rest.

[Meshooo]

OK. I will try to follow your suggestion and will post here my result. Indeed thank you very much.

[Image Analyst]

OK, I spent about 20 minutes putting together this for you. It gets you close but there are still somethings you have to do and I'll let you figure those out yourself. See attached test.m file.

[Meshooo]

Thank you very much. I am also working on it. I was wondering why two different color codes were generated in the regions I marked with RED and PINK although they have same curvature convex direction?

[Image Analyst]

You'd have to use the debugger to find the curvatures there and the colors assigned to them.

[Meshooo]

OK, I will try. Thank you.

[Meshooo]

I was thinking to use your code in different approach.

I = imread('conv.tif');  %same input image
binaryImage = im2bw(I);
%figure, imshow(binaryImage, []);
%
boundaries = bwboundaries(binaryImage);
%
x = boundaries{1}(:, 2);
y = boundaries{1}(:, 1);
windowSize = 19;
halfWidth = floor(windowSize/2);   %half window size
curvatures = zeros(size(x));
%
%k =10
for k = halfWidth+1 : length(x) - halfWidth
  theseX = x(k-halfWidth:k+halfWidth);   % 19 element of x
  theseY = y(k-halfWidth:k+halfWidth);   % 19 elemebt of y
  % Get a fit.
  coefficients = polyfit(theseX, theseY, 2);
  % Get the curvature
  curvatures(k) = coefficients(1);
  xc(k) = x(k);
  yc(k) = y(k);  

end

% Get rid of ridiculous curvatures (straight line segments).
curvatures(abs(curvatures) > 20) = 0;

Now reconstruct the original image from the boundaries using bound2im from (DIPUM)

bound = bwboundaries(binaryImage);
bound = bound{1};
[M N] = size(binaryImage);

B = bound2im(bound, M, N, min(bound(:, 1)), min(bound(:, 2)));

B_gray = uint8(B);  % convert the image to 8 bit intensity
[L V] = size(bound);

pix_Value = []

for k = 1:L
Value = B_gray(y(k), x(k)) % Y then X values
pix_Value(k,1) = Value
end

curv10 = curvatures*1000;  % increase the intensity 1000 times

%but the values of the contour back to the image
for k = 1:L
B_gray(y(k), x(k)) = curv10(k);
end

figure, imshow(B_gray, []) % this will show the image

Now how to color code properly this intensity image?

[Meshooo]

Tried this but not working well

imagesc(B_gray)
axis square
colormap('hot')

[Image Analyst]

No, that definitely won't do it. colormap() just affects the image itself, not the color of any graphics, like lines or markers, that are overlaid onto the image. First see if the curvatures that are computed seem reasonable. Then, if they are, make up a color map that is reasonable. The colormap may be non-linear since it looks like many of the curvatures are the same but there's some huge spikes that are outliers.

[Meshooo]

OK, thank you. Still couldn't figure out what was wrong with the curvature...

[Meshooo]

Any hints?

[Image Analyst]

First see if the curvatures that are computed seem reasonable. If not, you might try widening the window to include more points. Then, if the curvatures are reasonable, make up a color map that is reasonable.

[Meshooo]

I tried too much but it seems there is something wrong with the method. I tried to test a simple shape that has one curvature and changed the window size as well but no curvature was detected

result is:

[Meshooo]

However, I tried to measure the curvature based on the distance from the centroid and I got the following result:

But I think this is different from the real curvature. Right?

My code is below:

% this software calculates the curvature of a colsed shape based on the
% distance of contour pixels from the centroid using signature

A = imread('curv_circle.tif');
A = im2bw(A);
%
binaryImage = A;
%
s  = regionprops(A, 'centroid');
centroids = cat(1, s.Centroid);
imshow(A)
hold on 
plot(centroids(:,1), centroids(:,2), 'r*')
hold off
X0 = centroids(:,1)
Y0 = centroids(:,2)
%*************************************************************************
% *********** Curvature measurement based on contour signature ***********
%*************************************************************************
%I = imread('curv_circle2.tif');
%
%binaryImage = im2bw(I);
%figure, imshow(binaryImage)
%
bSq = bwboundaries(binaryImage, 'noholes');
%[curvatures, angleSq] = signature(bSq{1}, X0, Y0);
[curvatures, angleSq] = signature(bSq{1})
figure
plot(angleSq, curvatures)
axis ([0 500 0 600]);
grid on;
boundaries = bwboundaries(binaryImage);
%
x = boundaries{1}(:, 2);
y = boundaries{1}(:, 1);
bound = bwboundaries(binaryImage);
bound = bound{1};
[M N] = size(binaryImage);
B = bound2im(bound, M, N, min(bound(:, 1)), min(bound(:, 2)));
B_gray = uint8(B);  % convert the image to 8 bit intensity
[L V] = size(bound);
pix_Value = []
for k = 1:L
Value = B_gray(y(k), x(k)) % Y then X values
pix_Value(k,1) = Value
end
%%adjust the size of curvatures
[C R] = size(pix_Value);
curvatures = imresize(curvatures, [C 1])
%but the values of the contour back to the image
for k = 1:L
B_gray(y(k), x(k)) = curvatures(k);
end
figure, imshow(B_gray, []) % this will show the image
%
SE = strel('square', 5);
%
Dil = imdilate(B_gray, SE);
%
RGB = label2rgb(Dil, 'cool');
figure, imshow(RGB)

[Meshooo]

If I use the jet color coding I get this:

[Meshooo]

I also tried to correct your code by adding small part:

I = imread('curv_circle.tif');
binaryImage = im2bw(I);
%figure, imshow(binaryImage, []);
%
boundaries = bwboundaries(binaryImage);
%
x = boundaries{1}(:, 2);
y = boundaries{1}(:, 1);
windowSize = 10;
halfWidth = floor(windowSize/2);   %half window size
curvatures = zeros(size(x));
%
%k =10
for k = halfWidth+1 : length(x) - halfWidth
  theseX = x(k-halfWidth:k+halfWidth);   % 19 element of x
  theseY = y(k-halfWidth:k+halfWidth);   % 19 elemebt of y
  % Get a fit.
  coefficients = polyfit(theseX, theseY, 2);
  % Get the curvature
  curvatures(k) = coefficients(1);
  xc(k) = x(k);
  yc(k) = y(k);  

end
clc

% Get rid of ridiculous curvatures (straight line segments).
curvatures(abs(curvatures) > 20) = 0;

%curvatures = flip(curvatures)
%
plot(x, curvatures, 'b-');
hold on 
plot(y, curvatures, 'r-');

min_Cur = min(curvatures)
%
max_Cur = max(curvatures)
[L W] = size(curvatures);
%
CurCur = curvatures;
for i = 1 :L
    for j = 1 : W
        if CurCur (i,j) < 0  %negative values
           CurCur (i,j) =  CurCur (i,j) +(abs(min_Cur))
        elseif CurCur (i,j) > 0  %positive values
            CurCur (i,j) =  CurCur (i,j) + 100
        end
    end
end

curvatures = CurCur;

Make up a colormap

minC = min(curvatures)
maxC = max(curvatures)
range = ceil(maxC - minC)
myColorMap = jet(range);
% Display the image again.
imshow(binaryImage);
hold on;
%
for k = halfWidth+1 : length(x) - halfWidth
  % Get the index in the color map.
  thisIndex = round(size(myColorMap, 1) * (curvatures(k) - minC) / range)
  
    fprintf('For point #%d, the colormap index is %d\n', k, thisIndex);
  if thisIndex <= 0
    thisIndex = 1;
  end
  if isnan(thisIndex)
    thisIndex = 1;
  end
  % Extract out the RGB triplet for this particular row in the color map.
  thisColor = myColorMap(thisIndex, :);
  plot(x(k), y(k), '.', 'MarkerSize', 25, 'Color', thisColor);
end

the result is:

[sana3 sal]

Hello Mr @Image Analyst, I am wondering how to read the curvature? or how to analyse it?

[Image Analyst]

A large radius of curvature means a "flat" part of the curve while a small radius of curvature means "pointy" parts of the curve.

Whether it's positive or negative just says which side of the curve it's bending to. You might want to just take the absolute value of the curvature if you don't care which side it bends toward.

[Francesco Pignatelli]

@Image Analyst I have followed this post and I saw your comments and you mentioned that "Knowing the derivative of the curve model you used, and the coefficients polyfit() gave you, you can get the curvature at each point.". My question is, how is it possible to get the derivatives?

At the current stage, I am at:

boundaries = bwboundaries(binaryImage);

x = boundaries{1}(:, 2);

y = boundaries{1}(:, 1);

So I have x and y now.

[Image Analyst]

@Francesco Pignatelli see attached demo that does what I said. Adapt as needed. If you need more help, start your own discussion thread and attach your (x,y) data.

Answer 2

[Preetham Manjunatha]

Here this function measures the curvature and other shape properties. In addition, it displays the curvature.

Here is the T-section with sharp corners:

Answer 3

[linda cummings [e-mail]]

Very nice. Have you worked in gray image the this code ? Have you tried?

Linda Cummings

cummingslinda14 at gmail.com

Answer 4

[sana3 sal]

Hello Mr @Image Analyst, I am wondering how to read the curvature? or how to analyse it?

[cyril grandjean]

Hello @Image Analyst

among the results generated by the script which of the results represents the average curvature or a score of general curvature for all the shape? if not generated did you know how to get it?

thank you in advance

[Image Analyst]

Did you see the code he gave for computing curvature in this link to his comment.

[cyril grandjean]

Hello I read the code well but I was not sure that I understood correctly. So the coefficient provided in this code corresponds to the overall curvature of my form or a trend is that right?

thank you in advance