opencv4 cpp曲面拟合填补之最小二乘法拟合三次曲面（含代码）

void surfaceFitting3(IN cv::Mat mask, IN cv::Mat height, OUT cv::Mat & factor) const
{
	if (!(mask.depth() == CV_8UC1 && height.depth() == CV_16SC1))return;
	Mat left_hand = Mat::zeros(Size(10, 10), CV_64FC1);
	Mat right_hand = Mat::zeros(Size(1, 10), CV_64FC1);
	short z;
	for (int i = 0;i < mask.rows;++i) {
		for (int j = 0;j < mask.cols;++j) {
			if (mask.ptr<uchar>(i)[j] == 255) {
				z = height.ptr<short>(i)[j];
				left_hand.ptr<double>(0)[0] += pow(i, 6);
				left_hand.ptr<double>(1)[0] = left_hand.ptr<double>(0)[1] += pow(i, 3)*pow(j, 3);
				left_hand.ptr<double>(2)[0] = left_hand.ptr<double>(0)[2] += pow(i, 5)*j;
				left_hand.ptr<double>(3)[0] = left_hand.ptr<double>(0)[3] += pow(i, 4)*pow(j, 2);
				left_hand.ptr<double>(4)[0] = left_hand.ptr<double>(0)[4] += pow(i, 5);
				left_hand.ptr<double>(5)[0] = left_hand.ptr<double>(0)[5] += pow(i, 3)*pow(j, 2);
				left_hand.ptr<double>(6)[0] = left_hand.ptr<double>(0)[6] += pow(i, 4)*j;
				left_hand.ptr<double>(7)[0] = left_hand.ptr<double>(0)[7] += pow(i, 4);
				left_hand.ptr<double>(8)[0] = left_hand.ptr<double>(0)[8] += pow(i, 3)*j;
				left_hand.ptr<double>(9)[0] = left_hand.ptr<double>(0)[9] += pow(i, 3);
				left_hand.ptr<double>(1)[1] += pow(j, 6);
				left_hand.ptr<double>(2)[1] = left_hand.ptr<double>(1)[2] += pow(i, 2)*pow(j, 4);
				left_hand.ptr<double>(3)[1] = left_hand.ptr<double>(1)[3] += i*pow(j, 5);
				left_hand.ptr<double>(4)[1] = left_hand.ptr<double>(1)[4] += pow(i, 2)*pow(j, 3);
				left_hand.ptr<double>(5)[1] = left_hand.ptr<double>(1)[5] += pow(j, 5);
				left_hand.ptr<double>(6)[1] = left_hand.ptr<double>(1)[6] += i*pow(j, 4);
				left_hand.ptr<double>(7)[1] = left_hand.ptr<double>(1)[7] += i*pow(j, 3);
				left_hand.ptr<double>(8)[1] = left_hand.ptr<double>(1)[8] += pow(j, 4);
				left_hand.ptr<double>(9)[1] = left_hand.ptr<double>(1)[9] += pow(j, 3);
				left_hand.ptr<double>(2)[2] += pow(i, 4)*pow(j, 2);
				left_hand.ptr<double>(3)[2] = left_hand.ptr<double>(2)[3] += pow(i, 3)*pow(j, 3);
				left_hand.ptr<double>(4)[2] = left_hand.ptr<double>(2)[4] += pow(i, 4)*j;
				left_hand.ptr<double>(5)[2] = left_hand.ptr<double>(2)[5] += pow(i, 2)*pow(j, 3);
				left_hand.ptr<double>(6)[2] = left_hand.ptr<double>(2)[6] += pow(i, 3)*pow(j, 2);
				left_hand.ptr<double>(7)[2] = left_hand.ptr<double>(2)[7] += pow(i, 3)*j;
				left_hand.ptr<double>(8)[2] = left_hand.ptr<double>(2)[8] += pow(i, 2)*pow(j, 2);
				left_hand.ptr<double>(9)[2] = left_hand.ptr<double>(2)[9] += pow(i, 2)*j;
				left_hand.ptr<double>(3)[3] += pow(i, 2)*pow(j, 4);
				left_hand.ptr<double>(4)[3] = left_hand.ptr<double>(3)[4] += pow(i, 3)*pow(j, 2);
				left_hand.ptr<double>(5)[3] = left_hand.ptr<double>(3)[5] += i*pow(j, 4);
				left_hand.ptr<double>(6)[3] = left_hand.ptr<double>(3)[6] += pow(i, 2)*pow(j, 3);
				left_hand.ptr<double>(7)[3] = left_hand.ptr<double>(3)[7] += pow(i, 2)*pow(j, 2);
				left_hand.ptr<double>(8)[3] = left_hand.ptr<double>(3)[8] += i*pow(j, 3);
				left_hand.ptr<double>(9)[3] = left_hand.ptr<double>(3)[9] += i*pow(j, 2);
				left_hand.ptr<double>(4)[4] += pow(i, 4);
				left_hand.ptr<double>(5)[4] = left_hand.ptr<double>(4)[5] += pow(i, 2)*pow(j, 2);
				left_hand.ptr<double>(6)[4] = left_hand.ptr<double>(4)[6] += pow(i, 3)*j;
				left_hand.ptr<double>(7)[4] = left_hand.ptr<double>(4)[7] += pow(i, 3);
				left_hand.ptr<double>(8)[4] = left_hand.ptr<double>(4)[8] += pow(i, 2)*j;
				left_hand.ptr<double>(9)[4] = left_hand.ptr<double>(4)[9] += pow(i, 2);
				left_hand.ptr<double>(5)[5] += pow(j, 4);
				left_hand.ptr<double>(6)[5] = left_hand.ptr<double>(5)[6] += i*pow(j, 3);
				left_hand.ptr<double>(7)[5] = left_hand.ptr<double>(5)[7] += i*pow(j, 2);
				left_hand.ptr<double>(8)[5] = left_hand.ptr<double>(5)[8] += pow(j, 3);
				left_hand.ptr<double>(9)[5] = left_hand.ptr<double>(5)[9] += pow(j, 2);
				left_hand.ptr<double>(6)[6] += pow(i, 2)*pow(j, 2);
				left_hand.ptr<double>(7)[6] = left_hand.ptr<double>(6)[7] += pow(i, 2)*j;
				left_hand.ptr<double>(8)[6] = left_hand.ptr<double>(6)[8] += i*pow(j, 2);
				left_hand.ptr<double>(9)[6] = left_hand.ptr<double>(6)[9] += i*j;
				left_hand.ptr<double>(7)[7] += pow(i, 2);
				left_hand.ptr<double>(8)[7] = left_hand.ptr<double>(7)[8] += i*j;
				left_hand.ptr<double>(9)[7] = left_hand.ptr<double>(7)[9] += i;
				left_hand.ptr<double>(8)[8] += pow(j, 2);
				left_hand.ptr<double>(9)[8] = left_hand.ptr<double>(8)[9] += j;
				left_hand.ptr<double>(9)[9] += 1;
				right_hand.ptr<double>(0)[0] += pow(i, 3)*z;
				right_hand.ptr<double>(1)[0] += pow(j, 3)*z;
				right_hand.ptr<double>(2)[0] += pow(i, 2)*j*z;
				right_hand.ptr<double>(3)[0] += i*pow(j, 2)*z;
				right_hand.ptr<double>(4)[0] += pow(i, 2)*z;
				right_hand.ptr<double>(5)[0] += pow(j, 2)*z;
				right_hand.ptr<double>(6)[0] += i*j*z;
				right_hand.ptr<double>(7)[0] += i*z;
				right_hand.ptr<double>(8)[0] += j*z;
				right_hand.ptr<double>(9)[0] += z;
			}
		}
	}
	solve(left_hand, right_hand, factor);
}

#define DEF_DGM_I(x,y,a,b,c,d,e,f,g,h,i,j) a*x*x*x+b*y*y*y+c*x*x*y+d*x*y*y+e*x*x+f*y*y+g*x*y+h*x+i*y+j;
void mendSurface3(IN cv::Mat mask, IN cv::Mat factor, OUT cv::Mat & height) const
{
	if (!(mask.depth() == CV_8UC1 && factor.depth() == CV_64FC1 && height.depth() == CV_16SC1
		&& factor.size() == Size(1, 10) && mask.size() == height.size()))return;
	double a = factor.ptr<double>(0)[0], b = factor.ptr<double>(1)[0],
		c = factor.ptr<double>(2)[0], d = factor.ptr<double>(3)[0],
		e = factor.ptr<double>(4)[0], f = factor.ptr<double>(5)[0],
		g = factor.ptr<double>(6)[0], h = factor.ptr<double>(7)[0],
		i = factor.ptr<double>(8)[0], j = factor.ptr<double>(9)[0];
	for (int _row = 0;_row < mask.rows;++_row) {
		for (int _col = 0;_col < mask.cols;++_col) {
			if (mask.ptr<uchar>(_row)[_col] == 255) {
				double hT = DEF_DGM_I(_row, _col, a, b, c, d, e, f, g, h, i, j);
				height.ptr<short>(_row)[_col] = static_cast<short>(hT);
			}
		}
	}
}