向量点乘与叉乘

Mar 5, 2024 · 5 min read ·

Share on:

Overview

向量点乘，叉乘的概念，以及他们的应用及相关C++代码的实现。

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. 共线：$\vec{A}$=k*$\vec{B}$，其中k是一个标量；叉积是零向量（仅适用于三维空间）；对应坐标的比率相等。
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}$
如何验证C是投影点：
1. 验证其共线性：$\vec{AC} = k \vec{AB}$或$\vec{AC} \times \vec{AB}=\vec{0}$
2. 验证垂直：$\vec{PC}\cdot \vec{AB} = 0$
如何计算点P到线AB的距离d
1. 叉积的范数是由两个向量张成的平行四边形的面积$（\vec{AB} \times \vec{AP}）$
2. 基于几何原理，这个面积也等于距离（高）乘以边长 $d * |\vec{AB}|$
3. 所以$d = |\vec{AB} \times \vec{AP}|/|\vec{AB}|$

3 代码实现

第一个版本代码，不用额外的库，手搓一些Utility函数，透彻了解原理：

  1#include<iostream>
  2#include<cmath>
  3using namespace std;
  4
  5struct Point
  6{
  7    double x, y, z;
  8    // Overloading the multiplication operator
  9    Point operator*(double k) const 
 10    {
 11        return {k*x, k*y, k*z};
 12    }
 13    Point operator+(Point A) const
 14    {
 15        return {A.x + x, A.y + y, A.z + z};
 16    }
 17    bool operator==(Point A) const
 18    {
 19        if (A.x == x and A.y == y and A.z == z)
 20        {
 21            return true;
 22        }
 23        else
 24        {
 25            return false;
 26        }
 27    }
 28};
 29
 30double dotProduct(Point A, Point B)
 31{
 32    return A.x * B.x + A.y * B.y + A.z * B.z;
 33}
 34
 35Point crossProduct(Point A, Point B)
 36{
 37    return {A.y * B.z - A.z * B.y,
 38            A.x * B.z - A.z * B.x,
 39            A.x * B.y - A.y * B.x};
 40}
 41
 42float calcNorm(Point A)
 43{
 44    return sqrt(A.x * A.x + A.y * A.y + A.z * A.z);
 45}
 46
 47Point calcProjection(Point A, Point B, Point P)
 48{
 49    Point AB = {B.x - A.x, B.y - A.y, B.z - A.z};
 50    Point AP = {P.x - A.x, P.y - A.y, P.z - A.z};
 51    double dot_product = dotProduct(AB, AP);
 52    double k = dot_product / dotProduct(AB, AB);
 53    Point C = A + (AB * k);
 54    return C;
 55}
 56
 57bool verifyProjection(Point A, Point B, Point P, Point C)
 58{
 59    Point AC = {C.x - A.x, C.y - A.y, C.z - A.z};
 60    Point AB = {B.x - A.x, B.y - A.y, B.z - A.z};
 61    Point PC = {C.x - P.x, C.y - P.y, C.z - P.z};
 62    double dot_product = dotProduct(PC, AB);
 63    Point cross_product = crossProduct(AC, AB);
 64    Point zero_vec = {0, 0, 0};
 65    if (dot_product == 0 and cross_product == zero_vec)
 66    {
 67        return true;
 68    }
 69    else
 70    {
 71        return false;
 72    }
 73}
 74
 75float calcDistance(Point A, Point B, Point P)
 76{
 77    Point AB = {B.x - A.x, B.y - A.y, B.z - A.z};
 78    Point AP = {P.x - A.x, P.y - A.y, P.z - A.z};
 79    Point cross_product = crossProduct(AB, AP);
 80    float area_parallelogram = calcNorm(cross_product);
 81    return (area_parallelogram / calcNorm(AB));
 82}
 83
 84int main()
 85{
 86    // Line segment AB
 87    Point A = {0, 0, 0};
 88    Point B = {4, 0, 0};
 89    // Point P
 90    Point P = {5, 8, 0};
 91    // Project P to AB and get point C
 92    Point C = calcProjection(A, B, P);
 93    cout << "Projection Point C: (" << C.x << ", " << C.y << ", " << C.z << ")" << endl;
 94    if (verifyProjection(A, B, P, C))
 95    {
 96        cout << "Correct!" << endl;
 97    }
 98    else
 99    {
100        cout << "Incorrect." << endl;
101    }
102    cout << "Distance from P to AB is: " << calcDistance(A, B, P) << endl;
103    return 0;
104}

另外一个版本，通过使用Eigen库来避免自己写Utility函数，行数大大减少（君子生非异也，善假于物也。）：

 1#include <iostream>
 2#include <Eigen/Dense>
 3using namespace std;
 4using namespace Eigen;
 5
 6Vector3d calcProjection(Vector3d A, Vector3d B, Vector3d P)
 7{
 8    Vector3d AB = B - A;
 9    Vector3d AP = P - A;
10    float AC_norm = AB.dot(AP) / AB.norm();
11    Vector3d C = A + AC_norm / AB.norm() * AB;
12    return C;
13}
14
15bool verifyProjection(Vector3d A, Vector3d B, Vector3d P, Vector3d C)
16{
17    Vector3d AB = B - A;
18    Vector3d PC = P - C;
19    Vector3d AC = C - A;
20    Vector3d zero_vec = {0, 0, 0};
21    Vector3d cross_product = AB.cross(AC);
22    float dot_product = PC.dot(AB);
23    if (dot_product == 0 and cross_product == zero_vec)
24    {
25        return true;
26    }
27    else
28    {
29        return false;
30    }
31}
32
33float calcDistance(const Vector3d A, const Vector3d B, const Vector3d P)
34{
35    Vector3d AB = B - A;
36    Vector3d AP = P - A;
37    Vector3d cross_product = AB.cross(AP);
38    float area_parallelogram = cross_product.norm();
39    return area_parallelogram / AB.norm();
40}
41
42int main()
43{
44    Eigen::Vector3d A = {0, 0, 0};
45    Eigen::Vector3d B = {2, 0, 0};
46    Eigen::Vector3d P = {1, 3, 0};
47
48    Vector3d C = calcProjection(A, B, P);
49    cout << "Projection point is: " << C.x() << ", " << C.y() << ", " << C.z() << endl;
50    if (verifyProjection(A, B, P, C))
51    {
52        cout << "Verification passes!" << endl;
53    }
54    else
55    {
56        cout << "Verification failed." << endl;
57    }
58    cout << "Distance from P to AB is: " << calcDistance(A, B, P) << endl;
59}

有问题欢迎一起评论交流，我会回复或更新文章，谢谢！

上述步骤的代码源文件：

comments powered by Disqus