Fit curve edge and determine curvature

50 views (last 30 days)
Lai Jiang
Lai Jiang on 8 Mar 2022
Answered: yanqi liu on 10 Mar 2022
Dear Community,
I recently working on a project which involves determining the curvature of some bent thin foils. I collected the bent side edge of the foil using optical microscope like this attached.
I did some processings like binarizing, edge detection, dilate and smoothing and then was able to obtain something like this with the central line overlayed:
My question is what will be a good way to curve fit the central line and smooth it to be ready for the curvature determination. I found a curvature tool that may work but it is very sensitive to the noise:
Noted that the way I find the central line is by scanning line by line to find two pixel points from each sides, that means it is not complete symmetric regarding the curve shape as the bent foil edge is not align perfectly horizontal or vertical. Here is how it looks like with central line dilated and overlap on the original images:
I think there must be some better way to do this for curve fitting and curvature determination. Looking forward to discuss and learn from the community regarding this. Thanks in advance!
  1 Comment
Simon Chan
Simon Chan on 8 Mar 2022
Try function bwskel
data = imread('');
BW1 = imfill(rgb2gray(data)>50,'holes');
SE = strel('disk',10);
BW2 = imclose(BW1,SE);
BW3 = bwskel(BW2);
Vertices = [c,r];

Sign in to comment.

Answers (2)

Simon Chan
Simon Chan on 8 Mar 2022
Use function polyfit with higher degree.
data = imread('');
BW1 = imfill(rgb2gray(data)>50,'holes');
SE = strel('disk',10);
BW2 = imclose(BW1,SE);
BW3 = transpose(bwskel(BW2));
[y,x] = find(BW3);
[p,S,mu] = polyfit(x,y,8);
x1 = min(x):1:max(x);
[y1,delta] = polyval(p,x1,S,mu);
[L,R,k1] = curvature([x1' y1']); % Use function curvature
ax.YDir = 'normal';
hold on
grid on
function [L,R,k] = curvature(X)
% Radius of curvature and curvature vector for 2D or 3D curve
% [L,R,k] = curvature(X)
% X: 2 or 3 column array of x, y (and possibly z) coordiates
% L: Cumulative arc length
% R: Radius of curvature
% k: Curvature vector
% The scalar curvature value is 1./R
% Version 2.6: Calculates end point values for closed curve
N = size(X,1);
dims = size(X,2);
if dims == 2
X = [X,zeros(N,1)]; % Use 3D expressions for 2D as well
L = zeros(N,1);
R = NaN(N,1);
k = NaN(N,3);
for i = 2:N-1
[R(i),~,k(i,:)] = circumcenter(X(i,:)',X(i-1,:)',X(i+1,:)');
L(i) = L(i-1)+norm(X(i,:)-X(i-1,:));
if norm(X(1,:)-X(end,:)) < 1e-10 % Closed curve.
[R(1),~,k(1,:)] = circumcenter(X(end-1,:)',X(1,:)',X(2,:)');
R(end) = R(1);
k(end,:) = k(1,:);
L(end) = L(end-1) + norm(X(end,:)-X(end-1,:));
i = N;
L(i) = L(i-1)+norm(X(i,:)-X(i-1,:));
if dims == 2
k = k(:,1:2);
function [R,M,k] = circumcenter(A,B,C)
% Center and radius of the circumscribed circle for the triangle ABC
% A,B,C 3D coordinate vectors for the triangle corners
% R Radius
% M 3D coordinate vector for the center
% k Vector of length 1/R in the direction from A towards M
% (Curvature vector)
D = cross(B-A,C-A);
b = norm(A-C);
c = norm(A-B);
if nargout == 1
a = norm(B-C); % slightly faster if only R is required
R = a*b*c/2/norm(D);
if norm(D) == 0
R = Inf;
E = cross(D,B-A);
F = cross(D,C-A);
G = (b^2*E-c^2*F)/norm(D)^2/2;
M = A + G;
R = norm(G); % Radius of curvature
if R == 0
k = G;
elseif norm(D) == 0
R = Inf;
k = D;
k = G'/R^2; % Curvature vector

yanqi liu
yanqi liu on 10 Mar 2022
yes,sir,may be use thin can get image shrink effect,such as
warning off all
im = imread('');
bw = imclose(im2bw(rgb2gray(im), 0.2), strel('disk', 9));
bw2 = bwmorph(bw, 'thin', inf);
[y,x] = find(bw2);
p = polyfit(y,x,9);
yt = linspace(min(y), max(y), 1e3);
xt = polyval(p, yt);
figure; imshow(im); hold on;
h = imshow(label2rgb(bwlabel(bw)));
hold on; plot(xt, yt, 'r-', 'LineWidth',2);

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!