Unity3D绘制两圆柱体相贯线

1 相贯线

  项目中需要用代码去绘制两个圆柱体相贯线，花了两天时间可算整明白怎么画了。国内有关这方面的文章不多，所以花了很多时间，我这里总结一下，如果能够帮助到需要的同学，也算善莫大焉。
  相贯线就是两个圆柱体相交表面所形成的曲线。如图。

  最终使用Unity计算并绘制出的相贯线，如图。

  我们知道点构成线，线再构成面。特殊点可以构成特殊的曲线，特殊的曲线可以构成特殊的面。这个特殊是指，这些顶点的xyz坐标之间满足一定的关系，或者说满足一个方程。比方说顶点的$x$、$y$坐标满足$y=x^2$，那么当$z=0$时在空间中画出来就是一条抛物线。

  如果$z$从-1到1之间变化，那画出来就是个抛物面。

  所以我们只需把相贯线的方程算出来，然后就能够用代码把相贯线画出来。
  下面我将详细解释如何根据两个圆柱体的参数(位置、半径)去求解两个圆柱体的相贯线（需要一些矩阵方面的知识才能深刻理解背后的原理，不懂的强烈建议去学习一下，推荐《Unity Shader入门精要》第4章），然后会给出代码和工程去实现绘制相贯线。

2 求解相贯线

2.1 圆柱体方程

  我们知道，要列方程，必须先选定一个参考坐标系。只有选择了合适的参考系，我们才能很快地列出方程。如果参考系没选好，列方程将会变得很困难。对于圆柱体，我们一般把圆柱体轴线上的一个点当作原点O，然后把与圆柱底面平行的平面当作Oxz平面，圆柱的轴线当作y轴（这里的xyz我用的是Unity中的方向，大家可以随意选择的，不一定要按这个来设置）。
  我们知道二维坐标系(Oxz平面)下圆的方程为：$$x^2+z^2=R^2$$
  然后这个圆沿着y轴不断移动，就形成圆柱体，所以这个圆柱体的方程如下

  现在我们另外选一个参考坐标系，这个参考坐标系的原心在圆柱轴线“右”边2个单位，坐标系变成O’x’y’z’，那么以Ox’y’z’为参考坐标系列出来的圆柱体方程就变为

  那么现在，请思考一个问题，Oxyz坐标系下的xyz值与O’x’y’z’坐标系下的值x’y’z’的数学关系是什么？
  在上面的例子中，我们能知道O’x’y’z’坐标可以理解为Oxyz坐标系往右移动了2个单位，所以就有方程
$$\{\begin{array}{l}x+2=x^{\prime }\\ y=y^{\prime }\\ z=z^{\prime }\end{array}$$
  这一点很简单，但却特别重要（就是说两个参考坐标系的坐标，一定是存在一种转化关系的），下面部分将会深入一点探讨两个坐标系地相互转换问题。接下来我们看看怎么列相贯线的方程，同时看看如何求解。

2.2 求解相贯线方程

2.2.1 选择参考坐标系

  先设定一些参数。如图，两个圆柱体相交。大圆柱体的半径为R，小圆柱体的半径为r，然后两个圆柱体轴心的夹角为$\alpha$，蓝色圆柱体偏移灰色圆柱体的距离为e。

  第一步，我们先把两个立方体的方程列出来，要列方程，就得先选择参考坐标系。这里选择参考坐标系就有讲究了，不能随便选。在2.1节中提到，我们一般把圆柱体轴线上的一个点当作原点O，然后把与圆柱底面平行的平面当作Oxz平面，圆柱的轴线当作y轴。但是在这里，我们不这样选择参考系，而是把参考坐标系$xyz$轴的方向选为和Unity世界坐标的$xyz$方向一样，为什么？主要是为了方便我们之后生成线的顶点（理论上参考坐标系无论怎么选都是可以的，但是这样选是最方便的）。
  灰色圆柱体的坐标系怎么选？
  先确定原心$O$点。确定$O$点分两个步骤，①蓝色圆柱体的轴线与灰色圆柱体圆心处水平切面会有一个交点，②把这个交点再平移到灰色圆柱体的轴线上，即得到$O$点。

  然后取灰色圆柱体的轴线为$z$轴，两圆柱体的公垂线为$x$轴，垂直于$Oxz$平面同时朝向蓝色圆柱体一边的直线为$y$轴。最终，灰色圆柱体的参考坐标系如下图。记为$Oxyz$。

  这里强调一下，我们选择的参考坐标系 $Oxyz$ 与Unity的世界坐标系的方向是一模一样的，只是为了方便后面设置曲线的网格用。除此之外，两者没有任何关系! 这一点大家一定要记住并且理解。
   再来看蓝色圆柱体的参考坐标系怎么选。其实蓝色圆柱体的参考坐标系选择特别特别简单，就是把灰色圆柱体的参考坐标系旋转绕着x轴旋转一定角度（这里的角度大小是$180°-\alpha$，后面部分我们还会看到旋转方向是有正负之分的，由于我们选择的参考坐标系是个左手坐标系，所以后面我们把这个角度带入矩阵计算时使用的是$-(180°-\alpha )$，这一点后面会详细说到），使得z轴与灰色圆柱体的轴线平行，然后将旋转后的坐标系平移到蓝色圆柱体的轴线与灰色圆柱体圆心处水平切面的交点(看上上图)处，就得到蓝色圆柱体的参考坐标系。具体过程如下图。

  灰色圆柱体最终的参考坐标系如下。我们记为$O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$。

  两个坐标系放在一起如下。

2.2.2 两个参考坐标系之间的变换

  那么问题来了，两个参考坐标系之间如何变换或者说如何将坐标系$O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$中的一个点$(x^{\prime },y^{\prime },z^{\prime })$转换为坐标系$Oxyz$中$(x,y,z)$点呢？你可能会问，我们为什么要知道这两个坐标系的关系呢，没啥原因，因为后面求解相贯线方程的时候需要用到这个转换关系。
  这里我就直接说结论了，详细的解释请参考冯乐乐的《Unity Shader入门精要》P69。

如果我们已知坐标系B的3个坐标轴在A坐标系下的表示${\overrightarrow{x}}_B$、${\overrightarrow{y}}_B$、${\overrightarrow{z}}_B$，以及原点位置$O_B$，那么坐标系B到坐标系A的变换矩阵为
$${\vec{M}}_{B\to A}=\left[\begin{array}{cccc}|& |& |& |\\ {\overrightarrow{x}}_B& {\overrightarrow{y}}_B& {\overrightarrow{z}}_B& O_B\\ |& |& |& |\\ 0& 0& 0& 1\end{array}\right]$$
其中，| 表示按列展开。

  根据上面的结论，我们需要先要求出O^’x’y’z’几个坐标轴 $\overrightarrow{O^{\prime }{x}^{\prime }}$、 $\overrightarrow{O^{\prime }{y}^{\prime }}$、 $\overrightarrow{O^{\prime }{z}^{\prime }}$以及$O^{\prime }$在$Oxyz$坐标系下的表示，然后才能写出将$O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$坐标系中的点变换到$Oxyz$坐标系下的变换矩阵。
  这一点我们很容易就能办到，因为从上一节中，我们知道$O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$实质是$Oxyz$先绕x轴旋转-(180°-α)，然后再沿着$Oxyz$坐标系的x轴平移-e个单位得到的。
  注意我的表述，先绕x轴旋转$-(180°-\alpha )$，为什么是$-(180°-\alpha )$，而不是$180°-\alpha$呢？因为我们选择的参考坐标系$Oxyz$ 是左手坐标系，左手坐标系应遵循左手法则，即绕x轴旋转时，左手大拇指朝向x轴正方向，然后其他4个手指握拳弯曲的方向即是绕x轴旋转的正方向 。很显然，虽然$\overrightarrow{O^{\prime }{z}^{\prime }}$与$\overrightarrow{Oz}$之间的夹角大小为$180°-\alpha$，但是我们在使用变换矩阵计算时应该代入$-(180°-\alpha )$。
  一个向量，绕着它的参考坐标系的$x$轴旋转$\theta$的变换矩阵为
$${\vec{M}}_{绕x轴旋转\theta 角}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \theta & -\sin \theta \\ 0& \sin \theta & \cos \theta \end{array}\right]$$
  将$-(180°-\alpha )$代入$\theta$得到绕$Oxyz$的x轴旋转-(180°-α)的变换矩阵为
$${\vec{M}}_{绕Oxyz的x轴旋转-(180-\alpha )度}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos (-180+\alpha )& -\sin (-180+\alpha )\\ 0& \sin (-180+\alpha ))& \cos (-180+\alpha )\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& -\cos \alpha & \sin \alpha \\ 0& -\sin \alpha & -\cos \alpha \end{array}\right]$$

  因为$\overrightarrow{O^{\prime }{y}^{\prime }}$为$\overrightarrow{Oy}$(=$[0,1,0]$)绕x轴旋转$-(180°-\alpha )$后得到的，所以$\overrightarrow{O^{\prime }{y}^{\prime }}$在$Oxyz$坐标系下的表示为
$$\overrightarrow{O^{\prime }{y}^{\prime }}=\left[\begin{array}{ccc}1& 0& 0\\ 0& -\cos \alpha & \sin \alpha \\ 0& -\sin \alpha & -\cos \alpha \end{array}\right]\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ -\cos \alpha \\ -\sin \alpha \end{array}\right]$$
  同理可得$\overrightarrow{O^{\prime }{z}^{\prime }}$在$Oxyz$坐标系下的表示为
$$\overrightarrow{O^{\prime }{z}^{\prime }}=\left[\begin{array}{ccc}1& 0& 0\\ 0& -\cos \alpha & \sin \alpha \\ 0& -\sin \alpha & -\cos \alpha \end{array}\right]\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]=\left[\begin{array}{c}0\\ \sin \alpha \\ -\cos \alpha \end{array}\right]$$
  $\overrightarrow{O^{\prime }{x}^{\prime }}$在$Oxyz$坐标系下的表示与$\overrightarrow{Ox}$相同，仍然为
$$\overrightarrow{O^{\prime }{x}^{\prime }}=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]$$
  另外，$O’$是$O$点沿着$Oxyz$的x轴平移-e个单位得到的，所以$O^{\prime }$在$Oxyz$坐标系下的表示为
$$O^{\prime }=\left[\begin{array}{c}-e\\ 0\\ 0\end{array}\right]$$
  至此，我们便可根据冯乐乐书中给出的结论，直接写出从$O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$到$Oxyz$的变换矩阵
$${\vec{M}}_{O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }\to Oxyz}=\left[\begin{array}{cccc}|& |& |& |\\ \overrightarrow{O^{\prime }{x}^{\prime }}& \overrightarrow{O^{\prime }{y}^{\prime }}& \overrightarrow{O^{\prime }{z}^{\prime }}& \overrightarrow{O^{\prime }}\\ |& |& |& |\\ 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& -e\\ 0& -\cos \alpha & \sin \alpha & 0\\ 0& -\sin \alpha & -\cos \alpha & 0\\ 0& 0& 0& 1\end{array}\right]$$
  即，坐标系$O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$中的点$(x^{\prime },y^{\prime },z^{\prime })$变换为坐标系$Oxyz$中的$(x,y,z)$点，满足的等式为
$$\left[\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& -e\\ 0& -\cos \alpha & \sin \alpha & 0\\ 0& -\sin \alpha & -\cos \alpha & 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\\ 1\end{array}\right]$$
  即
$$\left[\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right]=\left[\begin{array}{c}{x}^{\prime }-e\\ -y^{\prime }\cos \alpha +z^{\prime }\sin \alpha \\ -y^{\prime }\sin \alpha -z^{\prime }\cos \alpha \\ 1\end{array}\right]$$
  亦即
$$\{\begin{array}{l}x=x^{\prime }-e\\ y=-y^{\prime }\cos \alpha +z^{\prime }\sin \alpha \\ z=-y^{\prime }\sin \alpha -z^{\prime }\cos \alpha \end{array}\text{\hspace{1em}(1)}$$

2.2.3 相贯线方程

  根据2.1节，我们可知灰色圆柱体的方程为
$$\{\begin{array}{l}{x}^2+y^2=R^2\\ z=z\end{array}$$
  蓝色圆柱体的方程为
$$\{\begin{array}{l}{y}^{\prime 2}+z^{\prime 2}=r^2\\ x^{\prime }=x^{\prime }\end{array}$$
  相贯线的方程其实很简单，相贯线上的点只要同时满足以上两个圆柱体的方程就行了，即
$$\{\begin{array}{l}{x}^2+y^2=R^2\\ x^{\prime 2}+y^{\prime 2}=r^2\end{array}\text{\hspace{1em}(2)}$$
  由(1)式，我们可以解出
$$\{\begin{array}{l}{x}^{\prime }=x+e\\ y^{\prime }=-y\cos \alpha -z\sin \alpha \end{array}$$
  将$x^{\prime }$、$y^{\prime }$代入(2)式解得
$$\{\begin{array}{l}y=\pm \sqrt{R^2-x^2}\\ z=\frac{1}{\sin \alpha }(\pm \sqrt{r^2-(x+e{)}^2}-y\cos \alpha )\end{array}\text{\hspace{1em}(3)}$$
  其中，$x\in [-e-r，-e+r]$，$\alpha \in (0°，180°)$。
  现在我们已经知道了$x$的范围，只需将$x$在此范围内不断变化，并代入公式(3)即可以算出$y$、$z$，也就求出了相贯线上所有点的坐标。
  你可能会问，为什么$x$的范围范围是$[-e-r，-e+r]$呢？其实很简单，是根据正视图来的，这个相信大家一看图就明白，就不多说了。

3 代码实现

  历经千辛万苦，终于把相贯线方程给列出来了。接下来就是写代码了。这个就不多说了，直接上码。
  项目链接：https://pan.baidu.com/s/17lWqUNDaRHnklZUWgKvB2Q
  提取码：d2u7

// 把画线的代码单独提出来看看，如果大家没使用过，去查查API吧。
mesh.vertices = vertices;
mesh.SetIndices(indices, MeshTopology.Lines, 0);

  完整代码如下。

using UnityEngine;

public class IntersectingLine : MonoBehaviour
{
    private Transform momPipe;                                                // 母管(灰色圆柱体).
    private Transform sonPipe;                                                // 子管(蓝色圆柱体).

    [Header("大圆柱体半径")]
    public float R;
    [Header("小圆柱体半径")]
    public float r;
    [Space(5)]
    public float e;
    public float alpha;

    private Mesh mesh;
    private MeshFilter mf;
    private Material lineMat;

    #region Unity_Method

    private void Start()
    {
        Init();

        GetParameter();
        SetOPos();
        CreateIntersectingLine(R, r, e, alpha, mesh);
    }

    private void Update()
    {
        if (Input.GetKeyDown(KeyCode.Space))
        {
            GetParameter();
            SetOPos();
            CreateIntersectingLine(R, r, e, alpha, mesh);
            Debug.Log("Recalculate intersecting line over.");
        }
    }

    #endregion

    #region Private_Method

    private void Init()
    {
        momPipe = GameObject.Find("MomPipe").transform;
        sonPipe = GameObject.Find("SonPipe").transform;

        mesh = new Mesh();
        MeshRenderer mr = gameObject.AddComponent<MeshRenderer>();
        mf = gameObject.AddComponent<MeshFilter>();
        mf.mesh = mesh;
        lineMat = new Material(Shader.Find("Unlit/Color"));
        lineMat.SetColor("_Color", Color.red);
        mr.material = lineMat;
    }

    /// <summary>
    /// 获取参数.
    /// </summary>
    private void GetParameter()
    {
        R = momPipe.localScale.x * 0.5f;
        r = sonPipe.localScale.x * 0.5f;
        e = Mathf.Abs(momPipe.position.x - sonPipe.position.x);
        alpha = Vector3.Angle(-momPipe.forward, sonPipe.forward);
    }

    /// <summary>
    /// 将自身挪动到O点.
    /// </summary>
    private void SetOPos()
    {
        // 把自身挪动至O点. 这一部分只适用于当前项目，如果灰色圆柱体有旋转需要另外处理。
        Vector3 dir = sonPipe.forward;
        // 根据蓝色圆柱体的轴线求出O'点对应的z值
        float z = -dir.z / dir.y * sonPipe.position.y + sonPipe.position.z;
        transform.position = new Vector3(momPipe.position.x, 0, z);
    }

    /// <summary>
    /// 创建相贯线.
    /// </summary>
    /// <param name="R">大圆柱体的半径.</param>
    /// <param name="r">小圆柱体的半径.</param>
    /// <param name="e">两圆柱体的偏移距离.</param>
    /// <param name="alpha">两圆柱体的夹角，角度值，范围为(0°, 180°).</param>
    /// <param name="mesh">相贯线的网格.</param>
    private bool CreateIntersectingLine(float R, float r, float e, float alpha, Mesh mesh)
    {
        if (R < r || (r + Mathf.Abs(e) > R) || alpha >= 180 || alpha <= 0)
        {
            Debug.LogError("Parameter error, please check." + " R: " + R.ToString("F1") + " r: " + r.ToString("F1") + " e: " + e.ToString("F1") + " alpha: " + alpha.ToString("F1"));
            return false;
        }

        int vertexCount = 100;                                                // 顶点数.
        Vector3[] vertices = new Vector3[vertexCount];
        int[] indices = new int[vertexCount * 2];

        float deltaRad = 2 * Mathf.PI / vertexCount;

        float alphaRad = alpha * Mathf.Deg2Rad;
        float sinAlpha = Mathf.Sin(alphaRad);
        float cosAlpha = Mathf.Cos(alphaRad);

        for (int i = 0; i < vertexCount; i++)
        {
            float rad = deltaRad * i;
            float x = -e + r * Mathf.Cos(rad);
            // 两个圆柱体相贯有两条相贯线，这里我们只取y为正的那一条
            float y = Mathf.Sqrt(R * R - x * x);
            // 取一下绝对值。由于float精度问题，当r = x+e 时，可能会使temp=-0.00000001，下面开方导致NaN
            float temp = Mathf.Abs(r * r - (x + e) * (x + e));
            float z = 0f;
            if (rad > Mathf.PI)
            {
                z = 1.0f / sinAlpha * (Mathf.Sqrt(temp) - y * cosAlpha);
            }
            else
            {
                z = 1.0f / sinAlpha * (-Mathf.Sqrt(temp) - y * cosAlpha);
            }

            // 注意，我们这里算出的x, y, z是以Oxyz坐标系为参考的，不是世界坐标更不是圆柱体的局部坐标.
            vertices[i] = new Vector3(x, y, z);

            indices[2 * i] = i;
            if (i == vertexCount - 1)
            {
                indices[2 * i + 1] = 0;
            }
            else
            {
                indices[2 * i + 1] = i + 1;
            }
        }

        mesh.vertices = vertices;
        mesh.SetIndices(indices, MeshTopology.Lines, 0);
        return true;
    }

    #endregion
}

  ps：这篇文章花了我周末两天时间，真心觉得写文章非常不容易，但是也发现写文章能够帮助自己大大加深对相关知识的理解，可能这也是大佬们坚持写技术博客的原因之一吧。希望能一直坚持下去，期待未来的自己。