(来点有用的)MATLAB中轴变换求图像骨架
by HPC_ZY
主要介绍中轴变换中的迭代删除法,并与MATLAB自带函数进行对比。
MAT算法
中轴变换(Medial Axis Transfrom),对二值(所以你的图只有0和1两个值)图像进行操作,获取目标的中轴(也称作骨架)。
常见方法有:迭代删除法,大圆法,距离函数法,模拟退火法等。
今天分享迭代删除法的实现原理,以及简单的优化。
我理想中的结果
迭代删除算法
- 算法实现
遍历所有点,判断其八领域0,1分布情况,并按照一定规则进行删除。
关于删除规则可参考:https://blog.csdn.net/qq826309057/article/details/73008608
% 用户提供二值图像bw
bw = double(bw); % 转为非逻辑类型,方便后面求和
[M,N] = size(bw);
skeleton = bw;
deleCount = 1; % 用于记录一次迭代删除的像素个数
while deleCount>0
% 删(右,下)
deleCount = 0;
label = zeros(M,N,'logical');
for i = 2:M-1
for j = 2:N-1
if skeleton(i,j) % 如果该点为1,判断周围像素
p = [skeleton(i-1,j),skeleton(i-1,j+1),skeleton(i,j+1),skeleton(i+1,j+1),...
skeleton(i+1,j),skeleton(i+1,j-1),skeleton(i,j-1),skeleton(i-1,j-1)];
Np = sum(p);
Tp = length(find(([p(2:end),p(1)]-p) == 1));
if Np>1 && Np<7 && Tp==1 && p(1)*p(3)*p(5)==0 && p(7)*p(3)*p(5)==0
deleCount = deleCount+1;
label(i,j) = 1;
end
end
end
end
skeleton(label) = 0;
% 删(左,上)
label = zeros(M,N,'logical');
for i = 2:M-1
for j = 2:N-1
if skeleton(i,j) % 如果该点为1,判断周围像素
p = [skeleton(i-1,j),skeleton(i-1,j+1),skeleton(i,j+1),skeleton(i+1,j+1),...
skeleton(i+1,j),skeleton(i+1,j-1),skeleton(i,j-1),skeleton(i-1,j-1)];
Np = sum(p);
Tp = length(find(([p(2:end),p(1)]-p) == 1));
if Np>1 && Np<7 && Tp==1 && p(1)*p(3)*p(7)==0 && p(1)*p(5)*p(7)==0
deleCount = deleCount+1;
label(i,j) = 1;
end
end
end
end
skeleton(label) = 0;
end
- 算法优化
为避免每次迭代都要进行大量判断,而根据删除规则,对所有情况建立查找表。
关于查找表的建立可参考:https://www.cnblogs.com/xianglan/archive/2011/01/01/1923779.html
% 用户提供二值图像bw
bw = double(bw);
[M,N] = size(bw);
skeleton = bw;
% 生成查找表
List1 = zeros(256,1,'logical');
List2 = zeros(256,1,'logical');
for n = 0:255
p = bitget(n,1:8);
Np = sum(p);
Tp = length(find(([p(2:end),p(1)]-p) == 1));
if Np>1 && Np<7 && Tp==1
% 右,下
if p(1)*p(3)*p(5)==0 && p(7)*p(3)*p(5)==0
List1(n+1) = true;
end
% 左,上
if p(1)*p(3)*p(7)==0 && p(1)*p(5)*p(7)==0
List2(n+1) = true;
end
end
end
% 压缩搜索范围
[x,y] = find(bw);
a = min(x)+1;
b = max(x)-1;
c = min(y)+1;
d = max(y)-1;
% 开始删除
deleCount = 1;
mat = [1,2,4,8,16,32,64,128];
while deleCount>0
deleCount = 0;
% 删(右,下)
label1 = zeros(M,N,'logical');
for i = a:b
for j = c:d
if skeleton(i,j)
p = [skeleton(i-1,j),skeleton(i-1,j+1),skeleton(i,j+1),skeleton(i+1,j+1),...
skeleton(i+1,j),skeleton(i+1,j-1),skeleton(i,j-1),skeleton(i-1,j-1)];
idx = sum(p.*mat)+1;
if List1(idx)
label1(i,j) = true;
deleCount = deleCount+1;
end
end
end
end
skeleton(label1) = 0;
% 删(左,上)
label2 = zeros(M,N,'logical');
for i = a:b
for j = c:d
if skeleton(i,j)
p = [skeleton(i-1,j),skeleton(i-1,j+1),skeleton(i,j+1),skeleton(i+1,j+1),...
skeleton(i+1,j),skeleton(i+1,j-1),skeleton(i,j-1),skeleton(i-1,j-1)];
idx = sum(p.*mat)+1;
if List2(idx)
label2(i,j) = true;
deleCount = deleCount+1;
end
end
end
end
skeleton(label2) = 0;
end
- 算法简化
有人会问,为什么这个算法每次迭代要分成两次删除,一次删除(右,下)的点,一次删除(左,上)的点。如果一次性完成(上,下,左,右)点的删除是否可以。我在一开始也有这个疑问,于是做了测试。(我个人认为不能简化)
将上述算法简化为一次性删除后,代码如下:
% 用户提供二值图像bw
bw = double(bw);
[M,N] = size(bw);
skeleton = bw;
% 生成查找表
List = zeros(256,1,'logical');
for n = 0:255
p = bitget(n,1:8);
Np = sum(p);
Tp = length(find(([p(2:end),p(1)]-p) == 1));
if Np>1 && Np<7 && Tp==1
if (p(1)*p(3)*p(5)==0 && p(7)*p(3)*p(5)==0)||...
(p(1)*p(3)*p(7)==0 && p(1)*p(5)*p(7)==0)
List(n+1) = true;
end
end
end
% 收缩搜索范围
[x,y] = find(bw);
a = min(x);
b = max(x);
c = min(y);
d = max(y);
% 开始迭代
deleCount = 1;
mat = [1,2,4,8,16,32,64,128];
while deleCount>0
deleCount = 0;
label = zeros(M,N,'logical');
for i = a:b
for j = c:d
if skeleton(i,j)
p = [skeleton(i-1,j),skeleton(i-1,j+1),skeleton(i,j+1),skeleton(i+1,j+1),...
skeleton(i+1,j),skeleton(i+1,j-1),skeleton(i,j-1),skeleton(i-1,j-1)];
idx = sum(p.*mat)+1;
if List(idx)
label(i,j) = true;
deleCount = deleCount+1;
end
end
end
end
skeleton(label) = 0;
end
- 实验结果
1)利用查找表,计算速度大大提升;
2)简化后结果发生变化。其中有好的变化,如左上角往外延伸;也有坏的,如部分地方断开(我判定为不可简化)。
MATLAB函数实现
代码非常简单,如下:
% 求骨架
skeleton = bwmorph(im,'skeleton',Inf);
% 去毛刺
skeleton = bwmorph(skeleton,'spur',16);
注:去毛刺,将骨架中的非闭合线进行删除,第三个参数为删除的像素个数
更多实验结果
其他
- 通过实验结果可以发现,我们简化后的算法与MATLAB自带的函数计算结果最为相似。
- 目前还没有深入研究MATLAB提取骨架使用的是什么算法。
- 以上几种结果都没到达我想要的样子。
.
有任何问题欢迎讨论,最后还是把测试代码上传
https://download.csdn.net/download/xsz591541060/11463945
由于很简单,不推荐下载,除非你买了年VIP。