直角点的自动增补与直角化处理（数字摄影测量）

直角点的自动增补与直角化处理

数字摄影测量关于直角点的自动增补与直角化处理基于MFC的实现
核心代码：

std::vector<CPoint> temp;

void CzhijiaodianView::OnLButtonDown(UINT nFlags, CPoint point)
{
	// TODO: 在此添加消息处理程序代码和/或调用默认值


	CDC* pDC = GetDC();
	CString str;
	str.Format(TEXT("%d,%d"), point.x, point.y);
	pDC->FillSolidRect(0, 0, 100, 100, GetSysColor(COLOR_WINDOW));
	pDC->TextOut(1, 0, str);
	
	temp.push_back(point);

	CView::OnLButtonDown(nFlags, point);

	//CDC* pDC = GetDC();//此时利用的是CWnd的成员函数GetDC
	for (size_t i = 0; i < temp.size()-1; i++)
	{
		pDC->MoveTo(temp[i]);
		pDC->LineTo(temp[i + 1]);
	}
	ReleaseDC(pDC);
}


void CzhijiaodianView::OnAzengbu()
{
	// TODO: 在此添加命令处理程序代码
	CPoint Xn_Yn;
	CDC* pDC = GetDC();//此时利用的是CWnd的成员函数GetDC
	double x1, y1, x2, y2, xn_, yn_, dx, dy, D;
	double Xn, Yn;
	x1 = temp[0].x;
	y1 = temp[0].y;
	x2 = temp[1].x;
	y2 = temp[1].y;
	xn_ = temp[temp.size() - 1].x;
	yn_ = temp[temp.size() - 1].y;
	dx = x2 - x1;
	dy = y2 - y1;
	D = pow(dx, 2) + pow(dy, 2);
	Xn = ((x1 * dx + y1 * dy) * dx - (yn_ * dx - xn_ * dy) * dy) / D;
	Yn = ((yn_ * dx - xn_ * dy) * dx + (x1 * dx + y1 * dy) * dy) / D;
	Xn_Yn.x = long(Xn);
	Xn_Yn.y = long(Yn);
	temp.push_back(Xn_Yn);
	pDC->MoveTo(temp[temp.size() - 2]);
	pDC->LineTo(temp[temp.size() - 1]);

	pDC->MoveTo(temp[temp.size() - 1]);
	pDC->LineTo(temp[0]);

	ReleaseDC(pDC);
}


void CzhijiaodianView::OnAzhijiaohua()
{
	// TODO: 在此添加命令处理程序代码
	CDC* pDC = GetDC();

	double xi=temp[0].x, yi=temp[0].y, xj = temp[1].x, yj = temp[1].y,
		xk = temp[2].x, yk = temp[2].y, xn = temp[3].x, yn = temp[3].y;
	double A[24], W[3], P[64] = {};

	double AT[24] = {}, AQ[24] = {}, AQAT[9] = {}, K[3] = {}, P_AT[24] = {}, V[8] = {};

	A[0] = xk-xi ;      A[1] = yk-yi;       A[2] = 2*xj-xi-xk;  A[3] = 2*yj-yi-yk; A[4] = xi-xj;
	A[5] = yi-yj;       A[6] = 0.0;         A[7] = 0.0;         A[8] = 0.0;        A[9] = 0.0;A[10] = xn-xk;
	A[11] = yn-yk;      A[12] = 2*xk-xj-xn; A[13] = 2*yk-yi-yn; A[14] = xj-xk;     A[15] = yj-yk;
	A[16] = xk - xn;    A[17] = yk - yn;    A[18] = 0.0;        A[19] = 0.0;       A[20] = xi-xn;A[21] = yi-yn;
	A[22] = 2*xn-xk-xi; A[23] = 2*yn-yk-yi;

	W[0] = -((xi - xj) * (xk - xj) + (yi - yj) * (yk - yi));
	W[1] = -((xj - xk) * (xn - xk) + (yj - yk) * (yn - yk));
	W[2] = -((xk - xn) * (xi - xn) + (yk - yn) * (yi - yn));

	for (int i = 0; i < 64; i += 9)
	{
		P[i] = 1.0;
		
	}

	MatrixTra(3, 8, A, AT);//求AT(8*3)
	MatrixInv(3, P);//求P的逆
	MatrixMul(A, P, 3, 8, 8, AQ);//求A*Q(3*8)
	MatrixMul(AQ, AT, 3, 8, 3, AQAT);//求AQAT(3*3)
	MatrixInv(3, AQAT);//求AQAT的逆矩阵(3*3)
	MatrixMul(AQAT, W, 3, 3, 1, K);//求K
	MatrixMul(P, AT, 8, 8, 3, P_AT);//求QAT
	MatrixMul(P_AT, K, 8, 3, 1, V);//求V

	std::vector<CPoint>temp_new;
	CPoint var;
	for (size_t i = 0; i < temp.size(); i++)
	{
		var.x = temp[i].x + V[2 * i];
		var.y = temp[i].y + V[2 * i + 1];
		temp_new.push_back(var);
	}
	
	for (size_t i = 0; i < temp_new.size() - 1; i++)
	{
		pDC->MoveTo(temp_new[i]);
		pDC->LineTo(temp_new[i + 1]);
	}
	pDC->MoveTo(temp_new[temp_new.size() - 1]);
	pDC->LineTo(temp_new[0]);
	ReleaseDC(pDC);


}