知识点

知识点 - 计算几何基础

讲义

点

我们把点 (mathbf r) 看成从 (mathbf 0) 到 (mathbf r)的向量 (vec{mathbf r})

#define ftype long double
struct point2d {
    ftype x, y;
    point2d() {}
    point2d(ftype x, ftype y): x(x), y(y) {}
    point2d& operator+=(const point2d &t) {
        x += t.x;
        y += t.y;
        return *this;
    }
    point2d& operator-=(const point2d &t) {
        x -= t.x;
        y -= t.y;
        return *this;
    }
    point2d& operator*=(ftype t) {
        x *= t;
        y *= t;
        return *this;
    }
    point2d& operator/=(ftype t) {
        x /= t;
        y /= t;
        return *this;
    }
    point2d operator+(const point2d &t) const {
        return point2d(*this) += t;
    }
    point2d operator-(const point2d &t) const {
        return point2d(*this) -= t;
    }
    point2d operator*(ftype t) const {
        return point2d(*this) *= t;
    }
    point2d operator/(ftype t) const {
        return point2d(*this) /= t;
    }
};
point2d operator*(ftype a, point2d b) {
    return b * a;
}

点乘

1. (mathbf a cdot mathbf b = mathbf b cdot mathbf a)
2. ((alpha cdot mathbf a)cdot mathbf b = alpha cdot (mathbf a cdot mathbf b))
3. ((mathbf a + mathbf b)cdot mathbf c = mathbf a cdot mathbf c + mathbf b cdot mathbf c)

若有单位向量：

[mathbf e_x = egin{pmatrix} 1 \ 0 \ 0 end{pmatrix}, mathbf e_y = egin{pmatrix} 0 \ 1 \ 0 end{pmatrix}, mathbf e_z = egin{pmatrix} 0 \ 0 \ 1 end{pmatrix}. ]

我们定义 (mathbf r = (x;y;z)) 表示 (r = x cdot mathbf e_x + y cdot mathbf e_y + z cdot mathbf e_z).
因为

[mathbf e_xcdot mathbf e_x = mathbf e_ycdot mathbf e_y = mathbf e_zcdot mathbf e_z = 1,\ mathbf e_xcdot mathbf e_y = mathbf e_ycdot mathbf e_z = mathbf e_zcdot mathbf e_x = 0]

所以对 (mathbf a = (x_1;y_1;z_1)) 和 (mathbf b = (x_2;y_2;z_2)) 有

[mathbf acdot mathbf b = (x_1 cdot mathbf e_x + y_1 cdotmathbf e_y + z_1 cdotmathbf e_z)cdot( x_2 cdotmathbf e_x + y_2 cdotmathbf e_y + z_2 cdotmathbf e_z) = x_1 x_2 + y_1 y_2 + z_1 z_2 ]

ftype dot(point2d a, point2d b) {
    return a.x * b.x + a.y * b.y;
}
ftype dot(point3d a, point3d b) {
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

​ 一些定义：

1. Norm of (mathbf a) (长度的平方): (|mathbf a|^2 = mathbf acdot mathbf a)
2. Length of (mathbf a): (|mathbf a| = sqrt{mathbf acdot mathbf a})
3. Projection of (mathbf a) onto (mathbf b)（投影）: (dfrac{mathbf acdotmathbf b}{|mathbf b|})
4. Angle between vectors（夹角）: (arccos left(dfrac{mathbf acdot mathbf b}{|mathbf a| cdot |mathbf b|} ight))
5. 从上一点说明点乘的正负可以用来判断锐角(acute)钝角(obtuse)直角(orthogonal).

ftype norm(point2d a) {
    return dot(a, a);
}
double abs(point2d a) {
    return sqrt(norm(a));
}
double proj(point2d a, point2d b) {
    return dot(a, b) / abs(b);
}
double angle(point2d a, point2d b) {
    return acos(dot(a, b) / abs(a) / abs(b));
}

叉乘

定义：

先定义三重积triple product (mathbf acdot(mathbf b imes mathbf c)) 为上面这个平行六面体的体积，于是我们可以得到(mathbf b imes mathbf c)的大小和方向。

性质：

1. (mathbf a imes mathbf b = -mathbf b imes mathbf a)
2. ((alpha cdot mathbf a) imes mathbf b = alpha cdot (mathbf a imes mathbf b))
3. (mathbf acdot (mathbf b imes mathbf c) = mathbf bcdot (mathbf c imes mathbf a) = -mathbf acdot( mathbf c imes mathbf b))
4. ((mathbf a + mathbf b) imes mathbf c = mathbf a imes mathbf c + mathbf b imes mathbf c).
   对任意的 (mathbf r) 有:
   [mathbf rcdot( (mathbf a + mathbf b) imes mathbf c) = (mathbf a + mathbf b) cdot (mathbf c imes mathbf r) = mathbf a cdot(mathbf c imes mathbf r) + mathbf bcdot(mathbf c imes mathbf r) = mathbf rcdot (mathbf a imes mathbf c) + mathbf rcdot(mathbf b imes mathbf c) = mathbf rcdot(mathbf a imes mathbf c + mathbf b imes mathbf c) ]
   这证明了第三点 ((mathbf a + mathbf b) imes mathbf c = mathbf a imes mathbf c + mathbf b imes mathbf c)
5. (|mathbf a imes mathbf b|=|mathbf a| cdot |mathbf b| sin heta)

因为

[mathbf e_x imes mathbf e_x = mathbf e_y imes mathbf e_y = mathbf e_z imes mathbf e_z = mathbf 0,\ mathbf e_x imes mathbf e_y = mathbf e_z,~mathbf e_y imes mathbf e_z = mathbf e_x,~mathbf e_z imes mathbf e_x = mathbf e_y]

于是我们可以算出 (mathbf a = (x_1;y_1;z_1)) 和 (mathbf b = (x_2;y_2;z_2)) 的叉乘结果:

[mathbf a imes mathbf b = (x_1 cdot mathbf e_x + y_1 cdot mathbf e_y + z_1 cdot mathbf e_z) imes (x_2 cdot mathbf e_x + y_2 cdot mathbf e_y + z_2 cdot mathbf e_z) = ]

[(y_1 z_2 - z_1 y_2)mathbf e_x + (z_1 x_2 - x_1 z_2)mathbf e_y + (x_1 y_2 - y_1 x_2) ]

用行列式表达的话：

[mathbf a imes mathbf b = egin{vmatrix}mathbf e_x & mathbf e_y & mathbf e_z \ x_1 & y_1 & z_1 \ x_2 & y_2 & z_2 end{vmatrix},~acdot(b imes c) = egin{vmatrix} x_1 & y_1 & z_1 \ x_2 & y_2 & z_2 \ x_3 & y_3 & z_3 end{vmatrix} ]

二维的叉乘 (namely the pseudo-scalar product)可以被定义为

[|mathbf e_zcdot(mathbf a imes mathbf b)| = |x_1 y_2 - y_1 x_2| ]

一个直观理解方式是为了计算(|mathbf a| cdot |mathbf b| sin heta) 将 (mathbf a)转 (90^circ)得到(widehat{mathbf a}=(-y_1;x_1)).于是(|widehat{mathbf a}cdotmathbf b|=|x_1y_2 - y_1 x_2|).

point3d cross(point3d a, point3d b) {
    return point3d(a.y * b.z - a.z * b.y,
                   a.z * b.x - a.x * b.z,
                   a.x * b.y - a.y * b.x);
}
ftype triple(point3d a, point3d b, point3d c) {
    return dot(a, cross(b, c));
}
ftype cross(point2d a, point2d b) {
    return a.x * b.y - a.y * b.x;
}

直线与平面

​ 一个直线可以被表示为一个起始点(mathbf r_0) 和一个方向向量(mathbf d) ，或者两个点(mathbf a) , (mathbf b).对应的方程为

[(mathbf r - mathbf r_0) imesmathbf d=0 \ (mathbf r - mathbf a) imes (mathbf b - mathbf a) = 0. ]

​ 一个平面可以被三个点确定： (mathbf a), (mathbf b) , (mathbf c)。或者一个初始点(mathbf r_0)和一组在这个平面里的向量(mathbf d_1) , (mathbf d_2)确定：

[(mathbf r - mathbf a)cdot((mathbf b - mathbf a) imes (mathbf c - mathbf a))=0\ (mathbf r - mathbf r_0)cdot(mathbf d_1 imes mathbf d_2)=0 ]

直线交点

(l_1:mathbf r = mathbf a_1 + t cdot mathbf d_1) 带入 (l_2:(mathbf r - mathbf a_2) imes mathbf d_2=0)

[(mathbf a_1 + t cdot mathbf d_1 - mathbf a_2) imes mathbf d_2=0 quadRightarrowquad t = dfrac{(mathbf a_2 - mathbf a_1) imesmathbf d_2}{mathbf d_1 imes mathbf d_2} ]

point2d intersect(point2d a1, point2d d1, point2d a2, point2d d2) {
    return a1 + cross(a2 - a1, d2) / cross(d1, d2) * d1;
}

三个平面交点

给你三个平面的初始点 (mathbf a_i) 和法向量 (mathbf n_i) 于是得到方程：

[egin{cases}mathbf rcdot mathbf n_1 = mathbf a_1cdot mathbf n_1, \ mathbf rcdot mathbf n_2 = mathbf a_2cdot mathbf n_2, \ mathbf rcdot mathbf n_3 = mathbf a_3cdot mathbf n_3end{cases} ]

用克拉默法则解：

point3d intersect(point3d a1, point3d n1, point3d a2, point3d n2, point3d a3, point3d n3) {
    point3d x(n1.x, n2.x, n3.x);
    point3d y(n1.y, n2.y, n3.y);
    point3d z(n1.z, n2.z, n3.z); 
    point3d d(dot(a1, n1), dot(a2, n2), dot(a3, n3));
    return point3d(triple(d, y, z),
                   triple(x, d, z),
                   triple(x, y, d)) / triple(n1, n2, n3);
}

模板

两个流派，一个是向量表示直线即两个点(mathbf a) , (mathbf b)，另一个是直线方程即(a_1 x + b_1 y + c_1 = 0)