VTK解决曲面拟合问题

1. 将问题简化为拟合二次曲面函数
   $$z=a{x}^2+b{y}^2+cxy+dx+ey+f$$
   写成矩阵乘法形式
   $$z=\left[\begin{array}{cccccc}{x}^2& y^2& xy& x& y& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\\ d\\ e\end{array}\right]$$
   记$M=\left[\begin{array}{c}a\\ b\\ c\\ d\\ e\end{array}\right]$, $X=\left[\begin{array}{cccccc}{x}^2& y^2& xy& x& y& 1\end{array}\right]$,$Y=z$，
   则$M=Y{X}^{-1}$
2. 最小二乘法求解线性方程组

#include <vtkMath.h>

namespace
{
	/* allocate memory for an nrow x ncol matrix */
	template <class TReal> TReal** create_matrix(long nrow, long ncol)
	{
		typedef TReal* TRealPointer;
		TReal** m = new TRealPointer[nrow];

		TReal* block = (TReal*)calloc(nrow * ncol, sizeof(TReal));
		m[0] = block;
		for (int row = 1; row < nrow; ++row)
		{
			m[row] = &block[row * ncol];
		}
		return m;
	}

	/* free a TReal matrix allocated with matrix() */
	template <class TReal> void free_matrix(TReal** m)
	{
		free(m[0]);
		delete[] m;
	}
} // namespace


int main(int, char*[])
{
	srand(42);
	// 拟合平面
	// z = ax^2 + by^2 + cxy + dx + ey + f
	// 等价于
	// Solve XM = Y;
	// [x^2 y^2 xy x y 1][a b c d e f].T = z
	int numberOfSamples = 50;
	int numberOfVariables = 6;
	double** x = create_matrix<double>(numberOfSamples, numberOfVariables);
	double** y = create_matrix<double>(numberOfSamples, 1);
	for (int i = 0; i < 50; i++)
	{
		// [a b c d e f] = [1,2,3,4,5,6]
		auto rx = rand() % 100;
		auto ry = rand() % 100;
		double noise = 1.0 / sqrt(6.28)*exp(-0.5*(rx*rx));
		x[i][0] = rx * rx;
		x[i][1] = ry * ry;
		x[i][2] = rx * ry;
		x[i][3] = rx;
		x[i][4] = ry;
		x[i][5] = 1;

		y[i][0] = rx * rx + 2 * ry * ry + 3 * rx * ry + 4 * rx + 5 * ry + 6 + noise;
	}


	double** m = create_matrix<double>(numberOfVariables, 1);
	vtkMath::SolveLeastSquares(numberOfSamples, x, numberOfVariables, y, 1, m);

	std::cout << "Solution is: "
		<< m[0][0] << " " << m[1][0] << " "
		<< m[2][0] << " " << m[3][0] << " "
		<< m[4][0] << " " << m[5][0] << std::endl;

	free_matrix(x);
	free_matrix(m);
	free_matrix(y);

	return EXIT_SUCCESS;
}

运行结果

[1]. https://kitware.github.io/vtk-examples/site/Cxx/Math/LeastSquares/
[2]. https://kitware.github.io/vtk-examples/site/Cxx/Math/HomogeneousLeastSquares/