|
个人想分享一些在大学中编写的一些程序，在进行坐标转换的时候，我们经常涉及到四参数与七参数的计算，在文章中，采用C#语言来进行编程，方便计算。 （1）四参数的计算： 在转换范围较小内不同的平面坐标转换通常采用二维四参数模型转换，二维四参数的转换模型的公式如下：  式中的x1，y1与x2，y2是两个坐标系下的坐标点；  是平移参数，单位为米； α是旋转参数，单位为弧度； m是尺度参数，无单位。 （2）七参数的计算： 两个空间直角坐标系进行转换计算就需要用到七个参数，其中包括：三个平移参数，（ΔX， ΔY ，ΔZ），三个旋转角度参数（εX，εY，εZ）以及一个尺度参数dK 公式如下：  式中的XT与X是用来表示P点在图中的两个坐标系O-XYZ与OT-XTYTZT的坐标向量  ∆X0是原点的平移向量，  R（ε）是一个旋转参数矩阵  为了便于计算，我们需要简化公式，因此假设当旋转角的值很小时，可以得到与得到公式的最终形式：  （3）四参数转换主要代码： publicstaticvoid Canshu4(Point2d[] p1, Point2d[] p2, int PointCount, refdouble rota, refdouble scale, refdouble dx, refdouble dy) { double u = 1.0, v = 0, Dx = 0.0, Dy = 0.0; int intCount = PointCount; //Matrix dx1 ;//误差方程改正数 Matrix B;//误差方程系数矩阵 // Matrix W ;//误差方程常数项 double[,] dx1 = newdouble[4, 1]; double[,] B1 = newdouble[2 * intCount, 4]; double[,] W1 = newdouble[2 * intCount, 1]; // Matrix BT, N, InvN, BTW; double[,] BT = newdouble[4, 2 * intCount]; double[,] N = newdouble[4, 4]; double[,] InvN = newdouble[4, 4]; double[,] BTW = newdouble[4, 1]; for (int i = 0; i < intCount; i++) { //计算误差方程系数 B1[2 * i, 0] = 1; B1[2 * i, 1] = 0; B1[2 * i, 2] = p1[i].X; B1[2 * i, 3] = -p1[i].Y; B1[2 * i + 1, 0] = 0; B1[2 * i + 1, 1] = 1; B1[2 * i + 1, 2] = p1[i].Y; B1[2 * i + 1, 3] = p1[i].X; } B = newMatrix(B1); for (int i = 0; i < intCount; i++) { //计算误差方程系常数 W1[2 * i, 0] = p2[i].X - u * p1[i].X + v * p1[i].Y - Dx; W1[2 * i + 1, 0] = p2[i].Y - u * p1[i].Y - v * p1[i].X - Dy; } //最小二乘求解 B.MatrixInver(B1, ref BT);//转置 B.MatrixMultiply(BT, B1, ref N); InvN = B.MatrixOpp(N); B.MatrixMultiply(BT, W1, ref BTW); B.MatrixMultiply(InvN, BTW, ref dx1); Dx = Dx + dx1[0, 0]; Dy = Dy + dx1[1, 0]; u = u + dx1[2, 0]; v = v + dx1[3, 0]; dx = Dx; dy = Dy; rota = Math.Atan(v / u); scale = u / Math.Cos(rota); dx = Math.Round(dx, 6); dy = Math.Round(dy, 6); rota = Math.Round(rota, 6); scale = Math.Round(scale, 6); } （4）七参数转换主要代码： publicstaticvoid Canshu7(Point3d[] p1, Point3d[] p2, int PointCount, refdouble rotax, refdouble rotay, refdouble rotaz, refdouble scale, refdouble dx, refdouble dy, refdouble dz) { double[,] B1 = newdouble[PointCount * 3, 7]; Matrix B = newMatrix(B1); double[,] dx1 = newdouble[7, 1];//V=B*X-L double[,] L = newdouble[PointCount * 3, 1]; double[,] BT = newdouble[7, PointCount * 3]; double[,] N = newdouble[7, 7]; double[,] InvN = newdouble[7, 7]; double[,] BTL = newdouble[7, 1]; //初始化L矩阵 for(int i=0;i<pointcount*3;i++)< p=""> { if(i%3==0) { L[i, 0] = p2[i / 3].X; } elseif(i%3==1) { L[i, 0] = p2[i / 3].Y; } elseif (i % 3 == 2) { L[i, 0] = p2[i / 3].Z; } } //初始化B矩阵 for (int i = 0; i < PointCount * 3; i++) { if(i%3==0) { B1[i, 0] = 1; B1[i, 1] = 0; B1[i, 2] = 0; B1[i, 3] = p1[i/3].X; B1[i, 4] = 0; B1[i, 5] = -p1[i/3].Z; B1[i, 6] = p1[i/3].Y; } elseif(i%3==1) { B1[i, 0] = 0; B1[i, 1] = 1; B1[i, 2] = 0; B1[i, 3] = p1[i / 3].Y; B1[i, 4] = p1[i/3].Z; B1[i, 5] = 0; B1[i, 6] = -p1[i / 3].X; } elseif (i % 3 == 2) { B1[i, 0] = 0; B1[i, 1] = 0; B1[i, 2] = 1; B1[i, 3] = p1[i / 3].Z; B1[i, 4] = -p1[i / 3].Y; B1[i, 5] = p1[i/3].X; B1[i, 6] = 0; } } //转置 B.MatrixInver(B1,ref BT); //法方程矩阵 //N=BT*B B.MatrixMultiply(BT,B1,ref N); //求逆 InvN=B.MatrixOpp(N); //BTL=BT*L B.MatrixMultiply(BT,L,ref BTL); //dx1=invN*BTL; B.MatrixMultiply(InvN,BTL,ref dx1); // dx = Math.Round(dx1[0, 0], 6); dy = Math.Round(dx1[1, 0], 6); dz = Math.Round(dx1[2, 0], 6); scale = Math.Round(dx1[3, 0], 6); rotax = Math.Round(dx1[4, 0] / dx1[3, 0], 6); rotay = Math.Round(dx1[5, 0] / dx1[3, 0], 6); rotaz = Math.Round(dx1[6, 0] / dx1[3, 0], 6); } |
|
|
