我们把点 \(\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;
}
若有单位向量:
\[\mathbf e_x = \begin{pmatrix} 1 \\\ 0 \\\ 0 \end{pmatrix}, \mathbf e_y = \begin{pmatrix} 0 \\\ 1 \\\ 0 \end{pmatrix}, \mathbf e_z = \begin{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_x\cdot \mathbf e_x = \mathbf e_y\cdot \mathbf e_y = \mathbf e_z\cdot \mathbf e_z = 1,\\\mathbf e_x\cdot \mathbf e_y = \mathbf e_y\cdot \mathbf e_z = \mathbf e_z\cdot \mathbf e_x = 0\]
所以对 \(\mathbf a = (x_1;y_1;z_1)\) 和 \(\mathbf b = (x_2;y_2;z_2)\) 有
\[\mathbf a\cdot \mathbf b = (x_1 \cdot \mathbf e_x + y_1 \cdot\mathbf e_y + z_1 \cdot\mathbf e_z)\cdot( x_2 \cdot\mathbf e_x + y_2 \cdot\mathbf e_y + z_2 \cdot\mathbf 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;
}
? 一些定义:
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 a\cdot(\mathbf b\times \mathbf c)\) 为上面这个平行六面体的体积,于是我们可以得到\(\mathbf b\times \mathbf c\)的大小和方向。
因为
\[\mathbf e_x\times \mathbf e_x = \mathbf e_y\times \mathbf e_y = \mathbf e_z\times \mathbf e_z = \mathbf 0,\\\mathbf e_x\times \mathbf e_y = \mathbf e_z,~\mathbf e_y\times \mathbf e_z = \mathbf e_x,~\mathbf e_z\times \mathbf e_x = \mathbf e_y\]
于是我们可以算出 \(\mathbf a = (x_1;y_1;z_1)\) 和 \(\mathbf b = (x_2;y_2;z_2)\) 的叉乘结果:
\[\mathbf a\times \mathbf b = (x_1 \cdot \mathbf e_x + y_1 \cdot \mathbf e_y + z_1 \cdot \mathbf e_z)\times (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\times \mathbf b = \begin{vmatrix}\mathbf e_x & \mathbf e_y & \mathbf e_z \\\ x_1 & y_1 & z_1 \\\ x_2 & y_2 & z_2 \end{vmatrix},~a\cdot(b\times c) = \begin{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_z\cdot(\mathbf a\times \mathbf b)| = |x_1 y_2 - y_1 x_2|
\]
一个直观理解方式是为了计算\(|\mathbf a| \cdot |\mathbf b| \sin \theta\) 将 \(\mathbf a\)转 \(90^\circ\)得到\(\widehat{\mathbf a}=(-y_1;x_1)\).于是\(|\widehat{\mathbf a}\cdot\mathbf 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)\times\mathbf d=0 \\ (\mathbf r - \mathbf a)\times (\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)\times (\mathbf c - \mathbf a))=0\(\mathbf r - \mathbf r_0)\cdot(\mathbf d_1\times \mathbf d_2)=0
\]
\(l_1:\mathbf r = \mathbf a_1 + t \cdot \mathbf d_1\) 带入 \(l_2:(\mathbf r - \mathbf a_2)\times \mathbf d_2=0\)
\[
(\mathbf a_1 + t \cdot \mathbf d_1 - \mathbf a_2)\times \mathbf d_2=0 \quad\Rightarrow\quad t = \dfrac{(\mathbf a_2 - \mathbf a_1)\times\mathbf d_2}{\mathbf d_1\times \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\) 于是得到方程:
\[
\begin{cases}\mathbf r\cdot \mathbf n_1 = \mathbf a_1\cdot \mathbf n_1, \\\ \mathbf r\cdot \mathbf n_2 = \mathbf a_2\cdot \mathbf n_2, \\\ \mathbf r\cdot \mathbf n_3 = \mathbf a_3\cdot \mathbf n_3\end{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\)
原文:https://www.cnblogs.com/SuuT/p/11370374.html