2017-04-10 22:58:50 qq826309057 阅读数 13179
• ###### MATLAB图像处理

全面系统的学习MATLAB在图像处理中的应用

20781 人正在学习 去看看 魏伟

## 数学中的矩

In mathematics, a moment is a specific quantitative measure, used in both mechanics and statistics, of the shape of a set of points. If the points represent mass, then the zeroth moment is the total mass, the first moment divided by the total mass is the center of mass, and the second moment is the rotational inertia. If the points represent probability density, then the zeroth moment is the total probability (i.e. one), the first moment is the mean, the second central moment is the variance, the third moment is the skewness, and the fourth moment (with normalization and shift) is the kurtosis. The mathematical concept is closely related to the concept of moment in physics.For a bounded distribution of mass or probability, the collection of all the moments (of all orders, from 0 to ∞) uniquely determines the distribution.

• c = 0 时， 称为 k 阶原点矩；
• c = E(x) 时，称为 k 阶中心矩。

## 图像的矩

### 零阶矩：

${M}_{00}=\sum _{I}\sum _{J}V\left(i,j\right)$

### 一阶矩：

${M}_{10}=\sum _{I}\sum _{J}i\cdot V\left(i,j\right)\phantom{\rule{0ex}{0ex}}{M}_{01}=\sum _{I}\sum _{J}j\cdot V\left(i,j\right)$

${x}_{c}=\frac{{M}_{10}}{{M}_{00}},{y}_{c}=\frac{{M}_{01}}{{M}_{00}}$

### 二阶矩

${M}_{20}=\sum _{I}\sum _{J}{i}^{2}\cdot V\left(i,j\right)\phantom{\rule{0ex}{0ex}}{M}_{02}=\sum _{I}\sum _{J}{j}^{2}\cdot V\left(i,j\right)\phantom{\rule{0ex}{0ex}}{M}_{11}=\sum _{I}\sum _{J}i\cdot j\cdot V\left(i,j\right)$

$\theta =\frac{1}{2}fastAtan2\left(2b,a-c\right)$

fastAtan2()为opencv的函数，输入向量，返回一个0-360的角度。

$\theta =\frac{1}{2}arctan\left(\frac{2b}{a-c}\right)$

## 利用矩求图像的重心，方向

Mat image = imread(imagename, 0);//读入灰度图
Mat binary;
//二值，椭圆是黑色的，所以反色下
threshold(image, binary, 200, 255, CV_THRESH_BINARY_INV);
Moments m = moments(binary, true);//moments()函数计算出三阶及一下的矩
Point2d center(m.m10 / m.m00, m.m01 / m.m00);//此为重心
//计算方向
double a = m.m20 / m.m00 - center.x*center.x;
double b = m.m11 / m.m00 - center.x*center.y;
double c = m.m02 / m.m00 - center.y*center.y;
double theta = fastAtan2(2*b,(a - c))/2;//此为形状的方向

## 总结

2016-12-08 11:41:15 m_buddy 阅读数 4103
• ###### MATLAB图像处理

全面系统的学习MATLAB在图像处理中的应用

20781 人正在学习 去看看 魏伟

## 2. 实现

//************************************************************************
// 函数名称:    	ClacCenterPointXY
// 访问权限:    	public
// 创建日期:		2016/11/30
// 创 建 人:
// 函数说明:		计算图像的型心（重心）
// 函数参数: 	cv::Mat & src_img		输入图像数据
// 函数参数: 	int & centerpoint_x		计算得到的图像的X轴型心或是重心
// 函数参数: 	int & centerpoint_y		计算得到的图像的Y轴型心或是重心
// 返 回 值:   	BOOL
//************************************************************************
bool CCalcMutualInfo::ClacCenterPointXY(cv::Mat& src_img, int& centerpoint_x, int& centerpoint_y)
{
if (!src_img.data)
{
cout << "LaplacianEnhance 需要计算直方图的图像没有图像数据" << endl;
return false;
}
if (CV_8UC1 != src_img.type())
{
cout << "LaplacianEnhance 输入的图像不满足的8 bit 灰度图像要求" << endl;
return false;
}

int rows(src_img.rows);	//height
int cols(src_img.cols);	//width
uchar* data = nullptr;

double m_sumx(0.0), m_sumy(0.0);
long m_sum(0);
for (int i=0; i<rows; i++)
{
data = src_img.ptr<uchar>(i);
for (int j=0; j<cols; j++)
{
m_sum += data[j];
m_sumx += i*data[j];
m_sumy += j*data[j];
}
}

centerpoint_x = static_cast<int>(m_sumx / m_sum + 0.5);
centerpoint_y = static_cast<int>(m_sumy / m_sum + 0.5);

return true;
}

//************************************************************************
// 函数名称:    	ClacCovarianceMatrix
// 访问权限:    	public
// 创建日期:		2016/11/30
// 创 建 人:
// 函数说明:		计算图像的协方差矩阵
// 函数参数: 	cv::Mat & src_img				输入图像数据
// 函数参数: 	cv::Mat & CovarianceMatrix		输出的写图像协方差矩阵
// 函数参数: 	int centerpoint_x				输入图像X轴的重心
// 函数参数: 	int centerpoint_y				输入图像Y轴的重心
// 返 回 值:   	BOOL
//************************************************************************
bool CCalcMutualInfo::ClacCovarianceMatrix(cv::Mat& src_img, cv::Mat& CovarianceMatrix, int centerpoint_x, int centerpoint_y)
{
if (!src_img.data)
{
cout << "LaplacianEnhance 需要计算直方图的图像没有图像数据" << endl;
return false;
}
if (CV_8UC1 != src_img.type())
{
cout << "LaplacianEnhance 输入的图像不满足的8 bit 灰度图像要求" << endl;
return false;
}

int rows(src_img.rows);	//height
int cols(src_img.cols);	//width
unsigned char* data = nullptr;

double C11(0.0), C12(0.0), C21(0.0), C22(0.0);
for (int i = 0; i < rows; i++)
{
data = src_img.ptr<uchar>(i);
for (int j = 0; j < cols; j++)
{
C11 += std::pow((double)(i - centerpoint_x), 2.0)*data[j];
C12 += (double)(i-centerpoint_x) * (double)(j - centerpoint_y) * data[j];
C21 = C12;
C22 += std::pow((double)(j - centerpoint_y), 2.0)*data[j];
}
}
cv::Mat temp = (cv::Mat_<double>(2,2) << C11, C12, C21, C22);
temp.copyTo(CovarianceMatrix);

return true;
}

//************************************************************************
// 函数名称:    	GetCosAngle
// 访问权限:    	public
// 创建日期:		2016/12/05
// 创 建 人:
// 函数说明:		获得模板协方差矩阵与待测试图像协方差矩阵之间的夹角
// 函数参数: 	cv::Mat & main_matrix_value	模板协方差矩阵的特征值
// 函数参数: 	cv::Mat & main_matrix_vec	模板协方差矩阵的特征向量
// 函数参数: 	cv::Mat & main_value		待测试图像协方差矩阵的特征值
// 函数参数: 	cv::Mat & main_vec			待测试图像协方差矩阵的特征向量
// 返 回 值:   	double
//************************************************************************
double CCalcMutualInfo::GetCosAngle(cv::Mat& main_matrix_value, cv::Mat& main_matrix_vec, cv::Mat& main_value, cv::Mat& main_vec)
{
double angle_cos(0.0);
double* main_vec1 = new double[2];
double* main_vec2 = new double[2];
this->GetMainVec(main_vec1, main_matrix_value, main_matrix_vec);
this->GetMainVec(main_vec2, main_value, main_vec);

angle_cos = std::acos(main_vec1[0]*main_vec2[0]+main_vec1[1]*main_vec2[1] / (std::sqrt(main_vec1[0]*main_vec1[0]+main_vec1[1]*main_vec1[1])*
std::sqrt(main_vec2[0]*main_vec2[0]+main_vec2[1]*main_vec2[1])));
//angle_cos = angle_cos*180/3.1415926f;

delete[] main_vec1;
main_vec1 = nullptr;
delete[] main_vec2;
main_vec2 = nullptr;

return angle_cos;
}

//************************************************************************
// 函数名称:    	GetMainVec
// 访问权限:    	public
// 创建日期:		2016/12/05
// 创 建 人:
// 函数说明:		获得主向量
// 函数参数: 	double * dst_vec		输出的图像的主轴方向
// 函数参数: 	cv::Mat & matrix_value	输入图像的协方差矩阵的特征值
// 函数参数: 	cv::Mat & matrix_vec	输入图像的协方差矩阵的特征向量
// 返 回 值:   	bool
//************************************************************************
bool CCalcMutualInfo::GetMainVec(double* dst_vec, cv::Mat& matrix_value, cv::Mat& matrix_vec)
{
if (!matrix_vec.data || !matrix_value.data)
{
cout << "no input data" <<  endl;
return false;
}
if (dst_vec == nullptr)
{
cout << "input pointer is nullptr" << endl;
return false;
}

double* data1 = nullptr;
double* data2 = nullptr;
data1 = matrix_value.ptr<double>(0);
data2 = matrix_value.ptr<double>(1);
int max_pos = data1[0] > data2[1] ? 0 : 1;

data1 = matrix_vec.ptr<double>(0);
data2 = matrix_vec.ptr<double>(1);
dst_vec[0] = data1[max_pos];
dst_vec[1] = data2[max_pos];

return true;
}

//************************************************************************
// 函数名称:    	ReposImg
// 访问权限:    	public
// 创建日期:		2016/12/08
// 创 建 人:
// 函数说明:		基于集合变换的简单粗暴配准
// 函数参数: 	std::vector<cv::Mat> & img_set	输入的图像序列
// 函数参数: 	cv::Mat & main_evalue			模板图像的协方差特征值
// 函数参数: 	cv::Mat & main_evec				模板图像的协方差特征向量
// 函数参数: 	int center_x					模板图像的X轴重心
// 函数参数: 	int center_y					模板图像的Y轴重心
// 返 回 值:   	bool
//************************************************************************
bool CCalcMutualInfo::ReposImg(std::vector<cv::Mat>& img_set, cv::Mat& main_evalue, cv::Mat& main_evec,
int center_x, int center_y)
{
if (img_set.size()<=0 || !main_evalue.data || !main_evec.data)
{
cout << "no input data " << endl;
}

for (int i=0; i<(int)img_set.size(); i++)
{
int center_x1, center_y1, center_x1_diff, center_y1_diff;
cv::Mat matrix_temp1, matrix_temp1_evalue, matrix_temp1_evec;
ClacCenterPointXY(img_set[i], center_x1, center_y1);
this->StretchImg(img_set[i], img_set[i], 0.488281 / 0.449219, 0.488281 / 0.449219, center_x1, center_y1);	//伸缩和截取图像

ClacCovarianceMatrix(img_set[i], matrix_temp1, center_x1, center_y1);
center_x1_diff = center_x - center_x1;		//计算重心点的差值
center_y1_diff = center_y - center_y1;
cout << "待测试图像"<< i+1 << "的重心差值：X = " << center_x1_diff << " Y=：" << center_y1_diff << endl;
cv::eigen(matrix_temp1, matrix_temp1_evalue, matrix_temp1_evec);
double angle = GetCosAngle(main_evalue, main_evec, matrix_temp1_evalue, matrix_temp1_evec);
cout << "待测试图像"<< i+1 << "与模板图像的夹角为：" << angle << endl << endl;
cv::Mat img_temp = img_set[i];
this->ImgMove(img_temp, img_temp, center_x1_diff, center_y1_diff);
this->ImgRotate(img_temp, img_temp, -angle);
img_temp.copyTo(img_set[i]);
}

//放缩图像到相同的尺度
/*for (int i = 0; i < (int)img_set.size(); i++)
{
this->StretchImg(img_set[i], img_set[i], 0.488281 / 0.449219, 0.488281 / 0.449219, center_x1, center_y1);
}*/

return true;
}

bool CCalcMutualInfo::ShowRectificationResult(cv::Mat& img1, cv::Mat& img2, std::string win_name)
{
if (!img1.data || !img2.data)
{
cout << "addweight input img is null" <<  endl;
return false;
}
cv::Mat img1_temp, img2_temp;
cv::threshold(img1, img1_temp, 10, 255, CV_THRESH_BINARY);
cv::threshold(img2, img2_temp, 10, 255, CV_THRESH_BINARY);

cv::Mat temp1(img1.rows, img1.cols, CV_8UC3, cv::Scalar::all(0));
cv::Mat temp2(img1.rows, img1.cols, CV_8UC3, cv::Scalar::all(0));
unsigned char *data1 =nullptr, *data2=nullptr;

for (int i=0; i<img1.rows; i++)
{
data1 = img1_temp.ptr<unsigned char>(i);
data2 = img2_temp.ptr<unsigned char>(i);
for (int j=0; j<img1.cols; j++)
{
temp1.at<cv::Vec3b>(i, j)[2] = data1[j];
temp2.at<cv::Vec3b>(i, j)[1] = data2[j];
}
}

cv::Mat result_img;
cv::addWeighted(temp1, 0.4, temp2, 0.5, 0, result_img);
cv::imshow(win_name, result_img);

return true;
}

int center_x, center_y;								//参考图像的中心点 X轴Y轴
cv::Mat matrix_main, matrix_evalue, matrix_evec;	//参考图像的协方差矩阵
ClacCenterPointXY(src_img, center_x, center_y);
ClacCovarianceMatrix(src_img, matrix_main, center_x, center_y);
cv::eigen(matrix_main, matrix_evalue, matrix_evec);

//待配准的图像序列
//////////////////////////////////////////////////////////////////////////
this->ReposImg(img_vec, matrix_evalue, matrix_evec, center_x, center_y);
this->ShowRectificationResult(src_img, img_vec[0], "result2");

## 4. 总结

2010-08-24 11:41:00 asmemgsd 阅读数 376
• ###### MATLAB图像处理

全面系统的学习MATLAB在图像处理中的应用

20781 人正在学习 去看看 魏伟

1.C++语言和面向对象的理解,数据结构和算法;

2.过去所从事工作的内容与重心(与当前招人条件是否符合);

3.职业规划,个人发展方向

1.总结过去工作用到的知识(MFC,WTL的消息处理,映射,循环等;图像处理的算法,控件的自绘,网络传输的I/O模型,数据结构,算法,);

2.加强基础知识的学习(基础知识包括:C++(虚函数,多态...),面向对象,GDI,GDI+,多线程,网络,数据结构与算法);

3.职业规划,确定个人的发展方向(个性优势,比较适合做哪种开发工作与适合什么职位,客户端,服务器端,管理...等等);

2018-08-29 20:57:20 aaron19890330 阅读数 1512
• ###### MATLAB图像处理

全面系统的学习MATLAB在图像处理中的应用

20781 人正在学习 去看看 魏伟

（2）将最前点投影在2D平面上；

（3）中值滤波和平滑处理；

（4）得到2D点集进行线性插值；

（5）以重心为中心判断是否闭环；

（6）2D点集进行椭圆拟合，判断是否椭圆、椭圆半径和椭圆度范围；

(7)  2D点集进行圆拟合，计算半径均值和标准差；

（8）判断2D点集方向：顺/逆时针。

1、2D点集线性插值算法

2、计算重心

3、闭环判断

4、椭圆拟合

５、圆拟合

６、顺/逆时针方向判断

# 1、2D点集线性插值算法

//************************************
// Description: 插值,根据两点间距离线性插值,两点间距离越大插值越多,距离小于min_dis则不插值返回原值
// Method:    InterpCurve
// FullName:  InterpCurve
// Access:    private
// Parameter: 插值前二维点集 std::vector<cv::Point2f> &data
// Parameter: 插值后二维点集 std::vector<cv::Point2f> &result
// Parameter: 最小插值距离 float min_dis
// Returns:
// Author:
// Date:      2018/08/28
// History:
//************************************
void InterpCurve(std::vector<cv::Point2f> &data, std::vector<cv::Point2f> &result, float min_dis) {
result.clear();
result.reserve(data.size());
result.push_back(data.at(0));
for (size_t i = 1; i < data.size(); i++) {
float dis = Distance(data.at(i - 1), data.at(i));
float num = dis / min_dis + 1;
float step = 1.0f / num;
float s = 0;
while (s <= 1.0) {
result.emplace_back(
cv::Point2f((1 - s) * data.at(i - 1).x + s * data.at(i).x,
(1.0f - s) * data.at(i - 1).y + s * data.at(i).y));
s += step;
}
}
//  vecp::print(result, "result");
}

# 2、计算重心

3D点计算方法类似。

//************************************
// Description: 将输入的point2d点集求和平均求重心
// Method:    CalcCentre
// FullName:  CalcCentre
// Access:    private
// Parameter: 二维点集 const std::vector<cv::Point2f> &point2ds
// Returns:   重心坐标 cv::Point2f
// Author:
// Date:      2018/08/28
// History:
//************************************
cv::Point2f CalcCentre(const std::vector<cv::Point2f> &point2ds) {
float x = 0, y = 0;
for (auto &pt: point2ds) {
x += pt.x;
y += pt.y;
}
return cv::Point2f(x / point2ds.size(), y / point2ds.size());
}

# 3、闭环判断

//************************************
// Description: 将输入的point2d点集判断是否闭环
// Method:    FindOneCircleFromUnformInCircle
// FullName:  FindOneCircleFromUnformInCircle
// Access:    private
// Parameter: 二维点集 const std::vector<cv::Point2f> &point2ds
// Parameter: 判断闭环坐标中心 const cv::Point2f &centre
// Parameter: 最大无点角度0-360 float max_angle
// Parameter: 划分区域数量 int part_count
// Returns:   bool
// Author:
// Date:      2018/08/28
// History:
//************************************
bool FindOneCircleFromUnformInCircle(const std::vector<cv::Point2f> &point2ds,
const cv::Point2f &centre,
float max_angle, int part_count) {
double each_path_dog = 2.0f * M_PI / part_count;
//cout<<"each_path_dog:"<<each_path_dog<<"  ";
std::vector<int> counts(part_count, 0);

int full_part_counts = 0;
for (auto &pt: point2ds) {
float alpha = atan2Ex(pt.x - centre.x, pt.y - centre.y);
auto idx = static_cast<size_t>(alpha / each_path_dog);
counts.at(idx)++;
}
//每个区域大于5个点算作填满
for (auto &count: counts) {
if (count > 5)
full_part_counts++;
}

double dao_each_path_dog = 1.0 / each_path_dog;
//填满区域大于需要填满的区域则返回true
return full_part_counts  >= int(DegToRad(360.0 - max_angle) * dao_each_path_dog);
}

# 4、椭圆拟合

//************************************
// Description: 将输入的point2d点集做椭圆拟合
// Method:    IsCircleFromEllipse
// FullName:  IsCircleFromEllipse
// Access:    private
// Parameter: 二维点集 const std::vector<cv::Point2f> &point2ds
// Parameter: 最大椭圆度（两半径比值） float min_ovality
// Parameter: 最大非拟合点 float min_error
// Returns:   bool
// Author:
// Date:      2018/08/28
// History:
//************************************
bool IsCircleFromEllipse(const std::vector<cv::Point2f> &point2ds,
float min_radius, float min_ovality, float min_error,
cv::RotatedRect &ellipse) {
//OpenCV椭圆拟合,要求至少有六个点输入
float error = fitEllipseEx(point2ds, ellipse);
float a = ellipse.size.width;
float b = ellipse.size.height;
//printf("a %f, b %f, ovality %f, error %f\n", a, b, a/b, error);
return !(a < min_radius || b < min_radius || a / b < min_ovality || a / b > 1.0 / min_ovality ||
error > min_error);
}

# ５、圆拟合

 /************************************
Description: 将输入的point2d点集做圆拟合,根据圆心计算半径均值和标准差
Method:    IsCircelFromDev
FullName:  IsCircelFromDev
Access:    private
Parameter: 最大半径标准差 float max_dev
Returns:   bool
Author:
Date:      2018/10/17
History:
************************************/
bool CircleAction::IsCircelFromDev(const std::vector<cv::Point2f> &point2ds, float min_radius, float max_dev) {
cv::Point2f centre = CalcCentre(point2ds);
for (size_t i = 0; i < point2ds.size(); i++) {
}

cv::Mat mean, dev;

}

# ６、顺/逆时针方向判断

//************************************
// Description: 将输入的point2d点集判断顺/逆时针方向
// Method:    IsClockwiseFromCross
// FullName:  IsClockwiseFromCross
// Access:    private
// Parameter: 二维点集 const std::vector<cv::Point2f> &point2ds
// Returns:   bool
// Author:
// Date:      2018/08/28
// History:
//************************************
bool IsClockwiseFromCross(const std::vector<cv::Point2f> &point2ds) {
int positive = 0, negative = 0;
for (size_t i = 1; i < point2ds.size() - 1; i++) {
//由i,i-1,i+1得到局部方向
int status = WhichClockWise(point2ds.at(i - 1), point2ds.at(i), point2ds.at(i + 1));
if (status == -1)
negative++;
else if (status == 1)
positive++;
}

//  if (abs(positive - negative) < 5)
//      std::cout << "[warning]: the clockwise may be wrong" << std::endl;
//按正、反方向数量判断返回多的反向
return (positive >= negative);
}
//************************************
// Description: 将输入的point2d点判断局部方向
// Method:    WhichClockWise
// FullName:  WhichClockWise
// Access:    private
// Parameter: 二维点a cv::Point2f a
// Parameter: 二维点b cv::Point2f b
// Parameter: 二维点c cv::Point2f c
// Returns:   int
// Author:
// Date:      2018/08/28
// History:
//************************************
int WhichClockWise(cv::Point2f a, cv::Point2f b, cv::Point2f c) {
double triangle_area = a.x * b.y - a.y * b.x + a.y * c.x - a.x * c.y + b.x * c.y - c.x * b.y;
if (triangle_area < 0) return -1;
else if (triangle_area > 0) return 1;
return 0;
}

2018-01-20 20:14:17 KYJL888 阅读数 2374
• ###### MATLAB图像处理

全面系统的学习MATLAB在图像处理中的应用

20781 人正在学习 去看看 魏伟

class Moments { public: ...... // 空间矩 double m00, m10, m01, m20, m11, m02, m30, m21, m12, m03;// 中心矩double mu20, mu11, mu02, mu30, mu21, mu12, mu03;// 中心归一化矩 double nu20, nu11, nu02, nu30, nu21, nu12, nu03; }空间矩的公式为：

A1,一阶矩就是 E(X),即样本均值。具体说来就是A1=(西格玛Xi)/n ----（1）
A2,二阶矩就是 E(X^2)即样本平方均值 ,具体说来就是 A2=(西格玛Xi^2)/n-----(2)
Ak,K阶矩就是 E(X^k)即样本K次方的均值,具体说来就是 Ak=(西格玛Xi^k)/n，-----(3)

1、在数学领域，矩 非常的常见
2、在计算机视觉中，使用2维离散形式的矩计算方法
3、使用矩，可以计算物体的面积，物体的质心等。
4、中心矩的计算方法是：某个矩除以0阶矩
5、高阶矩具有旋转不变性，尺度不变性，变换不变性等。

M20和M02为二阶矩，接下来计算物体形状的方向

OpenCV代码如下所示：
[cpp] view plain
1. #include <opencv2/opencv.hpp>
2. using namespace cv;
3. using namespace std;
4. Mat src; Mat src_gray;
5. int thresh = 30;
6. int max_thresh = 255;
7. RNG rng(12345);
8. int main(){
10.     cvtColor( src, src_gray, CV_BGR2GRAY );//灰度化
11.     GaussianBlur( src, src, Size(3,3), 0.1, 0, BORDER_DEFAULT );
12.     blur( src_gray, src_gray, Size(3,3) ); //滤波
13.     namedWindow( "image", CV_WINDOW_AUTOSIZE );
14.     imshow( "image", src );
15.     moveWindow("image",20,20);
16.     //定义Canny边缘检测图像
17.     Mat canny_output;
18.     vector<vector<Point> > contours;
19.     vector<Vec4i> hierarchy;
20.     //利用canny算法检测边缘
21.     Canny( src_gray, canny_output, thresh, thresh*3, 3 );
22.     namedWindow( "canny", CV_WINDOW_AUTOSIZE );
23.     imshow( "canny", canny_output );
24.     moveWindow("canny",550,20);
25.     //查找轮廓
26.     findContours( canny_output, contours, hierarchy, CV_RETR_TREE, CV_CHAIN_APPROX_SIMPLE, Point(0, 0) );
27.     //计算轮廓矩
28.     vector<Moments> mu(contours.size() );
29.     forint i = 0; i < contours.size(); i++ )
30.     { mu[i] = moments( contours[i], false ); }
31.     //计算轮廓的质心
32.     vector<Point2f> mc( contours.size() );
33.     forint i = 0; i < contours.size(); i++ )
34.     { mc[i] = Point2d( mu[i].m10/mu[i].m00 , mu[i].m01/mu[i].m00 ); }
35.
36.     //画轮廓及其质心并显示
37.     Mat drawing = Mat::zeros( canny_output.size(), CV_8UC3 );
38.     printf("\t\t 几何特性\n");
39.     forint i = 0; i< contours.size(); i++ )
40.     {
41.         Scalar color = Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
42.         drawContours( drawing, contours, i, color, 2, 8, hierarchy, 0, Point() );
43.         circle( drawing, mc[i], 4, color, -1, 8, 0 );
44.         rectangle(drawing, boundingRect(contours.at(i)), cvScalar(0,255,0));
45.         printf("目标%d - 面积：%.2f - 周长: %.2f ", i, mu[i].m00, contourArea(contours[i]), arcLength( contours[i], true ) );
46.         RotatedRect r = fitEllipse(contours.at(i));
47.         double majorAxis = r.size.height > r.size.width ? r.size.height : r.size.width;//长轴大小
48.         double minorAxis = r.size.height > r.size.width ? r.size.width : r.size.height;//短轴大小
49.         double area = mu[i].m00;//面积
50.         int perimeter = arcLength(contours.at(i), true);
51.         double orientation = r.angle;
53.         printf("- 偏移角度: %.1f\n\n", orientation);
54.         double diameter = sqrt((4*area)/3.1416);//直径
55.         double eccentricity = sqrt(1-pow(minorAxis/majorAxis,2));//离心率
56.         double roundness = pow(perimeter, 2)/(2*3.1416*area);//圆滑度
58.         char tam[100];
59.         sprintf(tam, "%.2f", orientation);
60.         putText(drawing, tam, Point(mc[i].x, mc[i].y), FONT_HERSHEY_SIMPLEX, 0.5, cvScalar(0,220,120),1.5);
61.     }
62.     namedWindow( "Contours", CV_WINDOW_AUTOSIZE );
63.     imshow( "Contours", drawing );
64.     moveWindow("Contours",1100,20);
65.     waitKey(0);
66.     src.release();
67.     src_gray.release();
68.     return 0;
69. }

## 矩、中心矩、质心、patch方向

author@jason_ql
http://blog.csdn.net/lql0716

# １、几何矩理论

## 1.1 矩与数学期望

• 数学期望

E(X)=i=1xiP(xi)

E(X)=+xf(x)dx

E(X)=ixiPX(xi)=jixiP(xi,yj)

E(Y)=jyjPY(yj)=jiyjP(xi,yj)

fX(x)=+f(x,y)dy

E(X)=+xfX(x)dx=+x(+f(x,y)dy)dx=++xf(x,y)dxdy

E(Y)=+yfY(y)dy=++yf(x,y)dxdy

• 原点矩

X是离散型随机变量，则：

νk(X)=ixkip(xi)

X是连续型随机变量，则：

νk(X)=+xkf(x)dx

• 中心距

μk(X)=E(XE(X))k
Xk中心距

X是离散型随机变量，则：

μk(X)=i(xiE(X))kp(xi)

X是连续型随机变量，则：

μk(X)=+(xE(X))kf(x)dx

• 原点矩与中心距

μ2=ν2ν21

μ3=ν33ν2ν1+2ν31

μ4=ν44ν3ν1+6ν2ν213ν41

# 1.2 图像的矩

• 原点矩：

mpq=x1My1Nxpyqf(x,y)

• 中心距：

μpq=x1My1N(xx0)p(yy0)qf(x,y)

• 归一化中心距：

ηpq=μpqμr00

其中r=p+q+22,p+q=2,3,...

• 一阶矩：
见下面第２节．

• 二阶矩：

M20=xyx2I(x,y)

M02=xyy2I(x,y)

M11=xyxyI(x,y)

M20M02分别表示图像围绕通过重心的垂直和水平轴线的惯性矩。
M30M03可以度量图像对于垂直和水平轴线的对称性等。

θ=arctan(b,(ac))2=arctan(b/(ac))2,θ[9090]

a=M20M00C20

b=2(M11M00C0C1)

c=M02M00C21

# 2、质心原理

• moments of a patch(矩):

mpq=x=r,y=rrxpyqI(x,y)(1)

• 角点为中心：

m00=x=r,y=rrx0y0I(x,y)=x=r,y=rrI(x,y)(1-1)

• 一阶矩m01

m01=x=r,y=rrx0y1I(x,y)=x=r,y=rryI(x,y)(1-2)

• 一阶矩m10

m10=x=r,y=rrx1y0I(x,y)=x=r,y=rrxI(x,y)(1-3)

• centroid（质心，亦可称为重心）：

C=(m10m00,m01m00)(2)

计算质心的优势：对噪声不敏感。当有外部噪声干扰时，计算出的质心不会偏离太大。从数学的角度来看，这种方法是计算一个连通域的质心（或一个团块儿blob的质心）。

• 构造一个向量OC,从角点中心O到质心C

• orientation of patch（方向）：

θ=atan2(m01,m10)(3)

建立以角点为圆心的坐标系，如图
在图中，P为角点，园内为取点区域，每个方格代表一个像素。
则质心Q可根据式（2）求得。

# 3、中心距函数moments()

Calculates all of the moments up to the third order of a polygon or rasterized shape.

C++: Moments moments(InputArray array, bool binaryImage=false )

Python: cv2.moments(array[, binaryImage]) → retval

C: void cvMoments(const CvArr* arr, CvMoments* moments, int binary=0 )

Python: cv.Moments(arr, binary=0) → moments

Parameters:
array – Raster image (single-channel, 8-bit or floating-point 2D array) or an array ( 1 \times N or N \times 1 ) of 2D points (Point or Point2f ).
binaryImage – If it is true, all non-zero image pixels are treated as 1’s. The parameter is used for images only.
moments – Output moments.

• The function computes moments, up to the 3rd order, of a vector shape or a rasterized shape. The results are returned in the structure Moments defined as:
class Moments
{
public:
Moments();
Moments(double m00, double m10, double m01, double m20, double m11,
double m02, double m30, double m21, double m12, double m03 );
Moments( const CvMoments& moments );
operator CvMoments() const;

// spatial moments
double  m00, m10, m01, m20, m11, m02, m30, m21, m12, m03;
// central moments
double  mu20, mu11, mu02, mu30, mu21, mu12, mu03;
// central normalized moments
double  nu20, nu11, nu02, nu30, nu21, nu12, nu03;
}12345678910111213141516

# 4、中心矩示例代码

• opencv2.4.13

## 4.1 C++版代码

#include <QCoreApplication>
#include <opencv2/opencv.hpp>
//  Qt Creator 4.2.0(Based on Qt 5.7.1)
//  OpenCV 2.4.13
using namespace cv;
using namespace std;

#define name1 "原图"
#define name2 "效果图"

cv::Mat img, gray;
int nThresh = 100;
int nMaxThresh = 255;
cv::RNG rng(12345);  //产生一个随机数
cv::Mat cannyImg;
std::vector<std::vector<cv::Point>> contours;
std::vector<cv::Vec4i> hierarchy;

//void on_ThreshChange( int, void* ){
//    //canny边缘检测
//    cv::Canny( gray, cannyImg, nThresh, nThresh*2, 3 );
//    //找轮廓
//    cv::findContours( cannyImg, contours, hierarchy, cv::RETR_TREE, cv::CHAIN_APPROX_SIMPLE, cv::Point( 0, 0 ) );
//    //计算矩
//    std::vector<cv::Moments> mu( contours.size() );
//    for(unsigned int i = 0; i < contours.size(); i++){
//        mu[i] = cv::moments( contours[i], false);
//    }
//    //计算中心矩
//    std::vector<cv::Point2f> mc( contours.size() );
//    for( unsigned int i = 0; i < contours.size(); i++ ){
//        mc[i] = cv::Point2f( static_cast<float>(mu[i].m10 / mu[i].m00), static_cast<float>(mu[i].m01 / mu[i].m00));
//    }
//    //画轮廓
//    cv::Mat drawing = cv::Mat::zeros( cannyImg.size(), CV_8UC3);
//    for( unsigned int i = 0; i < contours.size(); i++ ){
//        cv::Scalar color = cv::Scalar( rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255) );
//        cv::drawContours( drawing, contours, i, color, 2, 8, hierarchy, 0, cv::Point() );
//        cv::circle( drawing, mc[i], 4, color, -1, 8, 0 );
//    }

//    cv::namedWindow( name2, cv::WINDOW_NORMAL);
//    cv::imshow( name2, drawing );

//    std::cout << "输出内容：　面积和轮廓长度　\n" << std::endl;
//    for(unsigned int i = 0; i < contours.size(); i++ ){
//        std::cout << ">通过ｍ00计算出轮廓[" << i << "]的面积：(M_00) =" << mu[i].m00 << "\n OpenCV 函数计算出的面积　＝　" << cv::contourArea(contours[i]) << "长度：" << cv::arcLength( contours[i], true) <<  "\n\n" << std::endl;
//        cv::Scalar color = cv::Scalar( rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255) );
//        cv::drawContours(drawing, contours, i, color, 2, 8, hierarchy, 0, cv::Point() );
//        cv::circle( drawing, mc[i], 4, color, -1, 8, 0 );
//     }
//}

int main(){
cv::cvtColor( img, gray, cv::COLOR_BGR2GRAY );
cv::blur( gray, gray, cv::Size(3, 3) );

cv::namedWindow( name1, cv::WINDOW_NORMAL );
cv::imshow( name1, img );

//    cv::createTrackbar( "阈值", name1, &nThresh, nMaxThresh, on_ThreshChange );
//    on_ThreshChange( 0, 0 );

//canny边缘检测
cv::Canny( gray, cannyImg, nThresh, nThresh*2, 3 );
//找轮廓
cv::findContours( cannyImg, contours, hierarchy, cv::RETR_TREE, cv::CHAIN_APPROX_SIMPLE, cv::Point( 0, 0 ) );
//计算矩
std::vector<cv::Moments> mu( contours.size() );
for(unsigned int i = 0; i < contours.size(); i++){
mu[i] = cv::moments( contours[i], false);
}
//计算中心矩
std::vector<cv::Point2f> mc( contours.size() );
for( unsigned int i = 0; i < contours.size(); i++ ){
mc[i] = cv::Point2f( static_cast<float>(mu[i].m10 / mu[i].m00), static_cast<float>(mu[i].m01 / mu[i].m00));
}
//画轮廓
cv::Mat drawing = cv::Mat::zeros( cannyImg.size(), CV_8UC3);
for( unsigned int i = 0; i < contours.size(); i++ ){
cv::Scalar color = cv::Scalar( rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255) );
cv::drawContours( drawing, contours, i, color, 2, 8, hierarchy, 0, cv::Point() );
cv::circle( drawing, mc[i], 4, color, -1, 8, 0 );
}

cv::namedWindow( name2, cv::WINDOW_NORMAL);
cv::imshow( name2, drawing );

std::cout << "输出内容：　面积和轮廓长度　\n" << std::endl;
for(unsigned int i = 0; i < contours.size(); i++ ){
std::cout << ">通过ｍ00计算出轮廓[" << i << "]的面积：(M_00) =" << mu[i].m00 << "\n OpenCV 函数计算出的面积　＝　" << cv::contourArea(contours[i]) << "长度：" << cv::arcLength( contours[i], true) <<  "\n\n" << std::endl;
cv::Scalar color = cv::Scalar( rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255) );
cv::drawContours(drawing, contours, i, color, 2, 8, hierarchy, 0, cv::Point() );
cv::circle( drawing, mc[i], 4, color, -1, 8, 0 );
}
cv::waitKey(0);

return 0;
}

原图： 

## 4.2 Python版代码

# -*- coding: utf-8 -*-
"""
Created on Sun Mar 26 18:36:19 2017

@author: lql0716
"""

import cv2
import numpy as np

nThresh = 100
nMaxThresh = 255

gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
gray = cv2.blur(gray, (3,3))
cv2.namedWindow('img', cv2.WINDOW_NORMAL)
cv2.imshow('img', img)

cannyImg = cv2.Canny(gray, nThresh, nThresh*2, 3)
contours, hierarchy = cv2.findContours(cannyImg, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)

mu = []
mc = []
retval = np.array([])
for i in range(0, np.array(contours).shape[0]):
retval = cv2.moments(contours[i], False)
mu.append(retval)
mu = np.array(mu)
#print mu[0]['m10']
thetas = []
for i in range(0, np.array(contours).shape[0]):
if mu[i]['m00'] == 0.0:
a=0
b=0
else:
a = mu[i]['m10'] / mu[i]['m00']  #质心x坐标
b = mu[i]['m01'] / mu[i]['m00']  #质心y坐标

#根据二阶矩计算物体形状的方向
r1 = mu[i]['m20'] / mu[i]['m00'] - a*a
r2 = 2.0*(mu[i]['m11'] / mu[i]['m00'] - a*b)
r3 = mu[i]['m02'] / mu[i]['m00'] - b*b
#        print r1-r3
if r1-r3==0:
theta = np.pi / 2
else:
theta = np.arctan(r2/(r1-r3)) / 2
thetas.append(theta)
mc.append([a,b])
mc = np.array(mc)
drawing = np.zeros(img.shape, dtype = np.uint8)
for i in range(0, mc.shape[0]):
c1 = np.random.randint(0, 256)
c2 = np.random.randint(0, 256)
c3 = np.random.randint(0, 256)
cv2.drawContours(drawing, contours, i, (c1, c2, c3), 2, 8)
cv2.circle(drawing, (int(round(mc[i][0])), int(round(mc[i][1]))), 4, (c1, c2, c3), -1, 8, 0)
cv2.namedWindow('img2', cv2.WINDOW_NORMAL)
cv2.imshow('img2', drawing)
cv2.waitKey(0) 


# 5、Hu矩HuMoments()

mpq=x1My1Nxpyqf(x,y)

μpq=x1My1N(xx0)p(yy0)qf(x,y)

ηpq=μpqμr00

Hu矩利用二阶、三阶中心距构造了7个不变矩，它们在连续图像条件下可保持平移、旋转、伸缩不变，公式如下：
M1=η20+η02
M2=(η20η02)2+4η211
M3=(η303η12)2+(3η21η03)2
M4=(η30+η12)2+(η21+η03)2
M5=(η303η12)(η30+η12)((η30+η12)23(η21+η03)2)+(3η21η03)(η21+η03)(3(η30+η12)2(η21+η03