# Fit curve edge and determine curvature

50 views (last 30 days)
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: https://www.mathworks.com/matlabcentral/fileexchange/69452-curvature-of-a-1d-curve-in-a-2d-or-3d-space?s_tid=prof_contriblnk
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!
Simon Chan on 8 Mar 2022
Try function bwskel
BW1 = imfill(rgb2gray(data)>50,'holes');
SE = strel('disk',10);
BW2 = imclose(BW1,SE);
BW3 = bwskel(BW2);
[r,c]=find(BW3);
Vertices = [c,r];
plot(Vertices(:,1),Vertices(:,2),'r.'); Simon Chan on 8 Mar 2022
Use function polyfit with higher degree.
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
%
figure(1)
ax=gca;
imagesc(ax,pagetranspose(data));
colormap(ax,'gray');
axis(ax,'image');
ax.YDir = 'normal';
hold(ax,'on')
plot(ax,x,y,'b:');
hold on
plot(ax,x1,y1,'r--') %
figure(2)
plot(k1);
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
% 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
end
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,:));
end
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,:));
end
i = N;
L(i) = L(i-1)+norm(X(i,:)-X(i-1,:));
if dims == 2
k = k(:,1:2);
end
end
%%
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
% 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;
end
return
end
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;
else
k = G'/R^2; % Curvature vector
end
end

yanqi liu on 10 Mar 2022
yes，sir，may be use thin can get image shrink effect，such as
warning off all
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); 