yukicoder - No.622
No.622 点と三角柱の内外判定
外積を使う
点が平面のどちら側に存在するか判定する
まずは二次元で
点 (1, 0.5) と 点 (-1, 1) が直線 y = 2x + 1 のどちら側にあるかという問題で,
これは y = 2x + 1 を変形した 2x - y + 1 = 0 の左辺の二変数関数
f(x, y) = 2x - y + 1 の符号を考えればいい.
f(x, y) は直線上で 0 になり,直線を隔てて符号が入れ替わることが少し実験すればわかる.
で,プラマイがどっちだったか自分もよく忘れるので備忘録も兼ねてメモしておくと,
直線の法線ベクトルの向きがプラスで,反対側がマイナスになる.
法線ベクトルは f(x, y) の (x, y) の係数 (2, -1) を並べたものなのでわかりやすい.(図では正規化してある)
ここで,図の直線と平行な直線の方程式が,$\vec{x} = (x, y)$ として,
$\vec{n} \cdot \vec{x} = k$
と表せることにも留意しておく.
同様のことは三次元でも考えることができて,
こんな感じになる.適当に察してください.
やっと問題ですが,
3次元上の4点A,B,C,Dが与えられる。点Dが、三角形ABCを、三角形ABCが存在する平面(平面ABC)と垂直な方向に伸ばした三角柱の内部に存在するかどうか判定せよ。勿論、平面ABCがXY平面と平行であるとは限らない。
ということで,図にします.
$\vec{n} = \vec{AB} \times \vec{AC}$ とすると,これは平面 ABC の法線ベクトルになる.右ねじの向きとか適当に思い出してください.
で,D が三角柱の内部にあるかどうかの判定ですが,これは
- AB を通り,平面 ABC に垂直な平面 $P_{1}$
- BC を通り,平面 ABC に垂直な平面 $P_{2}$
- CA を通り,平面 ABC に垂直な平面 $P_{3}$
の 3 つの平面に囲まれていれば YES,そうでなければ NO になる.
まず平面 $P_{1}$ の方程式を求める.
$\vec{m_{1}} = \vec{n} \times \vec{AB}$
とおけば,これは点 $C$ を含む (図の奥の) 方向への平面 $P_{1}$ の法線ベクトルである.
ゆえに平面 $P_{1}$ に平行な平面の方程式は
$\vec{m_{1}} \cdot \vec{x} = k$
と表すことができて,(さっきどこかで一瞬説明した)
これが点 A を通ることから,
$k = \vec{m_{1}} \cdot \vec{a}$
と k が求まる.(ここで $\vec{a}$ は A の位置ベクトルを表す)
これをもとの k に代入して,
$\vec{m_{1}} \cdot \vec{x} = \vec{m_{1}} \cdot \vec{a} \quad \Leftrightarrow \quad \vec{m_{1}} \cdot ( \vec{x} - \vec{a} ) = 0$
が求める方程式になる.だいぶ綺麗になった.
同様のことを他の平面に対しても行うと,3 つの平面の方程式
$\vec{m_{1}} \cdot ( \vec{x} - \vec{a} ) = 0$
$\vec{m_{2}} \cdot ( \vec{x} - \vec{b} ) = 0$
$\vec{m_{3}} \cdot ( \vec{x} - \vec{c} ) = 0$
が得られる.
点 D がこの 3 つの平面に囲まれているための条件は,法線 $\vec{m_{1}}, \vec{m_{2}}, \vec{m_{3}}$ が全て内側に向いていることを考えると,
$\vec{m_{1}} \cdot ( \vec{d} - \vec{a} ) > 0$
$\vec{m_{2}} \cdot ( \vec{d} - \vec{b} ) > 0$
$\vec{m_{3}} \cdot ( \vec{d} - \vec{c} ) > 0$
と表すことができる.
あとはこれをプログラムに書けばいい.
#include "bits/stdc++.h" using namespace std; namespace NVec3 { struct Vec3 { double v[3]; Vec3() { v[0] = v[1] = v[2] = 0; } Vec3(double x, double y, double z) { v[0] = x; v[1] = y; v[2] = z; } double& operator[] (size_t n) & { return v[n]; } const double& operator[] (size_t n) const& { return v[n]; } double operator[] (size_t n) const&& { return move(v[n]); } Vec3& operator+=(const Vec3& w) { v[0] += w[0]; v[1] += w[1]; v[2] += w[2]; return *this; } Vec3& operator-=(const Vec3& w) { v[0] -= w[0]; v[1] -= w[1]; v[2] -= w[2]; return *this; } double operator*=(const Vec3& w) { return v[0] * w[0] + v[1] * w[1] + v[2] * w[2]; } template<typename T> Vec3& operator/=(const T& x) { v[0] /= x; v[1] /= x; v[2] /= x; return *this; } friend istream& operator >> (istream& is, Vec3& v) { is >> v[0] >> v[1] >> v[2]; return is; } friend ostream& operator << (ostream& os, Vec3& v) { os << "{" << v[0] << ", " << v[1] << ", " << v[2] << "}"; return os; } double length() { return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); } }; Vec3 operator+(const Vec3& t1, const Vec3& t2) { return Vec3(t1) += t2; } Vec3 operator-(const Vec3& t1, const Vec3& t2) { return Vec3(t1) -= t2; } double operator*(const Vec3& t1, const Vec3& t2) { return Vec3(t1) *= t2; } template<typename T> Vec3 operator/(const Vec3& t, const T& x) { return Vec3(t) /= x; } Vec3 cross_product(const Vec3& v1, const Vec3& v2) { return Vec3( v1[1] * v2[2] - v1[2] * v2[1], v1[2] * v2[0] - v1[0] * v2[2], v1[0] * v2[1] - v1[1] * v2[0] ); } } int main() { //clock_t start, end; //start = clock(); using namespace NVec3; Vec3 A, B, C, D; cin >> A >> B >> C >> D; Vec3 n = cross_product(B - A, C - A); Vec3 m1 = cross_product(n, B - A); Vec3 m2 = cross_product(n, C - B); Vec3 m3 = cross_product(n, A - C); if (m1 * (D - A) > 0 && m2 * (D - B) > 0 && m3 * (D - C) > 0) cout << "YES" << endl; else cout << "NO" << endl; //end = clock(); //printf("%d msec.\n", end - start); return 0; }
謎のこだわりでベクトルのクラスを作ってしまった.でもおかげで main は超スッキリしてる