SP7 BULK - The Bulk!

· · 题解

题目传送门:BULK - The Bulk!

解释: 

首先,找到已知为外部面的一个面,然后选择该面上的任意一条边。将面和属于它的边的组合称为“面边”。然后,从初始面边开始进行搜索(BFS或DFS或其他方式),为我们遇到的所有面边分配一个方向。这样我们就可以将所有外部面定向。然后,我们可以计算由这些面界定的多面体的体积。

多面体可能包含孔,因此只要在第一步之后仍然存在未访问的面边,就重复选择一个面边,从它开始搜索,并计算孔的体积。为了得到多面体的最终体积,需要减去所有孔的体积。注意,问题描述告诉我们多面体是连通的,因此孔本身不能包含岛屿。这保证了在第一段之后剩下的任何面边都属于孔(它们的多面体体积必须被减去)。

一旦所有面都定向了,计算体积就很容易了:我们选择一个任意的参考点,它不需要在多面体内部(例如,原点),对每个面进行三角化,然后求和所有有符号的四面体体积。我们可以在完成后取绝对值,或者确保第一个面边被定向“正确”(最好只是做前者)。

困难的部分在于分配所有方向。显然,一旦为特定的面边选择了一个方向,我们就可以很容易地为同一面上的所有其他边选择方向;但问题是如何到达共享相同边的其他面。如果我们只需要找到精确匹配(即,一个面描述中有连续的顶点 AB,另一个面描述也是如此),那么使用哈希表就很容易实现。但我们还必须考虑以下情况:

+--+--+

| | |

| +--+

| | |

+--+--+

在这里,本应该是一个面被任意地细分为3个面(问题描述允许这样做)。较小的两个面只与彼此完全匹配。因此,为了让搜索确定这两个面与大面(和多面体中的其他面)“连接”,我们需要更加小心。当我们为特定的面边 E 分配一个方向时,让 UV 成为两个端点。寻找属于其他面的边,这些边与点 U 相交并且具有与 E 相同的单位方向向量;这些是我们需要探索的边。然后,对 V 也进行同样的操作。

#include<bits/stdc++.h>
using namespace std;

struct point { int x, y, z; };
using face = vector<point>;
// Face-edges are oriented.
struct face_edge { int f, v, d; /* d = +/- 1 */ };

// unit vector codes:
// +/- 1 = +/- \hat{x}
// +/- 2 = +/- \hat{y}
// +/- 3 = +/- \hat{z}

struct ray {
    point pt;
    int u;  // unit vector code, defined above
    friend bool operator==(const ray& r1, const ray& r2) {
        return r1.pt.x == r2.pt.x && r1.pt.y == r2.pt.y &&
               r1.pt.z == r2.pt.z && r1.u == r2.u;

    }
};
namespace std {
template <> struct hash<ray> {
    size_t operator()(const ray& r) const {
        const long long w = r.pt.x + (r.pt.y << 10) + (r.pt.z << 20);
        return std::hash<long long>{}((w << 3) + r.u);
    }
};
}

using lookup_table = unordered_multimap<ray, face_edge>;

int get_uvec(point p1, point p2) {
    if (p1.x < p2.x) return 1;
    if (p1.x > p2.x) return -1;
    if (p1.y < p2.y) return 2;
    if (p1.y > p2.y) return -2;
    if (p1.z < p2.z) return 3;
    return -3;
}

long long volume(point p1, point p2, point p3) {
    long long r1 = p1.x * (p2.y * p3.z - p2.z * p3.y);
    long long r2 = p2.x * (p3.y * p1.z - p3.z * p1.y);
    long long r3 = p3.x * (p1.y * p2.z - p1.z * p2.y);
    return r1 + r2 + r3;
}

// the indicated face-edge is given positive orientation
long long dfs(const vector<face>& faces, vector<int>& orientation,
              face_edge fe, lookup_table& M, int& unvisited) {
    if (orientation[fe.f]) {
        if (orientation[fe.f] != fe.d) {
            cerr << "inconsistent orientations for face " << fe.f << '\n';
            abort();
        }
        return 0;
    }
    orientation[fe.f] = fe.d;
    while (unvisited < faces.size() && orientation[unvisited]) ++unvisited;
    const auto& face = faces[fe.f];
    long long result = 0;
    for (int i = 0; i < face.size(); i++) {
        for (int d = -1; d <= 1; d += 2) {
            const auto j = (i + d + face.size()) % face.size();
            const auto u = get_uvec(face[i], face[j]);
            const auto range = M.equal_range({face[i], u});
            if (range.first == range.second) {
                cerr << "current edge not found in lookup table\n";
                abort();
            }
            for (auto it = range.first; it != range.second; ++it) {
                const auto& new_fe = it->second;
                if (new_fe.f == fe.f) continue;
                result += dfs(faces, orientation,
                              {new_fe.f, new_fe.v, -new_fe.d * fe.d * d},
                              M, unvisited);
            }
        }
        if (i > 0 && i + 1 < face.size()) {
            // add a triangle to the result
            result += fe.d * volume(face[0], face[i], face[i+1]);
        }
    }
    return result;
}

void do_testcase() {
    int F; cin >> F;
    vector<face> faces;
    vector<int> orientation;  // 0 = unassigned, +1/-1 = assigned
    face_edge initial_fe = {-1, -1, 1};
    int least_z = 2000;
    for (int i = 0; i < F; i++) {
        face face;
        int P; cin >> P;
        while (P--) {
            int x, y, z; cin >> x >> y >> z;
            face.push_back(point{x, y, z});
            if (z < least_z) {
                least_z = z;
                initial_fe.f = i;
                initial_fe.v = face.size() - 1;
            }
        }
        faces.push_back(move(face));
        orientation.push_back(0);
    }
    // build lookup table
    lookup_table M;
    for (int i = 0; i < faces.size(); i++) {
        const auto& face = faces[i];
        for (int j = 0; j < face.size(); j++) {
            const int k = (j + 1) % face.size();
            const auto u = get_uvec(face[j], face[k]);
            M.emplace(ray{face[j], u}, face_edge{i, j, +1});
            M.emplace(ray{face[k], -u}, face_edge{i, k, -1});
        }
    }
    int unvisited = 0;
    auto vol = abs(dfs(faces, orientation, initial_fe, M, unvisited));
    while (unvisited < faces.size()) {
        const face_edge fe = {unvisited, 0, 1};
        vol -= abs(dfs(faces, orientation, fe, M, unvisited));
    }
    cout << "The bulk is composed of " << vol/6 << " units.\n";
}
int main() {
    ios::sync_with_stdio(false);
    int T; cin >> T;
    while (T--) {
        do_testcase();
    }
}