diff options
Diffstat (limited to 'source/blender/blenlib/tests/BLI_mesh_intersect_test.cc')
| -rw-r--r-- | source/blender/blenlib/tests/BLI_mesh_intersect_test.cc | 1074 |
1 files changed, 1074 insertions, 0 deletions
diff --git a/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc b/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc new file mode 100644 index 00000000000..237905a1bf6 --- /dev/null +++ b/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc @@ -0,0 +1,1074 @@ +/* Apache License, Version 2.0 */ + +#include "testing/testing.h" + +#include <algorithm> +#include <fstream> +#include <iostream> + +#include "PIL_time.h" + +#include "BLI_array.hh" +#include "BLI_math_mpq.hh" +#include "BLI_mesh_intersect.hh" +#include "BLI_mpq3.hh" +#include "BLI_task.h" +#include "BLI_vector.hh" + +#define DO_REGULAR_TESTS 1 +#define DO_PERF_TESTS 0 + +#ifdef WITH_GMP +namespace blender::meshintersect::tests { + +constexpr bool DO_OBJ = false; + +/* Build and hold an IMesh from a string spec. Also hold and own resources used by IMesh. */ +class IMeshBuilder { + public: + IMesh imesh; + IMeshArena arena; + + /* "Edge orig" indices are an encoding of <input face#, position in face of start of edge>. */ + static constexpr int MAX_FACE_LEN = 1000; /* Used for forming "orig edge" indices only. */ + + static int edge_index(int face_index, int facepos) + { + return face_index * MAX_FACE_LEN + facepos; + } + + static std::pair<int, int> face_and_pos_for_edge_index(int e_index) + { + return std::pair<int, int>(e_index / MAX_FACE_LEN, e_index % MAX_FACE_LEN); + } + + /* + * Spec should have form: + * #verts #faces + * mpq_class mpq_class mpq_clas [#verts lines] + * int int int ... [#faces lines; indices into verts for given face] + */ + IMeshBuilder(const char *spec) + { + std::istringstream ss(spec); + std::string line; + getline(ss, line); + std::istringstream hdrss(line); + int nv, nf; + hdrss >> nv >> nf; + if (nv == 0 || nf == 0) { + return; + } + arena.reserve(nv, nf); + Vector<const Vert *> verts; + Vector<Face *> faces; + bool spec_ok = true; + int v_index = 0; + while (v_index < nv && spec_ok && getline(ss, line)) { + std::istringstream iss(line); + mpq_class p0; + mpq_class p1; + mpq_class p2; + iss >> p0 >> p1 >> p2; + spec_ok = !iss.fail(); + if (spec_ok) { + verts.append(arena.add_or_find_vert(mpq3(p0, p1, p2), v_index)); + } + ++v_index; + } + if (v_index != nv) { + spec_ok = false; + } + int f_index = 0; + while (f_index < nf && spec_ok && getline(ss, line)) { + std::istringstream fss(line); + Vector<const Vert *> face_verts; + Vector<int> edge_orig; + int fpos = 0; + while (spec_ok && fss >> v_index) { + if (v_index < 0 || v_index >= nv) { + spec_ok = false; + continue; + } + face_verts.append(verts[v_index]); + edge_orig.append(edge_index(f_index, fpos)); + ++fpos; + } + if (fpos < 3) { + spec_ok = false; + } + if (spec_ok) { + Face *facep = arena.add_face(face_verts, f_index, edge_orig); + faces.append(facep); + } + ++f_index; + } + if (f_index != nf) { + spec_ok = false; + } + if (!spec_ok) { + std::cout << "Bad spec: " << spec; + return; + } + imesh = IMesh(faces); + } +}; + +/* Return a const Face * in mesh with verts equal to v0, v1, and v2, in + * some cyclic order; return nullptr if not found. + */ +static const Face *find_tri_with_verts(const IMesh &mesh, + const Vert *v0, + const Vert *v1, + const Vert *v2) +{ + Face f_arg({v0, v1, v2}, 0, NO_INDEX); + for (const Face *f : mesh.faces()) { + if (f->cyclic_equal(f_arg)) { + return f; + } + } + return nullptr; +} + +/* How many instances of a triangle with v0, v1, v2 are in the mesh? */ +static int count_tris_with_verts(const IMesh &mesh, const Vert *v0, const Vert *v1, const Vert *v2) +{ + Face f_arg({v0, v1, v2}, 0, NO_INDEX); + int ans = 0; + for (const Face *f : mesh.faces()) { + if (f->cyclic_equal(f_arg)) { + ++ans; + } + } + return ans; +} + +/* What is the starting position, if any, of the edge (v0, v1), in either order, in f? -1 if none. + */ +static int find_edge_pos_in_tri(const Vert *v0, const Vert *v1, const Face *f) +{ + for (int pos : f->index_range()) { + int nextpos = f->next_pos(pos); + if (((*f)[pos] == v0 && (*f)[nextpos] == v1) || ((*f)[pos] == v1 && (*f)[nextpos] == v0)) { + return static_cast<int>(pos); + } + } + return -1; +} + +# if DO_REGULAR_TESTS +TEST(mesh_intersect, Mesh) +{ + Vector<const Vert *> verts; + Vector<Face *> faces; + IMeshArena arena; + + verts.append(arena.add_or_find_vert(mpq3(0, 0, 1), 0)); + verts.append(arena.add_or_find_vert(mpq3(1, 0, 1), 1)); + verts.append(arena.add_or_find_vert(mpq3(0.5, 1, 1), 2)); + faces.append(arena.add_face(verts, 0, {10, 11, 12})); + + IMesh mesh(faces); + const Face *f = mesh.face(0); + EXPECT_TRUE(f->is_tri()); +} + +TEST(mesh_intersect, OneTri) +{ + const char *spec = R"(3 1 + 0 0 0 + 1 0 0 + 1/2 1 0 + 0 1 2 + )"; + + IMeshBuilder mb(spec); + IMesh imesh = trimesh_self_intersect(mb.imesh, &mb.arena); + imesh.populate_vert(); + EXPECT_EQ(imesh.vert_size(), 3); + EXPECT_EQ(imesh.face_size(), 1); + const Face &f_in = *mb.imesh.face(0); + const Face &f_out = *imesh.face(0); + EXPECT_EQ(f_in.orig, f_out.orig); + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(f_in[i], f_out[i]); + EXPECT_EQ(f_in.edge_orig[i], f_out.edge_orig[i]); + } +} + +TEST(mesh_intersect, TriTri) +{ + const char *spec = R"(6 2 + 0 0 0 + 4 0 0 + 0 4 0 + 1 0 0 + 2 0 0 + 1 1 0 + 0 1 2 + 3 4 5 + )"; + + /* Second triangle is smaller and congruent to first, resting on same base, partway along. */ + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 6); + EXPECT_EQ(out.face_size(), 6); + if (out.vert_size() == 6 && out.face_size() == 6) { + const Vert *v0 = mb.arena.find_vert(mpq3(0, 0, 0)); + const Vert *v1 = mb.arena.find_vert(mpq3(4, 0, 0)); + const Vert *v2 = mb.arena.find_vert(mpq3(0, 4, 0)); + const Vert *v3 = mb.arena.find_vert(mpq3(1, 0, 0)); + const Vert *v4 = mb.arena.find_vert(mpq3(2, 0, 0)); + const Vert *v5 = mb.arena.find_vert(mpq3(1, 1, 0)); + EXPECT_TRUE(v0 != nullptr && v1 != nullptr && v2 != nullptr); + EXPECT_TRUE(v3 != nullptr && v4 != nullptr && v5 != nullptr); + if (v0 != nullptr && v1 != nullptr && v2 != nullptr && v3 != nullptr && v4 != nullptr && + v5 != nullptr) { + EXPECT_EQ(v0->orig, 0); + EXPECT_EQ(v1->orig, 1); + const Face *f0 = find_tri_with_verts(out, v4, v1, v5); + const Face *f1 = find_tri_with_verts(out, v3, v4, v5); + const Face *f2 = find_tri_with_verts(out, v0, v3, v5); + const Face *f3 = find_tri_with_verts(out, v0, v5, v2); + const Face *f4 = find_tri_with_verts(out, v5, v1, v2); + EXPECT_TRUE(f0 != nullptr && f1 != nullptr && f2 != nullptr && f3 != nullptr && + f4 != nullptr); + /* For boolean to work right, there need to be two copies of the smaller triangle in the + * output. */ + EXPECT_EQ(count_tris_with_verts(out, v3, v4, v5), 2); + if (f0 != nullptr && f1 != nullptr && f2 != nullptr && f3 != nullptr && f4 != nullptr) { + EXPECT_EQ(f0->orig, 0); + EXPECT_TRUE(f1->orig == 0 || f1->orig == 1); + EXPECT_EQ(f2->orig, 0); + EXPECT_EQ(f3->orig, 0); + EXPECT_EQ(f4->orig, 0); + } + int e03 = find_edge_pos_in_tri(v0, v3, f2); + int e34 = find_edge_pos_in_tri(v3, v4, f1); + int e45 = find_edge_pos_in_tri(v4, v5, f1); + int e05 = find_edge_pos_in_tri(v0, v5, f3); + int e15 = find_edge_pos_in_tri(v1, v5, f0); + EXPECT_TRUE(e03 != -1 && e34 != -1 && e45 != -1 && e05 != -1 && e15 != -1); + if (e03 != -1 && e34 != -1 && e45 != -1 && e05 != -1 && e15 != -1) { + EXPECT_EQ(f2->edge_orig[e03], 0); + EXPECT_TRUE(f1->edge_orig[e34] == 0 || + f1->edge_orig[e34] == 1 * IMeshBuilder::MAX_FACE_LEN); + EXPECT_EQ(f1->edge_orig[e45], 1 * IMeshBuilder::MAX_FACE_LEN + 1); + EXPECT_EQ(f3->edge_orig[e05], NO_INDEX); + EXPECT_EQ(f0->edge_orig[e15], NO_INDEX); + } + } + } + if (DO_OBJ) { + write_obj_mesh(out, "tritri"); + } +} + +TEST(mesh_intersect, TriTriReversed) +{ + /* Like TriTri but with triangles of opposite orientation. + * This matters because projection to 2D will now need reversed triangles. */ + const char *spec = R"(6 2 + 0 0 0 + 4 0 0 + 0 4 0 + 1 0 0 + 2 0 0 + 1 1 0 + 0 2 1 + 3 5 4 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 6); + EXPECT_EQ(out.face_size(), 6); + if (out.vert_size() == 6 && out.face_size() == 6) { + const Vert *v0 = mb.arena.find_vert(mpq3(0, 0, 0)); + const Vert *v1 = mb.arena.find_vert(mpq3(4, 0, 0)); + const Vert *v2 = mb.arena.find_vert(mpq3(0, 4, 0)); + const Vert *v3 = mb.arena.find_vert(mpq3(1, 0, 0)); + const Vert *v4 = mb.arena.find_vert(mpq3(2, 0, 0)); + const Vert *v5 = mb.arena.find_vert(mpq3(1, 1, 0)); + EXPECT_TRUE(v0 != nullptr && v1 != nullptr && v2 != nullptr); + EXPECT_TRUE(v3 != nullptr && v4 != nullptr && v5 != nullptr); + if (v0 != nullptr && v1 != nullptr && v2 != nullptr && v3 != nullptr && v4 != nullptr && + v5 != nullptr) { + EXPECT_EQ(v0->orig, 0); + EXPECT_EQ(v1->orig, 1); + const Face *f0 = find_tri_with_verts(out, v4, v5, v1); + const Face *f1 = find_tri_with_verts(out, v3, v5, v4); + const Face *f2 = find_tri_with_verts(out, v0, v5, v3); + const Face *f3 = find_tri_with_verts(out, v0, v2, v5); + const Face *f4 = find_tri_with_verts(out, v5, v2, v1); + EXPECT_TRUE(f0 != nullptr && f1 != nullptr && f2 != nullptr && f3 != nullptr && + f4 != nullptr); + /* For boolean to work right, there need to be two copies of the smaller triangle in the + * output. */ + EXPECT_EQ(count_tris_with_verts(out, v3, v5, v4), 2); + if (f0 != nullptr && f1 != nullptr && f2 != nullptr && f3 != nullptr && f4 != nullptr) { + EXPECT_EQ(f0->orig, 0); + EXPECT_TRUE(f1->orig == 0 || f1->orig == 1); + EXPECT_EQ(f2->orig, 0); + EXPECT_EQ(f3->orig, 0); + EXPECT_EQ(f4->orig, 0); + } + int e03 = find_edge_pos_in_tri(v0, v3, f2); + int e34 = find_edge_pos_in_tri(v3, v4, f1); + int e45 = find_edge_pos_in_tri(v4, v5, f1); + int e05 = find_edge_pos_in_tri(v0, v5, f3); + int e15 = find_edge_pos_in_tri(v1, v5, f0); + EXPECT_TRUE(e03 != -1 && e34 != -1 && e45 != -1 && e05 != -1 && e15 != -1); + if (e03 != -1 && e34 != -1 && e45 != -1 && e05 != -1 && e15 != -1) { + EXPECT_EQ(f2->edge_orig[e03], 2); + EXPECT_TRUE(f1->edge_orig[e34] == 2 || + f1->edge_orig[e34] == 1 * IMeshBuilder::MAX_FACE_LEN + 2); + EXPECT_EQ(f1->edge_orig[e45], 1 * IMeshBuilder::MAX_FACE_LEN + 1); + EXPECT_EQ(f3->edge_orig[e05], NO_INDEX); + EXPECT_EQ(f0->edge_orig[e15], NO_INDEX); + } + } + } + if (DO_OBJ) { + write_obj_mesh(out, "tritrirev"); + } +} + +TEST(mesh_intersect, TwoTris) +{ + Array<mpq3> verts = { + mpq3(1, 1, 1), mpq3(1, 4, 1), mpq3(1, 1, 4), /* T0 */ + mpq3(2, 2, 2), mpq3(-3, 3, 2), mpq3(-4, 1, 3), /* T1 */ + mpq3(2, 2, 2), mpq3(-3, 3, 2), mpq3(0, 3, 5), /* T2 */ + mpq3(2, 2, 2), mpq3(-3, 3, 2), mpq3(0, 3, 3), /* T3 */ + mpq3(1, 0, 0), mpq3(2, 4, 1), mpq3(-3, 2, 2), /* T4 */ + mpq3(0, 2, 1), mpq3(-2, 3, 3), mpq3(0, 1, 3), /* T5 */ + mpq3(1.5, 2, 0.5), mpq3(-2, 3, 3), mpq3(0, 1, 3), /* T6 */ + mpq3(1, 0, 0), mpq3(-2, 3, 3), mpq3(0, 1, 3), /* T7 */ + mpq3(1, 0, 0), mpq3(-3, 2, 2), mpq3(0, 1, 3), /* T8 */ + mpq3(1, 0, 0), mpq3(-1, 1, 1), mpq3(0, 1, 3), /* T9 */ + mpq3(3, -1, -1), mpq3(-1, 1, 1), mpq3(0, 1, 3), /* T10 */ + mpq3(0, 0.5, 0.5), mpq3(-1, 1, 1), mpq3(0, 1, 3), /* T11 */ + mpq3(2, 1, 1), mpq3(3, 5, 2), mpq3(-2, 3, 3), /* T12 */ + mpq3(2, 1, 1), mpq3(3, 5, 2), mpq3(-2, 3, 4), /* T13 */ + mpq3(2, 2, 5), mpq3(-3, 3, 5), mpq3(0, 3, 10), /* T14 */ + mpq3(0, 0, 0), mpq3(4, 4, 0), mpq3(-4, 2, 4), /* T15 */ + mpq3(0, 1.5, 1), mpq3(1, 2.5, 1), mpq3(-1, 2, 2), /* T16 */ + mpq3(3, 0, -2), mpq3(7, 4, -2), mpq3(-1, 2, 2), /* T17 */ + mpq3(3, 0, -2), mpq3(3, 6, 2), mpq3(-1, 2, 2), /* T18 */ + mpq3(7, 4, -2), mpq3(3, 6, 2), mpq3(-1, 2, 2), /* T19 */ + mpq3(5, 2, -2), mpq3(1, 4, 2), mpq3(-3, 0, 2), /* T20 */ + mpq3(2, 2, 0), mpq3(1, 4, 2), mpq3(-3, 0, 2), /* T21 */ + mpq3(0, 0, 0), mpq3(4, 4, 0), mpq3(-3, 0, 2), /* T22 */ + mpq3(0, 0, 0), mpq3(4, 4, 0), mpq3(-1, 2, 2), /* T23 */ + mpq3(2, 2, 0), mpq3(4, 4, 0), mpq3(0, 3, 2), /* T24 */ + mpq3(0, 0, 0), mpq3(-4, 2, 4), mpq3(4, 4, 0), /* T25 */ + }; + struct two_tri_test_spec { + int t0; + int t1; + int nv_out; + int nf_out; + }; + Array<two_tri_test_spec> test_tris = Span<two_tri_test_spec>{ + {0, 1, 8, 8}, /* 0: T1 pierces T0 inside at (1,11/6,13/6) and (1,11/5,2). */ + {0, 2, 8, 8}, /* 1: T2 intersects T0 inside (1,11/5,2) and edge (1,7/3,8/3). */ + {0, 3, 8, 7}, /* 2: T3 intersects T0 (1,11/5,2) and edge-edge (1,5/2,5/2). */ + {4, 5, 6, 4}, /* 3: T5 touches T4 inside (0,2,1). */ + {4, 6, 6, 3}, /* 4: T6 touches T4 on edge (3/2,2/1/2). */ + {4, 7, 5, 2}, /* 5: T7 touches T4 on vert (1,0,0). */ + {4, 8, 4, 2}, /* 6: T8 shared edge with T4 (1,0,0)(-3,2,2). */ + {4, 9, 5, 3}, /* 7: T9 edge (1,0,0)(-1,1,1) is subset of T4 edge. */ + {4, 10, 6, 4}, /* 8: T10 edge overlaps T4 edge with seg (-1,1,0)(1,0,0). */ + {4, 11, 6, 4}, /* 9: T11 edge (-1,1,1)(0,1/2,1/2) inside T4 edge. */ + {4, 12, 6, 2}, /* 10: parallel planes, not intersecting. */ + {4, 13, 6, 2}, /* 11: non-parallel planes, not intersecting, all one side. */ + {0, 14, 6, 2}, /* 12: non-paralel planes, not intersecting, alternate sides. */ + /* Following are all coplanar cases. */ + {15, 16, 6, 8}, /* 13: T16 inside T15. Note: dup'd tri is expected. */ + {15, 17, 8, 8}, /* 14: T17 intersects one edge of T15 at (1,1,0)(3,3,0). */ + {15, 18, 10, 12}, /* 15: T18 intersects T15 at (1,1,0)(3,3,0)(3,15/4,1/2)(0,3,2). */ + {15, 19, 8, 10}, /* 16: T19 intersects T15 at (3,3,0)(0,3,2). */ + {15, 20, 12, 14}, /* 17: T20 intersects T15 on three edges, six intersects. */ + {15, 21, 10, 11}, /* 18: T21 intersects T15 on three edges, touching one. */ + {15, 22, 5, 4}, /* 19: T22 shares edge T15, one other outside. */ + {15, 23, 4, 4}, /* 20: T23 shares edge T15, one other outside. */ + {15, 24, 5, 4}, /* 21: T24 shares two edges with T15. */ + {15, 25, 3, 2}, /* 22: T25 same T15, reverse orientation. */ + }; + static int perms[6][3] = {{0, 1, 2}, {0, 2, 1}, {1, 0, 2}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}}; + + const int do_only_test = -1; /* Make this negative to do all tests. */ + for (int test = 0; test < test_tris.size(); ++test) { + if (do_only_test >= 0 && test != do_only_test) { + continue; + } + int tri1_index = test_tris[test].t0; + int tri2_index = test_tris[test].t1; + int co1_i = 3 * tri1_index; + int co2_i = 3 * tri2_index; + + const bool verbose = false; + + if (verbose) { + std::cout << "\nTest " << test << ": T" << tri1_index << " intersect T" << tri2_index + << "\n"; + } + + const bool do_all_perms = true; + const int perm_limit = do_all_perms ? 3 : 1; + + for (int i = 0; i < perm_limit; ++i) { + for (int j = 0; j < perm_limit; ++j) { + if (do_all_perms && verbose) { + std::cout << "\nperms " << i << " " << j << "\n"; + } + IMeshArena arena; + arena.reserve(2 * 3, 2); + Array<const Vert *> f0_verts(3); + Array<const Vert *> f1_verts(3); + for (int k = 0; k < 3; ++k) { + f0_verts[k] = arena.add_or_find_vert(verts[co1_i + perms[i][k]], k); + } + for (int k = 0; k < 3; ++k) { + f1_verts[k] = arena.add_or_find_vert(verts[co2_i + perms[i][k]], k + 3); + } + Face *f0 = arena.add_face(f0_verts, 0, {0, 1, 2}); + Face *f1 = arena.add_face(f1_verts, 1, {3, 4, 5}); + IMesh in_mesh({f0, f1}); + IMesh out_mesh = trimesh_self_intersect(in_mesh, &arena); + out_mesh.populate_vert(); + EXPECT_EQ(out_mesh.vert_size(), test_tris[test].nv_out); + EXPECT_EQ(out_mesh.face_size(), test_tris[test].nf_out); + bool constexpr dump_input = true; + if (DO_OBJ && i == 0 && j == 0) { + if (dump_input) { + std::string name = "test_tt_in" + std::to_string(test); + write_obj_mesh(in_mesh, name); + } + std::string name = "test_tt" + std::to_string(test); + write_obj_mesh(out_mesh, name); + } + } + } + } +} + +TEST(mesh_intersect, OverlapCluster) +{ + /* Chain of 5 overlapping coplanar tris. + * Ordered so that clustering will make two separate clusters + * that it will have to merge into one cluster with everything. */ + const char *spec = R"(15 5 + 0 0 0 + 1 0 0 + 1/2 1 0 + 1/2 0 0 + 3/2 0 0 + 1 1 0 + 1 0 0 + 2 0 0 + 3/2 1 0 + 3/2 0 0 + 5/2 0 0 + 2 1 0 + 2 0 0 + 3 0 0 + 5/2 1 0 + 0 1 2 + 3 4 5 + 9 10 11 + 12 13 14 + 6 7 8 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 16); + EXPECT_EQ(out.face_size(), 18); + if (DO_OBJ) { + write_obj_mesh(out, "overlapcluster"); + } +} + +TEST(mesh_intersect, TriCornerCross1) +{ + /* A corner formed by 3 tris, and a 4th crossing two of them. */ + const char *spec = R"(12 4 + 0 0 0 + 1 0 0 + 0 0 1 + 0 0 0 + 0 1 0 + 0 0 1 + 0 0 0 + 1 0 0 + 0 1 0 + 1 1 1/2 + 1 -2 1/2 + -2 1 1/2 + 0 1 2 + 3 4 5 + 6 7 8 + 9 10 11 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 10); + EXPECT_EQ(out.face_size(), 14); + if (DO_OBJ) { + write_obj_mesh(out, "test_tc_1"); + } +} + +TEST(mesh_intersect, TriCornerCross2) +{ + /* A corner formed by 3 tris, and a 4th coplanar with base. */ + const char *spec = R"(12 4 + 0 0 0 + 1 0 0 + 0 0 1 + 0 0 0 + 0 1 0 + 0 0 1 + 0 0 0 + 1 0 0 + 0 1 0 + 1 1 0 + 1 -2 0 + -2 1 0 + 0 1 2 + 3 4 5 + 6 7 8 + 9 10 11 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 7); + EXPECT_EQ(out.face_size(), 8); + if (DO_OBJ) { + write_obj_mesh(out, "test_tc_2"); + } +} + +TEST(mesh_intersect, TriCornerCross3) +{ + /* A corner formed by 3 tris, and a 4th crossing all 3. */ + const char *spec = R"(12 4 + 0 0 0 + 1 0 0 + 0 0 1 + 0 0 0 + 0 1 0 + 0 0 1 + 0 0 0 + 1 0 0 + 0 1 0 + 3/2 -1/2 -1/4 + -1/2 3/2 -1/4 + -1/2 -1/2 3/4 + 0 1 2 + 3 4 5 + 6 7 8 + 9 10 11 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 10); + EXPECT_EQ(out.face_size(), 16); + if (DO_OBJ) { + write_obj_mesh(out, "test_tc_3"); + } +} + +TEST(mesh_intersect, TetTet) +{ + const char *spec = R"(8 8 + 0 0 0 + 2 0 0 + 1 2 0 + 1 1 2 + 0 0 1 + 2 0 1 + 1 2 1 + 1 1 3 + 0 1 2 + 0 3 1 + 1 3 2 + 2 3 0 + 4 5 6 + 4 7 5 + 5 7 6 + 6 7 4 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 11); + EXPECT_EQ(out.face_size(), 20); + /* Expect there to be a triangle with these three verts, oriented this way, with original face 1. + */ + const Vert *v1 = mb.arena.find_vert(mpq3(2, 0, 0)); + const Vert *v8 = mb.arena.find_vert(mpq3(0.5, 0.5, 1)); + const Vert *v9 = mb.arena.find_vert(mpq3(1.5, 0.5, 1)); + EXPECT_TRUE(v1 != nullptr && v8 != nullptr && v9 != nullptr); + const Face *f = mb.arena.find_face({v1, v8, v9}); + EXPECT_NE(f, nullptr); + EXPECT_EQ(f->orig, 1); + int v1pos = f->vert[0] == v1 ? 0 : (f->vert[1] == v1 ? 1 : 2); + EXPECT_EQ(f->edge_orig[v1pos], NO_INDEX); + EXPECT_EQ(f->edge_orig[(v1pos + 1) % 3], NO_INDEX); + EXPECT_EQ(f->edge_orig[(v1pos + 2) % 3], 1001); + EXPECT_EQ(f->is_intersect[v1pos], false); + EXPECT_EQ(f->is_intersect[(v1pos + 1) % 3], true); + EXPECT_EQ(f->is_intersect[(v1pos + 2) % 3], false); + if (DO_OBJ) { + write_obj_mesh(out, "test_tc_3"); + } +} + +TEST(mesh_intersect, CubeCubeStep) +{ + const char *spec = R"(16 24 + 0 -1 0 + 0 -1 2 + 0 1 0 + 0 1 2 + 2 -1 0 + 2 -1 2 + 2 1 0 + 2 1 2 + -1 -1 -1 + -1 -1 1 + -1 1 -1 + -1 1 1 + 1 -1 -1 + 1 -1 1 + 1 1 -1 + 1 1 1 + 0 1 3 + 0 3 2 + 2 3 7 + 2 7 6 + 6 7 5 + 6 5 4 + 4 5 1 + 4 1 0 + 2 6 4 + 2 4 0 + 7 3 1 + 7 1 5 + 8 9 11 + 8 11 10 + 10 11 15 + 10 15 14 + 14 15 13 + 14 13 12 + 12 13 9 + 12 9 8 + 10 14 12 + 10 12 8 + 15 11 9 + 15 9 13 + )"; + + IMeshBuilder mb(spec); + IMesh out = trimesh_self_intersect(mb.imesh, &mb.arena); + out.populate_vert(); + EXPECT_EQ(out.vert_size(), 22); + EXPECT_EQ(out.face_size(), 56); + if (DO_OBJ) { + write_obj_mesh(out, "test_cubecubestep"); + } + + IMeshBuilder mb2(spec); + IMesh out2 = trimesh_nary_intersect( + mb2.imesh, 2, [](int t) { return t < 12 ? 0 : 1; }, false, &mb2.arena); + out2.populate_vert(); + EXPECT_EQ(out2.vert_size(), 22); + EXPECT_EQ(out2.face_size(), 56); + if (DO_OBJ) { + write_obj_mesh(out2, "test_cubecubestep_nary"); + } +} +# endif + +# if DO_PERF_TESTS + +static void get_sphere_params( + int nrings, int nsegs, bool triangulate, int *r_num_verts, int *r_num_faces) +{ + *r_num_verts = nsegs * (nrings - 1) + 2; + if (triangulate) { + *r_num_faces = 2 * nsegs + 2 * nsegs * (nrings - 2); + } + else { + *r_num_faces = nsegs * nrings; + } +} + +static void fill_sphere_data(int nrings, + int nsegs, + const double3 ¢er, + double radius, + bool triangulate, + MutableSpan<Face *> face, + int vid_start, + int fid_start, + IMeshArena *arena) +{ + int num_verts; + int num_faces; + get_sphere_params(nrings, nsegs, triangulate, &num_verts, &num_faces); + BLI_assert(num_faces == face.size()); + Array<const Vert *> vert(num_verts); + const bool nrings_even = (nrings % 2 == 0); + int half_nrings = nrings / 2; + const bool nsegs_even = (nsegs % 2) == 0; + const bool nsegs_four_divisible = (nsegs % 4 == 0); + int half_nsegs = nrings; + int quarter_nsegs = half_nsegs / 2; + double delta_phi = 2 * M_PI / nsegs; + double delta_theta = M_PI / nrings; + int fid = fid_start; + int vid = vid_start; + auto vert_index_fn = [nrings, num_verts](int seg, int ring) { + if (ring == 0) { /* Top vert. */ + return num_verts - 2; + } + if (ring == nrings) { /* Bottom vert. */ + return num_verts - 1; + } + return seg * (nrings - 1) + (ring - 1); + }; + auto face_index_fn = [nrings](int seg, int ring) { return seg * nrings + ring; }; + auto tri_index_fn = [nrings, nsegs](int seg, int ring, int tri) { + if (ring == 0) { + return seg; + } + if (ring < nrings - 1) { + return nsegs + 2 * (ring - 1) * nsegs + 2 * seg + tri; + } + return nsegs + 2 * (nrings - 2) * nsegs + seg; + }; + Array<int> eid = {0, 0, 0, 0}; /* Don't care about edge ids. */ + /* + * (x, y , z) is given from inclination theta and azimuth phi, + * where 0 <= theta <= pi; 0 <= phi <= 2pi. + * x = radius * sin(theta) cos(phi) + * y = radius * sin(theta) sin(phi) + * z = radius * cos(theta) + */ + for (int s = 0; s < nsegs; ++s) { + double phi = s * delta_phi; + double sin_phi; + double cos_phi; + /* Avoid use of trig functions for pi/2 divisible angles. */ + if (s == 0) { + /* phi = 0. */ + sin_phi = 0.0; + cos_phi = 1.0; + } + else if (nsegs_even && s == half_nsegs) { + /* phi = pi. */ + sin_phi = 0.0; + cos_phi = -1.0; + } + else if (nsegs_four_divisible && s == quarter_nsegs) { + /* phi = pi/2. */ + sin_phi = 1.0; + cos_phi = 0.0; + } + else if (nsegs_four_divisible && s == 3 * quarter_nsegs) { + /* phi = 3pi/2. */ + sin_phi = -1.0; + cos_phi = 0.0; + } + else { + sin_phi = sin(phi); + cos_phi = cos(phi); + } + for (int r = 1; r < nrings; ++r) { + double theta = r * delta_theta; + double r_sin_theta; + double r_cos_theta; + if (nrings_even && r == half_nrings) { + /* theta = pi/2. */ + r_sin_theta = radius; + r_cos_theta = 0.0; + } + else { + r_sin_theta = radius * sin(theta); + r_cos_theta = radius * cos(theta); + } + double x = r_sin_theta * cos_phi + center[0]; + double y = r_sin_theta * sin_phi + center[1]; + double z = r_cos_theta + center[2]; + const Vert *v = arena->add_or_find_vert(mpq3(x, y, z), vid++); + vert[vert_index_fn(s, r)] = v; + } + } + const Vert *vtop = arena->add_or_find_vert(mpq3(center[0], center[1], center[2] + radius), + vid++); + const Vert *vbot = arena->add_or_find_vert(mpq3(center[0], center[1], center[2] - radius), + vid++); + vert[vert_index_fn(0, 0)] = vtop; + vert[vert_index_fn(0, nrings)] = vbot; + for (int s = 0; s < nsegs; ++s) { + int snext = (s + 1) % nsegs; + for (int r = 0; r < nrings; ++r) { + int rnext = r + 1; + int i0 = vert_index_fn(s, r); + int i1 = vert_index_fn(s, rnext); + int i2 = vert_index_fn(snext, rnext); + int i3 = vert_index_fn(snext, r); + Face *f; + Face *f2 = nullptr; + if (r == 0) { + f = arena->add_face({vert[i0], vert[i1], vert[i2]}, fid++, eid); + } + else if (r == nrings - 1) { + f = arena->add_face({vert[i0], vert[i1], vert[i3]}, fid++, eid); + } + else { + if (triangulate) { + f = arena->add_face({vert[i0], vert[i1], vert[i2]}, fid++, eid); + f2 = arena->add_face({vert[i2], vert[i3], vert[i0]}, fid++, eid); + } + else { + f = arena->add_face({vert[i0], vert[i1], vert[i2], vert[i3]}, fid++, eid); + } + } + if (triangulate) { + int f_index = tri_index_fn(s, r, 0); + face[f_index] = f; + if (r != 0 && r != nrings - 1) { + int f_index2 = tri_index_fn(s, r, 1); + face[f_index2] = f2; + } + } + else { + int f_index = face_index_fn(s, r); + face[f_index] = f; + } + } + } +} + +static void spheresphere_test(int nrings, double y_offset, bool use_self) +{ + /* Make two uvspheres with nrings rings ad 2*nrings segments. */ + if (nrings < 2) { + return; + } + BLI_task_scheduler_init(); /* Without this, no parallelism. */ + double time_start = PIL_check_seconds_timer(); + IMeshArena arena; + int nsegs = 2 * nrings; + int num_sphere_verts; + int num_sphere_tris; + get_sphere_params(nrings, nsegs, true, &num_sphere_verts, &num_sphere_tris); + Array<Face *> tris(2 * num_sphere_tris); + arena.reserve(2 * num_sphere_verts, 2 * num_sphere_tris); + double3 center1(0.0, 0.0, 0.0); + fill_sphere_data(nrings, + nsegs, + center1, + 1.0, + true, + MutableSpan<Face *>(tris.begin(), num_sphere_tris), + 0, + 0, + &arena); + double3 center2(0.0, y_offset, 0.0); + fill_sphere_data(nrings, + nsegs, + center2, + 1.0, + true, + MutableSpan<Face *>(tris.begin() + num_sphere_tris, num_sphere_tris), + num_sphere_verts, + num_sphere_verts, + &arena); + IMesh mesh(tris); + double time_create = PIL_check_seconds_timer(); + // write_obj_mesh(mesh, "spheresphere_in"); + IMesh out; + if (use_self) { + out = trimesh_self_intersect(mesh, &arena); + } + else { + int nf = num_sphere_tris; + out = trimesh_nary_intersect( + mesh, 2, [nf](int t) { return t < nf ? 0 : 1; }, false, &arena); + } + double time_intersect = PIL_check_seconds_timer(); + std::cout << "Create time: " << time_create - time_start << "\n"; + std::cout << "Intersect time: " << time_intersect - time_create << "\n"; + std::cout << "Total time: " << time_intersect - time_start << "\n"; + if (DO_OBJ) { + write_obj_mesh(out, "spheresphere"); + } + BLI_task_scheduler_exit(); +} + +static void get_grid_params( + int x_subdiv, int y_subdiv, bool triangulate, int *r_num_verts, int *r_num_faces) +{ + *r_num_verts = x_subdiv * y_subdiv; + if (triangulate) { + *r_num_faces = 2 * (x_subdiv - 1) * (y_subdiv - 1); + } + else { + *r_num_faces = (x_subdiv - 1) * (y_subdiv - 1); + } +} + +static void fill_grid_data(int x_subdiv, + int y_subdiv, + bool triangulate, + double size, + const double3 ¢er, + MutableSpan<Face *> face, + int vid_start, + int fid_start, + IMeshArena *arena) +{ + if (x_subdiv <= 1 || y_subdiv <= 1) { + return; + } + int num_verts; + int num_faces; + get_grid_params(x_subdiv, y_subdiv, triangulate, &num_verts, &num_faces); + BLI_assert(face.size() == num_faces); + Array<const Vert *> vert(num_verts); + auto vert_index_fn = [x_subdiv](int ix, int iy) { return iy * x_subdiv + ix; }; + auto face_index_fn = [x_subdiv](int ix, int iy) { return iy * (x_subdiv - 1) + ix; }; + auto tri_index_fn = [x_subdiv](int ix, int iy, int tri) { + return 2 * iy * (x_subdiv - 1) + 2 * ix + tri; + }; + Array<int> eid = {0, 0, 0, 0}; /* Don't care about edge ids. */ + double r = size / 2.0; + double delta_x = size / (x_subdiv - 1); + double delta_y = size / (y_subdiv - 1); + int vid = vid_start; + for (int iy = 0; iy < y_subdiv; ++iy) { + for (int ix = 0; ix < x_subdiv; ++ix) { + double x = center[0] - r + ix * delta_x; + double y = center[1] - r + iy * delta_y; + double z = center[2]; + const Vert *v = arena->add_or_find_vert(mpq3(x, y, z), vid++); + vert[vert_index_fn(ix, iy)] = v; + } + } + int fid = fid_start; + for (int iy = 0; iy < y_subdiv - 1; ++iy) { + for (int ix = 0; ix < x_subdiv - 1; ++ix) { + int i0 = vert_index_fn(ix, iy); + int i1 = vert_index_fn(ix, iy + 1); + int i2 = vert_index_fn(ix + 1, iy + 1); + int i3 = vert_index_fn(ix + 1, iy); + if (triangulate) { + Face *f = arena->add_face({vert[i0], vert[i1], vert[i2]}, fid++, eid); + Face *f2 = arena->add_face({vert[i2], vert[i3], vert[i0]}, fid++, eid); + face[tri_index_fn(ix, iy, 0)] = f; + face[tri_index_fn(ix, iy, 1)] = f2; + } + else { + Face *f = arena->add_face({vert[i0], vert[i1], vert[i2], vert[i3]}, fid++, eid); + face[face_index_fn(ix, iy)] = f; + } + } + } +} + +static void spheregrid_test(int nrings, int grid_level, double z_offset, bool use_self) +{ + /* Make a uvsphere and a grid. + * The sphere is radius 1, has nrings rings and 2 * nrings segs, + * and is centered at (0,0,z_offset). + * The plane is 4x4, has 2**grid_level subdivisions x and y, + * and is centered at the origin. */ + if (nrings < 2 || grid_level < 1) { + return; + } + BLI_task_scheduler_init(); /* Without this, no parallelism. */ + double time_start = PIL_check_seconds_timer(); + IMeshArena arena; + int num_sphere_verts; + int num_sphere_tris; + int nsegs = 2 * nrings; + int num_grid_verts; + int num_grid_tris; + int subdivs = 1 << grid_level; + get_sphere_params(nrings, nsegs, true, &num_sphere_verts, &num_sphere_tris); + get_grid_params(subdivs, subdivs, true, &num_grid_verts, &num_grid_tris); + Array<Face *> tris(num_sphere_tris + num_grid_tris); + arena.reserve(num_sphere_verts + num_grid_verts, num_sphere_tris + num_grid_tris); + double3 center(0.0, 0.0, z_offset); + fill_sphere_data(nrings, + nsegs, + center, + 1.0, + true, + MutableSpan<Face *>(tris.begin(), num_sphere_tris), + 0, + 0, + &arena); + fill_grid_data(subdivs, + subdivs, + true, + 4.0, + double3(0, 0, 0), + MutableSpan<Face *>(tris.begin() + num_sphere_tris, num_grid_tris), + num_sphere_verts, + num_sphere_tris, + &arena); + IMesh mesh(tris); + double time_create = PIL_check_seconds_timer(); + // write_obj_mesh(mesh, "spheregrid_in"); + IMesh out; + if (use_self) { + out = trimesh_self_intersect(mesh, &arena); + } + else { + int nf = num_sphere_tris; + out = trimesh_nary_intersect( + mesh, 2, [nf](int t) { return t < nf ? 0 : 1; }, false, &arena); + } + double time_intersect = PIL_check_seconds_timer(); + std::cout << "Create time: " << time_create - time_start << "\n"; + std::cout << "Intersect time: " << time_intersect - time_create << "\n"; + std::cout << "Total time: " << time_intersect - time_start << "\n"; + if (DO_OBJ) { + write_obj_mesh(out, "spheregrid"); + } + BLI_task_scheduler_exit(); +} + +TEST(mesh_intersect_perf, SphereSphere) +{ + spheresphere_test(512, 0.5, false); +} + +TEST(mesh_intersect_perf, SphereGrid) +{ + spheregrid_test(512, 4, 0.1, false); +} + +# endif + +} // namespace blender::meshintersect::tests +#endif |