【數學】向量點乘、叉乘的理論、應用及程式碼實現（C++）

MengWoods發表於2024-04-14
點乘C++

前言

我總結了一下向量點乘，叉乘的概念，以及他們的應用及相關C++程式碼的實現。blog
這類問題也是技術面試經常碰到的，一次研究透了會有收穫。

1 向量

向量具有大小和方向。
共線向量：兩個平行的向量為共線向量。

1.1 叉積 Cross Product

$$\vec{a}\times\vec{b}=|\vec{a}||\vec{b}|\sin{\theta}\vec{n}$$

$\theta$是兩個向量之間的角度，$\vec{n}$是與兩個向量都垂直的單位向量，方向遵循右手定則（右手食指從$\vec{a}$劃到$\vec{b}$，大拇指的方向）。

兩個向量的叉積結果是一個與兩者都垂直的向量。叉積的幅度值大小等於由這兩個向量為邊組成的平行四邊形的面積。當兩個向量垂直時，大小也達到最大，及矩形的面積。（這個特性決定了他可以用來計算空間中一個點到一個直線的距離，利用幾何中平行四邊形的面積同時等於底乘高，後面會介紹。）

三維空間中，叉積的結果也可以用3x3矩陣的行列式表示。

$$
\vec{a}\times \vec{b}=\det(\vec{i},\vec{j},\vec{k};a_1,a_2,a_3;b_1,b_2,b_3)\
=(a_2b_3 - a_3,b_2)\vec{i}+(a_3b_1 - a_1,b_3)\vec{j}+(a_1b_2 - a_2,b_1)\vec{k}$$

1.2 點積 Dot Product

叉積給出一個向量結果，但點積給出一個標量結果。
它將向量的相同方向投影的的長度相乘，因此使用$\cos{\theta}$將其中一個向量投影到另一個上。所以如果兩個向量成直角，那麼結果為零。點積更容易理解一些。

$$
\vec{a}\cdot \vec{b} = |\vec{a}||\vec{b}|\cos{\theta}
$$

2 實際應用

  1. 判斷兩個向量是否：

    1. 共線：$\vec{A}$=k*$\vec{B}$，其中k是一個標量；叉積是零向量（僅適用於三維空間）；對應座標的比率相等。
    2. 垂直：點積為零。

  2. 計算點P線上段AB上的投影點C座標。

    1. 向量$\vec{AB}$，$\vec{AP}$，點積結果為D。$\vec{AB} \cdot \vec{AP}=D$
    2. 推導一下：$\vec{AB} \cdot \vec{AP}=|\vec{AB}| |\vec{AP}| \cos{\theta}=|\vec{AC}| |\vec{AB}| = D$ -> $D/|\vec{AB}| = |\vec{AC}|$
    3. 求比率：$k = |\vec{AC}|/|\vec{AB}| = D/{|\vec{AB}|^2}$
    4. 最終座標根據$A$和$k$可以求得：$C = A + k\vec{AB}$

  3. 如何驗證C是投影點：

    1. 驗證其共線性：$\vec{AC} = k \vec{AB}$或$\vec{AC} \times \vec{AB}=\vec{0}$
    2. 驗證垂直：$\vec{PC}\cdot \vec{AB} = 0$

  4. 如何計算點P到線AB的距離d

    1. 叉積的範數是由兩個向量張成的平行四邊形的面積$（\vec{AB} \times \vec{AP}）$
    2. 基於幾何原理，這個面積也等於距離（高）乘以邊長 $d * |\vec{AB}|$
    3. 所以$d = |\vec{AB} \times \vec{AP}|/|\vec{AB}|$

3 程式碼實現

第一個版本程式碼，不用額外的庫，手搓一些Utility函式，透徹瞭解原理：

#include<iostream>
#include<cmath>
using namespace std;

struct Point
{
    double x, y, z;
    // Overloading the multiplication operator
    Point operator*(double k) const 
    {
        return {k*x, k*y, k*z};
    }
    Point operator+(Point A) const
    {
        return {A.x + x, A.y + y, A.z + z};
    }
    bool operator==(Point A) const
    {
        if (A.x == x and A.y == y and A.z == z)
        {
            return true;
        }
        else
        {
            return false;
        }
    }
};

double dotProduct(Point A, Point B)
{
    return A.x * B.x + A.y * B.y + A.z * B.z;
}

Point crossProduct(Point A, Point B)
{
    return {A.y * B.z - A.z * B.y,
            A.x * B.z - A.z * B.x,
            A.x * B.y - A.y * B.x};
}

float calcNorm(Point A)
{
    return sqrt(A.x * A.x + A.y * A.y + A.z * A.z);
}

Point calcProjection(Point A, Point B, Point P)
{
    Point AB = {B.x - A.x, B.y - A.y, B.z - A.z};
    Point AP = {P.x - A.x, P.y - A.y, P.z - A.z};
    double dot_product = dotProduct(AB, AP);
    double k = dot_product / dotProduct(AB, AB);
    Point C = A + (AB * k);
    return C;
}

bool verifyProjection(Point A, Point B, Point P, Point C)
{
    Point AC = {C.x - A.x, C.y - A.y, C.z - A.z};
    Point AB = {B.x - A.x, B.y - A.y, B.z - A.z};
    Point PC = {C.x - P.x, C.y - P.y, C.z - P.z};
    double dot_product = dotProduct(PC, AB);
    Point cross_product = crossProduct(AC, AB);
    Point zero_vec = {0, 0, 0};
    if (dot_product == 0 and cross_product == zero_vec)
    {
        return true;
    }
    else
    {
        return false;
    }
}

float calcDistance(Point A, Point B, Point P)
{
    Point AB = {B.x - A.x, B.y - A.y, B.z - A.z};
    Point AP = {P.x - A.x, P.y - A.y, P.z - A.z};
    Point cross_product = crossProduct(AB, AP);
    float area_parallelogram = calcNorm(cross_product);
    return (area_parallelogram / calcNorm(AB));
}

int main()
{
    // Line segment AB
    Point A = {0, 0, 0};
    Point B = {4, 0, 0};
    // Point P
    Point P = {5, 8, 0};
    // Project P to AB and get point C
    Point C = calcProjection(A, B, P);
    cout << "Projection Point C: (" << C.x << ", " << C.y << ", " << C.z << ")" << endl;
    if (verifyProjection(A, B, P, C))
    {
        cout << "Correct!" << endl;
    }
    else
    {
        cout << "Incorrect." << endl;
    }
    cout << "Distance from P to AB is: " << calcDistance(A, B, P) << endl;
    return 0;
}

另外一個版本，透過使用Eigen庫來避免自己寫Utility函式，行數大大減少（君子生非異也，善假於物也。）：

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

Vector3d calcProjection(Vector3d A, Vector3d B, Vector3d P)
{
    Vector3d AB = B - A;
    Vector3d AP = P - A;
    float AC_norm = AB.dot(AP) / AB.norm();
    Vector3d C = A + AC_norm / AB.norm() * AB;
    return C;
}

bool verifyProjection(Vector3d A, Vector3d B, Vector3d P, Vector3d C)
{
    Vector3d AB = B - A;
    Vector3d PC = P - C;
    Vector3d AC = C - A;
    Vector3d zero_vec = {0, 0, 0};
    Vector3d cross_product = AB.cross(AC);
    float dot_product = PC.dot(AB);
    if (dot_product == 0 and cross_product == zero_vec)
    {
        return true;
    }
    else
    {
        return false;
    }
}

float calcDistance(const Vector3d A, const Vector3d B, const Vector3d P)
{
    Vector3d AB = B - A;
    Vector3d AP = P - A;
    Vector3d cross_product = AB.cross(AP);
    float area_parallelogram = cross_product.norm();
    return area_parallelogram / AB.norm();
}

int main()
{
    Eigen::Vector3d A = {0, 0, 0};
    Eigen::Vector3d B = {2, 0, 0};
    Eigen::Vector3d P = {1, 3, 0};

    Vector3d C = calcProjection(A, B, P);
    cout << "Projection point is: " << C.x() << ", " << C.y() << ", " << C.z() << endl;
    if (verifyProjection(A, B, P, C))
    {
        cout << "Verification passes!" << endl;
    }
    else
    {
        cout << "Verification failed." << endl;
    }
    cout << "Distance from P to AB is: " << calcDistance(A, B, P) << endl;
}

有問題歡迎一起評論交流，我會回覆或更新文章，謝謝！

相關文章