http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=0115
三次元幾何
問題概要
三次元空間上で、線分 s と三角形領域 t が与えられる。
s と t とが交点をもつかどうか判定せよ。
解法
まず、sを延長した直線と、tを延長した平面との交点 p を求める (三次元幾何版の交点ライブラリが必要になる)
- p が線分 s 上にある
- p が三角形領域 t 上にある
をともに満たすかどうかを判定する。
2 については、面積計算を用いるのが簡単である。
コード
#include <bits/stdc++.h> using namespace std; typedef double DD; const DD EPS = 1e-10; struct Point3D { DD x, y, z; Point3D(DD x = 0.0, DD y = 0.0, DD z = 0.0) : x(x), y(y), z(z) {} }; Point3D operator + (const Point3D &p, const Point3D &q) {return Point3D(p.x + q.x, p.y + q.y, p.z + q.z);} Point3D operator - (const Point3D &p, const Point3D &q) {return Point3D(p.x - q.x, p.y - q.y, p.z - q.z);} Point3D operator * (const Point3D &p, DD a) {return Point3D(p.x * a, p.y * a, p.z * a);} Point3D operator * (DD a, const Point3D &p) {return Point3D(a * p.x, a * p.y, a * p.z);} Point3D operator / (const Point3D &p, DD a) {return Point3D(p.x / a, p.y / a), p.z / a;} Point3D cross (const Point3D &p, const Point3D &q) { return Point3D(p.y * q.z - p.z * q.y, p.z * q.x - p.x * q.z, p.x * q.y - p.y * q.x); } DD dot(const Point3D &p, const Point3D &q) {return p.x * q.x + p.y * q.y + p.z * q.z;} DD abs(const Point3D &p) {return sqrt(dot(p, p));} DD area(const Point3D &a, const Point3D &b, const Point3D &c) { return abs(cross(b - a, c - a)) / 2; } struct Line3D : vector<Point3D> { Line3D(const Point3D &a = Point3D(), const Point3D &b = Point3D()) { this->push_back(a); this->push_back(b); } }; struct Plane : vector<Point3D> { Plane(const Point3D &a = Point3D(), const Point3D &b = Point3D(), const Point3D &c = Point3D()) { this->push_back(a); this->push_back(b); this->push_back(c); } }; vector<Point3D> crosspointSPL(const Line3D &s, const Plane &pl) { vector<Point3D> res; Point3D ph = cross(pl[1] - pl[0], pl[2] - pl[0]); DD baseLength = dot(s[1] - s[0], ph); if (abs(baseLength) < EPS) return vector<Point3D>(); DD crossLength = dot(pl[0] - s[0], ph); DD ratio = crossLength / baseLength; if (ratio < -EPS || ratio > 1.0 + EPS) return vector<Point3D>(); res.push_back(s[0] + (s[1] - s[0]) * ratio); return res; } Point3D my, en, bar[3]; int main() { while (cin >> my.x >> my.y >> my.z >> en.x >> en.y >> en.z) { for (int i = 0; i < 3; ++i) { cin >> bar[i].x >> bar[i].y >> bar[i].z; } Line3D myen(my, en); Plane barier(bar[0], bar[1], bar[2]); vector<Point3D> cps = crosspointSPL(myen, barier); if (cps.empty()) cout << "HIT" << endl; else { Point3D cp = cps[0]; if (area(cp, bar[1], bar[2]) + area(cp, bar[2], bar[0]) + area(cp, bar[0], bar[1]) - area(bar[0], bar[1], bar[2]) > EPS) cout << "HIT" << endl; else cout << "MISS" << endl; } } }