矩阵变换实现从二维到三维坐标映射的方法探究

如果我们面前有一张二维平面图和一个三维模型，如何让平面图和模型上的物体产生空间对应关系呢？本文将针对这一问题进行探究。

---

问题背景

举个例子，下图左边是我家房子的平面图，右边是我家房子的三维模型。家里装修我买了一个沙发，我在平面图上把沙发的位置标出来，在客厅区域放置一个图标，然后我希望三维模型能在客厅的相应位置自动添加一个沙发。

简单来讲就是如何把二维平面上的一个二维坐标一一映射到三维空间中的一个三维坐标

解决思路

在这个问题中我们有两个输入和三个输出，输入是二维平面$A_1$中的二维坐标$(x_1,y_1)$，输出是三维空间$B_1$中的三维坐标$(x_2,y_2,z_2)$，从函数的角度看就是$F(x_1,y_1)=(x_2,y_2,z_2)$。看到这样一个公式是不是觉得很抽象，完全无从下手，所以就要考虑这个问题能不能简化。
###问题简化
幸运的是我们一般的平面图都是很“正”的，不是“斜切”的，就是说平面图所在的二维平面$A_1$是和三维空间$B_1$的$xOy$平面平行的，也可以说平面$A_1$的$z$轴高度就是$z_2$，这样问题用函数表示不就变成了$F(x_1,y_1,z_2)=(x_2,y_2,z_2)$。如果再考虑已知$z_2$值的情况，那问题就更加简化了$F(x_1,y_1)=(x_2,y_2)$，两个输入两个输出感觉还有可能解决。
从直观的几何角度我们可以这样考虑，把平面图看作一张有弹性的透明薄膜，你拿着这张膜往模型上盖，大了就压缩，小了就拉伸，歪了就平移或者转一转，直到这张薄膜上的所有点都能“严丝合缝”的和模型对上，那么对于薄膜上的点在模型上总有唯一的一个与之重合的点，反过来也是一样。但其实我们仔细想想，如果在拉伸薄膜的过程中，始终是均匀变化的，只要薄膜上有两个确定的点和模型上两个确定的点已经重合了，那么其他的点肯定也是重合的。这样的话只要我们能知道这张薄膜是如何伸缩、平移和旋转的，那么在变化前的薄膜上随便选一个点，就一定能计算出薄膜变化后这个点会与模型上的哪个点重合。而对于薄膜这样的变化，我们很自然的就想到了线形代数中的“线形变换”，因为薄膜均匀伸缩、平移、选择的过程符合线形变换的特征。
我们可以用几张图简单的描述下整个过程。

平面$A_1$上的点要和空间$B_1$上的点一一映射（反过来不行，不然就变成多对一了）。

注意，这里的坐标系是完全不同的，X轴、Y轴只是示意，两边的坐标原点和单位长度是不同的。这里的“平面”等价于"坐标系"的概念。

用$A_2$来表示空间$B_1$上的平面$z=z_2$，在平面$A_2$上分别选择三个不同的参照点$i(x_i,y_i),j(x_j,y_j),k(x_k,y_k)$，再从平面$A_1$上选择对应的三个参照点$a(x_a,y_a),b(x_b,y_b),c(x_c,y_c)$。通过变换$A_1$让点$a,b,c$分别与平面$A_2$上的点$i,j,k$重合，这样的话平面$A_1$和平面$A_2$上所有的点都会一一对应。注意这里的坐标系仍然是不同的，不仅是坐标原点不同，而且基底也就算单位向量和单位长度是不同的。在线形代数中，我们可以用矩阵来表示从一个坐标系到另外一个同维度坐标系的线形变换过程。但目前我们仍无法用矩阵来表示从$A_1$到$A_2$的变换过程，因为线形变换的要求是原点重合，否则就不能称为 线形变换。那么如何让两个平面的原点重合呢？我们将继续简化。

实际上两个平面上分别选取两个（共四个）点也是可以的，也能解决这个问题，但分别选取三个（共六个）更利于我们进行计算和理解。另外两个坐标系的单位长度不同并不会影响最终结果，因为单位长度不管是1m、1mm还是1px，都可以抽象的看作是1，后续的计算过程也会体现这一特性。

分别将平面$A_1$中的点$b$和平面$A_2$中的点$j$作为坐标原点，这样平面$A_1$中的点$a,b,c$和平面$A_2$中的点$i,j,k$的坐标就分别变成了$(x_a-x_b,y_a-y_b),(0,0),(x_c-x_b,y_c-y_b)$和$(x_i-x_j,y_i-y_j),(0,0),(x_k-x_j,y_k-y_j)$。

既然点$b$和点$j$都是原点，我们就让他们重合。

其实这里的平面$A_1$和平面$A_2$已经不是原来的平面了，因为原点已经发生了改变，但为了保持前后连贯性，仍然使用$A_1,A_2$来表示这两个平面。

然后围绕原点旋转、缩放$A_1$，使点$i$和点$a$重合，点$k$和点$c$重合，最终$A_1$和$A_2$重合。

计算过程

现在从$A_1$到$A_2$的变换就可以用矩阵来表示了。用公式表示就是下面的形式。
$$T{A}_1=A_2$$

其中$T$表示从$A_1$到$A_2$的变换矩阵，对$A_1$应用线形变换$T$就会得到$A_2$。另外在线形代数中，点和向量往往是可以互换的，比如$(x_a-x_b,y_a-y_b)$可以用来表示从点$b$到点$a$的向量$v_a$，对向量$v_a$应用线形变换$T$就可以得到一个向量$v_i$，而$v_i$又可以用$(x_i-x_j,y_i-y_j)$来表示。
其中$$A_1=\left[\begin{array}{cc}{x}_a-x_b& x_c-x_b\\ y_a-y_b& y_c-y_b\end{array}\right],A_2=\left[\begin{array}{cc}{x}_i-x_j& x_k-x_j\\ y_i-y_j& y_k-y_j\end{array}\right],T=\left[\begin{array}{cc}a& c\\ b& d\end{array}\right]$$
即$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\left[\begin{array}{cc}{x}_a-x_b& x_c-x_b\\ y_a-y_b& y_c-y_b\end{array}\right]=\left[\begin{array}{cc}{x}_i-x_j& x_k-x_j\\ y_i-y_j& y_k-y_j\end{array}\right]$$
只要能求出变换矩阵$T$，也就是分别算出$a,b,c,d$的值，那么对于$A_1$上任一点$m(x_m,y_m)$，其以点$b$为起点的向量$v_m=(x_m-x_b,y_m-y_b)$经过变换矩阵$T$作用后都能在$A_2$上得到一个唯一的向量$v_n$，其中$v_n=(x_n-x_j,y_n-y_j)$，点$n(x_n,y_n)$即为点$m$唯一对应的点。
首先我们要求矩阵$T$，可以通过逆矩阵求解$$T{A}_1=A_2$$
$$T{A}_1{A}_1^{-1}=A_2{A}_1^{-1}$$
$$T=A_2{A}_1^{-1}$$
又从逆矩阵的定义可以知道
$$A_1^{-1}=\frac{1}{|A_1|}{A}_1^{\ast }$$
其中$|A_1|$为矩阵$A_1$的行列式，$A_1^{\ast }$是矩阵$A_1$的伴随矩阵，又
$$|A_1|=\left|\begin{array}{cc}{x}_a-x_b& x_c-x_b\\ y_a-y_b& y_c-y_b\end{array}\right|=(x_a-x_b)(y_c-y_b)-(y_a-y_b)(x_c-x_b)$$
$$A_1^{\ast }=\left[\begin{array}{cc}(-1{)}^2(y_c-y_b)& (-1{)}^3(y_a-y_b)\\ (-1{)}^3(x_c-x_b)& (-1{)}^4(x_a-x_b)\end{array}\right]$$
由于上面的参数在实际运算中都是常数，所以矩阵$T$是很容易计算出来的，下面的公式就不再赘述了。有了矩阵$T$之后，下一步就是根据点$m(x_m,y_m)$算出点$n(x_n,y_n)$了。这个过程用公式表示就是：
$$A_1=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$$
$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\left[\begin{array}{c}{x}_m-x_b\\ y_m-y_b\end{array}\right]=\left[\begin{array}{c}{x}_n-x_j\\ y_n-y_j\end{array}\right]$$
其中$b(x_b,y_b),j(x_j,y_j)$分别为原点。根据上面公式可以得到
$$x_n=a{x}_m+b{y}_m-a{x}_b-b{y}_b+x_j,y_n=c{x}_m+d{y}_m-c{x}_b-d{y}_b+y_j$$

为了更好的理解上述过程，推荐观看线性代数系列视频——线性代数的本质

代码实现

上述过程的关键在于矩阵的运算，那么如何在高级语言中实现矩阵运算呢？对于C#语言，我们可以使用数学计算类库MathNet。在其对应的类库MathNet.Numerics中提供了一系列矩阵运算方法。
使用VS2015创建一个控制台程序MatrixApp，打开NuGet管理器，搜索并安装MathNet类库。

计算变换矩阵

示例代码如下：

        private static DenseMatrix GetTranformMatrix()
        {
            var point_b = new DenseVector(new[] { -1.0, 4.0 });
            var point_c = new DenseVector(new[] { 3.0, 6.0 });
            Debug.WriteLine(String.Format("平面图参照点坐标：a {0},b {1},c {2}", point_a.ToString(), point_b.ToString(), point_c.ToString()));

            var point_j = new DenseVector(new[] { 1.0, 3.0 });
            var point_k = new DenseVector(new[] { 2.0, 1.0 });
            Debug.WriteLine(String.Format("模型参照点坐标：i {0},j {1},k {2}", point_i.ToString(), point_j.ToString(), point_k.ToString()));

            var vector_ba = point_b - point_a;
            var vector_ca = point_c - point_a;

            var vector_ji = point_j - point_i;
            var vector_ki = point_k -point_i;
            var matrixA = DenseMatrix.OfColumnVectors(new Vector<double>[] { vector_ba, vector_ca });
            Debug.WriteLine("变换前矩阵1：" + matrixA.ToString());
            var matrixB = DenseMatrix.OfColumnVectors(new Vector<double>[] { vector_ji, vector_ki });
            Debug.WriteLine("变换后矩阵2：" + matrixB.ToString());

            var valueA = matrixA.Determinant();//矩阵A的行列式
            Debug.WriteLine("矩阵1的行列式：" + valueA.ToString());

            var matrixAdjoint = DenseMatrix.OfArray(new double[,] { { 2.0, 0 }, { 0, -4.0 } });
            Debug.WriteLine("矩阵1的伴随矩阵：" + matrixAdjoint.ToString());

            var matrixAnti = matrixAdjoint.Divide(valueA);
            Debug.WriteLine("矩阵1的逆矩阵：" + matrixAnti.ToString());

            var matrixT = (DenseMatrix)(matrixB * matrixAnti);
            Debug.WriteLine("变换矩阵T：" + matrixT.ToString());

            return matrixT;
        }

输入二维坐标&输出三维坐标

        private static void Get3DFrom2D(DenseVector point,double height)
        {
            var matrixT = GetTranformMatrix();
            double[] beforeArry = point.ToArray();
            var vector_before = point - Program.point_a;
            var vector_after = matrixT * vector_before;
            var vector_result = vector_after + point_i;
            double[] resultArry = vector_result.ToArray();

            Debug.WriteLine(string.Format("平面图中的坐标为：({0},{1})", beforeArry[0], beforeArry[1]));
            Debug.WriteLine("模型标高为："+height);
            Debug.WriteLine(string.Format("映射到模型中的坐标为：({0},{1},{2})",resultArry[0],resultArry[1], height));
        }

结果验证

在主程序中输入数据测试一下：

        static DenseVector point_a = new DenseVector(new[] { 3.0, 4.0 });
        static DenseVector point_i = new DenseVector(new[] { 1.0, 1.0 });
        static void Main(string[] args)
        {
            var point2D = new DenseVector(new[] { 1.0, 4.0 });

            Get3DFrom2D(point2D, 12.1);
        }

测试用参照点的坐标分别为$a(3,4),b(-1,4),c(3,6),i(1,1),j(1,3),k(2,1)$，其中点$a$和点$i$作为原点，标高为12.1。为了更直观的验证结果，可以动手将上述点在白纸上画出来。运行程序输出的结果为：

平面图参照点坐标：a DenseVector 2-Double
3
4
,b DenseVector 2-Double
-1
 4
,c DenseVector 2-Double
3
6

模型参照点坐标：i DenseVector 2-Double
1
1
,j DenseVector 2-Double
1
3
,k DenseVector 2-Double
2
1

变换前矩阵1：DenseMatrix 2x2-Double
-4  0
 0  2

变换后矩阵2：DenseMatrix 2x2-Double
0  1
2  0

矩阵1的行列式：-8
矩阵1的伴随矩阵：DenseMatrix 2x2-Double
2   0
0  -4

矩阵1的逆矩阵：DenseMatrix 2x2-Double
-0.25    0
    0  0.5

变换矩阵T：DenseMatrix 2x2-Double
   0  0.5
-0.5    0

平面图中的坐标为：(1,4)
模型标高为：12.1
映射到模型中的坐标为：(1,2,12.1)

以上过程已经解决并实现了理论问题，下一步就需要在项目中完善相应的程序了。

---

后续思考

1.参照点的选取有什么要求吗？

参照点的选取一定不是随便的。平面$A_1$中的参照点$a,b,c$和平面$A_2$中的参照点$i,j,k$一定是要在实际意义上是一一对应的，比如平面图$A_1$上代表沙发的点$a$一定是要和模型$B_1$上代表沙发的点$i$对应，而不能和代表餐桌的点$j$对应。虽然将点$a$与点$j$对应也能够求出一个变换矩阵$T$，也能实现$A_1$到$B_1$的一一映射，但那绝不是我们想要的结果。
另外参照点的选取要尽量分散，避免重合或者共线

2.原点的选取会对结果有影响吗？

从上面的过程我们可以推导出，只要原点的选取遵循对应的原则，即如果$A_1$选点$a$为原点，$A_2$就要选对应的点$i$为原点，如果$A_1$选择$b$为原点，$A_2$就要选对应的点$j$为原点，这样的话原点的变化并不会对结果有影响。

3.计算的结果会与实际情况有误差吗？

有误差是必然的。虽然理论推导没有问题，但在我们实际操作过程中参照点的选取就已经产生了误差。而且由于MathNet类在计算过程中使用了Double类型，也会产生一些误差。但是在目前的应用场景下，能够确认沙发在客厅而不是在厨房就已经足够了，这点误差是可以接受的。上面提到的参照点的选取尽量分散且避免重合及共线也是为了尽量减小误差。

4.能够推广到更高维度的空间吗？

答案是肯定的。线性变换并不是只能应用于二维平面，推广到三维、四维甚至更高维空间也是成立的，只不过变换矩阵$T$增加到了3阶、4阶…n阶，计算量也会更大。

---

代码示例