h_nosonの日記

競プロ、CTFなど

yukicoder No.132 点と平面との距離

問題
No.132 点と平面との距離 - yukicoder
3次元上の点PN個の点Q_1,Q_2,...,Q_Nが与えられる.
\text{dist}(Q_i,Q_j,Q_k)Pから点Q_i,Q_j,Q_kを通る平面への距離とした時,
\sum_{i=1}^{N-2}\sum_{j=i+1}^{N-1}\sum_{k=j+1}^N\text{dist}(Q_i,Q_j,Q_k)を求める問題.

解法
P\text{-}Q_iベクトル,P\text{-}Q_jベクトル,P\text{-}Q_kベクトルに囲まれた三角錐の体積を求め,三角形Q_iQ_jQ_kの面積で割って3をかけることで,高さつまりPから点Q_i,Q_j,Q_kを通る平面への距離を求めることができる.
Q_iの座標を(Q_{ix},Q_{iy},Q_{iz})とした時,三角錐の体積は,
\frac{1}{6}
\left|
\begin{array}{ccc}
Q_{ix}&Q_{iy}&Q_{iz}\\
Q_{jx}&Q_{jy}&Q_{jz}\\
Q_{kx}&Q_{ky}&Q_{kz}
\end{array}
\right|
となり,三角形Q_iQ_jQ_kの面積は,ベクトル

\left|
\begin{array}{ccc}
e_x&e_y&e_z\\
Q_{jx}-Q_{ix}&Q_{jy}-Q_{iy}&Q_{jz}-Q_{iz}\\
Q_{kx}-Q_{ix}&Q_{ky}-Q_{iy}&Q_{kz}-Q_{iz}
\end{array}
\right|
の大きさの半分になる.(e_x,e_y,e_zはx軸,y軸,z軸の単位ベクトル)

ソースコード

#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>
using namespace std;

#define RREP(i,s,e) for (i = s; i >= e; i--)
#define rrep(i,n) RREP(i,(int)(n)-1,0)
#define REP(i,s,e) for (i = s; i <= e; i++)
#define rep(i,n) REP(i,0,(int)(n)-1)
#define INF 100000000

typedef long long ll;

double det(double* p, double* q, double* r) {
    int i;
    double ret = 0;
    rep (i,3) ret += p[i]*q[(i+1)%3]*r[(i+2)%3];
    rep (i,3) ret -= p[i]*q[(i+2)%3]*r[(i+1)%3];
    return ret;
}

double calc_area(double* p, double* q) {
    int i;
    double ret = 0;
    ret += (p[1]*q[2]-p[2]*q[1]) * (p[1]*q[2]-p[2]*q[1]);
    ret += (p[2]*q[0]-p[0]*q[2]) * (p[2]*q[0]-p[0]*q[2]);
    ret += (p[0]*q[1]-p[1]*q[0]) * (p[0]*q[1]-p[1]*q[0]);
    return sqrt(ret);
}

int main() {
    int i, j, k, l, n;
    double p[3], a[3], b[3];
    double xyz[300][3];
    double ans = 0;
    cin >> n;
    rep (i,3) cin >> p[i];
    rep (i,n) {
        cin >> xyz[i][0] >> xyz[i][1] >> xyz[i][2];
        rep (j,3) xyz[i][j] -= p[j];
    }
    rep (i,n) REP (j,i+1,n-1) REP (k,j+1,n-1) {
        double vol = abs(det(xyz[i],xyz[j],xyz[k])/6);
        rep (l,3) a[l] = xyz[j][l] - xyz[i][l];
        rep (l,3) b[l] = xyz[k][l] - xyz[i][l];
        double area = calc_area(a,b) / 2;
        ans += vol * 3 / area;
    }
    cout << fixed << setprecision(10) << ans << endl;
    return 0;
}