Skip to content

Commit 6c1da13

Browse files
thewtexthewtex
committed
feat(mesh-filters): add slice-mesh pipeline
1 parent f5a2063 commit 6c1da13

File tree

2 files changed

+141
-94
lines changed

2 files changed

+141
-94
lines changed

‎packages/mesh-filters/CMakeLists.txt

Lines changed: 39 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,16 +66,31 @@ FetchContent_Declare(
6666
GIT_SHALLOW TRUE
6767
)
6868

69-
FetchContent_MakeAvailable(pmp geogram)
69+
set(MESH_PLANE_INTERSECTION_GIT_REPOSITORY "https://github.com/intents-software/mesh-plane-intersection.git")
70+
# 2024-12-14 master
71+
set(MESH_PLANE_INTERSECTION_GIT_TAG "b632819a0027e8c3a6e64036c4984fe0f5c2c0d4")
72+
FetchContent_Declare(
73+
mesh_plane_intersection
74+
GIT_REPOSITORY ${MESH_PLANE_INTERSECTION_GIT_REPOSITORY}
75+
GIT_TAG ${MESH_PLANE_INTERSECTION_GIT_TAG}
76+
GIT_SHALLOW TRUE
77+
)
78+
79+
FetchContent_MakeAvailable(pmp geogram mesh_plane_intersection)
7080
include_directories(${pmp_SOURCE_DIR}/src)
7181
include_directories(${geogram_SOURCE_DIR}/src/lib)
82+
include_directories(${mesh_plane_intersection_SOURCE_DIR}/src)
7283

7384
foreach(pipeline geogram-conversion repair smooth-remesh)
7485
add_executable(${pipeline} ${pipeline}.cxx)
7586
target_link_libraries(${pipeline} PUBLIC ${ITK_LIBRARIES} geogram)
7687
target_include_directories(${pipeline} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
7788
endforeach()
7889

90+
add_executable(slice-mesh slice-mesh.cxx)
91+
target_link_libraries(slice-mesh PUBLIC ${ITK_LIBRARIES})
92+
target_include_directories(slice-mesh PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
93+
7994
enable_testing()
8095

8196
add_test(NAME geogram-conversion
@@ -124,6 +139,29 @@ add_test(NAME smooth-remesh-brain
124139
${CMAKE_CURRENT_BINARY_DIR}/brain-smooth-remesh.off
125140
--number-points 30
126141
)
142+
143+
add_test(NAME slice-mesh-help
144+
COMMAND slice-mesh
145+
--help
146+
)
147+
add_test(NAME slice-mesh-suzanne
148+
COMMAND slice-mesh
149+
${CMAKE_CURRENT_SOURCE_DIR}/test/data/baseline/suzanne-repair.off
150+
${CMAKE_CURRENT_SOURCE_DIR}/test/data/input/slice-mesh-plane.json
151+
${CMAKE_CURRENT_BINARY_DIR}/slice-mesh-suzanne.vtk
152+
)
153+
add_test(NAME slice-mesh-cow
154+
COMMAND slice-mesh
155+
${CMAKE_CURRENT_SOURCE_DIR}/test/data/baseline/cow-repair.off
156+
${CMAKE_CURRENT_SOURCE_DIR}/test/data/input/slice-mesh-plane.json
157+
${CMAKE_CURRENT_BINARY_DIR}/slice-mesh-cow.vtk
158+
)
159+
add_test(NAME slice-mesh-brain
160+
COMMAND slice-mesh
161+
${CMAKE_CURRENT_SOURCE_DIR}/test/data/baseline/brain-repair.off
162+
${CMAKE_CURRENT_SOURCE_DIR}/test/data/input/slice-mesh-plane.json
163+
${CMAKE_CURRENT_BINARY_DIR}/slice-mesh-brain.vtk
164+
)
127165
# # Interesting backtrace on exit
128166
# if (NOT CMAKE_BUILD_TYPE STREQUAL "Debug")
129167
# add_test(NAME mesh-filters

‎packages/mesh-filters/slice-mesh.cxx

Lines changed: 102 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -20,127 +20,136 @@
2020
#include "itkInputMesh.h"
2121
#include "itkOutputMesh.h"
2222
#include "itkSupportInputMeshTypes.h"
23+
#include "itkInputTextStream.h"
24+
#include "itkPolyLineCell.h"
2325

24-
#include "itkmeshToGeogramMesh.h"
25-
#include "itkgeogramMeshToMesh.h"
26+
#include "MeshPlaneIntersect.hpp"
2627

27-
#include <geogram/mesh/mesh_remesh.h>
28-
#include <geogram/mesh/mesh_geometry.h>
29-
#include <geogram/basic/command_line.h>
30-
#include <geogram/basic/command_line_args.h>
31-
32-
#include <geogram/basic/common.h>
33-
#include <geogram/basic/process.h>
34-
#include <geogram/basic/logger.h>
35-
#include <geogram/basic/command_line.h>
28+
#include "glaze/glaze.hpp"
3629

3730
template <typename TMesh>
38-
int repairMesh(itk::wasm::Pipeline &pipeline, const TMesh *inputMesh)
31+
int sliceMesh(itk::wasm::Pipeline &pipeline, const TMesh *inputMesh)
3932
{
4033
using MeshType = TMesh;
4134
constexpr unsigned int Dimension = MeshType::PointDimension;
4235
using PixelType = typename MeshType::PixelType;
36+
using PointIdentifierType = typename MeshType::PointIdentifier;
37+
using CellIdentifierType = typename MeshType::CellIdentifier;
38+
39+
using FloatType = typename MeshType::CoordRepType;
40+
using IndexType = typename MeshType::CellIdentifier;
41+
using MeshPlaneIntersectType = MeshPlaneIntersect<FloatType, IndexType>;
42+
using Vec3D = typename MeshPlaneIntersectType::Vec3D;
43+
using Face = typename MeshPlaneIntersectType::Face;
44+
using Plane = typename MeshPlaneIntersectType::Plane;
45+
using Planes = std::vector<Plane>;
46+
47+
using OutputMeshType = itk::Mesh<uint16_t, Dimension>;
48+
using CellAutoPointer = typename OutputMeshType::CellAutoPointer;
49+
using PolyLineCellType = itk::PolyLineCell<typename OutputMeshType::CellType>;
50+
using PointDataContainer = typename OutputMeshType::PointDataContainer;
51+
using CellDataContainer = typename OutputMeshType::CellDataContainer;
4352

4453
pipeline.get_option("input-mesh")->required()->type_name("INPUT_MESH");
4554

46-
double numberPointsPercent = 75.0;
47-
pipeline.add_option("--number-points", numberPointsPercent, "Number of points as a percent of the bounding box diagonal. Output may have slightly more points.");
48-
49-
double triangleShapeAdaptation = 1.0;
50-
pipeline.add_option("--triangle-shape-adaptation", triangleShapeAdaptation, "Triangle shape adaptation factor. Use 0.0 to disable.");
51-
52-
double triangleSizeAdaptation = 0.0;
53-
pipeline.add_option("--triangle-size-adaptation", triangleSizeAdaptation, "Triangle size adaptation factor. Use 0.0 to disable.");
54-
55-
uint32_t normalIterations = 3;
56-
pipeline.add_option("--normal-iterations", normalIterations, "Number of normal smoothing iterations.");
57-
58-
uint32_t lloydIterations = 5;
59-
pipeline.add_option("--lloyd-iterations", lloydIterations, "Number of Lloyd relaxation iterations.");
60-
61-
uint32_t newtonIterations = 30;
62-
pipeline.add_option("--newton-iterations", newtonIterations, "Number of Newton iterations.");
63-
64-
uint32_t newtonM = 7;
65-
pipeline.add_option("--newton-m", newtonM, "Number of Newton evaluations per step for Hessian approximation.");
66-
67-
uint64_t lfsSamples = 10000;
68-
pipeline.add_option("--lfs-samples", lfsSamples, "Number of samples for size adaptation if triangle size adaptation is not 0.0.");
55+
itk::wasm::InputTextStream planesJson;
56+
pipeline.add_option("planes", planesJson, "An array of plane locations to slice the mesh. Each plane is defined by an array of 'origin' and 'spacing' values.")->required()->type_name("INPUT_JSON");
6957

70-
itk::wasm::OutputMesh<MeshType> outputMesh;
71-
pipeline.add_option("output-mesh", outputMesh, "The output repaired mesh.")->type_name("OUTPUT_MESH");
58+
itk::wasm::OutputMesh<OutputMeshType> outputMesh;
59+
pipeline.add_option("polylines", outputMesh, "The output mesh comprised of polylines. Cell data indicates whether part of a closed line. Point data indicates the slice index.")->type_name("OUTPUT_MESH");
7260

7361
ITK_WASM_PARSE(pipeline);
7462

75-
// GEO::initialize(GEO::GEOGRAM_INSTALL_ALL);
76-
77-
GEO::Logger::initialize(); GEO::Logger::instance()->set_quiet(true);
78-
GEO::CmdLine::initialize();
79-
const auto flags = GEO::GEOGRAM_INSTALL_ALL;
80-
GEO::Process::initialize(flags);
81-
GEO::CmdLine::import_arg_group("algo");
63+
std::vector<Vec3D> vertices;
64+
vertices.resize(inputMesh->GetNumberOfPoints());
65+
for (itk::SizeValueType i = 0; i < inputMesh->GetNumberOfPoints(); ++i)
66+
{
67+
typename MeshType::PointType point = inputMesh->GetPoint(i);
68+
Vec3D vertex;
69+
for (unsigned int d = 0; d < Dimension; d++)
70+
{
71+
vertex[d] = point[d];
72+
}
73+
vertices[i] = vertex;
74+
}
8275

83-
GEO::Mesh geoMesh(Dimension, false);
84-
itk::meshToGeogramMesh<MeshType>(inputMesh, geoMesh);
76+
std::vector<Face> faces;
77+
faces.resize(inputMesh->GetNumberOfCells());
78+
using CellIterator = typename MeshType::CellsContainer::ConstIterator;
79+
CellIterator cellIterator = inputMesh->GetCells()->Begin();
80+
CellIterator cellEnd = inputMesh->GetCells()->End();
8581

86-
if (geoMesh.facets.nb() == 0)
82+
while (cellIterator != cellEnd)
8783
{
88-
std::cerr << "The input mesh has no facets." << std::endl;
89-
return EXIT_FAILURE;
84+
const auto cell = cellIterator.Value();
85+
if (cell->GetNumberOfPoints() != 3)
86+
{
87+
std::cerr << "Only triangle cells are supported." << std::endl;
88+
return EXIT_FAILURE;
89+
}
90+
typename MeshType::CellType::PointIdIterator pointIdIterator = cell->PointIdsBegin();
91+
Face face;
92+
for (unsigned int j = 0; j < cell->GetNumberOfPoints(); ++j)
93+
{
94+
face[j] = *pointIdIterator;
95+
++pointIdIterator;
96+
}
97+
faces[cellIterator.Index()] = face;
98+
++cellIterator;
9099
}
91100

92-
if (!geoMesh.facets.are_simplices())
101+
const std::string planesJsonString(std::istreambuf_iterator<char>(planesJson.Get()), {});
102+
auto deserializedAttempt = glz::read_json<Planes>(planesJsonString);
103+
if (!deserializedAttempt)
93104
{
94-
std::cerr << "The mesh needs to be simplicial. Use repair." << std::endl;
105+
const std::string descriptiveError = glz::format_error(deserializedAttempt, planesJsonString);
106+
std::cerr << "Failed to deserialize planes: " << descriptiveError << std::endl;
95107
return EXIT_FAILURE;
96108
}
109+
const auto planes = deserializedAttempt.value();
97110

98-
uint64_t numberPoints = static_cast<uint64_t>(std::ceil(numberPointsPercent * 0.01 * geoMesh.vertices.nb()));
99-
100-
std::cout << "Remeshing with " << numberPoints << " points" << std::endl;
111+
typename OutputMeshType::Pointer outputMeshPtr = OutputMeshType::New();
101112

102-
if(triangleShapeAdaptation != 0.0)
103-
{
104-
triangleShapeAdaptation *= 0.02;
105-
GEO::compute_normals(geoMesh);
106-
if(normalIterations != 0)
107-
{
108-
GEO::simple_Laplacian_smooth(geoMesh, normalIterations, true);
109-
}
110-
GEO::set_anisotropy(geoMesh, triangleShapeAdaptation);
111-
}
112-
else
113-
{
114-
geoMesh.vertices.set_dimension(3);
115-
}
113+
PointIdentifierType pointId = 0;
114+
CellIdentifierType cellId = 0;
115+
size_t sliceIndex = 0;
116+
PointDataContainer *pointData = outputMeshPtr->GetPointData();
117+
CellDataContainer *cellData = outputMeshPtr->GetCellData();
116118

117-
if(triangleSizeAdaptation != 0.0)
118-
{
119-
GEO::compute_sizing_field(geoMesh, triangleSizeAdaptation, lfsSamples);
120-
}
121-
else
119+
typename MeshPlaneIntersectType::Mesh meshPlaneIntersect(vertices, faces);
120+
for (const auto &plane : planes)
122121
{
123-
GEO::AttributesManager& attributes = geoMesh.vertices.attributes();
124-
if(attributes.is_defined("weight"))
122+
const auto paths = meshPlaneIntersect.Intersect(plane);
123+
for (const auto &path : paths)
125124
{
126-
attributes.delete_attribute_store("weight");
125+
CellAutoPointer cell;
126+
cell.TakeOwnership(new PolyLineCellType);
127+
PointIdentifierType polyPointId = 0;
128+
const auto initialPointId = pointId;
129+
for (const auto &point : path.points)
130+
{
131+
typename OutputMeshType::PointType outputPoint;
132+
for (unsigned int d = 0; d < Dimension; d++)
133+
{
134+
outputPoint[d] = point[d];
135+
}
136+
outputMeshPtr->SetPoint(pointId, outputPoint);
137+
cell->SetPointId(polyPointId++, pointId);
138+
pointId++;
139+
pointData->push_back(static_cast<uint16_t>(sliceIndex));
140+
}
141+
if (path.isClosed)
142+
{
143+
cell->SetPointId(polyPointId++, initialPointId);
144+
}
145+
outputMeshPtr->SetCell(cellId++, cell);
146+
const uint16_t isClosed = path.isClosed ? 1 : 0;
147+
cellData->push_back(isClosed);
127148
}
149+
++sliceIndex;
128150
}
129151

130-
GEO::Mesh remesh;
131-
132-
GEO::remesh_smooth(geoMesh, remesh, numberPoints, 0, lloydIterations, newtonIterations, newtonM);
133-
134-
GEO::MeshElementsFlags what = GEO::MeshElementsFlags(
135-
GEO::MESH_VERTICES | GEO::MESH_EDGES | GEO::MESH_FACETS | GEO::MESH_CELLS
136-
);
137-
geoMesh.clear();
138-
geoMesh.copy(remesh, true, what);
139-
140-
typename MeshType::Pointer itkMesh = MeshType::New();
141-
itk::geogramMeshToMesh<MeshType>(geoMesh, itkMesh);
142-
143-
outputMesh.Set(itkMesh);
152+
outputMesh.Set(outputMeshPtr);
144153

145154
return EXIT_SUCCESS;
146155
}
@@ -154,18 +163,18 @@ class PipelineFunctor
154163
using MeshType = TMesh;
155164

156165
itk::wasm::InputMesh<MeshType> testMesh;
157-
pipeline.add_option("input-mesh", testMesh, "The input mesh")->type_name("INPUT_MESH");
166+
pipeline.add_option("input-mesh", testMesh, "The input triangle mesh.")->type_name("INPUT_MESH");
158167

159168
ITK_WASM_PRE_PARSE(pipeline);
160169

161170
typename MeshType::ConstPointer inputMeshRef = testMesh.Get();
162-
return repairMesh<MeshType>(pipeline, inputMeshRef);
171+
return sliceMesh<MeshType>(pipeline, inputMeshRef);
163172
}
164173
};
165174

166175
int main(int argc, char *argv[])
167176
{
168-
itk::wasm::Pipeline pipeline("remesh", "Smooth and remesh a mesh to improve quality.", argc, argv);
177+
itk::wasm::Pipeline pipeline("slice-mesh", "Slice a mesh along planes into polylines.", argc, argv);
169178

170179
return itk::wasm::SupportInputMeshTypes<PipelineFunctor,
171180
// uint8_t,

0 commit comments

Comments
 (0)