// LocalHistogramOfOrientedGradients.cpp : 定义控制台应用程序的入口点。
//局部方向梯度直方图
#include <iostream>
#include <opencv2\opencv.hpp>
using
namespace
std
;
using
namespace
cv
;
/*
计算给定像素位置上的加权方向梯度直方图(orientation gradient histogram)
img:原始图像
pt: 指定的像素点
radius: 统计半径
isSmoothed:是否平滑输出直方图
用高斯函数进行中心加权;
isWeighted,和 weighted_sigma: 是否加权和计算权重的标准差
hist: 要输出的直方图
n: 直方图的bin个数
返回值:直方图的峰值(最大值)
*/
static
float
calcOrientationHist
(
const
Mat
&
img
,
Point
pt
,
int
radius
,
float
*
hist
,
int
n
,
int
isSmoothed
,
int
isWeighted
,
float
weighted_sigma
)
{
//radius应该是以pt为中心的正方形的边长的一半
int
i
,
j
,
k
,
len
=
(
radius
*
2
+
1
)
*
(
radius
*
2
+
1
);
//使用高斯函数中心加权
float
expf_scale
=
-
1.f
/
(
2.f
*
weighted_sigma
*
weighted_sigma
);
//为什么加4呢,是为了给临时直方图开辟多余的4个存储位置,
//用来存放temphist[-1],temphist[-2],temphist[n],temphist[n+1]的
//为什么加n呢,这n个位置是留给temphist[0 ... n-1]的
//为什么len*4呢,这4个len长度的数组位置是留给X、Y、W以及方向Ori的
AutoBuffer
<
float
>
buf
(
len
*
4
+
n
+
4
);
//X是横向梯度,Y是纵向梯度,Mag是梯度幅值=sqrt(X^2+Y^2), Ori是梯度方向 = arctan(Y/X)
float
*
X
=
buf
,
*
Y
=
X
+
len
,
*
Mag
=
X
,
*
Ori
=
Y
+
len
,
*
W
=
Ori
+
len
;
float
*
temphist
=
W
+
len
+
2
;
//加2是用来存放temphist[-1],temphist[-2]的
//临时直方图清零
for
(
i
=
0
;
i
<
n
;
i
++
)
temphist
[
i
]
=
0.f
;
//从上往下,从左往右扫描求横向,纵向的梯度值以及对应的权值
for
(
i
=
-
radius
,
k
=
0
;
i
<=
radius
;
i
++
)
{
int
y
=
pt
.
y
+
i
;
//指向原图像img的第pt.y+i行
if
(
y
<=
0
||
y
>=
img
.
rows
-
1
)
//边界检查
continue
;
for
(
j
=
-
radius
;
j
<=
radius
;
j
++
)
{
int
x
=
pt
.
x
+
j
;
//指向原图像img的第pt.x+j列
if
(
x
<=
0
||
x
>=
img
.
cols
-
1
)
//边界检查
continue
;
//横向梯度
float
dx
=
(
float
)(
img
.
at
<
uchar
>
(
y
,
x
+
1
)
-
img
.
at
<
uchar
>
(
y
,
x
-
1
));
//纵向梯度
float
dy
=
(
float
)(
img
.
at
<
uchar
>
(
y
-
1
,
x
)
-
img
.
at
<
uchar
>
(
y
+
1
,
x
));
//保存纵向梯度和横向梯度
X
[
k
]
=
dx
;
Y
[
k
]
=
dy
;
//计算加权数组
if
(
isWeighted
)
W
[
k
]
=
(
i
*
i
+
j
*
j
)
*
expf_scale
;
else
W
[
k
]
=
1.f
;
//如果不加权,则每个统计点上的权重是一样的
k
++
;
}
}
//把实际的统计点数复制给len,由于矩形局部邻域有可能超出图像边界,
len
=
k
;
//所以实际的点数总是小于等于 (radius*2+1)*(radius*2+1)
//在指定像素点的邻域内 计算梯度幅值, 梯度方向 and 权重
exp
(
W
,
W
,
len
);
//权重
fastAtan2
(
Y
,
X
,
Ori
,
len
,
true
);
//梯度方向
magnitude
(
X
,
Y
,
Mag
,
len
);
//梯度幅值
//填充临时直方图,直方图的横轴是梯度方向方向角度[0,360),bin宽度为n/360;
//纵轴是梯度幅值乘以对应的权重
for
(
k
=
0
;
k
<
len
;
k
++
)
{
int
bin
=
cvRound
((
n
/
360.f
)
*
Ori
[
k
]);
//求出第k个角度Ori[k]的bin索引号
if
(
bin
>=
n
)
bin
-=
n
;
if
(
bin
<
0
)
bin
+=
n
;
temphist
[
bin
]
+=
W
[
k
]
*
Mag
[
k
];
}
if
(
isSmoothed
)
{
// 直方图平滑,平滑后放入输出直方图数组中
temphist
[
-
1
]
=
temphist
[
n
-
1
];
temphist
[
-
2
]
=
temphist
[
n
-
2
];
temphist
[
n
]
=
temphist
[
0
];
temphist
[
n
+
1
]
=
temphist
[
1
];
for
(
i
=
0
;
i
<
n
;
i
++
)
{
hist
[
i
]
=
(
temphist
[
i
-
2
]
+
temphist
[
i
+
2
])
*
(
1.f
/
16.f
)
+
(
temphist
[
i
-
1
]
+
temphist
[
i
+
1
])
*
(
4.f
/
16.f
)
+
temphist
[
i
]
*
(
6.f
/
16.f
);
}
}
else
//不平滑直方图
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
hist
[
i
]
=
temphist
[
i
];
}
}
//求直方图梯度的最大值
float
maxval
=
hist
[
0
];
for
(
i
=
1
;
i
<
n
;
i
++
)
maxval
=
std
::
max
(
maxval
,
hist
[
i
]);
return
maxval
;
}
void
DrawHist
(
Mat
&
hist
,
string
&
winname
)
{
Mat
drawHist
;
int
histSize
=
hist
.
rows
;
// 创建直方图画布
int
hist_w
=
360
;
int
hist_h
=
360
;
//直方图图像的宽度和高度
int
bin_w
=
cvRound
(
(
double
)
hist_w
/
histSize
);
//直方图中一个矩形条纹的宽度
Mat
histImage
(
hist_w
,
hist_h
,
CV_8UC3
,
Scalar
(
0
,
0
,
0
)
);
//创建画布图像
/// 将直方图归一化到范围 [ 0, histImage.rows ]
normalize
(
hist
,
drawHist
,
0
,
histImage
.
rows
,
NORM_MINMAX
,
-
1
,
Mat
()
);
/// 在直方图画布上画出直方图
for
(
int
i
=
1
;
i
<
histSize
;
i
++
)
{
//矩形图表示
rectangle
(
histImage
,
Point
((
i
-
1
)
*
bin_w
,
hist_h
),
Point
(
i
*
bin_w
,
hist_h
-
cvRound
(
drawHist
.
at
<
float
>
(
i
-
1
))),
Scalar
(
0
,
0
,
255
),
1
,
8
,
0
);
//折线表示
/* line( histImage, Point( bin_w*(i-1), hist_h - cvRound(hist.at<float>(i-1)) ) ,
Point( bin_w*(i), hist_h - cvRound(hist.at<float>(i)) ),
Scalar( 0, 0, 255), 1, 8, 0 );*/
}
/// 显示直方图
cv
::
namedWindow
(
winname
,
1
);
cv
::
imshow
(
winname
,
histImage
);
}
int
main
(
int
argc
,
char
**
argv
)
{
const
string
filename
=
"meinv2.jpg"
;
Mat
Image
=
imread
(
filename
,
1
);
if
(
Image
.
empty
())
{
std
::
cout
<<
"无法读取图像...."
<<
endl
;
getchar
();
}
Mat
grayImage
;
//灰度图像
cvtColor
(
Image
,
grayImage
,
CV_BGR2GRAY
);
//-------------------------------计算方向梯度直方图的参数------------------------
Point
center
(
grayImage
.
cols
/
2
,
grayImage
.
rows
/
2
);
cout
<<
"原图中心点: "
<<
center
<<
endl
;
int
radius
=
120
;
//局部邻域半径
Rect
StatisticArea
(
center
.
x
-
radius
,
center
.
y
-
radius
,
2
*
radius
,
2
*
radius
);
int
bins_count
=
60
;
//bin的数目
float
sigma
=
radius
*
0.5f
;
//加权函数的标准差,设为统计区域的半径的一半
bool
isSmoothed
=
true
;
//是否平滑直方图
bool
isWeighted
=
false
;
//不加权
//-------------------------------计算原图中心点的方向梯度直方图---------------------
//声明一个直方图数组
Mat
originHist
=
Mat
::
zeros
(
bins_count
,
1
,
CV_32FC1
);
float
*
oh
=
(
float
*
)
originHist
.
data
;
if
(
!
originHist
.
isContinuous
())
{
cout
<<
"直方图地址存储不连续"
<<
endl
;
getchar
();
}
calcOrientationHist
(
grayImage
,
center
,
radius
,
oh
,
bins_count
,
isSmoothed
,
isWeighted
,
sigma
);
//绘制直方图
string
winname
=
"Origin hist"
;
DrawHist
(
originHist
,
winname
);
//---------------计算旋转图像的方向梯度直方图-------------------------------------------
//step:1 计算绕图像中点顺时针旋转30度缩放因子为1的旋转矩阵
Point
rot_center
=
center
;
//旋转中心
double
angle
=
30.0
;
//旋转角度
double
scale
=
1.
;
//缩放因子
/// 通过上面的旋转细节信息求得旋转矩阵
Mat
rot_mat
=
getRotationMatrix2D
(
rot_center
,
angle
,
scale
);
cout
<<
"旋转矩阵:"
<<
rot_mat
<<
endl
;
/// 调用仿射变换函数旋转原始图像
Mat
rotate_dst
;
warpAffine
(
grayImage
,
rotate_dst
,
rot_mat
,
grayImage
.
size
()
);
//声明一个直方图数组
Mat_
<
float
>
rotateHist
=
Mat
::
zeros
(
bins_count
,
1
,
CV_32FC1
);
float
*
rh
=
(
float
*
)
rotateHist
.
data
;
calcOrientationHist
(
rotate_dst
,
center
,
radius
,
rh
,
bins_count
,
isSmoothed
,
isWeighted
,
sigma
);
//绘制直方图
string
winname2
=
"rotated hist"
;
DrawHist
(
rotateHist
,
winname2
);
//- -利用旋转前和旋转后的两个直方图的纵轴主峰对应的角度bin估算图像的旋转角度--------------------
cout
<<
"直方图bin的宽度: "
<<
(
360.f
/
bins_count
)
<<
"度"
<<
endl
;
//求出旋转前图像方向梯度直方图的主峰位置对应的bin角度值
Point
max_location1
;
//直方图主峰对应的bin位置
float
angle1
;
minMaxLoc
(
originHist
,
0
,
0
,
0
,
&
max_location1
);
angle1
=
max_location1
.
y
*
(
360.f
/
bins_count
);
cout
<<
"未旋转之前的方向梯度直方图的主峰位置:"
<<
max_location1
.
y
<<
endl
<<
" ,对应角度值: "
<<
angle1
<<
endl
;;
//求出旋转后图像方向梯度直方图的主峰位置对应的bin角度值
Point
max_location2
;
//直方图主峰对应的bin位置
float
angle2
;
minMaxLoc
(
rotateHist
,
0
,
0
,
0
,
&
max_location2
);
angle2
=
max_location2
.
y
*
(
360.f
/
bins_count
);
cout
<<
"旋转之后的方向梯度直方图的主峰位置:"
<<
max_location2
.
y
<<
endl
<<
" ,对应角度值: "
<<
angle2
<<
endl
;;
cout
<<
"真实旋转角度:"
<<
angle
<<
endl
;
cout
<<
"估计的旋转角度: "
<<
angle2
-
angle1
<<
endl
;
rectangle
(
grayImage
,
StatisticArea
,
Scalar
(
0
),
2
);
imshow
(
"origin img"
,
grayImage
);
rectangle
(
rotate_dst
,
StatisticArea
,
Scalar
(
0
),
2
);
imshow
(
"rotated img"
,
rotate_dst
);
waitKey
(
0
);
return
0
;
}