cgal
cgal copied to clipboard
CGAL
How to get outlines of a Polyhedron_3 so I can render the edges as 3D polylines instead of shading edges
Hello,
I'm trying to figure out a way to create a vector of 3D poly-lines for each defining feature of a polyhedron. I'm creating an application that will serve as a 3D conversational tubing notch generator for both laser and plasma cutting applications. User will define the host tube (tube to be cut) then the tubes on each end of the cut along with the orientations. Tubes can be round, square, or rectangular. I'll use boolean subtractions to generate the face notches into the host pipe. All that functionality works just fine, I just can't find a way to get the resulting 3D polyline from each side of the host pipe (before or after) the subtraction is done. I would expect 4 3D polylines that lie on the outside edges and inside edges of the pipe.
Issue Details
I need to do this for two reasons 1 . Id like to render the lines of these features around the polyhedron's so that the edges are visibly defined (Most CAD programs do this) without shading. 2. I need the 3D poly-lines in order to create toolpaths after a boolean subtraction occurs from the host pipe to the cut pipe. The host pipes and cut pipes can be round, square, and rectangular.
Source Code
My current attempt, although I've spent days trying everything I could possible think of to do this, is as follows.
typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef CGAL::Polyhedron_3<K> Polyhedron;
Polyhedron cube_poly;
this->make_cube(cube_poly, 100.0f); //Produces a cube at the 0,0,0 with size of 100.0f units
Mesh mesh1, mesh2;
CGAL::copy_face_graph(this->body_poly, mesh1); //This is a pipe at 0,0,0 with a diameter of 5, a length of 5, and a wall thickness of 1
CGAL::copy_face_graph(cube_poly, mesh2);
// Compute the intersection of the two meshes
Mesh intersection;
if(CGAL::Polygon_mesh_processing::corefine_and_compute_intersection(mesh1, mesh2, intersection)) {
LOG_INFO("Intersection computed successfully.");
} else {
LOG_ERROR("Intersection could not be computed.");
}
LOG_INFO("intersection.size_of_vertices(): {}", intersection.number_of_vertices());
// Extract the boundary edges of the intersection
std::vector<std::vector<K::Point_3>> contours;
for(auto edge : edges(intersection))
{
std::vector<K::Point_3> contour;
auto h = halfedge(edge, intersection);
contour.push_back(intersection.point(source(h, intersection)));
contour.push_back(intersection.point(target(h, intersection)));
contours.push_back(contour);
}
My idea here was to create a huge cube that completely covers the pipe then do an intersection in order to retrieve the edges of all the intersections. The trouble is, I can't figure out a way to only get the edges of the intersection, it currently creates "contours" that have two vertices that represent each intersection so the result is just the same as rendering the wire frame of the triangular mesh. I'm hoping someone here can help me come up with a solution or point me in a direction I have not tried.
My meshes will aways be water tight/manifold and I originally thought I might be able to get the connected components (which I thought was a way to retrieve each surface (multiple faces with the same normal) but connected components is always 1 for my meshes and for a pipe I thought it may be 4 (front surface, rear surface, and outside surface, and inside surface)
I'm new to CGAL and hoping I'm missing something simple to provide this kind of functionality because it seems like most CAD programs do this with ease.
Environment
- Operating system (Mac M10):
- Compiler: Clang 14.0.03 --std=c++17
- CGAL version: 5.6
- Boost version: 1.83.0
- Other libraries versions if used (GLM, OpenGL, wxWidgets):
I'm not sure to understand correctly what you want to get. Are you interested in getting in the output the edges that are common to both mesh1 and mesh2? In such a case, you should use the named parameter edge_is_constrained_map() in the output named parameters and filter output edges using that map. Such a map is created in this example.
I added a screenshot of what i'm talking about, see the black lines rendered around the polyhedron. I was able to figure out a way to do it (Not by doing the way I was trying to do it when I asked the question originally) but I'm curious if there's a better way. Right now I'm iterating the half-edges and determining which faces are more than 45 degrees from another face, then the edge gets pushed back to a vector of pairs. Then the duplicates get removed. Then I chain them. It works okay on simple geometry but it does break on my resulting tube notch. So in the snapshot, I end up with 678 chains instead of 4 like I need. Chaining the edges (maybe I should use the terminology "sharp edges") into complete closed contours is important so I can create toolpaths of those edges translated into Gcode.
My method is like as follows
class MeshVertex
{
public:
/*
Define epsilon for floating point comparisons downstream
*/
static constexpr double EPSILON = 0.0001; // Add this line
MeshVertex()
{
this->x = 0;
this->y = 0;
this->z = 0;
};
MeshVertex(double xi, double yi, double zi)
{
this->x = xi;
this->y = yi;
this->z = zi;
};
double x;
double y;
double z;
bool operator==(const MeshVertex& other) const {
return std::abs(x - other.x) < EPSILON
&& std::abs(y - other.y) < EPSILON
&& std::abs(z - other.z) < EPSILON;
};
bool operator!=(const MeshVertex& other) const {
return !(*this == other);
};
MeshVertex& operator=(const MeshVertex& other) {
if (this != &other) {
x = other.x;
y = other.y;
z = other.z;
}
return *this;
};
/*
Mathmatic function on MeshVertex
*/
void offset(MeshVertex offset)
{
this->x += offset.x;
this->y += offset.y;
this->z += offset.z;
};
void scale(double scale)
{
this->x *= scale;
this->y *= scale;
this->z *= scale;
};
void rotate(MeshVertex axis, MeshVertex angle)
{
// Create a quaternion for each axis of rotation
glm::quat rotX = glm::angleAxis(static_cast<float>(angle.x), glm::vec3(1.0f, 0.0f, 0.0f));
glm::quat rotY = glm::angleAxis(static_cast<float>(angle.y), glm::vec3(0.0f, 1.0f, 0.0f));
glm::quat rotZ = glm::angleAxis(static_cast<float>(angle.z), glm::vec3(0.0f, 0.0f, 1.0f));
// Combine the quaternions
glm::quat rot = rotX * rotY * rotZ;
// Apply the rotation to the point
glm::vec3 point = glm::vec3(this->x, this->y, this->z);
point = rot * point;
// Update the point coordinates
this->x = point.x;
this->y = point.y;
this->z = point.z;
};
MeshVertex to_radians()
{
return MeshVertex(this->x * M_PI / 180.0f, this->y * M_PI / 180.0f, this->z * M_PI / 180.0f);
};
};
std::vector<std::vector<MeshySpace::MeshVertex>> MeshySpace::Mesh::GetEdges()
{
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
typedef Polyhedron::Halfedge_handle Halfedge_handle;
Polyhedron polyhedron = this->body_poly;
CGAL::Polygon_mesh_processing::triangulate_faces(polyhedron);
std::map<Polyhedron::Face_handle, Kernel::Vector_3> face_normals_map;
auto face_normals = boost::make_assoc_property_map(face_normals_map);
CGAL::Polygon_mesh_processing::compute_face_normals(polyhedron, face_normals);
double threshold = 45;
std::vector<std::pair<MeshySpace::MeshVertex, MeshySpace::MeshVertex>> sharp_edges;
for(auto h = polyhedron.halfedges_begin(); h != polyhedron.halfedges_end(); ++h)
{
auto normal1 = face_normals[h->face()];
auto normal2 = face_normals[h->opposite()->face()];
auto dot_product = normal1 * normal2;
auto cos_angle = CGAL::to_double(dot_product);
auto angle_in_radians = std::acos(cos_angle);
auto angle_in_degrees = angle_in_radians * 180.0 / CGAL_PI;
if(angle_in_degrees > threshold) {
auto p1 = h->vertex()->point();
auto p2 = h->opposite()->vertex()->point();
MeshySpace::MeshVertex p1_vertex(CGAL::to_double(p1.x()), CGAL::to_double(p1.y()), CGAL::to_double(p1.z()));
MeshySpace::MeshVertex p2_vertex(CGAL::to_double(p2.x()), CGAL::to_double(p2.y()), CGAL::to_double(p2.z()));
sharp_edges.push_back(std::make_pair(p1_vertex, p2_vertex));
}
}
//LOG_INFO("sharp_edges.size = {}", sharp_edges.size());
// Merge duplicate or near-duplicate edges
std::vector<std::pair<MeshySpace::MeshVertex, MeshySpace::MeshVertex>> unique_edges;
for (const auto& edge : sharp_edges) {
bool is_duplicate = false;
for (const auto& unique_edge : unique_edges) {
if ((edge.first == unique_edge.first && edge.second == unique_edge.second) ||
(edge.first == unique_edge.second && edge.second == unique_edge.first)) {
is_duplicate = true;
break;
}
}
if (!is_duplicate) {
unique_edges.push_back(edge);
}
}
//LOG_INFO("unique_edges.size = {}", unique_edges.size());
sharp_edges = unique_edges;
std::vector<std::vector<MeshySpace::MeshVertex>> contours;
//Form contours from the sharp edges
while (!sharp_edges.empty()) {
std::vector<MeshySpace::MeshVertex> contour;
auto edge = sharp_edges.back();
sharp_edges.pop_back();
contour.push_back(edge.first);
contour.push_back(edge.second);
bool found;
do {
found = false;
for (auto it = sharp_edges.begin(); it != sharp_edges.end(); ++it) {
if (it->first == edge.second) {
edge = *it;
sharp_edges.erase(it);
contour.push_back(edge.second);
found = true;
break;
} else if (it->second == edge.second) {
edge = std::make_pair(it->second, it->first);
sharp_edges.erase(it);
contour.push_back(edge.second);
found = true;
break;
}
}
} while (found);
contours.push_back(contour);
}
//LOG_INFO("contours.size = {}", contours.size());
return contours;
}