Curvature estimation in contour extracted from gray image

Hi,

I have the edges/contours extacted from the gray image, now I intend to estimate the curvature.

For e.g., I have attached my sample contour here..

Could someone help me with this out ??

Answer 1

[KSSV]

You may use this file exchange function:

https://in.mathworks.com/matlabcentral/fileexchange/32696-2d-line-curvature-and-normals

Check the code:

load('Test image.mat') ;
[y,x] = find(I) ; 
 
idx = boundary(x,y) ; 
x = x(idx) ; y = y(idx) ;

Vertices = [x y] ; 
k=LineCurvature2D(Vertices)  ; 
N=LineNormals2D(Vertices) ;

figure(1)
scatter(x,y,20,k,'filled')
colorbar
title('curvature')

figure(2)
plot(x,y)
hold on
plot([Vertices(:,1) Vertices(:,1)+10*N(:,1)]',[Vertices(:,2) Vertices(:,2)+10*N(:,2)]');
title('Normals')

[Turbulence Analysis]

Hi,

Thanks for the info.

I do tried, but while exceuting the LineCurvature2D(my image);I am just getting all values as NaN...

[KSSV]

Did you check the input? Input has any NaN's?

[Turbulence Analysis]

I checked, no NaN's in input...

[KSSV]

Attach your contour coordinates..

[Turbulence Analysis]

Here, it is

[Turbulence Analysis]

I also tried with the solution provided in one of the previous threads https://www.mathworks.com/matlabcentral/answers/164349-how-to-calculate-the-curvature-of-a-boundaries-in-binary-images#comment_253244

% Read in a standard MATLAB color demo image.
folder = 'D:\curvature';
baseFileName = 'Test image.bmp';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
	% Didn't find it there.  Check the search path for it.
	fullFileName = baseFileName; % No path this time.
	if ~exist(fullFileName, 'file')
		% Still didn't find it.  Alert user.
		errorMessage = sprintf('Error: %s does not exist.', fullFileName);
		uiwait(warndlg(errorMessage));
		return;
	end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
	% It's not really gray scale like we expected - it's color.
	% Convert it to gray scale by taking only the green channel.
	grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the original gray scale image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 

% Get the boundaries, but first you need to binarieze the image into a logical image.
binaryImage = grayImage > 250;
binaryImage =~binaryImage ;
% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 2, 2);
set(gca, 'ydir', 'reverse');
title('Outlines, from bwboundaries()', 'FontSize', fontSize);
hold on;
boundaries = bwboundaries(binaryImage);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
	thisBoundary = boundaries{k};
	plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
end
hold off;

% Get the curvature at each point of the first curve.
x = boundaries{1}(:, 2);
y = boundaries{1}(:, 1);
windowSize = 19;
halfWidth = floor(windowSize/2);
curvatures = zeros(size(x));
for k = halfWidth+1 : length(x) - halfWidth
	theseX = x(k-halfWidth:k+halfWidth);
	theseY = y(k-halfWidth:k+halfWidth);
	% 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;
subplot(2, 2, 3);
plot(x, curvatures, 'b-');
hold on 
plot(y, curvatures, 'r-');
title('Curvatures', 'FontSize', fontSize);

% Make up a colormap
minC = min(curvatures)
maxC = max(curvatures)
range = ceil(maxC - minC)
myColorMap = jet(range);
% Display the image again.
subplot(2,2,4);
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
title('Curvatures Indicated by Colors', 'FontSize', fontSize);

However, I am not getting the desired curvature... I am expecting the curvature value should be maximum in the region shown with black box in Desired.bmp....

Thanks in advance !!

[KSSV]

It should work with the mentioned file exchange functions. Edited the answer, check it.

[Turbulence Analysis]

Hi,

May thanks for the sugegstion, only concern is, the input image and obtained curvature image looks different.. Please have at the attachments..

load('Test image.mat') ;
imshow(I);
[y,x] = find(I) ; 
idx = boundary(x,y) ; 
x = x(idx) ; y = y(idx) ;

Vertices = [x y] ; 
k=LineCurvature2D(Vertices)  ; 
N=LineNormals2D(Vertices) ;

figure(1)
scatter(x,y,50,(k),'filled');
colorbar;
title('curvature');

figure(2)
plot(x,y)
hold on
plot([Vertices(:,1) Vertices(:,1)+10*N(:,1)]',[Vertices(:,2) Vertices(:,2)+10*N(:,2)]');
title('Normals')

[KSSV]

You can check like this:

imshow(I)
hold on
plot(x,y,'r')

The obtained points x, y will not form a closed contour, so using boundary they are arranged. So boundary misses some points. If you can arrange the points (x,y), you can get exactly the same curve shown in image.

[Turbulence Analysis]

Thanks.. Now,I used, axis equal command, Now It's better.

But, I ddidn't understand why the curvature scatter plot always dipslays in inverted ??

scatter(x,y,50,k,'filled');
axis equal

[Turbulence Analysis]

While executing, N=LineNormals2D(Vertices); it genartes two columns , Is it normal vectors in x, y.. ??