// 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;
}
