移动最小二乘法MLS基于C++代码实现

这是MLS基于C++代码实现，矩阵类运算再我的另一个资源中。代码仅供学习的，不能直接运行，需要调用。

#include 
#include
#include
#include “Matrix.h“

#include
#include
#define M 8//采样点的个数
#define N 6//基数个数，线性基三个，二次基六个，三次基十个

using namespace std;
using namespace cv;
int max_x=0;
int min_x=0;

typedef struct point
{
	double x;
	double y;
}point;
typedef struct Lnode
{
	int data;
	struct Lnode *next;
}Lnode*linklist;
extern void Beizer（vectornklist>&p_M）;
Matrix_Calc m_c;//定义Matrix_Calc对象

//将m行N列数组转换为m行N列的矩阵
Matrix Trans_Matrix（float p[][N]int m）
{
	Matrix C（mN）;
	C.init_Matrix（）;
	for（int i=0;i		for（int j=0;j		{
			C.write（ijp[i][j]）;
		}
	return C;
}
//将1行N列数组转换为1行N列的矩阵
Matrix Trans_Matrix_One（float p[]int n）
{
	Matrix C（1n）;
	C.init_Matrix（）;
	for（int i=0;i		for（int j=0;j		{
			C.write（ijp[j]）;
		}
	return C;
}
//计算A矩阵个元素的函数
void f（float w[]float x[]float y[]float sumf[][N]float p[][N]）
{
	for（int i=0;i		for（int j=0;j		{
			sumf[i][0] = sumf[i][0]+w[j]*1*p[j][i];
			sumf[i][1] = sumf[i][1] + w[j] * x[j] * p[j][i];
			sumf[i][2] = sumf[i][2] + w[j] * y[j]* p[j][i];
			sumf[i][3] = sumf[i][3] + w[j] * x[j] * x[j]* p[j][i];
			sumf[i][4] = sumf[i][4] + w[j] * x[j] * y[j]* p[j][i];
			sumf[i][5] = sumf[i][5] + w[j] * y[j] * y[j]* p[j][i];
		}
}
//移动最小二乘法的具体计算过程，参照论文“基于移动最小二乘法的曲线曲面拟合”，AB矩阵参照论文“移动最小二乘法的研究”
int MLS_Calc（int x_valint y_valfloat x[]float y[]float z[]）
{
	int max_delta=max_x-min_x;//区域半径
	float p[M][N]={0};
	float sumf[N][N]={0};
	float w[M]={0};
    for（int j=0;j	{
        float s=fabs（（x[j]-x_val））/max_delta;
        if（s<=0.5）
            w[j]=2/3.0-4*s*s+4*s*s*s;
        else
		{
            if（s<=1）
                w[j]=4/3.0-4*s+4*s*s-4*s*s*s/3.0;
            else
                w[j]=0;
		}
		p[j][0]=1;//每个采样点计算基函数
		p[j][1]=x[j];
		p[j][2]=y[j];
		p[j][3]=x[j]*x[j];
		p[j][4]=x[j]*y[j];
		p[j][5]=y[j]*y[j];
	}
 	f（wxysumfp）;//计算得出A矩阵
	
	float p1[N];
	Matrix A=Trans_Matrix（sumfN）;
	Matrix A_1=m_c.Matrix_copy（&A）;
	m_c.Matrix_inv（&A_1）;//求A矩阵的逆A_1

	Matrix B（N1）;//求矩阵B，N行M列
	B.init_Matrix（）;
	for（int j=0;j	{
		p1[0]=1*w[j];
		p1[1]=x[j]*w[j];
		p1[2]=y[j]*w[j];
		p1[3]=x[j]*x[j]*w[j];
		p1[4]=x[j]*y[j]*w[j];
		p1[5]=y[j]*y[j]*w[j];
		Matrix P=Trans_Matrix_One（p1N）;//数组P1转成1行N列的P矩阵
		if（j==0）//第一列直接赋值
		{
			for（int i=0;i				B.write（i0p1[i]）;
		}
		else
		{
			m_c.Matrix_trans（&P）;//矩阵转置，P转为N行1列矩阵
			m_c.Matrix_addCols（&B&P）;//矩阵B列附加，形成N行M列矩阵
		}
		P.free_Matrix（）;
	}

	float D[N]={1x_valy_valx_val*x_valx_val*y_valy_val*y_val};
	Matrix D1=Trans_Matrix_One（DN）;//转成1行N列矩阵

	Matrix D_A1_mul（1N）;//定义矩阵并初始化相乘的结果矩阵，1行N列
	D_A1_mul.init_Matrix（）;
	if（m_c.Matrix_mul（&D1&A_1&D_A1_mul）==-1）
		cout<<“矩阵有误1！“;//1行N列矩阵乘以N行N列矩阵得到结果为1行N列
	
	Matrix D_A1_B_mul（1M）;//定义矩阵并初始化相乘的结果矩阵，1行M列
	D_A1_

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        6649  2017-05-02 14:21  MLS\MLS.cpp
     文件         268  2017-05-03 22:02  MLS\代码说明.txt