Matlab：图像轮廓的曲率计算

给定一个连通区域的图像，如下图所示，想要求其轮廓像素点的曲率。

理论上，下图红框中的轮廓像素的曲率应尽可能大，而蓝框的曲率应比较小。

1.首先对图像二值化，并通过任意算子的边缘提取，得到初步的轮廓（白色像素点）：

I=im2bw(I);
eg=edge(I,'canny');

2.但可能存在某处的轮廓的“厚度”大于1，这样会影响之后按照顺序的轮廓点像素提取，所以对轮廓图像eg使用bwmorph进行细化：

eg=bwmorph(eg,'thin',Inf);

得到细化后的轮廓：

3.细化后的轮廓，保证了某像素点的8邻域内，有且仅有两个其他轮廓点（一个是遍历的上一个像素点，一个是遍历的下一个像素点）。于是我们从任意一个轮廓点(sx,sy)出发，沿着其8邻域内且未遍历过的像素点前进，直至回到起始点。除初始点外，其他点可供选择的邻域像素有且仅有一个，从而保证了程序能够正确地从初始点出发，遍历所有点后回到初始点。

因此我们定义矩阵x和y，用来存放按照顺序遍历的所有像素点的x坐标和y坐标，sx和sy为当前遍历点的坐标。dir为8邻域的方向，对于(sx,sy)的第i个邻域，其坐标为（sx+dir(i,1),sy+dir(i,2)）。定义一个函数findpoint，用来查找当前邻域像素点是否已经遍历过。如果邻域像素点值为1（白色），并且未遍历过，那么sx,xy更新为这个邻域像素点的坐标，并存入x,y。为了防止进入无限循环，设置了状态f，如果8个邻域像素都不能满足条件，说明走到了断头路，那么直接break。

x=[];
y=[];[sx,sy]=find(eg==1,1);%find the first point
x=[x,sx];
y=[y,sy];
dir=[[-1,-1];[-1,0];[-1,1];[0,1];[1,1];[1,0];[1,-1];[0,-1]];
while(true)f=0;for i=1:8if eg(sx+dir(i,1),sy+dir(i,2))==1 && findpoint(x,y,sx+dir(i,1),sy+dir(i,2))sx=sx+dir(i,1);sy=sy+dir(i,2);f=1;breakendendif (f==0)%结束条件1：走到断头路  一般情况下不会满足breakendif(sx==x(1) && sy==y(1)) %结束条件2：走到起始点breakendx=[x,sx];y=[y,sy];
end

函数findpoint，仅当sx和sy同时出现在x y的同一个索引下的值时，才返回0，表明该点已经遍历过，否则返回1.

function f=findpoint(x,y,sx,sy)n=length(x);f=1;for i=1:nif x(i)==sx && y(i)==syf=0;break;endend
end

4.这样得到的x和y连个矩阵（列表/向量），每个x(i),y(i)坐标的点，都会和x(i-1),y(i-1)以及x(i+1),y(i+1)相邻，因此用相邻的x(i)和x(i+1)的差值来近似在x方向的导数，y同理。因此得到一阶差值（导数）矩阵dfx1、dfy1；二阶差值（导数）在一阶差值（导数）基础上继续操作即可，得到dfx2,dfy2。

dfx1=[diff(x,1),x(1)-x(end)];
dfy1=[diff(y,1),y(1)-y(end)];
dfx2=[diff(dfx1,1),dfx1(1)-dfx1(end)];
dfy2=[diff(dfy1,1),dfy1(1)-dfy1(end)];

相邻差值可以用diff直接求得，但是会丢失最后一个点的差值，用x(1)-x(end)进行补充即可。

5.现在得到的差值矩阵dfx1(i),dfy(i)表示x(i)，y(i)点的“导数”，为了直观地显示图像中每个点的曲率，我们将dfx1等映射到二维矩阵dx1等，dfx1(i,j),dfy(i,j)表示点i,j的“导数”。

dx1=zeros(m,n);
dx2=zeros(m,n);
dy1=zeros(m,n);
dy2=zeros(m,n);
for i=1:length(x)dx1(x(i),y(i))=dfx1(i);dy1(x(i),y(i))=dfy1(i);dx2(x(i),y(i))=dfx2(i);dy2(x(i),y(i))=dfy2(i);
end

6.按照曲率公式进行计算曲率矩阵K：

K=zeros(m,n);
K=abs(dx1.*dy2-dx2.*dy1)./((dx1.^2+dy1.^2+0.1).^(3/2));

为了防止分母为0，在分母加入0.1。

直接绘制K的图像：

颜色越白，代表曲率越大，但是这样的结果很不满意，究其原因，是因为用离散点的差值来近似连续曲线的导数，误差太大（比如用(sin(1)-sin(0))/(1-0)近似sinx在x=0的导数，得到0.017，这和实际导数为1的误差过大）。因此，通过相邻差值的方法计算导数或者曲率效果并不好。

因此，我们试图去使用x(i+1)-x(i-1)代替x(i+1)-x(i)来近似x(i)，y(i)在x方向的导数：

dfx1=([diff(x,1),x(1)-x(end)]+[x(1)-x(end),diff(x,1)])/2;
dfy1=([diff(y,1),y(1)-y(end)]+[y(1)-y(end),diff(y,1)])/2;
dfx2=([diff(dfx1,1),dfx1(1)-dfx1(end)]+[dfx1(1)-dfx1(end),diff(dfx1,1)])/2;
dfy2=([diff(dfy1,1),dfy1(1)-dfy1(end)]+[dfy1(1)-dfy1(end),diff(dfy1,1)])/2;

但是结果仍不可观，似乎准确地在离散的点中准确地近似导数是不太可能的事（或者尝试到附近n个点的差值？），特别是在特别小的图像而言。试想曲率的定义似乎和“直率”恰恰相反，也就是说这个像素点附近（这里指的是连续）的点呈现地越接近于直线，曲率越低，反而如果发生拐弯，曲率应该比较大，这不和线性回归有点像么！

因此，采用了线性回归的思路，对每个点的附近的点组成的线段的弯曲程度进行计算。假设取点x(i),y(i)之前共pre和点，之后共next个点，那么线段由pre+next+1个点组成，其坐标为x(i-pre:i+next),y(i-pre:i+next)。

比如当前点用红色表示，之前的点用绿色表示，之后的点用蓝色表示，对他们进行线性回归后，所得的不呈直线的概率，或者回归的误差等参数都可以作为这个点的值。在此以regress返回的stats(1)为例，直线越直，这个stats(1)越接近于1，反之接近于0.因此在这用1-stats(1)表示弯曲程度。依次计算后得到的K为：

白色亮点分布在拐弯处，还是比较符合我们的要求的，OVER。

完整代码：

（为了避免讨论临界值问题，在线性回归时对x和y的边界进行了前后扩充）

I=imread('平滑.jpg');
I=im2bw(I);
[m,n]=size(I);
eg=edge(I,'canny');
eg=bwmorph(eg,'thin',Inf);
x=[];
y=[];[sx,sy]=find(eg==1,1);%find the first point
x=[x,sx];
y=[y,sy];
dir=[[-1,-1];[-1,0];[-1,1];[0,1];[1,1];[1,0];[1,-1];[0,-1]];
while(true)f=0;for i=1:8if eg(sx+dir(i,1),sy+dir(i,2))==1 && findpoint(x,y,sx+dir(i,1),sy+dir(i,2))sx=sx+dir(i,1);sy=sy+dir(i,2);f=1;breakendendif (f==0)breakendif(sx==x(1) && sy==y(1))breakendx=[x,sx];y=[y,sy];
end
%用线性回归类比曲率
pre=5;
next=5;
x=[x(end-pre+1:end),x,x(1:next)];
y=[y(end-pre+1:end),y,y(1:next)];
K=zeros(m,n);
for i=pre+1:length(x)-nextX=x(i-pre:i+next);X=[ones(size(X')),X'];Y=y(i-pre:i+next);[b,bint,r,rint,stats]=regress(Y',X,0.05);K(x(i),y(i))=1-stats(1);
end
%传统方式计算曲率的代码
% dx1=zeros(m,n);
% dx2=zeros(m,n);
% dy1=zeros(m,n);
% dy2=zeros(m,n);% dfx1=[diff(x,1),x(1)-x(end)];
% dfy1=[diff(y,1),y(1)-y(end)];
% dfx2=[diff(dfx1,1),dfx1(1)-dfx1(end)];
% dfy2=[diff(dfy1,1),dfy1(1)-dfy1(end)];%
% dfx1=([diff(x,1),x(1)-x(end)]+[x(1)-x(end),diff(x,1)])/2;
% dfy1=([diff(y,1),y(1)-y(end)]+[y(1)-y(end),diff(y,1)])/2;
% dfx2=([diff(dfx1,1),dfx1(1)-dfx1(end)]+[dfx1(1)-dfx1(end),diff(dfx1,1)])/2;
% dfy2=([diff(dfy1,1),dfy1(1)-dfy1(end)]+[dfy1(1)-dfy1(end),diff(dfy1,1)])/2;
%
% for i=1:length(x)
%     dx1(x(i),y(i))=dfx1(i);
%     dy1(x(i),y(i))=dfy1(i);
%     dx2(x(i),y(i))=dfx2(i);
%     dy2(x(i),y(i))=dfy2(i);
% end
%
% K=zeros(m,n);
% K=abs(dx1.*dy2-dx2.*dy1)./((dx1.^2+dy1.^2+0.1).^(3/2));
function f=findpoint(x,y,sx,sy)n=length(x);f=1;for i=1:nif x(i)==sx && y(i)==syf=0;break;endend
end

Matlab：图像轮廓的曲率计算相关推荐

matlab 图像轮廓填充颜色,基于Matlab的图形轮廓提取及填充
计算机工程应用技术本栏目责任编辑: 贾薇薇电脑知识与技术基于 Matlab 的图形轮廓提取及填充井艾斌,柳青,孟祥增 (山东师范大学, 山东济南 250014) 摘要: 提取图形的形状特征是 ...
一种新型鱼眼图像轮廓提取算法
from: http://www.scimao.com/read/2307651 摘要:提取鱼眼图像轮廓是利用鱼眼图像的前提.传统提取鱼眼图像轮廓的扫描线逼近法对噪点抑制能力不强,精度差.本 ...
matlab 提取图像轮廓（图像边缘提取）
利用edge()函数提取图像轮廓,绘制出对象的边界和提取边界坐标信息,matlab实现代码如下: close all;clear all;clc; % 提取图像轮廓,提取图像边缘 I = imread ...
Python机器视觉--OpenCV进阶(核心)--图像轮廓查找识别,绘制图像轮廓与图像轮廓的面积周长计算
1.图像轮廓查找识别与绘制图像轮廓 1.1 什么是图像轮廓图像轮廓是具有相同颜色或灰度的连续点的曲线. 轮廓在形状分析和物体的检测和识别中很有用. 轮廓的作用: 用于图形分析物体的识别和检测注意 ...
基于hvs图像水印matlab和psnr nc的计算首先读取图像和水印，进行图像加印
基于hvs图像水印matlab和psnr nc的计算首先读取图像和水印,进行图像加印然后进行攻击攻击方式有白噪声,裁剪,旋转10度,压缩,和无攻击,然后最后还原水印. ID:3131461758 ...
图片缩放算法 matlab,图像放大算法总结及MATLAB源程序.doc
图像放大算法总结及MATLAB源程序 1,插值算法(3种): (1)最邻近插值(近邻取样法): 最近插值的的思想很简单就是把这个非整数坐标作一个四舍五入,取最近的整数点坐标处的点的颜色.可见,最邻近插 ...
MATLAB图像融合分割系统
摘要图像分割是一种重要的图像分析技术.对图像分割的研究一直是图像技术研究中的热点和焦点.图像分割是一个很关键的图像分析技术,是由图像处理进到图像分析的关键步骤.它的目的就是把图像中感兴趣的那部分分 ...
[MATLAB] 图像的插值算法1：MATLAB中的插值函数及其原理
MATLAB图像插值算法文章集: 插值函数及其原理 https://blog.csdn.net/Effend/article/details/82870144 最近邻插值 https://blog.c ...
基于中科院-CASIA-GaitDatasetB步态图像轮廓数据库的步态周期检测与步态角度特征
基于中科院-CASIA-GaitDatasetB步态图像轮廓数据库的步态周期检测与步态角度特征由于最近研究的需要,开始将万恶的双手伸向了人体下肢步态(如步态周期检测.步态特征提取.步态相的划分 ...

Matlab：图像轮廓的曲率计算

Matlab：图像轮廓的曲率计算相关推荐

最新文章

热门文章