2017-04-10 22:58:50 qq826309057 阅读数 13179

数学中的矩

维基百科关于 Moment (mathematics) 的介绍:

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.

设 X 和 Y 是随机变量,c 为常数,k 为正整数,
如果E(|Xc|k)存在,则称E(|Xc|k)为 X 关于点 c 的 k 阶矩。

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

如果E(|Xc1|p|Yc2|q)存在,则称其为 X,Y 关于 c 点 p+q 阶矩。


图像的矩

零阶矩:

M00=IJV(i,j)

这里的图像是单通道图像,Vij表示图像在(i,j)点上的灰度值。
我们可以发现,当图像为二值图时,M00就是这个图像上白色区域的总和,因此,M00可以用来求二值图像(轮廓,连通域)的面积。

一阶矩:

M10=IJiV(i,j)M01=IJjV(i,j)

当图像为二值图时,V(i,j) 只有0(黑),1(白)两个值。M10 就是图像上所以白色区域 x坐标值的累加。因此,一阶矩可以用来求二值图像的重心:
xc=M10M00,yc=M01M00

二阶矩

M20=IJi2V(i,j)M02=IJj2V(i,j)M11=IJijV(i,j)

二阶矩可以用来求物体形状的方向。
θ=12fastAtan2(2b,ac)

其中:a=M20M00xc2,b=M11M00xcycc=M02M00yc2

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


这里修改一下,我之前在看二阶矩求物体形状方向的时候,有公式是这么写的:

θ=12arctan(2bac)

我当时很奇怪,arctan的范围是(π2,π2),再除以2就只有(π4,π4),但是物体形状的方向肯定不止这个范围,而且我在测试的时候也感觉不对,因此就写成了上面那个用opencv函数表示的公式。
其实,这个公式是对的,只不过需要再根据分子分母的正负性来判断一下象限范围,而fastAtan2()已经加入了这个功能。
顺便说一下 math.h 中atan()求出来的是弧度,而fastAtan2()求出来的是角度。


利用矩求图像的重心,方向

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;//此为形状的方向

图片测试:
这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述
这里写图片描述
从上图可以发现,这样算出来的方向是以x轴正方向为0,第四象限为0π/2,第三象限为π/2π 的范围内。


总结

图像的矩通常描述了该图像的全局特征,并提供了大量的关于该图像不同类型的几何特性信息,比如大小,位置,方向及形状等。在模式识别,目标分类,目标识别等有着广泛的应用。

2016-12-08 11:41:15 m_buddy 阅读数 4103

1. 概述

最近在处理图像配准方面的任务,初步接触医学图像配准,这篇文章讲述基于图像原始像素信息的简单配准方法。通过计算得到图像的重心、主轴方向使得图像能够定位到统一坐标系下。这种方法简单粗暴,不能实现精细的配准
图像重心的计算公式:

图像主轴的计算方法:

图像的主轴方向定义为协方差矩阵最大特征值对应的特征向量

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;
}
调用:(src_img为参考的模板图像, img_vec为待配准的数据集;调用中使用到了前一篇博客中提到的基础集合变换方法,需要的朋友请查看之前的文章)
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");

3. 结果

左边为未配准直接叠加的图像,右边为经过处理之后得到的结果,红色的层为CT的图像(提供轮廓信息),绿色的层为MRI图像(提供脑部细节信息)

4. 总结

以上给出的实现方法,相当简单粗暴,自然效果也不好。配准仅仅针对像素信息还有许多的方法,互信息也是很不错的度量手段。


2010-08-24 11:41:00 asmemgsd 阅读数 376

今天在公司总结了下这:

面试比较注重的内容:

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

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

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

得到几点重要的东西:

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

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

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

2018-08-29 20:57:20 aaron19890330 阅读数 1512

场景:前装摄像头

检测动作:单手指画圈,需要判断画圈方向和圈数。

     

步骤: (1)取得3D图像序列最前点;

            (2)将最前点投影在2D平面上;

            (3)中值滤波和平滑处理;

            (4)得到2D点集进行线性插值;

            (5)以重心为中心判断是否闭环;

            (6)2D点集进行椭圆拟合,判断是否椭圆、椭圆半径和椭圆度范围;

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

            (8)判断2D点集方向:顺/逆时针。

其中手的识别使用深度学习进行,根据TensorFlow训练完的模型判断出图像中有手且为单手指,
去掉背景、挖出只含有手的图像;中值滤波用的OpenCV medianBlur函数,其他主要函数算法如下:


Table of Contents

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

2、计算重心

3、闭环判断

4、椭圆拟合

5、圆拟合

6、顺/逆时针方向判断


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、椭圆拟合

使用了OpenCV的fitEllipseEx函数。

//************************************
// Description: 将输入的point2d点集做椭圆拟合
// Method:    IsCircleFromEllipse
// FullName:  IsCircleFromEllipse
// Access:    private
// Parameter: 二维点集 const std::vector<cv::Point2f> &point2ds
// Parameter: 最小半径 float min_radius
// 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);
}

5、圆拟合

以重心为圆心设置半径均值和标准差阈值:

 /************************************
 Description: 将输入的point2d点集做圆拟合,根据圆心计算半径均值和标准差
 Method:    IsCircelFromDev
 FullName:  IsCircelFromDev
 Access:    private
 Parameter: 最小半径均值 float min_radius
 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);
    std::vector<float> radius(point2ds.size(), 0.0f);
    for (size_t i = 0; i < point2ds.size(); i++) {
        radius.at(i) = Distance(centre, point2ds.at(i));
    }

    cv::Mat mean, dev;
    cv::meanStdDev(radius, mean, dev);

    double mean_radius = mean.at<double>(0, 0), dev_radius = dev.at<double>(0, 0);
    //std::cout << "radius mean---dev: " <<mean_radius << "dev_radius: " << dev_radius << std::endl;
    return (mean_radius > min_radius && dev_radius < max_dev);
}

6、顺/逆时针方向判断

//************************************
// 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

参考 http://blog.csdn.net/yang6464158/article/details/42459595

特征矩的知识在概率论和数理统计中有介绍,空间矩的方法在图像应用中比较广泛,包括零阶矩求面积、一阶矩确定重心、二阶矩确定主方向、二阶矩和三阶矩可以推导出七个不变矩Hu不变矩,不变矩具有旋转,平移、缩放等不变性,因此在工业应用和模式识别中得到广泛的应用。

目标物体灰度函数特征矩的公式定义如下:

如果是二值图像,那么f(x,y)就变成

在OpenCV中,可以很方便的计算多边形区域的3阶特征矩,opencv中的矩主要包括以下几种:空间矩,中心矩和中心归一化矩。

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; }空间矩的公式为:

image

可以知道,对于01二值化的图像,m00即为轮廓的面积。中心矩的公式为:

image

其中:

image

归一化的中心矩公式为:

image
参考于链接

二阶中心距,也叫作方差,它告诉我们一个随机变量在它均值附近波动的大小,方差越大,波动性越大。方差也相当于机械运动中以重心为转轴的转动惯量。(The moment of inertia.)

三阶中心距告诉我们一个随机密度函数向左或向右偏斜的程度。

在均值不为零的情况下,原点距只有纯数学意义。

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)

不变矩的物理含义:
如果把图像看成是一块质量密度不均匀的薄板,其图像的灰度分布函数f(x,y)就是薄板的密度分布函数,则其各阶矩有着不同的含义,如零阶矩表示它的总质量;一阶矩表示它的质心;二阶矩又叫惯性矩,表示图像的大小和方向。事实上,如果仅考虑阶次为2的矩集,则原始图像等同于一个具有确定的大小、方向和离心率,以图像质心为中心且具有恒定辐射率的椭圆。由三阶矩以下矩构成的七个矩不变量具有平移、旋转和尺度不变性等等。当密度分布函数发生改变时,图像的实质没有改变,仍然可以看做一个薄板,只是密度分布有所改变。虽然此时各阶矩的值可能发生变化,但由各阶矩计算出的不变矩仍具有平移、旋转和尺度不变性。通过这个思想,可对图像进行简化处理,保留最能反映目标特性的信息,再用简化后的图像计算不变矩特征,可减少计算量。

研究表明,只有基于二阶矩的不变矩对二维物体的描述才是真正的与旋转、平移和尺度无关的。较高阶的矩对于成像过程中的误差,微小的变形等因素非常敏感,所以相应的不变矩基本上不能用于有效的物体识别。即使是基于二阶矩的不变矩也只能用来识别外形相差特别大的物理,否则他们的不变矩会因为很相似而不能识别。

在OpenCV中,还可以很方便的得到Hu不变距,Hu不变矩在图像旋转、缩放、平移等操作后,仍能保持矩的不变性,所以有时候用Hu不变距更能识别图像的特征。

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


我们很熟悉概率论中的一阶矩二阶矩高阶矩,但是很多人可能和我一样,不明白图像中矩是拿来干嘛的。

在计算机视觉的书中,虽然有提到矩,但是讲的很泛泛也很笼统。自然Google百度这些东西也是靠不牢的。在阅读了相关论文之后,我终于大致对矩在图像中的应用有了了解。
其实矩除了在概率论中有体现,在几何中也是学过的。比方说零阶矩是物体的质量,一阶矩和零阶矩可以算出物体的中心,而二阶矩是用来计算物体的方向的。
拿图像出来来说,图像可以看成是一个平板的物体,其一阶矩和零阶矩就可以拿来计算某个形状的重心,而二阶矩就可以拿来计算形状的方向。
光说不练假把式,下面给出他们的计算公式:

其中M00即零阶矩


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



OpenCV代码如下所示:
[cpp] view plain copy
  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(){  
  9.     src = imread( "opencv-logo.png" ,CV_LOAD_IMAGE_COLOR );  
  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;  
  52.         double orientation_rads = orientation*3.1416/180;  
  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);//圆滑度  
  57.         line(drawing, Point(mc[i].x, mc[i].y), Point(mc[i].x+30*cos(orientation_rads), mc[i].y+30*sin(orientation_rads)), cvScalar(0,0,200), 3);  
  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.1 矩与数学期望

  • 数学期望

定义(一维离散):设X[a,b],密度为f(x),数学期望为: 

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


定义(一维连续):设X为连续型随机变量,其概率密度为f(x),则X的数学期望为: 

E(X)=+xf(x)dx

注:假定广义积分绝对收敛,即+|x|f(x)dx存在


定义(二维离散):对于离散变量(X,Y)P(xi,yi)PX(xi)=jP(xi,yj)期望为:

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

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


定义(二维连续):连续变量(X,Y)f(x,y)

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

  • 原点矩

定义1:设X是随机变量,则称νk(X)=E(Xk)Xk原点矩

X是离散型随机变量,则:

νk(X)=ixkip(xi)

X是连续型随机变量,则:

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

  • 中心距

定义2:设X是随机变量,则称 

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

X是离散型随机变量,则: 

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

X是连续型随机变量,则: 

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

  • 原点矩与中心距

当中心距中的E(X)为0时,此时为k阶原点矩,即原点矩是中心距的特殊情况。

一阶原点矩就是数学期望,二阶中心距就是方差,在实际中常用低阶矩,高于四阶矩极少使用。

原点矩与中心距的关系式: 

μ2=ν2ν21

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

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

以上可对μr用组合数拆开得到。

1.2 图像的矩

把图像的像素看做密度函数f(x,y),对该像素点求期望E,即是图像的矩(原点矩)。具体的求解过程参看下面第2节。

一般来说,一阶矩零阶矩可以计算某个形状的重心,二阶矩可以计算形状的方向。

图像的矩主要表征了图像区域的几何特征,又称几何矩,由于具有旋转、平移、尺度等不变的特兴奋,所以又称为不变矩。 
利用不变矩可以计算出物体的圆形度(物体形状和园的接近程度)、物体的矩形度(物体形状和矩形的接近程度)、物体的水平和垂直对称性、物体的主轴方向、扁度等。

  • 原点矩: 

    mpq=x1My1Nxpyqf(x,y)

  • 中心距: 

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

  • 归一化中心距: 

    ηpq=μpqμr00

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

  • 一阶矩: 
    见下面第2节.

  • 二阶矩:

M20=xyx2I(x,y)

M02=xyy2I(x,y)

M11=xyxyI(x,y)

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

物体形状的方向: 

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

其中:

根据一阶矩的质心C=(M10M00,M01M00)

a=M20M00C20

b=2(M11M00C0C1)

c=M02M00C21

2、质心原理

在图像处理中,一阶矩与形状有关,二阶矩显示曲线围绕直线平均值到扩展程度,三阶矩是关于平均值到对称性到测量.由二阶矩和三阶矩可以导出一组共7个不变矩.而不变矩是图像到统计特性,满足平移,伸缩,旋转均不变到不变性.

  • 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;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16

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 << ">通过m00计算出轮廓[" << 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(){
    img = cv::imread( "/home/jason/jason2/photo/1.jpg" );
    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 << ">通过m00计算出轮廓[" << 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

img = cv2.imread('D:/photo/04.jpg')
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

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

当图像变化时,mpq也变化,而μpq具有平移不变形,但对旋转依然敏感。如果用归一化中心距,则可以保持平移不变性、伸缩不变性。

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