From f81382f604d58d1ac9b73584622996de15248360 Mon Sep 17 00:00:00 2001 From: Lukas Matena Date: Thu, 19 May 2022 16:59:05 +0200 Subject: [PATCH] Generic sphere shape is now created by recursive division of an icosahedron --- src/libslic3r/TriangleMesh.cpp | 151 ++++++++++++++++++-------- tests/fff_print/test_trianglemesh.cpp | 14 ++- 2 files changed, 117 insertions(+), 48 deletions(-) diff --git a/src/libslic3r/TriangleMesh.cpp b/src/libslic3r/TriangleMesh.cpp index d5fe83c26..82dc254ca 100644 --- a/src/libslic3r/TriangleMesh.cpp +++ b/src/libslic3r/TriangleMesh.cpp @@ -1013,58 +1013,121 @@ indexed_triangle_set its_make_pyramid(float base, float height) // Generates mesh for a sphere centered about the origin, using the generated angle // to determine the granularity. // Default angle is 1 degree. -//FIXME better to discretize an Icosahedron recursively http://www.songho.ca/opengl/gl_sphere.html indexed_triangle_set its_make_sphere(double radius, double fa) { - int sectorCount = int(ceil(2. * M_PI / fa)); - int stackCount = int(ceil(M_PI / fa)); - float sectorStep = float(2. * M_PI / sectorCount); - float stackStep = float(M_PI / stackCount); - + // First build an icosahedron (taken from http://www.songho.ca/opengl/gl_sphere.html) indexed_triangle_set mesh; + + const float PI = 3.1415926f; + const float H_ANGLE = PI / 180 * 72; // 72 degree = 360 / 5 + const float V_ANGLE = atanf(1.0f / 2); // elevation = 26.565 degree + auto& vertices = mesh.vertices; - vertices.reserve((stackCount - 1) * sectorCount + 2); - for (int i = 0; i <= stackCount; ++ i) { - // from pi/2 to -pi/2 - double stackAngle = 0.5 * M_PI - stackStep * i; - double xy = radius * cos(stackAngle); - double z = radius * sin(stackAngle); - if (i == 0 || i == stackCount) - vertices.emplace_back(Vec3f(float(xy), 0.f, float(z))); - else - for (int j = 0; j < sectorCount; ++ j) { - // from 0 to 2pi - double sectorAngle = sectorStep * j; - vertices.emplace_back(Vec3d(xy * std::cos(sectorAngle), xy * std::sin(sectorAngle), z).cast()); - } - } + auto& indices = mesh.indices; + vertices.resize(12); + indices.reserve(20); - auto& facets = mesh.indices; - facets.reserve(2 * (stackCount - 1) * sectorCount); - for (int i = 0; i < stackCount; ++ i) { - // Beginning of current stack. - int k1 = (i == 0) ? 0 : (1 + (i - 1) * sectorCount); - int k1_first = k1; - // Beginning of next stack. - int k2 = (i == 0) ? 1 : (k1 + sectorCount); - int k2_first = k2; - for (int j = 0; j < sectorCount; ++ j) { - // 2 triangles per sector excluding first and last stacks - int k1_next = k1; - int k2_next = k2; - if (i != 0) { - k1_next = (j + 1 == sectorCount) ? k1_first : (k1 + 1); - facets.emplace_back(k1, k2, k1_next); + float z, xy; + float hAngle1 = -PI / 2 - H_ANGLE / 2; + + vertices[0] = stl_vertex(0, 0, radius); // the first top vertex at (0, 0, r) + + for (int i = 1; i <= 5; ++i) { + z = radius * sinf(V_ANGLE); + xy = radius * cosf(V_ANGLE); + vertices[i] = stl_vertex(xy * cosf(hAngle1), xy * sinf(hAngle1), z); + vertices[i+5] = stl_vertex(xy * cosf(hAngle1 + H_ANGLE / 2), xy * sinf(hAngle1 + H_ANGLE / 2), -z); + hAngle1 += H_ANGLE; + + indices.emplace_back(stl_triangle_vertex_indices(i, i < 5 ? i+1 : 1, 0)); + indices.emplace_back(stl_triangle_vertex_indices(i, i+5, i < 5 ? i+1 : 1)); + indices.emplace_back(stl_triangle_vertex_indices(i+5, i+6 < 11 ? i+6 : 6, i+6 < 11 ? i+1 : 1)); + indices.emplace_back(stl_triangle_vertex_indices(i+5, 11, i+6 < 11 ? i+6 : 6)); + } + vertices[11] = stl_vertex(0, 0, -radius); // the last bottom vertex at (0, 0, -r) + + + // We have a beautiful icosahedron. Now subdivide the triangles. + std::vector neighbors = its_face_neighbors(mesh); // This is cheap, the mesh is small. + + const double side_len_limit = radius * fa; + const double side_len = (vertices[1] - vertices[0]).norm(); + const int iterations = std::ceil(std::log2(side_len / side_len_limit)); + + indices.reserve(indices.size() * std::pow(4, iterations)); + vertices.reserve(vertices.size() * std::pow(2, iterations)); + + struct DividedEdge { + int neighbor = -1; + int middle_vertex_idx; + std::pair children_idxs; + }; + + for (int iter=0; iter> divided_triangles(indices.size()); + std::vector new_neighbors(4*indices.size()); + + size_t orig_indices_size = indices.size(); + for (int i=0; i, 3> edge_children = { std::make_pair(i,last_triangle_idx + 2), + std::make_pair(last_triangle_idx + 2,last_triangle_idx + 3), + std::make_pair(last_triangle_idx + 3,i) }; + + std::array middle_vertices_idxs; + std::array, 3> new_neighbors_per_edge; + + for (int n=0; n<3; ++n) { // for all three edges + const int edge_neighbor = neighbors[i][n]; + + if (divided_triangles[edge_neighbor][0].neighbor == -1) { + // This n-th edge is not yet divided. Divide it now. + vertices.emplace_back(0.5 * (vertices[indices[i][n]] + vertices[indices[i][n == 2 ? 0 : n+1]])); + vertices.back() *= radius / vertices.back().norm(); + middle_vertices_idxs[n] = vertices.size()-1; + + // Save information about what we did. + int j = -1; + while (divided_triangles[i][++j].neighbor != -1); + + divided_triangles[i][j] = { edge_neighbor, int(vertices.size()-1), edge_children[n] }; + new_neighbors_per_edge[n] = std::make_pair(-1,-1); + } else { + // This edge is already divided. Get the index of the middle point. + int j = -1; + while (divided_triangles[edge_neighbor][++j].neighbor != i); + middle_vertices_idxs[n] = divided_triangles[edge_neighbor][j].middle_vertex_idx; + new_neighbors_per_edge[n] = divided_triangles[edge_neighbor][j].children_idxs; + std::swap(new_neighbors_per_edge[n].first, new_neighbors_per_edge[n].second); + + // We have saved the middle-point. We are looking for edges leading to/from it. + int idx = -1; while (indices[new_neighbors_per_edge[n].first][++idx] != middle_vertices_idxs[n]); + new_neighbors[new_neighbors_per_edge[n].first][idx] = edge_children[n].first; + new_neighbors[new_neighbors_per_edge[n].second][idx] = edge_children[n].second; + } } - if (i + 1 != stackCount) { - k2_next = (j + 1 == sectorCount) ? k2_first : (k2 + 1); - facets.emplace_back(k1_next, k2, k2_next); - } - k1 = k1_next; - k2 = k2_next; + + // Add three new triangles, reindex the old one. + const int last_index = indices.size() - 1; + indices.emplace_back(stl_triangle_vertex_indices(middle_vertices_idxs[0], middle_vertices_idxs[1], middle_vertices_idxs[2])); + new_neighbors[indices.size()-1] = Vec3i(last_index+2, last_index+3, i); + + indices.emplace_back(stl_triangle_vertex_indices(middle_vertices_idxs[0], indices[i][1], middle_vertices_idxs[1])); + new_neighbors[indices.size()-1] = Vec3i(new_neighbors_per_edge[0].second, new_neighbors_per_edge[1].first, last_index+1); + + indices.emplace_back(stl_triangle_vertex_indices(middle_vertices_idxs[2], middle_vertices_idxs[1], indices[i][2])); + new_neighbors[indices.size()-1] = Vec3i(last_index+1, new_neighbors_per_edge[1].second, new_neighbors_per_edge[2].first); + + indices[i][1] = middle_vertices_idxs[0]; + indices[i][2] = middle_vertices_idxs[2]; + new_neighbors[i] = Vec3i(new_neighbors_per_edge[0].first, last_index+1, new_neighbors_per_edge[2].second); + } + neighbors = std::move(new_neighbors); } - return mesh; } diff --git a/tests/fff_print/test_trianglemesh.cpp b/tests/fff_print/test_trianglemesh.cpp index 6faaf1584..eff39ed35 100644 --- a/tests/fff_print/test_trianglemesh.cpp +++ b/tests/fff_print/test_trianglemesh.cpp @@ -220,10 +220,16 @@ SCENARIO( "make_xxx functions produce meshes.") { GIVEN("make_sphere() function") { WHEN("make_sphere() is called with arguments 10, PI / 3") { TriangleMesh sph = make_sphere(10, PI / 243.0); - THEN("Resulting mesh has one point at 0,0,-10 and one at 0,0,10") { - const std::vector &verts = sph.its.vertices; - REQUIRE(std::count_if(verts.begin(), verts.end(), [](const Vec3f& t) { return is_approx(t, Vec3f(0.f, 0.f, 10.f)); } ) == 1); - REQUIRE(std::count_if(verts.begin(), verts.end(), [](const Vec3f& t) { return is_approx(t, Vec3f(0.f, 0.f, -10.f)); } ) == 1); + THEN( "Edge length is smaller than limit but not smaller than half of it") { + double len = (sph.its.vertices[sph.its.indices[0][0]] - sph.its.vertices[sph.its.indices[0][1]]).norm(); + double limit = 10*PI/243.; + REQUIRE(len <= limit); + REQUIRE(len >= limit/2.); + } + THEN( "Vertices are about the correct distance from the origin") { + bool all_vertices_ok = std::all_of(sph.its.vertices.begin(), sph.its.vertices.end(), + [](const stl_vertex& pt) { return is_approx(pt.squaredNorm(), 100.f); }); + REQUIRE(all_vertices_ok); } THEN( "The mesh volume is approximately 4/3 * pi * 10^3") { REQUIRE(abs(sph.volume() - (4.0/3.0 * M_PI * std::pow(10,3))) < 1); // 1% tolerance?