功能實現的數學模型選擇 -- (點與線段交點的例子)

　　在實現不少數學相關的功能的時候, 數學模型的選擇尤其重要, 一個是在準確性上, 一個是在擴展性上.

好比最近我要計算一個射線跟一個線段的交點, 初中的時候就學過, 直線和直線外一點的交點怎樣計算, 這裏

直接上已經寫完的兩段代碼.

簡圖

 

　　由於是計算機, 確定是經過兩點來肯定直線方程了, 首先是用點斜式的計算方法( 原始方程 y = kx + b, 推導過程略 ): 

        /// <summary>
        /// 點斜式推理出的交點方程
        /// </summary>
        /// <param name="ray"></param>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="point"></param>
        /// <returns></returns>
        [System.Obsolete("Incorrect", false)]
        public static bool IntersectRayAndLineSegment_XOZ_PointSlope(Ray ray, Vector3 p1, Vector3 p2, out Vector3 point)
        {
            var p3 = new Vector3(ray.origin.x, 0, ray.origin.z);

            float k1 = (p2.z - p1.z) / (p2.x - p1.x);
            float k2 = -1.0f / k1;
            float b1 = p1.z - (k1 * p1.x);
            float b2 = p3.z - (k2 * p3.x);

            float x = (b2 - b1) / (k1 - k2);
            float z = x * k1 + b1;

            point = new Vector3(x, 0, z);
            var rayDir = ray.direction;
            rayDir.y = 0;

            bool intersect = Vector3.Dot(point - p3, rayDir) > 0;
            if(intersect)
            {
                intersect = (x >= Mathf.Min(p1.x, p2.x) && x <= Mathf.Max(p1.x, p2.x)) && (z >= Mathf.Min(p1.z, p2.z) && z <= Mathf.Max(p1.z, p2.z));
            }
            return intersect;
        }

 

　　而後是兩點式的計算方法( 原始方程 (y-y1)/(x-x1) = (y-y2)/(x-x2), 推導過程略 ):

    /// <summary>
    /// 兩點式推出的交點方程
    /// </summary>
    /// <param name="ray"></param>
    /// <param name="p1"></param>
    /// <param name="p2"></param>
    /// <param name="point"></param>
    /// <returns></returns>
    public static bool IntersectRayAndLineSegment_XOZ_TwoPoint(Ray ray, Vector3 p1, Vector3 p2, out Vector3 point)
    {
        point = Vector3.zero;
        Vector3 p3 = new Vector3(ray.origin.x, 0, ray.origin.z);
        Vector3 p4 = ray.GetPoint(1.0f);
        p4.y = 0;

        float a1, b1, c1;
        PointToLineEquation_2D_XOZ(p1, p2, out a1, out b1, out c1);
        float a2, b2, c2;
        PointToLineEquation_2D_XOZ(p3, p4, out a2, out b2, out c2);

        float D = a1 * b2 - a2 * b1;
        if(D == 0)
        {
            return false;
        }
        float D1 = -c1 * b2 + c2 * b1;
        float D2 = -c2 * a1 + a2 * c1;
        point = new Vector3(D1 / D, 0, D2 / D);
        var rayDir = ray.direction;
        rayDir.y = 0;

        bool intersect = Vector3.Dot(point - p3, rayDir) > 0;
        if(intersect)
        {
            intersect = (point.x >= Mathf.Min(p1.x, p2.x) && point.x <= Mathf.Max(p1.x, p2.x)) && (point.z >= Mathf.Min(p1.z, p2.z) && point.z <= Mathf.Max(p1.z, p2.z));
        }
        return intersect;
    }
    public static void PointToLineEquation_2D_XOZ(Vector3 p1, Vector3 p2, out float a, out float b, out float c)
    {
        float x1 = p1.x;
        float y1 = p1.z;

        float x2 = p2.x;
        float y2 = p2.z;

        a = y2 - y1;
        b = x1 - x2;
        c = (y1 * x2) - (y2 * x1);
    }

 

　　在大部分狀況下它們的結果是正確的, 但是在線段是平行於座標軸的狀況下, 點斜式的結果就不對了, 往回看計算過程: 

            var p3 = new Vector3(ray.origin.x, 0, ray.origin.z);

            float k1 = (p2.z - p1.z) / (p2.x - p1.x);
            float k2 = -1.0f / k1;
            float b1 = p1.z - (k1 * p1.x);
            float b2 = p3.z - (k2 * p3.x);

            float x = (b2 - b1) / (k1 - k2);
            float z = x * k1 + b1;

 

　　點斜式在剛開始就計算了斜率k, 而後k被做爲了分母參與計算, 這樣在平行於x軸的時候斜率就是0, 被除以後獲得無窮大(或無窮小), 在平行y軸的時候( 兩點的x相等 ), k直接就是無窮大(或無窮小).

因此點斜式在計算平行於座標軸的狀況就不行了. 點斜式的問題就在於過早進行了除法計算, 並且分母還可能爲0, 這在使用計算機的系統裏面天生就存在重大隱患.

 

　　而後是兩點式, 它的過程是由兩點式推算出通常方程的各個變量( 通常方程 ax + by + c = 0 ), 而後聯立兩個通常方程來進行求解, 它的穩定性就體如今這裏:

        a = y2 - y1;
        b = x1 - x2;
        c = (y1 * x2) - (y2 * x1);

　　它的初始計算徹底不涉及除法, 第一步天生就是計算上穩定的.

　　而後是計算聯立方程的方法, 這裏直接用了二元一次線性方程組的求解:

        float D = a1 * b2 - a2 * b1;
        if(D == 0)
        {
            return false;
        }
        float D1 = -c1 * b2 + c2 * b1;
        float D2 = -c2 * a1 + a2 * c1;
        point = new Vector3(D1 / D, 0, D2 / D);

 

　　使用行列式計算的好處是很容易判斷出是否有解, D==0 時就是無解, 這也是計算穩定性的保證, 而後這裏是只計算x, z平面的, 也就是2D的,

若是要計算3D的或更高階的, 只須要把行列式和通常方程封裝成任意元的多元方法便可, 例如通常方程爲( ax + by + cz + d = 0 )甚至更高階的(  ax + by + cz + dw......+n = 0  ),

而後行列式求多元線性方程組的解就更簡單了, 這就提供了很強大的擴展性, 不只2維, 3維, 甚至N維都能提供解. 

　　因此在作功能前仔細想好怎樣作, 慎重選擇數學模型是很是重要的. 雖然這些活看似初中生也能作, 大學生也能作, 差距仍是一目瞭然的吧.

 

補充一下上面說過的多元線性方程組的解, 分紅行列式和線性方程組:

    public class Determinant
    {
        public int order { get; private set; }
        public float[,] matrix { get; private set; }

        private int m_index = 0;

        public Determinant(int order)
        {
            if(order < 2)
            {
                throw new System.Exception("order >= 2 !!");
            }
            this.order = order;
            CreateMatrix();
        }

        #region Main Funcs
        public void SetColumn(int column, params float[] values)
        {
            if(column >= 0 && column < order && values != null)
            {
                for(int i = 0, imax = Mathf.Min(order, values.Length); i < imax; i++)
                {
                    matrix[i, column] = values[i];
                }
            }
        }
        public void SetRow(int row, params float[] values)
        {
            if(row >= 0 && row < order && values != null)
            {
                for(int i = 0, imax = Mathf.Min(order, values.Length); i < imax; i++)
                {
                    matrix[row, i] = values[i];
                }
            }
        }
        public void SetRow(params float[] values)
        {
            SetRow(m_index, values);
            m_index++;
        }
        public Determinant Clone()
        {
            Determinant tag = new Determinant(this.order);
            System.Array.Copy(matrix, tag.matrix, matrix.Length);
            return tag;
        }
        public Determinant GetSubDeterminant(int row, int column)
        {
            Determinant tag = new Determinant(this.order - 1);
            for(int i = 0, tagRow = 0; i < order; i++)
            {
                if(i == row)
                {
                    continue;
                }
                for(int j = 0, tagColumn = 0; j < order; j++)
                {
                    if(j == column)
                    {
                        continue;
                    }
                    tag.matrix[tagRow, tagColumn] = this.matrix[i, j];
                    tagColumn++;
                }
                tagRow++;
            }
            return tag;
        }
        public float GetValue()
        {
            if(order == 2)
            {
                float value = matrix[0, 0] * matrix[1, 1] - matrix[0, 1] * matrix[1, 0];
                return value;
            }
            else
            {
                float value = 0.0f;
                int row = order - 1;

                for(int column = 0; column < order; column++)
                {
                    var aij = matrix[row, column] * ((row + column) % 2 > 0 ? -1 : 1);
                    var subMatrix = GetSubDeterminant(row, column);
                    var mij = subMatrix.GetValue();
                    float tempValue = (aij * mij);
                    value += tempValue;
                }

                return value;
            }
        }
        #endregion

        #region Help Funcs
        private void CreateMatrix()
        {
            matrix = new float[order, order];
        }
        #endregion
    }

    public class DeterminantEquation
    {
        public Determinant determinant { get; private set; }
        public int order
        {
            get
            {
                return determinant != null ? determinant.order : 0;
            }
        }

        public DeterminantEquation(int order)
        {
            determinant = new Determinant(order);
        }

        #region Main Funcs
        public float[] CaculateVariables(params float[] values)
        {
            var D = determinant.GetValue();
            if(D == 0)
            {
                Debug.LogWarning("This Determinant has no Result !!");
                return null;
            }

            float[] variables = new float[order];
            for(int i = 0; i < order; i++)
            {
                var subDeterminant = determinant.Clone();
                subDeterminant.SetColumn(i, values);
                var Di = subDeterminant.GetValue();
                var variable = Di / D;
                variables[i] = variable;
            }
            return variables;
        }
        #endregion
    }

　　由於多元方程組的通常方程就是 NxN 的等邊矩陣, 因此行列數量是同樣的, 用二維數組表示. 各個解經過克拉默法則計算得出.

而行列式值的計算就經過按行展開, 降階經過代數餘子式計算. 這樣就能夠計算任意階的方程組了.

　　那麼兩點式的計算代碼改一改, 就能看出擴展性了:

        /// <summary>
        /// 兩點式推出的交點方程
        /// </summary>
        /// <param name="ray"></param>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="point"></param>
        /// <returns></returns>
        public static bool IntersectRayAndLineSegment_XOZ_TwoPoint(Ray ray, Vector3 p1, Vector3 p2, out Vector3 point)
        {
            point = Vector3.zero;
            Vector3 p3 = new Vector3(ray.origin.x, 0, ray.origin.z);
            Vector3 p4 = ray.GetPoint(1.0f);
            p4.y = 0;

            float a1, b1, c1;
            PointToLineEquation_2D_XOZ(p1, p2, out a1, out b1, out c1);
            float a2, b2, c2;
            PointToLineEquation_2D_XOZ(p3, p4, out a2, out b2, out c2);

            //float D = a1 * b2 - a2 * b1;
            //if(D == 0)
            //{
            //    return false;
            //}
            //float D1 = -c1 * b2 + c2 * b1;
            //float D2 = -c2 * a1 + a2 * c1;
            //point = new Vector3(D1 / D, 0, D2 / D);

            DeterminantEquation determinantEquation = new DeterminantEquation(2);
            determinantEquation.determinant.SetRow(a1, b1);
            determinantEquation.determinant.SetRow(a2, b2);
            var variables = determinantEquation.CaculateVariables(-c1, -c2);
            if(variables == null)
            {
                return false;
            }
            point = new Vector3(variables[0], 0, variables[1]);

            var rayDir = ray.direction;
            rayDir.y = 0;

            bool intersect = Vector3.Dot(point - p3, rayDir) > 0;
            if(intersect)
            {
                intersect = (point.x >= Mathf.Min(p1.x, p2.x) && point.x <= Mathf.Max(p1.x, p2.x)) && (point.z >= Mathf.Min(p1.z, p2.z) && point.z <= Mathf.Max(p1.z, p2.z));
            }
            return intersect;
        }
        /// <summary>
        /// 通常方程變量計算
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="c"></param>
        public static void PointToLineEquation_2D_XOZ(Vector3 p1, Vector3 p2, out float a, out float b, out float c)
        {
            float x1 = p1.x;
            float y1 = p1.z;

            float x2 = p2.x;
            float y2 = p2.z;

            a = y2 - y1;
            b = x1 - x2;
            c = (y1 * x2) - (y2 * x1);
        }

　　把行列式的類( Determinant ) 改爲結構體應該在效率上會高效一些.

PS: 數學模型這裏有兩個, 第一個是經過各個點獲取聯立方程, 第二個是怎樣解聯立方程

　　點斜式獲取的是點斜式方程

　　y = k1*x + b1 

　　y = k2*x + b2 

　　

　　兩點式獲取的是標準方程

　　 a1*x + b1*y + c1 = 0

 　　a2*x + b2*y + c2 = 0

這兩個方程的計算差異上面已經說了.

解聯立方程的方法也用了兩種, 點斜式用的是通常推論, 就是普通的求解過程, 兩點式使用了線性代數(注意點斜式的聯立方程也能用線性方程解的, 只是做爲對比),

能夠預想到若是聯立方程的維數增長到不少的時候通常求解過程的難度會大大增長, 

好比

　　3元方程 -> 6個變量

　　4元方程 -> 24個變量

　　5元方程 -> 120個變量... 因此初中生最多就教到3元方程

而 DeterminantEquation 這個類就能很簡單的計算出多元方程式的解. 因此在擴展性上勝出.