Fit curve edge and determine curvature

17 次查看（过去 30 天）
Lai Jiang
Lai Jiang 2022-3-8
回答： yanqi liu 2022-3-10
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!
  1 个评论
Simon Chan
Simon Chan 2022-3-8
Ran in:
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.');

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Simon Chan
Simon Chan 2022-3-8
Ran in:
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
yanqi liu
yanqi liu 2022-3-10
Ran in:
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);

类别

在 Help Center 和 File Exchange 中查找有关 Biomedical Imaging 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by