Fit curve edge and determine curvature

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]

Try function bwskel

data = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/918279/50um-508um-1.jpg');
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.');

Answer 1

[Simon Chan]

Use function polyfit with higher degree.

data = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/918279/50um-508um-1.jpg');
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
%   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
  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
%  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;
    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

Answer 2

[yanqi liu]

yes，sir，may be use thin can get image shrink effect，such as

warning off all
im = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/918279/50um-508um-1.jpg');
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)));
set(h,'AlphaData',0.5)
hold on; plot(xt, yt, 'r-', 'LineWidth',2);