Steger(线条中心提取算法)
1、算法原理
Steger算法基于Hessian矩阵,能够实现对线条中心的亚像素提取。对于图像上任一点(x,y),Hessian矩阵可以表示为:
H
(
x
,
y
)
=
[
h
x
x
h
x
y
h
x
y
h
y
y
]
H(x,y)=\left[ \begin{matrix} h_{xx} & h_{xy}\\ h_{xy} & h_{yy}\\ \end{matrix} \right]
H(x,y)=[hxxhxyhxyhyy]
该矩阵由图像某点沿不同方向的四个偏导数组成,该矩阵可以确定z(x,y)的二阶方向导数最大绝对值的方向。
Steger提取线条中心的步骤如下:
①高斯平滑;
②基于Hessian矩阵得到该像素点的法向;
③法向上泰勒展开拟合亚像素坐标。
Hessian 矩阵最大特征值对应的特征向量对应于线条的法线方向,用
(
n
x
,
n
y
)
(n_x,n_y)
(nx,ny) 表示,以点
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)为基准点,
图像上每个点都可以计算 Hessian矩阵,由于Hessian矩阵的特征值就是形容其在该点附近特征向量方向的凹凸性,特征值越大,凸性越强。线条中心点的Hessian矩阵特征值相比其他位置更大,所以可以设置阈值进行筛选,线条中心的亚像素坐标为:
(
p
x
,
p
y
)
=
(
x
0
+
t
n
x
,
y
0
+
t
N
y
)
其
中
,
t
=
−
r
x
n
x
+
r
y
n
y
r
x
x
n
x
2
+
2
r
x
y
n
x
n
y
+
r
y
y
n
y
2
若
存
在
点
满
足
(
t
n
x
,
t
n
y
)
ϵ
[
−
0.5
,
0.5
]
∗
[
−
0.5
,
0.5
]
(p_x,p_y) = (x_0 + tn_x,y_0+tN_y) \\其中,\\ t=-\frac{r_xn_x+r_yn_y}{r_{xx}n_x^2+2r_{xy}n_xn_y+r_{yy}n_y^2} \\若存在点满足\\ (tn_x,tn_y)\epsilon[-0.5,0.5]*[-0.5,0.5]
(px,py)=(x0+tnx,y0+tNy)其中,t=−rxxnx2+2rxynxny+ryyny2rxnx+ryny若存在点满足(tnx,tny)ϵ[−0.5,0.5]∗[−0.5,0.5]
即一阶导数为零的点位于当前像素内(找拐点附近的点,因为实际中心截面不是规则的曲线),且方向的二阶导数大于指定的阈值(确定是极值点),则该点
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0)为线条的中心点,
(
p
x
,
p
y
)
(p_x,p_y)
(px,py)则为亚像素坐标。
(
t
n
x
,
t
n
y
)
(tn_x,tn_y)
(tnx,tny)是中心点与亚像素的偏移值,所以其绝对值不能超过0.5,超过0.5就超出该像素的范围。
2、参考代码
void StegerCV()
{
Mat srcImg = imread("srcImg.png", 1);
Mat grayImg;
cvtColor(srcImg, grayImg, CV_BGR2GRAY);
// 高斯滤波
grayImg.convertTo(grayImg, CV_32FC1);
GaussianBlur(grayImg, grayImg, Size(0, 0), 6, 6);
// 一阶偏导数
Mat m1, m2;
m1 = (Mat_<float>(1, 2) << 1, -1); // x偏导
m2 = (Mat_<float>(2, 1) << 1, -1); // y偏导
Mat dx, dy;
filter2D(grayImg, dx, CV_32FC1, m1);
filter2D(grayImg, dy, CV_32FC1, m2);
//二阶偏导数
Mat m3, m4, m5;
m3 = (Mat_<float>(1, 3) << 1, -2, 1); // 二阶x偏导
m4 = (Mat_<float>(3, 1) << 1, -2, 1); // 二阶y偏导
m5 = (Mat_<float>(2, 2) << 1, -1, -1, 1); // 二阶xy偏导
Mat dxx, dyy, dxy;
filter2D(grayImg, dxx, CV_32FC1, m3);
filter2D(grayImg, dyy, CV_32FC1, m4);
filter2D(grayImg, dxy, CV_32FC1, m5);
// hessian矩阵
double maxD = -1;
int imgcol = img.cols;
int imgrow = img.rows;
vector<double> Pt;
for (int i=0;i<imgcol;i++)
{
for (int j=0;j<imgrow;j++)
{
if (img0.at<uchar>(j,i)>200)
{
Mat hessian(2, 2, CV_32FC1);
hessian.at<float>(0, 0) = dxx.at<float>(j, i);
hessian.at<float>(0, 1) = dxy.at<float>(j, i);
hessian.at<float>(1, 0) = dxy.at<float>(j, i);
hessian.at<float>(1, 1) = dyy.at<float>(j, i);
Mat eValue;
Mat eVectors;
eigen(hessian, eValue, eVectors);
double nx, ny;
double fmaxD = 0;
if (fabs(eValue.at<float>(0,0))>= fabs(eValue.at<float>(1,0))) //求特征值最大时对应的特征向量
{
nx = eVectors.at<float>(0, 0);
ny = eVectors.at<float>(0, 1);
fmaxD = eValue.at<float>(0, 0);
}
else
{
nx = eVectors.at<float>(1, 0);
ny = eVectors.at<float>(1, 1);
fmaxD = eValue.at<float>(1, 0);
}
double t = -(nx*dx.at<float>(j, i) + ny*dy.at<float>(j, i)) / (nx*nx*dxx.at<float>(j,i)+2*nx*ny*dxy.at<float>(j,i)+ny*ny*dyy.at<float>(j,i));
if (fabs(t*nx) <= 0.5 && fabs(t*ny) <= 0.5)
{
double x_sub = i + t*nx;
double y_sub = j + t*ny;
Pt.push_back(x_sub);
Pt.push_back(x_sub);
}
}
}
}
for (int k = 0;k<Pt.size()/2;k++)
{
Point rpt;
rpt.x = Pt[2 * k + 0];
rpt.y = Pt[2 * k + 1];
circle(srcImg, rpt, 1, Scalar(0, 0, 255));
}
imshow("result", srcImg);
waitKey(0);
}