D3D坐标系统下3D世界坐标映射到2D屏幕坐标的平移矩阵

D3D中绘画3D模型基本上就是靠3个矩阵World, View, Projection来联合进行模型位置定位、视角定位及透视变形的，这与2D绘制一个图形只需要给出屏幕上的一个像素坐标就能进行定位有着非常大的不同。在某些场合，我们想根据屏幕上的像素坐标来绘制3D模型，一般可以通过用正交投影代替透视投影就能轻松进行绘制，但在一些极其特殊的情况下我们还想让绘制出的3D模型保持原先指定的3个矩阵所有变换，这就需要通过这3个矩阵及目标像素点反求出一个平移矩阵，这篇文章就是介绍该怎么做的。

---

先简单介绍一下在D3D里3D空间坐标在不做任何变化情况下是怎么转换成2D屏幕坐标。D3D的空间坐标在3个轴范围是[-1,1]，D3D绘画到2D屏幕上的区域叫做ViewPort，ViewPort的横坐标范围是[Left, Right]，纵坐标范围是[Top, Bottom]，ViewPort还有近点(near)和远点(far)2个值代表z轴，则D3D空间坐标里的一点sp(x,y,z)的x坐标Xsp对应到屏幕上坐标sc(x,y)的x坐标Xsc的对应公式就是（这个详细解释网上有，这里不再赘述）：

$$\frac{Xsc - Left}{Right - Left} = \frac{Xsp - (-1)}{1 - (-1)}$$
让Width = Right - Left，上面公式最后就变成：
$$Xsc = \frac{(Xsp + 1) * Width}{2} + Left$$
同样得出y轴：
$$Ysc = \frac{(Ysp + 1) * Height}{2} + Top$$
z轴（虽然屏幕上是没有z坐标的，但这里给出是为了后面推导的时候能统一处理）：
$$Zsc = Zsp * (far - near) + near$$

---

接下来简单说一下齐次坐标，D3D系统里用得坐标并不只是(x,y,z)，而是用4维坐标(x,y,z,w)来表示，这里齐次坐标与3维坐标的转换关系为：

$$\begin{array}{c} \begin{bmatrix} x\y\z\w\\end{bmatrix} = \begin{bmatrix} \frac{x}{w}\\frac{y}{w}\\frac{z}{w}\1\\end{bmatrix} \end{array}$$
也就是齐次坐标的每一项都除以w，这个过程叫做齐次坐标正常化(Normalize)也可以叫标准化，得到的坐标就是正常化坐标，当然要注意一点是当齐次坐标表示向量的时w=0，这时候则不用去除w，向量就代表了坐标，也就是当w=0时不用做任何处理。最后简单讲一下w的作用，w从绘画结果来看实质上就是一个非线性变形，透视投影就是利用了w来进行变换的。

齐次坐标是通过一个正常化3维坐标与变换矩阵相乘得到的，因为D3D是左手坐标系，是用变换矩阵左乘坐标，变换矩阵为4x4矩阵：

$$\begin{array}{ccc} \begin{bmatrix} x & y & z & 1 \end{bmatrix} * \begin{bmatrix} M11 & M12 & M13 & M14\M21 & M22 & M23 & M24\M31 & M32 & M33 & M34\M41 & M42 & M43 & M44\\end{bmatrix} \end{array}$$
这个过程被称为Transform，结果是
X = x * M11 + y * M21 + z * M31 + M41
Y = x * M12 + y * M22 + z * M32 + M42
Z = x * M13 + y * M23 + z * M33 + M43
W = x * M14 + y * M24 + z * M34 + M44
这个结果是等会推导要用到的。

---

接下来开始推导：
我们已经知道空间坐标一点怎么通过变换矩阵得到齐次坐标，然后只要把齐次坐标正常化后就能知道变换后正确的空间坐标，再通过空间坐标到屏幕坐标的转换就能得到3D空间坐标映射到2D屏幕上的一点了，这个过程叫做Project。我们现在已知2D坐标一点P(x,y)，可以简单通过反向刚才操作(Unproject)求出空间的一个坐标Q(x1,y1,z1)，我们只要能让系统把这个Q坐标当成0坐标Z来进行绘画就大功告成了，问题是这个Q坐标的各个分量不仅是包含了平移信息也包含了旋转视角变换和投影信息，所以让0坐标Z对V的分量进行简单平移是不行的。我们假设Z通过对V(x’,y’,z’)进行平移后进行Project就能得到正确的P：

已知转换矩阵 M，设f(v)为齐次坐标转换成屏幕坐标函数
则f( (Z + V) * M) = P
其中Z为(x, y, z, 1), V为(x’,y’,z’,0),
Z + V = (x + x’, y + y’, z + z’, 1)
把这个带入之前Transform后就有：
f(
X’ = (x + x’) * M11 + (y + y’) * M21 + (z + z’) * M31 + M41
Y’ = (x + x’) * M12 + (y + y’) * M22 + (z + z’) * M32 + M42
Z’ = (x + x’) * M13 + (y + y’) * M23 + (z + z’) * M33 + M43
W’ = (x + x’) * M14 + (y + y’) * M24 + (z + z’) * M34 + M44
) = P
简化后就是：
f(
X’ = x’ * M11 + y’ * M21 + z’ * M31 + X
Y’ = x’ * M12 + y’ * M22 + z’ * M32 + Y
Z’ = x’ * M13 + y’ * M23 + z’ * M33 + Z
W’ = x’ * M14 + y’ * M24 + z’ * M34 + W
) = P

而右边P可以表示为原先的0点坐标Z(x,y,z)通过Project操作后的平移，也就是
P = f(Z* M) + Offset(x,y)
化成Transform后的形式就是：
P = f((X, Y, Z, W)) + (Ox, Oy)

这样就得到：
f(
X’ = x’ * M11 + y’ * M21 + z’ * M31 + X
Y’ = x’ * M12 + y’ * M22 + z’ * M32 + Y
Z’ = x’ * M13 + y’ * M23 + z’ * M33 + Z
W’ = x’ * M14 + y’ * M24 + z’ * M34 + W
)
= f((X, Y, Z, W)) + (Ox, Oy)

然后让我们把f这个函数展开，就是带入3d坐标映射成2d坐标的转换，当然在转换前还需要进行齐次坐标正常化：

$$\begin{cases} \frac{(\frac{X‘}{W‘} + 1) * Width}{2} + Left = \frac{(\frac{X}{W} + 1) * Width}{2} + Left + Ox\\frac{(\frac{Y‘}{W‘} + 1) * Height}{2} + Top = \frac{(\frac{Y}{W} + 1) * Height}{2} + Top + Oy\\frac{Z‘}{W‘} * (far - near) + near = \frac{Z}{W} * (far - near) + near \end{cases}$$
然后用刚才X’,Y’,Z’,W’公式右边分别替换X’,Y’,Z’,W’，整理后得到关于x’,y’,z’的方程组：
$$\left\{ \begin{array}{ll} x‘*(M11-M14*(\frac{X}{W} + \frac{Ox * Width}{2})) + y‘*(M21-M24*(\frac{X}{W} + \frac{Ox * Width}{2})) + z‘*(M31-M34*(\frac{X}{W} + \frac{Ox * Width}{2})) &= W * (\frac{X}{W} + \frac{Ox * Width}{2}) - X\x‘*(M12+M14*(-\frac{Y}{W} + \frac{Oy * Height}{2})) + y‘*(M22+M24*(-\frac{Y}{W} + \frac{Oy * Height}{2})) + z‘*(M32+M34*(-\frac{Y}{W} + \frac{Oy * Height}{2})) &= -W * (-\frac{Y}{W} + \frac{Oy * Height}{2}) - Y\x‘*(M13-M14*(\frac{Z}{W}*(far-near) + near) + y‘*(M23-M24*(\frac{Z}{W}*(far-near) + near) + z‘*(M33-M34*(\frac{Z}{W}*(far-near) + near) &= W * (\frac{Z}{W}*(far-near) + near) - Z \end{array} \right.$$

这是个标准的三元一次方程组，可用矩阵表示为：

$$\begin{array} {cccc} \begin{bmatrix} L1 & L2 & L3 \M1 & M2 & M3 \N1 & N2 & N3 \end{bmatrix} & * & \begin{bmatrix} x‘ \y‘ \z‘ \end{bmatrix} & = & \begin{bmatrix} A1\\A2\\A3 \end{bmatrix} \end{array}$$

解这个方程组只要用高斯-若尔当消元法对增广矩阵（如下）消元就可以了：

$$\left[ \begin{array}{ccc|c} L1&L2&L3&A1\ M1&M2&M3&A2\ N1&N2&N3&A3 \end{array} \right]$$

最后直接给出消元结果：
z’ = ((A3 * L1 - A1 * N1) * (M2 * L1 - L2 * M1) - (A2 * L1 - A1 * M1) * (N2 * L1 - L2 * N1)) / ((N3 * L1 - L3 * N1) * (M2 * L1 - L2 * M1) - (M3 * L1 - L3 * M1) * (N2 * L1 - L2 * N1))
y’ = ((A2 * L1 - A1 * M1) - (M3 * L1 - L3 * M1) * z) / (M2 * L1 - L2 * M1)
x’ = A1 / L1 - z * L3 / L1 - y * L2 / L1

这就是3x3矩阵的解，特别的当z=0，就得到2x2矩阵解为：
y’ = (A2 * L1 - A1 * M1) / (M2 * L1 - L2 * M1)
x’ = (A1 / L1 - y * L2 / L1)

至此终于求出了V(x’,y’,z’)，我们只要对Z点平移V后就能正确Project到P点。

---

把以上过程转为编码（使用了SharpDX）：

        //pos为期望要绘画的屏幕坐标
        private Matrix Get2DTranslationMatrix(ViewportF viewPort, Vector2 pos, Matrix world, Matrix view, Matrix projection)
        {
            //先取空间0点坐标对应的屏幕坐标
            var screenZ = viewPort.Project(Vector3.Zero, projection, view, world);

            //相对于viewPort左上角求出期望坐标与0点坐标在屏幕坐标系的偏移值O(x,y)
            var diff = new Vector3(pos.X + viewPort.X, pos.Y + viewPort.Y, 0) - screenZ;

           //由于D3D是通过3个矩阵进行变换，而我推导只需要1个，所以就把3个矩阵相乘合成一个变换矩阵进行计算
            var projM = world * view * projection;

            //求出转换后的齐次坐标，并确保w不为0
            var transV = Vector3.Transform(Vector3.Zero, projM);
            if (MathUtil.IsZero(transV.W))
            {
                transV.W = 1.0f;
            }

            var w = viewPort.Width;
            var h = viewPort.Height;
            /*
            C1,C2在消元时已消去
            C3,C4也没必要，理由见下
            */
            var C1 = viewPort.X;
            var C2 = viewPort.Y;
            var C3 = viewPort.MaxDepth;
            var C4 = viewPort.MinDepth;

           //设定增广矩阵所有系数
            var L1 = projM.M11 - projM.M14 * (transV.X / transV.W + diff.X * 2 / w);
            var L2 = projM.M21 - projM.M24 * (transV.X / transV.W + diff.X * 2 / w);
            var L3 = projM.M31 - projM.M34 * (transV.X / transV.W + diff.X * 2 / w);

            var M1 = projM.M12 + projM.M14 * (-transV.Y / transV.W + diff.Y * 2 / h);
            var M2 = projM.M22 + projM.M24 * (-transV.Y / transV.W + diff.Y * 2 / h);
            var M3 = projM.M32 + projM.M34 * (-transV.Y / transV.W + diff.Y * 2 / h);

           /*
           z系数应该是
           var N1 = projM.M13 - projM.M14 * (transV.Z / transV.W * (C3 - C4) + C4);
            var N2 = projM.M23 - projM.M24 * (transV.Z / transV.W * (C3 - C4) + C4);
            var N3 = projM.M33 - projM.M34 * (transV.Z / transV.W * (C3 - C4) + C4);
            但一般情况下viewport始终是在0点创建所以near=0,far=1，可以简化
           */
            var N1 = projM.M13 - projM.M14 * (transV.Z / transV.W);
            var N2 = projM.M23 - projM.M24 * (transV.Z / transV.W);
            var N3 = projM.M33 - projM.M34 * (transV.Z / transV.W);

            var A1 = transV.W * (transV.X / transV.W + diff.X * 2 / w) - transV.X;
            var A2 = -transV.W * (-transV.Y / transV.W + diff.Y * 2 / h) - transV.Y;
            /*
            var A3 = transV.W * (transV.Z / transV.W * (C3 - C4) + C4) - transV.Z;
            由于near=0,far=1,A3可以简化为0
            */
            var A3 = 0;

            /*
            可以不考虑z坐标，直接用2x2消元结果代替
            var y1 = (A2 * L1 - A1 * M1) / (M2 * L1 - L2 * M1);
            var x1 = (A1 / L1 - y * L2 / L1);
            这样上面的系数L3,M3,N1,N2,N3,A3都不用计算，下面返回的时候z可以取0
            */
            var z = ((A3 * L1 - A1 * N1) * (M2 * L1 - L2 * M1) - (A2 * L1 - A1 * M1) * (N2 * L1 - L2 * N1)) / ((N3 * L1 - L3 * N1) * (M2 * L1 - L2 * M1) - (M3 * L1 - L3 * M1) * (N2 * L1 - L2 * N1));
            var y = ((A2 * L1 - A1 * M1) - (M3 * L1 - L3 * M1) * z) / (M2 * L1 - L2 * M1);
            var x = A1 / L1 - z * L3 / L1 - y * L2 / L1;

            //返回平移矩阵
            return Matrix.Translation(x, y, z);
        }

最后使用的时候把world左乘平移矩阵得到一个新world就行了：
world = Get2DTranslationMatrix(…) * world

至此大功告成。

简单介绍一下运用场景：
摄像头实时采集画面后在VR眼镜(oculus rift)里播放，就是在D3D里做一个二维贴图，现在又对贴图进行了图形分析，在图像内匹配到的图形（利用openCV)上画出一个3D模型（也就是AR)，而VR设备有着自己特殊的变换矩阵，所以通过这个可以轻松把模型画到侦测到的图形标记上，换句话说就是在没有AR库的辅助变换下完成了AR效果。