AOJ 0115 Starship UAZ Advance

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;
        }
    }
}