From 91e6f669ad1ba31cef628b62456d1991f1c1b166 Mon Sep 17 00:00:00 2001 From: Jirawat I Date: Sun, 8 Aug 2021 21:07:57 +0700 Subject: [PATCH] Release 1.5.2 --- CMakeLists.txt | 24 +- blender/SCAFFOLDER_OP_generate_mesh.py | 5 +- blender/SCAFFOLDER_settings.py | 4 +- cmake/FindLIBIGL.cmake | 2 +- cmake/FindSol2.cmake | 4 +- docs/cmd.md | 34 +- docs/installation.md | 2 +- docs/python.md | 44 +- docs/tutorial_1.md | 4 +- include/MaxHeap.hpp | 6 +- include/Mesh.h | 2 +- include/MeshOperation.h | 2 +- include/OptimalSlice.hpp | 2 +- include/QuadricSimplification.h | 2 +- include/Scaffolder.h | 63 ++ include/Scaffolder_2.h | 141 ---- include/dualmc/dualmc.h | 992 +++++++++++++++++++------ include/dualmc/dualmc.tpp | 474 ------------ include/implicit_function.h | 13 +- include/utils.h | 1 - setup.cfg | 3 +- setup.py | 2 +- src/Main.cpp | 73 +- src/Scaffolder.cpp | 111 +++ src/SliceTest.cpp | 1 - src/python/PyScaffolder.cpp | 16 +- src/python/PyScaffolder.hpp | 21 +- src/python/implements.cpp | 150 +++- src/python/test.py | 144 +++- 29 files changed, 1343 insertions(+), 999 deletions(-) create mode 100644 include/Scaffolder.h delete mode 100644 include/Scaffolder_2.h delete mode 100644 include/dualmc/dualmc.tpp create mode 100644 src/Scaffolder.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index a5c9c87..2a5f177 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -117,23 +117,17 @@ endif() if(BUILD_SCAFFOLDER) # Add your project files file(GLOB MAIN_SOURCES - "${PROJECT_SOURCE_DIR}/include/*.c" - "${PROJECT_SOURCE_DIR}/include/*.cpp" - "${PROJECT_SOURCE_DIR}/include/*.tpp" - "${PROJECT_SOURCE_DIR}/include/*.h" - "${PROJECT_SOURCE_DIR}/include/*.hpp" "${PROJECT_SOURCE_DIR}/include/*/*.c" "${PROJECT_SOURCE_DIR}/include/*/*.cpp" - "${PROJECT_SOURCE_DIR}/include/*/*.tpp" - "${PROJECT_SOURCE_DIR}/include/*/*.h" - "${PROJECT_SOURCE_DIR}/include/*/*.hpp" "${PROJECT_SOURCE_DIR}/src/utils.cpp" + "${PROJECT_SOURCE_DIR}/src/Scaffolder.cpp" "${VCG_INCLUDE_DIR}/wrap/ply/plylib.cpp" "${OpenCTM_INCLUDE_DIR}/*.c" "${OpenCTM_INCLUDE_DIR}/liblzma/*.c" ) include(update_deps_file) update_deps_file("main_sources" "${MAIN_SOURCES}") + # Build binary add_executable(${PROJECT_NAME} ${MAIN_SOURCES} ${PROJECT_SOURCE_DIR}/src/Main.cpp) @@ -145,6 +139,7 @@ if(BUILD_SCAFFOLDER) ${SOL2_INCLUDE_DIR} ${LUA_INCLUDE_DIR} ${OpenCTM_INCLUDE_DIR} + ${LIBIGL_INCLUDE_DIR} "${OpenCTM_INCLUDE_DIR}/liblzma/" ) target_link_libraries(${PROJECT_NAME} PRIVATE @@ -155,10 +150,15 @@ if(BUILD_SCAFFOLDER) ${LUA_LIBRARIES} ) - add_executable("${PROJECT_NAME}.SliceTest" ${MAIN_SOURCES} ${PROJECT_SOURCE_DIR}/src/SliceTest.cpp) + add_executable("${PROJECT_NAME}.SliceTest" + ${PROJECT_SOURCE_DIR}/src/SliceTest.cpp + ${PROJECT_SOURCE_DIR}/src/utils.cpp + ${VCG_INCLUDE_DIR}/wrap/ply/plylib.cpp + ) target_include_directories("${PROJECT_NAME}.SliceTest" PRIVATE "${PROJECT_SOURCE_DIR}/include" ${VCG_INCLUDE_DIR} + ${LIBIGL_INCLUDE_DIR} ${EIGEN3_INCLUDE_DIR} ${OpenCTM_INCLUDE_DIR} "${OpenCTM_INCLUDE_DIR}/liblzma/" @@ -177,8 +177,8 @@ if(BUILD_PYSCAFFOLDER) file(GLOB PYTHON_SOURCES "${PROJECT_SOURCE_DIR}/src/python/*.cpp" - "${PROJECT_SOURCE_DIR}/src/python/*.hpp" - "${PROJECT_SOURCE_DIR}/src/utils.cpp" + ${PROJECT_SOURCE_DIR}/src/Scaffolder.cpp + ${PROJECT_SOURCE_DIR}/src/utils.cpp ) include(update_deps_file) update_deps_file("python_sources" "${PYTHON_SOURCES}") @@ -186,10 +186,12 @@ if(BUILD_PYSCAFFOLDER) pybind11_add_module("Py${PROJECT_NAME}" ${PYTHON_SOURCES}) target_include_directories("Py${PROJECT_NAME}" PRIVATE "${PROJECT_SOURCE_DIR}/include" + "${PROJECT_SOURCE_DIR}/src/python/" ${TBB_INCLUDE_DIR} ${VCG_INCLUDE_DIR} ${EIGEN3_INCLUDE_DIR} ${SOL2_INCLUDE_DIR} + ${LIBIGL_INCLUDE_DIR} ${LUA_INCLUDE_DIR} ${PyBind11_INCLUDE_DIR} ${OpenCTM_INCLUDE_DIR} diff --git a/blender/SCAFFOLDER_OP_generate_mesh.py b/blender/SCAFFOLDER_OP_generate_mesh.py index 5f5eb73..ef3e4dd 100644 --- a/blender/SCAFFOLDER_OP_generate_mesh.py +++ b/blender/SCAFFOLDER_OP_generate_mesh.py @@ -79,7 +79,10 @@ def callback(v): f = read_faces(mesh) def gen(self): - mesh = PyScaffolder.generate_mesh(v, f, params, callback) + if ('__version__' in PyScaffolder.__dict__): + mesh = PyScaffolder.generate_scaffold(v, f, params, callback) + else: + mesh = PyScaffolder.generate_mesh(v, f, params, callback) template_str = "Porosity: %.3f\nSurface Area: %.3f\nSurface Area Ratio: %.3f" props.result1 = template_str % (mesh.porosity, mesh.surface_area, mesh.surface_area_ratio) diff --git a/blender/SCAFFOLDER_settings.py b/blender/SCAFFOLDER_settings.py index e5112e0..7bf8224 100644 --- a/blender/SCAFFOLDER_settings.py +++ b/blender/SCAFFOLDER_settings.py @@ -38,8 +38,8 @@ def set_coff(self, context): is_build_inverse: BoolProperty(name="Build inverse", default=False) grid_size: IntProperty(name="Grid size", default=100, min=10) - grid_offset: IntProperty(name="Grid offset", default=5, min=0) - smooth_step: IntProperty(name="smooth step", default=5, min=0) + grid_offset: IntProperty(name="Grid offset", default=3, min=0) + smooth_step: IntProperty(name="smooth step", default=0, min=0) k_slice: IntProperty(name="k slice", default=100, min=1) k_polygon: IntProperty(name="k polygon", default=4, min=1) diff --git a/cmake/FindLIBIGL.cmake b/cmake/FindLIBIGL.cmake index cf2068b..59c0b71 100644 --- a/cmake/FindLIBIGL.cmake +++ b/cmake/FindLIBIGL.cmake @@ -30,7 +30,7 @@ set(LIBIGL_USE_STATIC_LIBRARY OFF CACHE INTERNAL "prefer STATIC build") set(LIBIGL_SKIP_DOWNLOAD ON CACHE INTERNAL "Skip download libigl") # set(HUNTER_ENABLED ON CACHE INTERNAL "Enable Hunter package manager support") -message(STATUS "USE DIR: ${LIBIGL_INCLUDE_DIR}") +message(STATUS "USE IGL DIR: ${LIBIGL_INCLUDE_DIR}") list(APPEND CMAKE_MODULE_PATH "${LIBIGL_INCLUDE_DIR}/../cmake") if (NOT EIGEN3_INCLUDE_DIR) set(EIGEN3_INCLUDE_DIR "${LIBIGL_INCLUDE_DIR}/../external/eigen") diff --git a/cmake/FindSol2.cmake b/cmake/FindSol2.cmake index 3dc1472..54951aa 100644 --- a/cmake/FindSol2.cmake +++ b/cmake/FindSol2.cmake @@ -56,7 +56,7 @@ file(GLOB SOL2_HEADER_SOURCES ${SOL2_INCLUDE_DIR}/sol*.hpp) source_group(sol2 FILES ${SOL2_HEADER_SOURCES}) list(APPEND CMAKE_MODULE_PATH "${SOL2_INCLUDE_DIR}/../cmake/Modules") -message(STATUS "USE DIR: ${SOL2_INCLUDE_DIR}") +message(STATUS "USE SOL2 DIR: ${SOL2_INCLUDE_DIR}") # # # Libraries # Here, we pull in all the necessary libraries for building examples and tests @@ -143,4 +143,4 @@ endif() if (NOT LUA_FOUND AND NOT LUABUILD_FOUND) message(FATAL_ERROR "sol2 Lua \"${SOL2_LUA_VERSION}\" not found and could not be targeted for building") endif() -message(STATUS "USE DIR: ${LUA_LIBRARIES} ${LUA_INCLUDE_DIR}") +message(STATUS "USE LUA DIR: ${LUA_LIBRARIES} ${LUA_INCLUDE_DIR}") diff --git a/docs/cmd.md b/docs/cmd.md index 72b5c16..15c4315 100644 --- a/docs/cmd.md +++ b/docs/cmd.md @@ -5,9 +5,9 @@ Usage: Scaffolder [OPTION...] INPUT OUTPUT PARAMETERS -h, --help Print help - -i, --input INPUT Input file (STL) + -i, --input INPUT Input file (STL/PLY/OFF/OBJ/VMI) -o, --output OUTPUT Output filename with extension - stl,ply,obj,off [default: out] + stl,ply,obj,off,ctm --params PARAMETERS Combined parameters list: surface[,coff,isolevel,grid_size,k_slice,k_polygon] -q, --quiet Disable verbose output [default: false] @@ -17,7 +17,7 @@ Usage: -n, --surface NAME implicit surface: rectlinear, schwarzp, schwarzd, gyroid, double-p, double-d, double-gyroiod, lidinoid, schoen_iwp, neovius, bcc, - tubular_g_ab, tubular_g_c [default: schwarzp] + tubular_g_ab, tubular_g_c (default: bcc) -g, --grid_size INT (0..60000) Grid size [default: 100] -s, --shell INT (0..60000) Outer thickness (layers) [default:0] @@ -27,23 +27,31 @@ Usage: technique ( [default: false] --export_microstructure Analysis microstructure and export the 2D contours (for debugging) [default: false] + --export_jpeg [X|Y|Z],INT + Export 2D JPEG (Default: Z,100) --k_slice INT (0..60000) K_slice: the number of slicing layers in each direction (used in microstructure analysis) - [default: 100] + (default: 100) --k_polygon INT (>0) K_polygon: the number of closest outer - contour (used in microstructure analysis) [default: 4] + contour (used in microstructure analysis) (default: + 4) -z, --size_optimize DOUBLE (0..1) - Experimental Quadric simplification [default: 0] + Experimental Quadric simplification (default: + 0) --smooth_step INT (0..60000) - Smooth with laplacian (default: 5) - --dirty Disable autoclean [default false] + Smooth with laplacian (default: 0) + --dirty Disable autoclean --minimum_diameter DOUBLE (0..1) - used for removing small orphaned (between 0-1) [default: 0.25] + used for removing small orphaned (between + 0-1) (default: 0.25) --format FORMAT (default, csv) - Format of logging output [default: default] - --output_inverse additional output inverse scaffold [default: false] - --fix_self_intersect Experimental fix self-intersect faces - [default: false] + Format of logging output (default: default) + --output_inverse additional output inverse scaffold + --fix_self_intersect INT Experimental fix self-intersect faces + (default: 0) + --mean_curvature INT Size of mean curvature histogram (default: 0) + --no_intersect Generate 3D mesh without intersect with + original mesh (false) ``` ## Examples diff --git a/docs/installation.md b/docs/installation.md index 10d825f..f0b56e2 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -111,7 +111,7 @@ params = PyScaffolder.Parameter() params.coff = 1.0 # Test the mesh generation from v,f -a = PyScaffolder.[generate_mesh]()(v, f, params) +a = PyScaffolder.generate_scaffold(v, f, params) print(a.porosity) print(a.surface_area_ratio) print(len(a.v), len(a.f)) diff --git a/docs/python.md b/docs/python.md index 8636a31..277ac15 100644 --- a/docs/python.md +++ b/docs/python.md @@ -57,11 +57,11 @@ namespace PyScaffolder { } ``` -The main core functions are `generate_mesh` and `slice_test` which resemble the standalone program: `Scaffolder` and `Scaffolder.Slice`. +The main core functions are `generate_scaffold` and `slice_test` which resemble the standalone program: `Scaffolder` and `Scaffolder.Slice`. -## generate_mesh +## generate_scaffold ```python -def generate_mesh(vertices, faces, params=Parameter(), callback=None) +def generate_scaffold(vertices, faces, params=Parameter(), callback=None) ``` : Generate 3D mesh from input vertices, faces, and parameters ### Parameters @@ -113,7 +113,7 @@ def generate_mesh(vertices, faces, params=Parameter(), callback=None) # use callback to show % progression def progression(n): - print(x) + print(n) mesh_info = generate_mesh(v, f, callback=progression) ``` @@ -147,11 +147,45 @@ def slice_test(vertices, faces, k_slice=100, k_polygon=4, direction=0, callback= print(len(pore_size.minFeret), len(pore_size.maxFeret)) ``` +## marching_cubes +```python +def marching_cubes(f, grid_size=(100,100,100), v_min=(0,0,0), delta=0.01, clean=False, callback=None) +``` +: Slice input mesh into pore sizes + ### Parameters + * `f`: numpy.array representing a discrete isosurface where F(x,y,z) = 0 is boundary + * `grid_size`: *Optional* array, list, tuple, or int representing the number of voxels + * `v_min`: *Optional* array, list, tuple the coordinate of the corner of grid + * `delta`: *Optional* array, list, tuple or double representing the dimension of a voxel + * `clean` *Optional* boolean that enable mesh cleaning after marching cubes + * `callback`: *Optional* function with one `integer` parameter indicating the progression + ###Return + `#!python (v, f)` representing vertices and faces in `np.array` + + ###Example + ```python + import PyScaffolder + Fxyz = [ + 1,1,1,1, + 1,-1,-1,1, + 1,-1,-1,1, + 1,1,1,1 + ]*4 + (v, f) = PyScaffolder.marching_cubes( + Fxyz, + grid_size=4, + delta=0.25, + v_min=(-.5, -.5, -.5), + callback = lambda x: print(x) + ) + ``` + ## Parameter ```python class Parameter: is_build_inverse = False is_intersect = True + verbose = False shell = 0 grid_offset = 5 smooth_step = 5 @@ -190,4 +224,4 @@ class MeshInfo: surface_area = 0.0 surface_area_ratio = 0.0 ``` -: Collection of triangular mesh resulted from [generate_mesh](#generate_mesh) \ No newline at end of file +: Collection of triangular mesh resulted from [generate_scaffold](#generate_scaffold) \ No newline at end of file diff --git a/docs/tutorial_1.md b/docs/tutorial_1.md index 39dfbbc..3df6c55 100644 --- a/docs/tutorial_1.md +++ b/docs/tutorial_1.md @@ -7,7 +7,7 @@ angular frequency ($w$), and isolevel ($t$) are required (See [Command line](cmd ## Angular frequency According to [prior study](https://github.com/nodtem66/Scaffolder/blob/master/data/data_visualization.ipynb), -**$w$ is inversely proportational to the pore size** (see [Tips](tips.md)). That is when $w$ increases, +**$w$ is inversely proportational to the pore size** (see [Tips](tips.md)). That is, when $w$ increases, the pore size will be decreased. For gyroid, $w \approx \pi/\text{pore size}$. If we want pore size of 0.25 (unit as same as `sphere.stl`), we try to initially set $w = \pi/0.25 \approx 12$ @@ -107,7 +107,7 @@ We use binary search to adjust the new value of $t$ to find the optimal value wi | -- | -- | -- | -- | -- | -- | -- | -- | -- | | Porosity | 0.508 | 0.85 | 0.67 | 0.59 | 0.634 | 0.613 | 0.602 | 0.597 -We can also use the classic numerical method called [Secant method](https://en.wikipedia.org/wiki/Secant_method) to find the optimal $t$ for the target porosity. +We can also use the classical numerical method called [Secant method](https://en.wikipedia.org/wiki/Secant_method) to find the optimal $t$ for the target porosity. $$ F(t) = Porosity(t) - 0.6 = 0 diff --git a/include/MaxHeap.hpp b/include/MaxHeap.hpp index 7e44eb8..f87ee58 100644 --- a/include/MaxHeap.hpp +++ b/include/MaxHeap.hpp @@ -1,4 +1,5 @@ -#pragma once +#ifndef MAXHEAP_HPP_INCLUDED +#define MAXHEAP_HPP_INCLUDED #include #include @@ -42,4 +43,5 @@ class MaxHeap { (*arr)[0] = d; heapify(0); } -}; \ No newline at end of file +}; +#endif \ No newline at end of file diff --git a/include/Mesh.h b/include/Mesh.h index de57ba2..f815fe9 100644 --- a/include/Mesh.h +++ b/include/Mesh.h @@ -1,6 +1,6 @@ -#pragma once #ifndef MESH_INCLUDED #define MESH_INCLUDED + #include #include #include diff --git a/include/MeshOperation.h b/include/MeshOperation.h index 43f11e3..fe00767 100644 --- a/include/MeshOperation.h +++ b/include/MeshOperation.h @@ -1,6 +1,6 @@ -#pragma once #ifndef MESHOPERATION_INCLUDED #define MESHOPERATION_INCLUDED + #include "Mesh.h" #include "utils.h" diff --git a/include/OptimalSlice.hpp b/include/OptimalSlice.hpp index 3eb191c..5d01965 100644 --- a/include/OptimalSlice.hpp +++ b/include/OptimalSlice.hpp @@ -2,9 +2,9 @@ * This is the implement slice algorithm from Rodrigo, 2017 (An Optimal Algorithm for 3D Triangle Mesh Slicing) * which is claimed to be faster than slic3r and CGAL method. */ -#pragma once #ifndef OPTIMALSLICE_INCLUDED #define OPTIMALSLICE_INCLUDED + #ifndef SLICE_PRECISION // There's a problem of double hashing with the precision less than 1e-8 (e.g. 1e-10) // when performed contour constructing diff --git a/include/QuadricSimplification.h b/include/QuadricSimplification.h index 686f79b..d2a96e5 100644 --- a/include/QuadricSimplification.h +++ b/include/QuadricSimplification.h @@ -19,9 +19,9 @@ * for more details. * * * ****************************************************************************/ -#pragma once #ifndef QUADSIM_INCLUDED #define QUADSIM_INCLUDED + #include #include #include diff --git a/include/Scaffolder.h b/include/Scaffolder.h new file mode 100644 index 0000000..b905955 --- /dev/null +++ b/include/Scaffolder.h @@ -0,0 +1,63 @@ +#ifndef SCAFFOLDER_INCLUDED +#define SCAFFOLDER_INCLUDED + +#ifndef _XOPEN_SOURCE +#define _XOPEN_SOURCE 600 +#endif +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "Mesh.h" +#include "sol/sol.hpp" +#include "ProgressBar.hpp" + +#define VERSION "v1.5.2" +#define PROGRESS_BAR_COLUMN 40 + +#define SCAFFOLDER_FORMAT_DEFAULT 0 +#define SCAFFOLDER_FORMAT_CSV 1 + + +typedef struct index_type { + size_t x; size_t y; size_t z; +} index_type; +typedef std::map Queue_t; + +inline size_t indexFromIJK(size_t i, size_t j, size_t k, Eigen::RowVector3i grid_size); +inline void indexToIJK(size_t index, Eigen::RowVector3i grid_size, index_type& r); +inline bool MarkAndSweepNeighbor( + Eigen::VectorXd& W, + index_type& index, + Queue_t& queue, + Eigen::RowVector3i grid_size, + double value = 0.0, + bool findAbove = false +); + +void marching_cube( + TMesh& mesh, + Eigen::MatrixXd& Fxyz, + Eigen::RowVector3i grid_size, + Eigen::RowVector3d& Vmin, + double delta, + bool verbose = true, + bool dirty = false +); +bool null_callback(int pos, const char* str); + +bool qsim_callback(int pos, const char* str); +extern ProgressBar qsim_progress; + +#endif \ No newline at end of file diff --git a/include/Scaffolder_2.h b/include/Scaffolder_2.h deleted file mode 100644 index b50e7f7..0000000 --- a/include/Scaffolder_2.h +++ /dev/null @@ -1,141 +0,0 @@ -#pragma once -#ifndef SCAFFOLDER_INCLUDED -#define SCAFFOLDER_INCLUDED -#ifndef _XOPEN_SOURCE -#define _XOPEN_SOURCE 600 -#endif -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include "cxxopts.hpp" -#include "dualmc/dualmc.h" -#include "sol/sol.hpp" -#include "toojpeg/toojpeg.h" - -#include "ProgressBar.hpp" -#include "OptimalSlice.hpp" -#include "implicit_function.h" -#include "MeshOperation.h" -#include "utils.h" -#include "QuadricSimplification.h" - -#define VERSION "v1.5.1" -#define PROGRESS_BAR_COLUMN 40 - -#define SCAFFOLDER_FORMAT_DEFAULT 0 -#define SCAFFOLDER_FORMAT_CSV 1 - -#define LOG if (log_format == SCAFFOLDER_FORMAT_DEFAULT) log -#define CSV if (log_format == SCAFFOLDER_FORMAT_CSV) log - -typedef struct index_type { - size_t x; size_t y; size_t z; -} index_type; -typedef std::map Queue_t; - -// Flatten between 1D and 3D -// https://stackoverflow.com/questions/7367770/how-to-flatten-or-index-3d-array-in-1d-array -inline size_t indexFromIJK(size_t i, size_t j, size_t k, Eigen::RowVector3i grid_size) { - return i + grid_size(0) * (j + grid_size(1) * k); -} - -inline void indexToIJK(size_t index, Eigen::RowVector3i grid_size, index_type& r) { - r.z = index / grid_size(0) / grid_size(1); - index -= r.z * grid_size(0) * grid_size(1); - r.y = index / grid_size(0); - r.x = index % grid_size(0); -} - -inline bool MarkAndSweepNeighbor(Eigen::VectorXd& W, index_type& index, Queue_t& queue, Eigen::RowVector3i grid_size, double value = 0.0, bool findAbove = false) { - bool isBorder = false; - for (int8_t di = -1; di <= 1; di++) { - for (int8_t dj = -1; dj <= 1; dj++) { - for (int8_t dk = -1; dk <= 1; dk++) { - if (di == 0 && dj == 0 && dk == 0) continue; - const size_t id = indexFromIJK(index.x + di, index.y + dj, index.z + dk, grid_size); - //std::cout << value << " " << W(id) << std::endl; - if ((findAbove && W(id) >= value - eps2) || (!findAbove && W(id) <= value + eps2)) { - isBorder = true; - break; - } - } - if (isBorder) break; - } - if (isBorder) break; - } - if (isBorder) { - for (int8_t di = -1; di <= 1; di++) { - for (int8_t dj = -1; dj <= 1; dj++) { - for (int8_t dk = -1; dk <= 1; dk++) { - if (di == 0 && dj == 0 && dk == 0) continue; - const size_t id = indexFromIJK(index.x + di, index.y + dj, index.z + dk, grid_size); - if (W(id) >= 0.5 && W(id) < 1.1) { - queue.insert({ id, true }); - } - } - } - } - } - return isBorder; -} - -inline void marching_cube(TMesh &mesh, Eigen::MatrixXd &Fxyz, Eigen::RowVector3i grid_size, Eigen::RowVector3d &Vmin, double delta, bool verbose = true, bool dirty = false) { - { - if (verbose) std::cout << "[Marching Cube] " << std::endl; - dualmc::DualMC builder; - std::vector mc_vertices; - std::vector mc_quads; - builder.build((double const*)Fxyz.data(), grid_size(0), grid_size(1), grid_size(2), 0, true, false, mc_vertices, mc_quads); - TMesh::VertexIterator vi = vcg::tri::Allocator::AddVertices(mesh, mc_vertices.size()); - TMesh::FaceIterator fi = vcg::tri::Allocator::AddFaces(mesh, mc_quads.size() * 2); - std::vector vp(mc_vertices.size()); - for (size_t i = 0, len = mc_vertices.size(); i < len; i++, ++vi) { - vp[i] = &(*vi); - vi->P() = TMesh::CoordType( - Vmin(0) + mc_vertices[i].x * delta, - Vmin(1) + mc_vertices[i].y * delta, - Vmin(2) + mc_vertices[i].z * delta - ); - } - for (size_t i = 0, len = mc_quads.size(); i < len; i++, ++fi) { - fi->V(0) = vp[mc_quads[i].i0]; - fi->V(1) = vp[mc_quads[i].i1]; - fi->V(2) = vp[mc_quads[i].i2]; - ++fi; - fi->V(0) = vp[mc_quads[i].i2]; - fi->V(1) = vp[mc_quads[i].i3]; - fi->V(2) = vp[mc_quads[i].i0]; - } - if (!dirty) { - vcg::tri::Clean::RemoveDuplicateFace(mesh); - vcg::tri::Clean::RemoveDuplicateVertex(mesh); - vcg::tri::Clean::RemoveUnreferencedVertex(mesh); - } - } -} - -inline bool null_callback(int pos, const char* str) { - return true; -} - -inline void set_shorten_function(sol::state& lua) { - lua.script("abs, acos, asin, atan, atan2 = math.abs, math.acos, math.atan, math.atan2"); - lua.script("ceil, cos, deg, exp, floor = math.ceil, math.cos, math.deg, math.exp, math.floor"); - lua.script("log, log10, max, min, mod = math.log, math.log10, math.max, math.min, math.mod"); - lua.script("pow, rad, sin, sqrt, tan = math.pow, math.rad, math.sin, math.sqrt, math.tan"); - lua.script("frexp, ldexp, random, randomseed = math.frexp, math.ldexp, math.random, math.randomseed"); - lua.script("local pi = math.pi"); -} -#endif \ No newline at end of file diff --git a/include/dualmc/dualmc.h b/include/dualmc/dualmc.h index 74b653c..93f3cf9 100644 --- a/include/dualmc/dualmc.h +++ b/include/dualmc/dualmc.h @@ -17,238 +17,782 @@ #include #include #include "ProgressBar.hpp" +#include "oneapi/tbb.h" #ifndef PROGRESS_BAR_COLUMN #define PROGRESS_BAR_COLUMN 40 #endif +#define USE_PARALLEL namespace dualmc { - typedef double VertexComponentsType; - typedef int32_t QuadIndexType; - - /// vertex structure for dual points - struct Vertex { - /// non-initializing constructor - Vertex(); - - /// initializing constructor - Vertex(VertexComponentsType x, VertexComponentsType y, VertexComponentsType z); - - /// initializing constructor - Vertex(Vertex const& v); - - // components - VertexComponentsType x, y, z; - }; - - /// quad indices structure - struct Quad { - /// non-initializing constructor - Quad(); - - /// initializing constructor - Quad(QuadIndexType i0, QuadIndexType i1, QuadIndexType i2, QuadIndexType i3); - - // quad indices - QuadIndexType i0, i1, i2, i3; - - }; - - /// \class DualMC - /// \author Dominik Wodniok - /// \date 2009 - /// Class which implements the dual marching cubes algorithm from Gregory M. Nielson. - /// Faces and vertices of the standard marching cubes algorithm correspond to - /// vertices and faces in the dual algorithm. As a vertex in standard marching cubes - /// usually is shared by 4 faces, the dual mesh is entirely made from quadrangles. - /// Unfortunately, under rare circumstances the original algorithm can create - /// non-manifold meshes. See the remarks of the original paper on this. - /// The class optionally can guarantee manifold meshes by taking the Manifold - /// Dual Marching Cubes approach from Rephael Wenger as described in - /// chapter 3.3.5 of his book "Isosurfaces: Geometry, Topology, and Algorithms". - template class DualMC { - public: - // typedefs - typedef T VolumeDataType; - - /// Extracts the iso surface for a given volume and iso value. - /// Output is a list of vertices and a list of indices, which connect - /// vertices to quads. - /// The quad mesh either uses shared vertex indices or is a quad soup if - /// desired. - void build( - VolumeDataType const* data, - int32_t const dimX, int32_t const dimY, int32_t const dimZ, - VolumeDataType const iso, - bool const generateManifold, - bool const generateSoup, - std::vector& vertices, - std::vector& quads - ); - - private: - - /// Extract quad mesh with shared vertex indices. - void buildSharedVerticesQuads( - VolumeDataType const iso, - std::vector& vertices, - std::vector& quads - ); - - /// Extract quad soup. - void buildQuadSoup( - VolumeDataType const iso, - std::vector& vertices, - std::vector& quads - ); - - - private: - - /// enum with edge codes for a 12-bit voxel edge mask to indicate - /// grid edges which intersect the ISO surface of classic marching cubes - enum DMCEdgeCode { - EDGE0 = 1, - EDGE1 = 1 << 1, - EDGE2 = 1 << 2, - EDGE3 = 1 << 3, - EDGE4 = 1 << 4, - EDGE5 = 1 << 5, - EDGE6 = 1 << 6, - EDGE7 = 1 << 7, - EDGE8 = 1 << 8, - EDGE9 = 1 << 9, - EDGE10 = 1 << 10, - EDGE11 = 1 << 11, - FORCE_32BIT = 0xffffffff - }; - - /// get the 8-bit in-out mask for the voxel corners of the cell cube at (cx,cy,cz) - /// and the given iso value - int getCellCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso) const; - - /// Get the 12-bit dual point code mask, which encodes the traditional - /// marching cube vertices of the traditional marching cubes face which - /// corresponds to the dual point. - /// This is also where the manifold dual marching cubes algorithm is - /// implemented. - int getDualPointCode(int32_t const cx, int32_t const cy, int32_t const cz, - VolumeDataType const iso, DMCEdgeCode const edge) const; - - /// Given a dual point code and iso value, compute the dual point. - void calculateDualPoint(int32_t const cx, int32_t const cy, int32_t const cz, - VolumeDataType const iso, int const pointCode, Vertex& v) const; - - /// Get the shared index of a dual point which is uniquly identified by its - /// cell cube index and a cube edge. The dual point is computed, - /// if it has not been computed before. - QuadIndexType getSharedDualPointIndex(int32_t const cx, int32_t const cy, int32_t const cz, - VolumeDataType const iso, DMCEdgeCode const edge, - std::vector& vertices); - - /// Compute a linearized cell cube index. - int32_t gA(int32_t const x, int32_t const y, int32_t const z) const; - - private: - // static lookup tables needed for (manifold) dual marching cubes - - /// Dual Marching Cubes table - /// Encodes the edge vertices for the 256 marching cubes cases. - /// A marching cube case produces up to four faces and ,thus, up to four - /// dual points. - static int32_t const dualPointsList[256][4]; - - /// Table which encodes the ambiguous face of cube configurations, which - /// can cause non-manifold meshes. - /// Needed for manifold dual marching cubes. - static uint8_t const problematicConfigs[256]; - - private: - - /// convenience volume extent array for x-,y-, and z-dimension - int32_t dims[3]; - - /// convenience volume data point - VolumeDataType const* data; - - /// store whether the manifold dual marching cubes algorithm should be - /// applied. - bool generateManifold; - - /// Dual point key structure for hashing of shared vertices - struct DualPointKey { - // a dual point can be uniquely identified by ite linearized volume cell - // id and point code - int32_t linearizedCellID; - int pointCode; - /// Equal operator for unordered map - bool operator==(DualPointKey const& other) const; - }; - - /// Functor for dual point key hash generation - struct DualPointKeyHash { - size_t operator()(DualPointKey const& k) const { - return size_t(k.linearizedCellID) | (size_t(k.pointCode) << 32u); - } - }; - - /// Hash map for shared vertex index computations - std::unordered_map pointToIndex; - }; - - // inline function definitions - - //------------------------------------------------------------------------------ - - inline - Vertex::Vertex() {} - - //------------------------------------------------------------------------------ - - inline - Vertex::Vertex( - VertexComponentsType x, - VertexComponentsType y, - VertexComponentsType z - ) : x(x), y(y), z(z) {} - - //------------------------------------------------------------------------------ - - inline - Vertex::Vertex(Vertex const& v) : x(v.x), y(v.y), z(v.z) {} - - //------------------------------------------------------------------------------ - - inline - Quad::Quad() {} - - //------------------------------------------------------------------------------ - - inline - Quad::Quad( - QuadIndexType i0, - QuadIndexType i1, - QuadIndexType i2, - QuadIndexType i3 - ) : i0(i0), i1(i1), i2(i2), i3(i3) {} - - //------------------------------------------------------------------------------ - - template inline - int32_t DualMC::gA(int32_t const x, int32_t const y, int32_t const z) const { - return x + dims[0] * (y + dims[1] * z); - } - - //------------------------------------------------------------------------------ - template inline - bool DualMC::DualPointKey::operator==(typename DualMC::DualPointKey const& other) const { - return linearizedCellID == other.linearizedCellID && pointCode == other.pointCode; - } - -#include "dualmc.tpp" + typedef double VertexComponentsType; + typedef int32_t QuadIndexType; + + /// vertex structure for dual points + struct Vertex { + /// non-initializing constructor + Vertex(); + + /// initializing constructor + Vertex(VertexComponentsType x, VertexComponentsType y, VertexComponentsType z); + + /// initializing constructor + Vertex(Vertex const& v); + + // components + VertexComponentsType x, y, z; + }; + + /// quad indices structure + struct Quad { + /// non-initializing constructor + Quad(); + + /// initializing constructor + Quad(QuadIndexType i0, QuadIndexType i1, QuadIndexType i2, QuadIndexType i3); + + // quad indices + QuadIndexType i0, i1, i2, i3; + + }; + + /// \class DualMC + /// \author Dominik Wodniok + /// \date 2009 + /// Class which implements the dual marching cubes algorithm from Gregory M. Nielson. + /// Faces and vertices of the standard marching cubes algorithm correspond to + /// vertices and faces in the dual algorithm. As a vertex in standard marching cubes + /// usually is shared by 4 faces, the dual mesh is entirely made from quadrangles. + /// Unfortunately, under rare circumstances the original algorithm can create + /// non-manifold meshes. See the remarks of the original paper on this. + /// The class optionally can guarantee manifold meshes by taking the Manifold + /// Dual Marching Cubes approach from Rephael Wenger as described in + /// chapter 3.3.5 of his book "Isosurfaces: Geometry, Topology, and Algorithms". + template class DualMC { + + public: + // typedefs + typedef T VolumeDataType; + + // Callback function to send the progression + std::function callback = NULL; + + /// Extracts the iso surface for a given volume and iso value. + /// Output is a list of vertices and a list of indices, which connect + /// vertices to quads. + /// The quad mesh either uses shared vertex indices or is a quad soup if + /// desired. + void build( + VolumeDataType const* data, + int32_t const dimX, int32_t const dimY, int32_t const dimZ, + VolumeDataType const iso, + bool const generateManifold, + bool const generateSoup, + std::vector& vertices, + std::vector& quads + ); + + private: + + /// Extract quad mesh with shared vertex indices. + void buildSharedVerticesQuads( + VolumeDataType const iso, + std::vector& vertices, + std::vector& quads + ); + + /// Extract quad soup. + void buildQuadSoup( + VolumeDataType const iso, + std::vector& vertices, + std::vector& quads + ); + + + private: + + /// enum with edge codes for a 12-bit voxel edge mask to indicate + /// grid edges which intersect the ISO surface of classic marching cubes + enum DMCEdgeCode { + EDGE0 = 1, + EDGE1 = 1 << 1, + EDGE2 = 1 << 2, + EDGE3 = 1 << 3, + EDGE4 = 1 << 4, + EDGE5 = 1 << 5, + EDGE6 = 1 << 6, + EDGE7 = 1 << 7, + EDGE8 = 1 << 8, + EDGE9 = 1 << 9, + EDGE10 = 1 << 10, + EDGE11 = 1 << 11, + FORCE_32BIT = 0xffffffff + }; + + /// get the 8-bit in-out mask for the voxel corners of the cell cube at (cx,cy,cz) + /// and the given iso value + int getCellCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso) const; + + /// Get the 12-bit dual point code mask, which encodes the traditional + /// marching cube vertices of the traditional marching cubes face which + /// corresponds to the dual point. + /// This is also where the manifold dual marching cubes algorithm is + /// implemented. + int getDualPointCode(int32_t const cx, int32_t const cy, int32_t const cz, + VolumeDataType const iso, DMCEdgeCode const edge) const; + + /// Given a dual point code and iso value, compute the dual point. + void calculateDualPoint(int32_t const cx, int32_t const cy, int32_t const cz, + VolumeDataType const iso, int const pointCode, Vertex& v) const; + + /// Compute a linearized cell cube index. + int32_t gA(int32_t const x, int32_t const y, int32_t const z) const; + + private: + // static lookup tables needed for (manifold) dual marching cubes + + /// Dual Marching Cubes table + /// Encodes the edge vertices for the 256 marching cubes cases. + /// A marching cube case produces up to four faces and ,thus, up to four + /// dual points. + static int32_t const dualPointsList[256][4]; + + /// Table which encodes the ambiguous face of cube configurations, which + /// can cause non-manifold meshes. + /// Needed for manifold dual marching cubes. + static uint8_t const problematicConfigs[256]; + + private: + + /// convenience volume extent array for x-,y-, and z-dimension + int32_t dims[3]; + + /// convenience volume data point + VolumeDataType const* data; + + /// store whether the manifold dual marching cubes algorithm should be + /// applied. + bool generateManifold; + + /// Dual point key structure for hashing of shared vertices + struct DualPointKey { + // a dual point can be uniquely identified by ite linearized volume cell + // id and point code + int32_t linearizedCellID; + int pointCode; + /// Equal operator for unordered map + bool operator==(DualPointKey const& other) const; + }; + + /// Functor for dual point key hash generation + struct DualPointKeyHash { + size_t operator()(DualPointKey const& k) const { + return size_t(k.linearizedCellID) | (size_t(k.pointCode) << 32u); + } + }; + + // TBB Mutex type + typedef tbb::queuing_mutex MutexType; + typedef tbb::queuing_mutex CallbackMutexType; + + + /// Hash map for shared vertex index computations + typedef std::unordered_map PointToIndexMap; + //std::unordered_map pointToIndex; + + /// Get the shared index of a dual point which is uniquly identified by its + /// cell cube index and a cube edge. The dual point is computed, + /// if it has not been computed before. + QuadIndexType getSharedDualPointIndex(int32_t const cx, int32_t const cy, int32_t const cz, + VolumeDataType const iso, DMCEdgeCode const edge, + std::vector& vertices, PointToIndexMap& p); + }; + + // inline function definitions + + //------------------------------------------------------------------------------ + + inline + Vertex::Vertex() {} + + //------------------------------------------------------------------------------ + + inline + Vertex::Vertex( + VertexComponentsType x, + VertexComponentsType y, + VertexComponentsType z + ) : x(x), y(y), z(z) {} + + //------------------------------------------------------------------------------ + + inline + Vertex::Vertex(Vertex const& v) : x(v.x), y(v.y), z(v.z) {} + + //------------------------------------------------------------------------------ + + inline + Quad::Quad() {} + + //------------------------------------------------------------------------------ + + inline + Quad::Quad( + QuadIndexType i0, + QuadIndexType i1, + QuadIndexType i2, + QuadIndexType i3 + ) : i0(i0), i1(i1), i2(i2), i3(i3) {} + + //------------------------------------------------------------------------------ + + template inline + int32_t DualMC::gA(int32_t const x, int32_t const y, int32_t const z) const { + return x + dims[0] * (y + dims[1] * z); + } + + //------------------------------------------------------------------------------ + template inline + bool DualMC::DualPointKey::operator==(typename DualMC::DualPointKey const& other) const { + return linearizedCellID == other.linearizedCellID && pointCode == other.pointCode; + } + + template inline + int DualMC::getCellCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso) const { + // determine for each cube corner if it is outside or inside + int code = 0; + if (data[gA(cx, cy, cz)] >= iso) + code |= 1; + if (data[gA(cx + 1, cy, cz)] >= iso) + code |= 2; + if (data[gA(cx, cy + 1, cz)] >= iso) + code |= 4; + if (data[gA(cx + 1, cy + 1, cz)] >= iso) + code |= 8; + if (data[gA(cx, cy, cz + 1)] >= iso) + code |= 16; + if (data[gA(cx + 1, cy, cz + 1)] >= iso) + code |= 32; + if (data[gA(cx, cy + 1, cz + 1)] >= iso) + code |= 64; + if (data[gA(cx + 1, cy + 1, cz + 1)] >= iso) + code |= 128; + return code; + } + + //------------------------------------------------------------------------------ + + template inline + int DualMC::getDualPointCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso, DMCEdgeCode const edge) const { + int cubeCode = getCellCode(cx, cy, cz, iso); + + // is manifold dual marching cubes desired? + if (generateManifold) { + // The Manifold Dual Marching Cubes approach from Rephael Wenger as described in + // chapter 3.3.5 of his book "Isosurfaces: Geometry, Topology, and Algorithms" + // is implemente here. + // If a problematic C16 or C19 configuration shares the ambiguous face + // with another C16 or C19 configuration we simply invert the cube code + // before looking up dual points. Doing this for these pairs ensures + // manifold meshes. + // But this removes the dualism to marching cubes. + + // check if we have a potentially problematic configuration + uint8_t const direction = problematicConfigs[uint8_t(cubeCode)]; + // If the direction code is in {0,...,5} we have a C16 or C19 configuration. + if (direction != 255) { + // We have to check the neighboring cube, which shares the ambiguous + // face. For this we decode the direction. This could also be done + // with another lookup table. + // copy current cube coordinates into an array. + int32_t neighborCoords[] = { cx,cy,cz }; + // get the dimension of the non-zero coordinate axis + unsigned int const component = direction >> 1; + // get the sign of the direction + int32_t delta = (direction & 1) == 1 ? 1 : -1; + // modify the correspong cube coordinate + neighborCoords[component] += delta; + // have we left the volume in this direction? + if (neighborCoords[component] >= 0 && neighborCoords[component] < (dims[component] - 1)) { + // get the cube configuration of the relevant neighbor + int neighborCubeCode = getCellCode(neighborCoords[0], neighborCoords[1], neighborCoords[2], iso); + // Look up the neighbor configuration ambiguous face direction. + // If the direction is valid we have a C16 or C19 neighbor. + // As C16 and C19 have exactly one ambiguous face this face is + // guaranteed to be shared for the pair. + if (problematicConfigs[uint8_t(neighborCubeCode)] != 255) { + // replace the cube configuration with its inverse. + cubeCode ^= 0xff; + } + } + } + } + for (int i = 0; i < 4; ++i) + if (dualPointsList[cubeCode][i] & edge) { + return dualPointsList[cubeCode][i]; + } + return 0; + } + + //------------------------------------------------------------------------------ + + template inline + void DualMC::calculateDualPoint(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso, int const pointCode, Vertex& v) const { + // initialize the point with lower voxel coordinates + v.x = cx; + v.y = cy; + v.z = cz; + + // compute the dual point as the mean of the face vertices belonging to the + // original marching cubes face + Vertex p; + p.x = 0; + p.y = 0; + p.z = 0; + int points = 0; + + // sum edge intersection vertices using the point code + if (pointCode & EDGE0) { + p.x += ((float)iso - (float)data[gA(cx, cy, cz)]) / ((float)data[gA(cx + 1, cy, cz)] - (float)data[gA(cx, cy, cz)]); + points++; + } + + if (pointCode & EDGE1) { + p.x += 1.0f; + p.z += ((float)iso - (float)data[gA(cx + 1, cy, cz)]) / ((float)data[gA(cx + 1, cy, cz + 1)] - (float)data[gA(cx + 1, cy, cz)]); + points++; + } + + if (pointCode & EDGE2) { + p.x += ((float)iso - (float)data[gA(cx, cy, cz + 1)]) / ((float)data[gA(cx + 1, cy, cz + 1)] - (float)data[gA(cx, cy, cz + 1)]); + p.z += 1.0f; + points++; + } + + if (pointCode & EDGE3) { + p.z += ((float)iso - (float)data[gA(cx, cy, cz)]) / ((float)data[gA(cx, cy, cz + 1)] - (float)data[gA(cx, cy, cz)]); + points++; + } + + if (pointCode & EDGE4) { + p.x += ((float)iso - (float)data[gA(cx, cy + 1, cz)]) / ((float)data[gA(cx + 1, cy + 1, cz)] - (float)data[gA(cx, cy + 1, cz)]); + p.y += 1.0f; + points++; + } + + if (pointCode & EDGE5) { + p.x += 1.0f; + p.z += ((float)iso - (float)data[gA(cx + 1, cy + 1, cz)]) / ((float)data[gA(cx + 1, cy + 1, cz + 1)] - (float)data[gA(cx + 1, cy + 1, cz)]); + p.y += 1.0f; + points++; + } + + if (pointCode & EDGE6) { + p.x += ((float)iso - (float)data[gA(cx, cy + 1, cz + 1)]) / ((float)data[gA(cx + 1, cy + 1, cz + 1)] - (float)data[gA(cx, cy + 1, cz + 1)]); + p.z += 1.0f; + p.y += 1.0f; + points++; + } + + if (pointCode & EDGE7) { + p.z += ((float)iso - (float)data[gA(cx, cy + 1, cz)]) / ((float)data[gA(cx, cy + 1, cz + 1)] - (float)data[gA(cx, cy + 1, cz)]); + p.y += 1.0f; + points++; + } + + if (pointCode & EDGE8) { + p.y += ((float)iso - (float)data[gA(cx, cy, cz)]) / ((float)data[gA(cx, cy + 1, cz)] - (float)data[gA(cx, cy, cz)]); + points++; + } + + if (pointCode & EDGE9) { + p.x += 1.0f; + p.y += ((float)iso - (float)data[gA(cx + 1, cy, cz)]) / ((float)data[gA(cx + 1, cy + 1, cz)] - (float)data[gA(cx + 1, cy, cz)]); + points++; + } + + if (pointCode & EDGE10) { + p.x += 1.0f; + p.y += ((float)iso - (float)data[gA(cx + 1, cy, cz + 1)]) / ((float)data[gA(cx + 1, cy + 1, cz + 1)] - (float)data[gA(cx + 1, cy, cz + 1)]); + p.z += 1.0f; + points++; + } + + if (pointCode & EDGE11) { + p.z += 1.0f; + p.y += ((float)iso - (float)data[gA(cx, cy, cz + 1)]) / ((float)data[gA(cx, cy + 1, cz + 1)] - (float)data[gA(cx, cy, cz + 1)]); + points++; + } + + // divide by number of accumulated points + float invPoints = 1.0f / (float)points; + p.x *= invPoints; + p.y *= invPoints; + p.z *= invPoints; + + // offset point by voxel coordinates + v.x += p.x; + v.y += p.y; + v.z += p.z; + } + + //------------------------------------------------------------------------------ + + template inline + QuadIndexType DualMC::getSharedDualPointIndex( + int32_t const cx, int32_t const cy, int32_t const cz, + VolumeDataType const iso, DMCEdgeCode const edge, + std::vector& vertices, + PointToIndexMap& pointToIndex + ) { + // create a key for the dual point from its linearized cell ID and point code + DualPointKey key; + key.linearizedCellID = gA(cx, cy, cz); + key.pointCode = getDualPointCode(cx, cy, cz, iso, edge); + + // have we already computed the dual point? + auto iterator = pointToIndex.find(key); + if (iterator != pointToIndex.end()) { + // just return the dual point index + return iterator->second; + } + else { + // create new vertex and vertex id + QuadIndexType newVertexId = vertices.size(); + vertices.emplace_back(); + calculateDualPoint(cx, cy, cz, iso, key.pointCode, vertices.back()); + // insert vertex ID into map and also return it + pointToIndex[key] = newVertexId; + return newVertexId; + } + } + + //------------------------------------------------------------------------------ + + template inline + void DualMC::build( + VolumeDataType const* data, + int32_t const dimX, int32_t const dimY, int32_t const dimZ, + VolumeDataType const iso, + bool const generateManifold, + bool const generateSoup, + std::vector& vertices, + std::vector& quads + ) { + + // set members + this->dims[0] = dimX; + this->dims[1] = dimY; + this->dims[2] = dimZ; + this->data = data; + this->generateManifold = generateManifold; + + // clear vertices and quad indices + vertices.clear(); + quads.clear(); + + // generate quad soup or shared vertices quad list + //if (generateSoup) { + buildQuadSoup(iso, vertices, quads); + //} + ///else { + //buildSharedVerticesQuads(iso, vertices, quads); + //} + } + + //------------------------------------------------------------------------------ + + template inline + void DualMC::buildQuadSoup( + VolumeDataType const iso, + std::vector& vertices, + std::vector& quads + ) { + + int32_t const reducedX = dims[0] - 2; + int32_t const reducedY = dims[1] - 2; + int32_t const reducedZ = dims[2] - 2; + + if (callback != NULL) callback(0); + size_t total_progress = reducedZ; + size_t progress = 0; + + // iterate voxels +#ifndef USE_PARALLEL + for (int32_t z = 0; z < reducedZ; ++z) { +#define local_vertices vertices +#else + static MutexType writeMutex; + static MutexType callbackMutex; + static tbb::affinity_partitioner ap; + tbb::parallel_for( + tbb::blocked_range(0, reducedZ), + [&](const tbb::blocked_range r) { + std::vector local_vertices; + Vertex vertex0; + Vertex vertex1; + Vertex vertex2; + Vertex vertex3; + int pointCode; + for (int32_t z = r.begin(); z < r.end(); z++) { +#endif + for (int32_t y = 0; y < reducedY; ++y) + for (int32_t x = 0; x < reducedX; ++x) { + // construct quad for x edge + if (z > 0 && y > 0) { + // is edge intersected? + bool const entering = data[gA(x, y, z)] < iso && data[gA(x + 1, y, z)] >= iso; + bool const exiting = data[gA(x, y, z)] >= iso && data[gA(x + 1, y, z)] < iso; + if (entering || exiting) { + // generate quad + pointCode = getDualPointCode(x, y, z, iso, EDGE0); + calculateDualPoint(x, y, z, iso, pointCode, vertex0); + + pointCode = getDualPointCode(x, y, z - 1, iso, EDGE2); + calculateDualPoint(x, y, z - 1, iso, pointCode, vertex1); + + pointCode = getDualPointCode(x, y - 1, z - 1, iso, EDGE6); + calculateDualPoint(x, y - 1, z - 1, iso, pointCode, vertex2); + + pointCode = getDualPointCode(x, y - 1, z, iso, EDGE4); + calculateDualPoint(x, y - 1, z, iso, pointCode, vertex3); + + if (entering) { + local_vertices.emplace_back(vertex0); + local_vertices.emplace_back(vertex1); + local_vertices.emplace_back(vertex2); + local_vertices.emplace_back(vertex3); + } + else { + local_vertices.emplace_back(vertex0); + local_vertices.emplace_back(vertex3); + local_vertices.emplace_back(vertex2); + local_vertices.emplace_back(vertex1); + } + } + } + + // construct quad for y edge + if (z > 0 && x > 0) { + // is edge intersected? + bool const entering = data[gA(x, y, z)] < iso && data[gA(x, y + 1, z)] >= iso; + bool const exiting = data[gA(x, y, z)] >= iso && data[gA(x, y + 1, z)] < iso; + if (entering || exiting) { + // generate quad + pointCode = getDualPointCode(x, y, z, iso, EDGE8); + calculateDualPoint(x, y, z, iso, pointCode, vertex0); + + pointCode = getDualPointCode(x, y, z - 1, iso, EDGE11); + calculateDualPoint(x, y, z - 1, iso, pointCode, vertex1); + + pointCode = getDualPointCode(x - 1, y, z - 1, iso, EDGE10); + calculateDualPoint(x - 1, y, z - 1, iso, pointCode, vertex2); + + pointCode = getDualPointCode(x - 1, y, z, iso, EDGE9); + calculateDualPoint(x - 1, y, z, iso, pointCode, vertex3); + + if (exiting) { + local_vertices.emplace_back(vertex0); + local_vertices.emplace_back(vertex1); + local_vertices.emplace_back(vertex2); + local_vertices.emplace_back(vertex3); + } + else { + local_vertices.emplace_back(vertex0); + local_vertices.emplace_back(vertex3); + local_vertices.emplace_back(vertex2); + local_vertices.emplace_back(vertex1); + } + } + } + + // construct quad for z edge + if (x > 0 && y > 0) { + // is edge intersected? + bool const entering = data[gA(x, y, z)] < iso && data[gA(x, y, z + 1)] >= iso; + bool const exiting = data[gA(x, y, z)] >= iso && data[gA(x, y, z + 1)] < iso; + if (entering || exiting) { + // generate quad + pointCode = getDualPointCode(x, y, z, iso, EDGE3); + calculateDualPoint(x, y, z, iso, pointCode, vertex0); + + pointCode = getDualPointCode(x - 1, y, z, iso, EDGE1); + calculateDualPoint(x - 1, y, z, iso, pointCode, vertex1); + + pointCode = getDualPointCode(x - 1, y - 1, z, iso, EDGE5); + calculateDualPoint(x - 1, y - 1, z, iso, pointCode, vertex2); + + pointCode = getDualPointCode(x, y - 1, z, iso, EDGE7); + calculateDualPoint(x, y - 1, z, iso, pointCode, vertex3); + + if (exiting) { + local_vertices.emplace_back(vertex0); + local_vertices.emplace_back(vertex1); + local_vertices.emplace_back(vertex2); + local_vertices.emplace_back(vertex3); + } + else { + local_vertices.emplace_back(vertex0); + local_vertices.emplace_back(vertex3); + local_vertices.emplace_back(vertex2); + local_vertices.emplace_back(vertex1); + } + } + } + //++progress; + } + } +#ifdef USE_PARALLEL + { + MutexType::scoped_lock lock(writeMutex); + vertices.insert(vertices.end(), local_vertices.begin(), local_vertices.end()); + progress += (r.end() - r.begin()); + } + if (callback != NULL) { + int p = progress * 100 / (total_progress); + CallbackMutexType::scoped_lock lock(callbackMutex); + callback(p >= 100 ? 99 : p); + } + }, ap); +#else +#undef local_vertices + progress += (r.end() - r.begin()); + callback(progress * 100 / total_progress); +#endif + // generate triangle soup quads + size_t const numQuads = vertices.size() / 4; + quads.reserve(numQuads); + for (size_t i = 0; i < numQuads; ++i) { + quads.emplace_back(i * 4, i * 4 + 1, i * 4 + 2, i * 4 + 3); + } + if (callback != NULL) { + callback(100); + } + } + + /*/------------------------------------------------------------------------------ + + template inline + void DualMC::buildSharedVerticesQuads( + VolumeDataType const iso, + std::vector& vertices, + std::vector& quads + ) { + + + int32_t const reducedX = dims[0] - 2; + int32_t const reducedY = dims[1] - 2; + int32_t const reducedZ = dims[2] - 2; + + QuadIndexType i0, i1, i2, i3; + + //pointToIndex.clear(); + + //ProgressBar progress(reducedX * reducedY * reducedZ, PROGRESS_BAR_COLUMN); + + // iterate voxels + #ifndef USE_PARALLEL + for (int32_t z = 0; z < reducedZ; ++z) { + //#define local_quads quads + //#define local_vertices vertices + #else + tbb::spin_mutex writeMutex; + static tbb::affinity_partitioner ap; + tbb::parallel_for( + tbb::blocked_range(0, reducedZ), + [&](const tbb::blocked_range r) { + std::vector local_vertices; + std::vector local_quads; + PointToIndexMap pointToIndex; + for (int32_t z = r.begin(); z < r.end(); z++) { + #endif + for (int32_t y = 0; y < reducedY; ++y) + for (int32_t x = 0; x < reducedX; ++x) { + // construct quads for x edge + if (z > 0 && y > 0) { + bool const entering = data[gA(x, y, z)] < iso && data[gA(x + 1, y, z)] >= iso; + bool const exiting = data[gA(x, y, z)] >= iso && data[gA(x + 1, y, z)] < iso; + if (entering || exiting) { + // generate quad + i0 = getSharedDualPointIndex(x, y, z, iso, EDGE0, local_vertices, pointToIndex); + i1 = getSharedDualPointIndex(x, y, z - 1, iso, EDGE2, local_vertices, pointToIndex); + i2 = getSharedDualPointIndex(x, y - 1, z - 1, iso, EDGE6, local_vertices, pointToIndex); + i3 = getSharedDualPointIndex(x, y - 1, z, iso, EDGE4, local_vertices, pointToIndex); + + if (entering) { + local_quads.emplace_back(i0, i1, i2, i3); + } + else { + local_quads.emplace_back(i0, i3, i2, i1); + } + } + } + + // construct quads for y edge + if (z > 0 && x > 0) { + bool const entering = data[gA(x, y, z)] < iso && data[gA(x, y + 1, z)] >= iso; + bool const exiting = data[gA(x, y, z)] >= iso && data[gA(x, y + 1, z)] < iso; + if (entering || exiting) { + // generate quad + i0 = getSharedDualPointIndex(x, y, z, iso, EDGE8, local_vertices, pointToIndex); + i1 = getSharedDualPointIndex(x, y, z - 1, iso, EDGE11, local_vertices, pointToIndex); + i2 = getSharedDualPointIndex(x - 1, y, z - 1, iso, EDGE10, local_vertices, pointToIndex); + i3 = getSharedDualPointIndex(x - 1, y, z, iso, EDGE9, local_vertices, pointToIndex); + + if (exiting) { + local_quads.emplace_back(i0, i1, i2, i3); + } + else { + local_quads.emplace_back(i0, i3, i2, i1); + } + } + } + + // construct quads for z edge + if (x > 0 && y > 0) { + bool const entering = data[gA(x, y, z)] < iso && data[gA(x, y, z + 1)] >= iso; + bool const exiting = data[gA(x, y, z)] >= iso && data[gA(x, y, z + 1)] < iso; + if (entering || exiting) { + // generate quad + i0 = getSharedDualPointIndex(x, y, z, iso, EDGE3, local_vertices, pointToIndex); + i1 = getSharedDualPointIndex(x - 1, y, z, iso, EDGE1, local_vertices, pointToIndex); + i2 = getSharedDualPointIndex(x - 1, y - 1, z, iso, EDGE5, local_vertices, pointToIndex); + i3 = getSharedDualPointIndex(x, y - 1, z, iso, EDGE7, local_vertices, pointToIndex); + + if (exiting) { + local_quads.emplace_back(i0, i1, i2, i3); + } + else { + local_quads.emplace_back(i0, i3, i2, i1); + } + } + } + // progress bar + //++progress; + } + //progress.display(); + } + #ifdef USE_PARALLEL + { + tbb::spin_mutex::scoped_lock lock(writeMutex); + QuadIndexType last_index = vertices.size(); + for (Quad& q : local_quads) { + q.i0 += last_index; + q.i1 += last_index; + q.i2 += last_index; + q.i3 += last_index; + } + vertices.insert(vertices.end(), local_vertices.begin(), local_vertices.end()); + quads.insert(quads.end(), local_quads.begin(), local_quads.end()); + } + }, ap); + std::cout << vertices.size() << "," << quads.size() << std::endl; + #endif + //progress.done(); + }*/ #include "dualmc_tables.tpp" -} // END: namespace dualmc + } // END: namespace dualmc #endif // DUALMC_H_INCLUDED \ No newline at end of file diff --git a/include/dualmc/dualmc.tpp b/include/dualmc/dualmc.tpp deleted file mode 100644 index a6c614d..0000000 --- a/include/dualmc/dualmc.tpp +++ /dev/null @@ -1,474 +0,0 @@ -// Copyright (C) 2017, Dominik Wodniok -// This software may be modified and distributed under the terms -// of the BSD 3-Clause license. See the LICENSE.txt file for details. - -/// \file dualmc.tpp -/// \author Dominik Wodniok -/// \date 2009 - -//------------------------------------------------------------------------------ - -template inline -int DualMC::getCellCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso) const { - // determine for each cube corner if it is outside or inside - int code = 0; - if(data[gA(cx,cy,cz)] >= iso) - code |= 1; - if(data[gA(cx+1,cy,cz)] >= iso) - code |= 2; - if(data[gA(cx,cy+1,cz)] >= iso) - code |= 4; - if(data[gA(cx+1,cy+1,cz)] >= iso) - code |= 8; - if(data[gA(cx,cy,cz+1)] >= iso) - code |= 16; - if(data[gA(cx+1,cy,cz+1)] >= iso) - code |= 32; - if(data[gA(cx,cy+1,cz+1)] >= iso) - code |= 64; - if(data[gA(cx+1,cy+1,cz+1)] >= iso) - code |= 128; - return code; -} - -//------------------------------------------------------------------------------ - -template inline -int DualMC::getDualPointCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso, DMCEdgeCode const edge) const { - int cubeCode = getCellCode(cx, cy, cz, iso); - - // is manifold dual marching cubes desired? - if(generateManifold) { - // The Manifold Dual Marching Cubes approach from Rephael Wenger as described in - // chapter 3.3.5 of his book "Isosurfaces: Geometry, Topology, and Algorithms" - // is implemente here. - // If a problematic C16 or C19 configuration shares the ambiguous face - // with another C16 or C19 configuration we simply invert the cube code - // before looking up dual points. Doing this for these pairs ensures - // manifold meshes. - // But this removes the dualism to marching cubes. - - // check if we have a potentially problematic configuration - uint8_t const direction = problematicConfigs[uint8_t(cubeCode)]; - // If the direction code is in {0,...,5} we have a C16 or C19 configuration. - if(direction != 255) { - // We have to check the neighboring cube, which shares the ambiguous - // face. For this we decode the direction. This could also be done - // with another lookup table. - // copy current cube coordinates into an array. - int32_t neighborCoords[] = {cx,cy,cz}; - // get the dimension of the non-zero coordinate axis - unsigned int const component = direction >> 1; - // get the sign of the direction - int32_t delta = (direction & 1) == 1 ? 1 : -1; - // modify the correspong cube coordinate - neighborCoords[component] += delta; - // have we left the volume in this direction? - if(neighborCoords[component] >= 0 && neighborCoords[component] < (dims[component]-1)) { - // get the cube configuration of the relevant neighbor - int neighborCubeCode = getCellCode(neighborCoords[0], neighborCoords[1], neighborCoords[2], iso); - // Look up the neighbor configuration ambiguous face direction. - // If the direction is valid we have a C16 or C19 neighbor. - // As C16 and C19 have exactly one ambiguous face this face is - // guaranteed to be shared for the pair. - if(problematicConfigs[uint8_t(neighborCubeCode)] != 255) { - // replace the cube configuration with its inverse. - cubeCode ^= 0xff; - } - } - } - } - for(int i = 0; i < 4; ++i) - if(dualPointsList[cubeCode][i] & edge) { - return dualPointsList[cubeCode][i]; - } - return 0; -} - -//------------------------------------------------------------------------------ - -template inline -void DualMC::calculateDualPoint(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso, int const pointCode, Vertex & v) const { - // initialize the point with lower voxel coordinates - v.x = cx; - v.y = cy; - v.z = cz; - - // compute the dual point as the mean of the face vertices belonging to the - // original marching cubes face - Vertex p; - p.x=0; - p.y=0; - p.z=0; - int points = 0; - - // sum edge intersection vertices using the point code - if(pointCode & EDGE0) { - p.x += ((float)iso - (float)data[gA(cx,cy,cz)])/((float)data[gA(cx+1,cy,cz)]-(float)data[gA(cx,cy,cz)]); - points++; - } - - if(pointCode & EDGE1) { - p.x += 1.0f; - p.z += ((float)iso - (float)data[gA(cx+1,cy,cz)])/((float)data[gA(cx+1,cy,cz+1)]-(float)data[gA(cx+1,cy,cz)]); - points++; - } - - if(pointCode & EDGE2) { - p.x += ((float)iso - (float)data[gA(cx,cy,cz+1)])/((float)data[gA(cx+1,cy,cz+1)]-(float)data[gA(cx,cy,cz+1)]); - p.z += 1.0f; - points++; - } - - if(pointCode & EDGE3) { - p.z += ((float)iso - (float)data[gA(cx,cy,cz)])/((float)data[gA(cx,cy,cz+1)]-(float)data[gA(cx,cy,cz)]); - points++; - } - - if(pointCode & EDGE4) { - p.x += ((float)iso - (float)data[gA(cx,cy+1,cz)])/((float)data[gA(cx+1,cy+1,cz)]-(float)data[gA(cx,cy+1,cz)]); - p.y += 1.0f; - points++; - } - - if(pointCode & EDGE5) { - p.x += 1.0f; - p.z += ((float)iso - (float)data[gA(cx+1,cy+1,cz)])/((float)data[gA(cx+1,cy+1,cz+1)]-(float)data[gA(cx+1,cy+1,cz)]); - p.y += 1.0f; - points++; - } - - if(pointCode & EDGE6) { - p.x += ((float)iso - (float)data[gA(cx,cy+1,cz+1)])/((float)data[gA(cx+1,cy+1,cz+1)]-(float)data[gA(cx,cy+1,cz+1)]); - p.z += 1.0f; - p.y += 1.0f; - points++; - } - - if(pointCode & EDGE7) { - p.z += ((float)iso - (float)data[gA(cx,cy+1,cz)])/((float)data[gA(cx,cy+1,cz+1)]-(float)data[gA(cx,cy+1,cz)]); - p.y += 1.0f; - points++; - } - - if(pointCode & EDGE8) { - p.y += ((float)iso - (float)data[gA(cx,cy,cz)])/((float)data[gA(cx,cy+1,cz)]-(float)data[gA(cx,cy,cz)]); - points++; - } - - if(pointCode & EDGE9) { - p.x += 1.0f; - p.y += ((float)iso - (float)data[gA(cx+1,cy,cz)])/((float)data[gA(cx+1,cy+1,cz)]-(float)data[gA(cx+1,cy,cz)]); - points++; - } - - if(pointCode & EDGE10) { - p.x += 1.0f; - p.y += ((float)iso - (float)data[gA(cx+1,cy,cz+1)])/((float)data[gA(cx+1,cy+1,cz+1)]-(float)data[gA(cx+1,cy,cz+1)]); - p.z += 1.0f; - points++; - } - - if(pointCode & EDGE11) { - p.z += 1.0f; - p.y += ((float)iso - (float)data[gA(cx,cy,cz+1)])/((float)data[gA(cx,cy+1,cz+1)]-(float)data[gA(cx,cy,cz+1)]); - points++; - } - - // divide by number of accumulated points - float invPoints = 1.0f / (float)points; - p.x*= invPoints; - p.y*= invPoints; - p.z*= invPoints; - - // offset point by voxel coordinates - v.x += p.x; - v.y += p.y; - v.z += p.z; -} - -//------------------------------------------------------------------------------ - -template inline -QuadIndexType DualMC::getSharedDualPointIndex( - int32_t const cx, int32_t const cy, int32_t const cz, - VolumeDataType const iso, DMCEdgeCode const edge, - std::vector & vertices - ) { - // create a key for the dual point from its linearized cell ID and point code - DualPointKey key; - key.linearizedCellID = gA(cx,cy,cz); - key.pointCode = getDualPointCode(cx,cy,cz,iso,edge); - - // have we already computed the dual point? - auto iterator = pointToIndex.find(key); - if(iterator != pointToIndex.end()) { - // just return the dual point index - return iterator->second; - } else { - // create new vertex and vertex id - QuadIndexType newVertexId = vertices.size(); - vertices.emplace_back(); - calculateDualPoint(cx,cy,cz,iso,key.pointCode, vertices.back()); - // insert vertex ID into map and also return it - pointToIndex[key] = newVertexId; - return newVertexId; - } -} - -//------------------------------------------------------------------------------ - -template inline -void DualMC::build( - VolumeDataType const * data, - int32_t const dimX, int32_t const dimY, int32_t const dimZ, - VolumeDataType const iso, - bool const generateManifold, - bool const generateSoup, - std::vector & vertices, - std::vector & quads - ) { - - // set members - this->dims[0] = dimX; - this->dims[1] = dimY; - this->dims[2] = dimZ; - this->data = data; - this->generateManifold = generateManifold; - - // clear vertices and quad indices - vertices.clear(); - quads.clear(); - - // generate quad soup or shared vertices quad list - if(generateSoup) { - buildQuadSoup(iso,vertices,quads); - } else { - buildSharedVerticesQuads(iso,vertices,quads); - } -} - -//------------------------------------------------------------------------------ - -template inline -void DualMC::buildQuadSoup( - VolumeDataType const iso, - std::vector & vertices, - std::vector & quads - ) { - - int32_t const reducedX = dims[0] - 2; - int32_t const reducedY = dims[1] - 2; - int32_t const reducedZ = dims[2] - 2; - - Vertex vertex0; - Vertex vertex1; - Vertex vertex2; - Vertex vertex3; - int pointCode; - - ProgressBar progress(reducedX*reducedY*reducedZ, PROGRESS_BAR_COLUMN); - - // iterate voxels - for(int32_t z = 0; z < reducedZ; ++z) { - for(int32_t y = 0; y < reducedY; ++y) - for(int32_t x = 0; x < reducedX; ++x) { - // construct quad for x edge - if(z > 0 && y > 0) { - // is edge intersected? - bool const entering = data[gA(x,y,z)] < iso && data[gA(x+1,y,z)] >= iso; - bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x+1,y,z)] < iso; - if(entering || exiting){ - // generate quad - pointCode = getDualPointCode(x,y,z,iso,EDGE0); - calculateDualPoint(x,y,z,iso,pointCode, vertex0); - - pointCode = getDualPointCode(x,y,z-1,iso,EDGE2); - calculateDualPoint(x,y,z-1,iso,pointCode, vertex1); - - pointCode = getDualPointCode(x,y-1,z-1,iso,EDGE6); - calculateDualPoint(x,y-1,z-1,iso,pointCode, vertex2); - - pointCode = getDualPointCode(x,y-1,z,iso,EDGE4); - calculateDualPoint(x,y-1,z,iso,pointCode, vertex3); - - if(entering) { - vertices.emplace_back(vertex0); - vertices.emplace_back(vertex1); - vertices.emplace_back(vertex2); - vertices.emplace_back(vertex3); - } else { - vertices.emplace_back(vertex0); - vertices.emplace_back(vertex3); - vertices.emplace_back(vertex2); - vertices.emplace_back(vertex1); - } - } - } - - // construct quad for y edge - if(z > 0 && x > 0) { - // is edge intersected? - bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y+1,z)] >= iso; - bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y+1,z)] < iso; - if(entering || exiting){ - // generate quad - pointCode = getDualPointCode(x,y,z,iso,EDGE8); - calculateDualPoint(x,y,z,iso,pointCode, vertex0); - - pointCode = getDualPointCode(x,y,z-1,iso,EDGE11); - calculateDualPoint(x,y,z-1,iso,pointCode, vertex1); - - pointCode = getDualPointCode(x-1,y,z-1,iso,EDGE10); - calculateDualPoint(x-1,y,z-1,iso,pointCode, vertex2); - - pointCode = getDualPointCode(x-1,y,z,iso,EDGE9); - calculateDualPoint(x-1,y,z,iso,pointCode, vertex3); - - if(exiting) { - vertices.emplace_back(vertex0); - vertices.emplace_back(vertex1); - vertices.emplace_back(vertex2); - vertices.emplace_back(vertex3); - } else { - vertices.emplace_back(vertex0); - vertices.emplace_back(vertex3); - vertices.emplace_back(vertex2); - vertices.emplace_back(vertex1); - } - } - } - - // construct quad for z edge - if(x > 0 && y > 0) { - // is edge intersected? - bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y,z+1)] >= iso; - bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y,z+1)] < iso; - if(entering || exiting){ - // generate quad - pointCode = getDualPointCode(x,y,z,iso,EDGE3); - calculateDualPoint(x,y,z,iso,pointCode, vertex0); - - pointCode = getDualPointCode(x-1,y,z,iso,EDGE1); - calculateDualPoint(x-1,y,z,iso,pointCode, vertex1); - - pointCode = getDualPointCode(x-1,y-1,z,iso,EDGE5); - calculateDualPoint(x-1,y-1,z,iso,pointCode, vertex2); - - pointCode = getDualPointCode(x,y-1,z,iso,EDGE7); - calculateDualPoint(x,y-1,z,iso,pointCode, vertex3); - - if(exiting) { - vertices.emplace_back(vertex0); - vertices.emplace_back(vertex1); - vertices.emplace_back(vertex2); - vertices.emplace_back(vertex3); - } else { - vertices.emplace_back(vertex0); - vertices.emplace_back(vertex3); - vertices.emplace_back(vertex2); - vertices.emplace_back(vertex1); - } - } - } - ++progress; - } - progress.display(); - } - - // generate triangle soup quads - size_t const numQuads = vertices.size() / 4; - quads.reserve(numQuads); - for (size_t i = 0; i < numQuads; ++i) { - quads.emplace_back(i * 4, i * 4 + 1, i * 4 + 2, i * 4 + 3); - } - progress.done(); -} - -//------------------------------------------------------------------------------ - -template inline -void DualMC::buildSharedVerticesQuads( - VolumeDataType const iso, - std::vector & vertices, - std::vector & quads - ) { - - - int32_t const reducedX = dims[0] - 2; - int32_t const reducedY = dims[1] - 2; - int32_t const reducedZ = dims[2] - 2; - - QuadIndexType i0,i1,i2,i3; - - pointToIndex.clear(); - - ProgressBar progress(reducedX*reducedY*reducedZ, PROGRESS_BAR_COLUMN); - - // iterate voxels - for(int32_t z = 0; z < reducedZ; ++z) { - for(int32_t y = 0; y < reducedY; ++y) - for(int32_t x = 0; x < reducedX; ++x) { - // construct quads for x edge - if(z > 0 && y > 0) { - bool const entering = data[gA(x,y,z)] < iso && data[gA(x+1,y,z)] >= iso; - bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x+1,y,z)] < iso; - if(entering || exiting){ - // generate quad - i0 = getSharedDualPointIndex(x,y,z,iso,EDGE0,vertices); - i1 = getSharedDualPointIndex(x,y,z-1,iso,EDGE2,vertices); - i2 = getSharedDualPointIndex(x,y-1,z-1,iso,EDGE6,vertices); - i3 = getSharedDualPointIndex(x,y-1,z,iso,EDGE4,vertices); - - if(entering) { - quads.emplace_back(i0,i1,i2,i3); - } else { - quads.emplace_back(i0,i3,i2,i1); - } - } - } - - // construct quads for y edge - if(z > 0 && x > 0) { - bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y+1,z)] >= iso; - bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y+1,z)] < iso; - if(entering || exiting){ - // generate quad - i0 = getSharedDualPointIndex(x,y,z,iso,EDGE8,vertices); - i1 = getSharedDualPointIndex(x,y,z-1,iso,EDGE11,vertices); - i2 = getSharedDualPointIndex(x-1,y,z-1,iso,EDGE10,vertices); - i3 = getSharedDualPointIndex(x-1,y,z,iso,EDGE9,vertices); - - if(exiting) { - quads.emplace_back(i0,i1,i2,i3); - } else { - quads.emplace_back(i0,i3,i2,i1); - } - } - } - - // construct quads for z edge - if(x > 0 && y > 0) { - bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y,z+1)] >= iso; - bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y,z+1)] < iso; - if(entering || exiting){ - // generate quad - i0 = getSharedDualPointIndex(x,y,z,iso,EDGE3,vertices); - i1 = getSharedDualPointIndex(x-1,y,z,iso,EDGE1,vertices); - i2 = getSharedDualPointIndex(x-1,y-1,z,iso,EDGE5,vertices); - i3 = getSharedDualPointIndex(x,y-1,z,iso,EDGE7,vertices); - - if(exiting) { - quads.emplace_back(i0,i1,i2,i3); - } else { - quads.emplace_back(i0,i3,i2,i1); - } - } - } - // progress bar - ++progress; - } - progress.display(); - } - progress.done(); -} \ No newline at end of file diff --git a/include/implicit_function.h b/include/implicit_function.h index 4660b4a..06ce18b 100644 --- a/include/implicit_function.h +++ b/include/implicit_function.h @@ -9,7 +9,7 @@ #include #include #include -#include +#include "sol/sol.hpp" typedef double FT; const FT pi = 3.141592653589793238462643383279502884L; @@ -242,7 +242,7 @@ class Implicit_function : public Function { FT operator ()(FT x, FT y, FT z) { // Since we use FREP where F(x,y,z) >= 0 defined as a solid // then we inversed the implicit function to match FREP - // For example, f(x,y,z) = x^2+y^2+z^2 - 1, We want the solid region of f(x, y, z) <= 0 + // For example, f(x,y,z) = x^2+y^2+z^2 - 1, We want the solid region inside the circle where f(x, y, z) <= 0 // so that F(x,y,z) = -f(x,y,z) >= 0 if (thickness <= eps) { if (isLuaFunction) { @@ -299,4 +299,13 @@ Function* isosurface(std::string name, FT t) { } throw std::runtime_error("Implicit surface called `" + name + "` doesn't exist"); } + +void set_shorten_function(sol::state& lua) { + lua.script("abs, acos, asin, atan, atan2 = math.abs, math.acos, math.atan, math.atan2"); + lua.script("ceil, cos, deg, exp, floor = math.ceil, math.cos, math.deg, math.exp, math.floor"); + lua.script("log, log10, max, min, mod = math.log, math.log10, math.max, math.min, math.mod"); + lua.script("pow, rad, sin, sqrt, tan = math.pow, math.rad, math.sin, math.sqrt, math.tan"); + lua.script("frexp, ldexp, random, randomseed = math.frexp, math.ldexp, math.random, math.randomseed"); + lua.script("local pi = math.pi"); +} #endif \ No newline at end of file diff --git a/include/utils.h b/include/utils.h index 31f00d5..b33cb5e 100644 --- a/include/utils.h +++ b/include/utils.h @@ -1,4 +1,3 @@ -#pragma once #ifndef UTILS_INCLUDED #define UTILS_INCLUDED #include diff --git a/setup.cfg b/setup.cfg index d46e7d8..b2055ca 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,5 +1,4 @@ [flake8] # these error codes interfere with Black ignore = E203, E231, E501, W503, B950 -select = C,E,F,W,B,B9 -version = 0.0.1 \ No newline at end of file +select = C,E,F,W,B,B9 \ No newline at end of file diff --git a/setup.py b/setup.py index 3676c01..490ced9 100644 --- a/setup.py +++ b/setup.py @@ -109,7 +109,7 @@ def build_extension(self, ext): # logic and declaration, and simpler if you include description/version in a file. setup( name="PyScaffolder", - version="1.5.1.dev1", + version="1.5.2", author="Jirawat Iamsamang", author_email="nodtem66@gmail.com", description="Python wrapper for scaffolder program", diff --git a/src/Main.cpp b/src/Main.cpp index 8f80239..7d8ba1d 100644 --- a/src/Main.cpp +++ b/src/Main.cpp @@ -1,16 +1,14 @@ -#include "Scaffolder_2.h" +#include "Scaffolder.h" +#include "OptimalSlice.hpp" +#include "MeshOperation.h" +#include "utils.h" #include "QuadricSimplification.h" +#include "cxxopts.hpp" +#include "toojpeg/toojpeg.h" +#include "implicit_function.h" -ProgressBar qsim_progress(100, 40); -inline bool qsim_callback(int pos, const char* str) { - if (pos >= 0 && pos <= 100) { - qsim_progress.update(pos); - qsim_progress.display(); - } - if (pos >= 100) - qsim_progress.done(); - return true; -} +#define LOG if (log_format == SCAFFOLDER_FORMAT_DEFAULT) log +#define CSV if (log_format == SCAFFOLDER_FORMAT_CSV) log int main(int argc, char* argv[]) { @@ -29,7 +27,7 @@ int main(int argc, char* argv[]) uint8_t export_axis = 'X'; uint16_t grid_offset = 3; uint16_t shell = 0; - uint16_t smooth_step = 5; + uint16_t smooth_step = 0; uint16_t k_slice = 100; uint16_t k_polygon = 4; size_t minimum_grid_size = 100; @@ -41,21 +39,23 @@ int main(int argc, char* argv[]) double coff = pi; double minimum_diameter = 0.25; - // format Eigen object to readable text - Eigen::IOFormat CleanFmt(4, Eigen::DontAlignCols, ", ", "\n", "[", "]"); - Eigen::IOFormat CSVFmt(-1, Eigen::DontAlignCols, ", ", ", "); - // file parameters std::string filename = "", dir = ""; - std::string surface_name = "schwarzp"; + std::string surface_name = "bcc"; std::string input_file = ""; std::string output_format = ""; std::uint8_t log_format = SCAFFOLDER_FORMAT_DEFAULT; + // Eigen format + // format Eigen object to readable text + Eigen::IOFormat CleanFmt(4, Eigen::DontAlignCols, ", ", "\n", "[", "]"); + Eigen::IOFormat CSVFmt(-1, Eigen::DontAlignCols, ", ", ", "); + // TRY-CATCH block // Parse the program options try { - cxxopts::Options options("Scaffolder", "Scaffolder - generate 3D scaffold from STL file based on implicit surface"); + std::string description("Scaffolder (%VERSION%) - generate 3D scaffold from STL file based on implicit surface"); + cxxopts::Options options("Scaffolder", description.replace(description.find("%VERSION%"), sizeof("%VERSION%") - 1, VERSION)); options.positional_help("INPUT OUTPUT PARAMETERS").show_positional_help(); options.add_options() ("h,help", "Print help") @@ -65,24 +65,24 @@ int main(int argc, char* argv[]) ("q,quiet", "Disable verbose output [default: false]") ("c,coff", "Angular frequency (pore size adjustment) default:PI", cxxopts::value(), "DOUBLE") ("t,isolevel", "isolevel (porosity adjustment) [default: 0]", cxxopts::value(), "DOUBLE") - ("n,surface", "implicit surface: rectlinear, schwarzp, schwarzd, gyroid, double-p, double-d, double-gyroiod, lidinoid, schoen_iwp, neovius, bcc, tubular_g_ab, tubular_g_c [default: schwarzp]", cxxopts::value(), "NAME") + ("n,surface", "implicit surface: rectlinear, schwarzp, schwarzd, gyroid, double-p, double-d, double-gyroiod, lidinoid, schoen_iwp, neovius, bcc, tubular_g_ab, tubular_g_c (default: bcc)", cxxopts::value(), "NAME") ("g,grid_size", "Grid size [default: 100]", cxxopts::value(), "INT (0..60000)") ("s,shell", "Outer thickness (layers) [default:0]", cxxopts::value(), "INT (0..60000)") ("grid_offset", "[default:3]", cxxopts::value(), "INT (0..60000)") ("m,microstructure", "Analysis microstructure with Slice contour technique ( [default: false]") ("export_microstructure", "Analysis microstructure and export the 2D contours (for debugging) [default: false]") - ("export_jpeg", "Export 2D JPEG [Default: Z,100]", cxxopts::value>(), "[X|Y|Z],INT") - ("k_slice", "K_slice: the number of slicing layers in each direction (used in microstructure analysis) [default: 100]", cxxopts::value(), "INT (0..60000)") - ("k_polygon", "K_polygon: the number of closest outer contour (used in microstructure analysis) [default: 4]", cxxopts::value(), "INT (>0)") - ("z,size_optimize", "Experimental Quadric simplification [default: 0]", cxxopts::value(), "DOUBLE (0..1)") - ("smooth_step", "Smooth with laplacian (default: 5)", cxxopts::value(), "INT (0..60000)") - ("dirty", "Disable autoclean [default false]") - ("minimum_diameter", "used for removing small orphaned (between 0-1) [default: 0.25]", cxxopts::value(), "DOUBLE (0..1)") - ("format", "Format of logging output [default: default]", cxxopts::value(), "FORMAT (default, csv)") - ("output_inverse", "additional output inverse scaffold [default: false]") - ("fix_self_intersect", "Experimental fix self-intersect faces [default: 0]", cxxopts::value()->default_value("0"), "INT") + ("export_jpeg", "Export 2D JPEG (Default: Z,100)", cxxopts::value>(), "[X|Y|Z],INT") + ("k_slice", "K_slice: the number of slicing layers in each direction (used in microstructure analysis) (default: 100)", cxxopts::value(), "INT (0..60000)") + ("k_polygon", "K_polygon: the number of closest outer contour (used in microstructure analysis) (default: 4)", cxxopts::value(), "INT (>0)") + ("z,size_optimize", "Experimental Quadric simplification (default: 0)", cxxopts::value(), "DOUBLE (0..1)") + ("smooth_step", "Smooth with laplacian", cxxopts::value()->default_value("0"), "INT (0..60000)") + ("dirty", "Disable autoclean", cxxopts::value()->default_value("false"), "") + ("minimum_diameter", "used for removing small orphaned (between 0-1)", cxxopts::value()->default_value("0.25"), "DOUBLE (0..1)") + ("format", "Format of logging output", cxxopts::value()->default_value("default"), "FORMAT (default, csv)") + ("output_inverse", "additional output inverse scaffold", cxxopts::value()->default_value("false"), "") + ("fix_self_intersect", "Experimental fix self-intersect faces", cxxopts::value()->default_value("0"), "INT") ("mean_curvature", "Size of mean curvature histogram", cxxopts::value()->default_value("0"), "INT") - ("intersect", "Generate 3D mesh without intersect with original mesh [default: true]", cxxopts::value()->default_value("true"), ""); + ("no_intersect", "Generate 3D mesh without intersect with original mesh (false)", cxxopts::value(), ""); options.parse_positional({ "input", "output", "params" }); bool isEmptyOption = (argc == 1); cxxopts::ParseResult result = options.parse(argc, argv); @@ -128,7 +128,7 @@ int main(int argc, char* argv[]) } } if (result.count("dirty")) dirty = result["dirty"].as(); - if (result.count("intersect")) is_intersect = result["intersect"].as(); + if (result.count("no_intersect")) is_intersect = !result["no_intersect"].as(); if (result.count("microstructure")) { is_analysis_microstructure = result["microstructure"].as(); no_output = true; @@ -250,7 +250,7 @@ int main(int argc, char* argv[]) << "-- Mean curvature: " << (is_mean_curvature ? "True" : "False") << std::endl << "-- Export JPEG: " << (is_export_jpeg ? "Yes" : "No") << std::endl << "-- Axis: " << export_axis << std::endl - << "-- Build: " << (no_output ? "No" : "Yes") << (is_build_inverse ? " (Inverse)" : "") << (is_intersect ? "(Intersect)" : "") << std::endl; + << "-- Build: " << (no_output ? "No" : "Yes") << (is_build_inverse ? " (Inverse)" : "") << (is_intersect ? " (Intersect)" : "") << std::endl; else CSV << "Surface,coff,shell,thickness,grid_size,grid_offset,smooth_step,input_file,"; @@ -538,9 +538,7 @@ int main(int argc, char* argv[]) if (qsim_percent > 0) { LOG << "[Quadric Simplificaion: " << qsim_percent << "]" << std::endl; - qsim_progress.reset(); vcg::tri::mesh_quad_simplification(mesh, qsim_percent, qsim_callback); - qsim_progress.done(); LOG << "-- Info: " << mesh.VN() << " vertices " << mesh.FN() << " faces" << std::endl; } @@ -749,7 +747,12 @@ int main(int argc, char* argv[]) << mesh.VN() << ',' << mesh.FN() << std::endl; } else { - LOG << "[Warning] The scaffold isn't a manifold. The grid_offset should be increased or enable --fix_self_intersect option" << std::endl; + LOG << "[Warning] The scaffold isn't a manifold, so the suggestions are here:\n" + << "* Increase --grid_offset \n" + << "* Enable --fix_self_intersect option \n" + << "* Adjust --smooth_step \n" + << "* OR use the external 3D mesh editor to fix the problem, which I hope you can do it" + << std::endl; else CSV << ",,,," << mesh.VN() << ',' << mesh.FN() << std::endl; } diff --git a/src/Scaffolder.cpp b/src/Scaffolder.cpp new file mode 100644 index 0000000..cc7450d --- /dev/null +++ b/src/Scaffolder.cpp @@ -0,0 +1,111 @@ +#include "Scaffolder.h" +#include "dualmc/dualmc.h" + +ProgressBar qsim_progress(100, 40); +bool qsim_callback(int pos, const char* str) { + if (pos >= 0 && pos <= 100) { + qsim_progress.update(pos); + qsim_progress.display(); + } + if (pos >= 100) + qsim_progress.done(); + return true; +} + +ProgressBar marching_progress(100, 40); +void marching_callback(int pos) { + if (pos >= 0 && pos <= 100) { + marching_progress.update(pos); + marching_progress.display(); + } + if (pos >= 100) + marching_progress.done(); +} + +// Flatten between 1D and 3D +// https://stackoverflow.com/questions/7367770/how-to-flatten-or-index-3d-array-in-1d-array +inline size_t indexFromIJK(size_t i, size_t j, size_t k, Eigen::RowVector3i grid_size) { + return i + grid_size(0) * (j + grid_size(1) * k); +} + +inline void indexToIJK(size_t index, Eigen::RowVector3i grid_size, index_type& r) { + r.z = index / grid_size(0) / grid_size(1); + index -= r.z * grid_size(0) * grid_size(1); + r.y = index / grid_size(0); + r.x = index % grid_size(0); +} + +inline bool MarkAndSweepNeighbor(Eigen::VectorXd& W, index_type& index, Queue_t& queue, Eigen::RowVector3i grid_size, double value, bool findAbove) { + bool isBorder = false; + for (int8_t di = -1; di <= 1; di++) { + for (int8_t dj = -1; dj <= 1; dj++) { + for (int8_t dk = -1; dk <= 1; dk++) { + if (di == 0 && dj == 0 && dk == 0) continue; + const size_t id = indexFromIJK(index.x + di, index.y + dj, index.z + dk, grid_size); + //std::cout << value << " " << W(id) << std::endl; + if ((findAbove && W(id) >= value - EPSIL) || (!findAbove && W(id) <= value + EPSIL)) { + isBorder = true; + break; + } + } + if (isBorder) break; + } + if (isBorder) break; + } + if (isBorder) { + for (int8_t di = -1; di <= 1; di++) { + for (int8_t dj = -1; dj <= 1; dj++) { + for (int8_t dk = -1; dk <= 1; dk++) { + if (di == 0 && dj == 0 && dk == 0) continue; + const size_t id = indexFromIJK(index.x + di, index.y + dj, index.z + dk, grid_size); + if (W(id) >= 0.5 && W(id) < 1.1) { + queue.insert({ id, true }); + } + } + } + } + } + return isBorder; +} + +void marching_cube(TMesh& mesh, Eigen::MatrixXd& Fxyz, Eigen::RowVector3i grid_size, Eigen::RowVector3d& Vmin, double delta, bool verbose, bool dirty) { + + dualmc::DualMC builder; + std::vector mc_vertices; + std::vector mc_quads; + if (verbose) { + std::cout << "[Marching Cube] " << std::endl; + builder.callback = marching_callback; + } + builder.build((double const*)Fxyz.data(), grid_size(0), grid_size(1), grid_size(2), 0, true, true, mc_vertices, mc_quads); + TMesh::VertexIterator vi = vcg::tri::Allocator::AddVertices(mesh, mc_vertices.size()); + TMesh::FaceIterator fi = vcg::tri::Allocator::AddFaces(mesh, mc_quads.size() * 2); + std::vector vp(mc_vertices.size()); + for (size_t i = 0, len = mc_vertices.size(); i < len; i++, ++vi) { + vp[i] = &(*vi); + vi->P() = TMesh::CoordType( + Vmin(0) + mc_vertices[i].x * delta, + Vmin(1) + mc_vertices[i].y * delta, + Vmin(2) + mc_vertices[i].z * delta + ); + } + for (size_t i = 0, len = mc_quads.size(); i < len; i++, ++fi) { + fi->V(0) = vp[mc_quads[i].i0]; + fi->V(1) = vp[mc_quads[i].i1]; + fi->V(2) = vp[mc_quads[i].i2]; + ++fi; + fi->V(0) = vp[mc_quads[i].i2]; + fi->V(1) = vp[mc_quads[i].i3]; + fi->V(2) = vp[mc_quads[i].i0]; + } + if (!dirty) { + vcg::tri::Clean::RemoveDuplicateFace(mesh); + vcg::tri::Clean::RemoveDuplicateVertex(mesh); + vcg::tri::Clean::RemoveUnreferencedVertex(mesh); + } + +} + +bool null_callback(int pos, const char* str) { + return true; +} \ No newline at end of file diff --git a/src/SliceTest.cpp b/src/SliceTest.cpp index 9c3fd40..c0f25d5 100644 --- a/src/SliceTest.cpp +++ b/src/SliceTest.cpp @@ -10,7 +10,6 @@ #include "MeshOperation.h" ProgressBar import_progress(100, 40); -ProgressBar qsim_progress(100, 40); int _pos = 0; bool import_callback(int pos, const char* str) { diff --git a/src/python/PyScaffolder.cpp b/src/python/PyScaffolder.cpp index 7f20b11..2fcac5f 100644 --- a/src/python/PyScaffolder.cpp +++ b/src/python/PyScaffolder.cpp @@ -5,6 +5,8 @@ namespace py = pybind11; PYBIND11_MODULE(PyScaffolder, m) { m.doc() = "PyScaffolder generate isosurface from implicit function"; + m.attr("__version__") = VERSION; + py::class_(m, "PoreSize", py::dynamic_attr()) .def(py::init<>()) .def_readwrite("minFeret", &PyScaffolder::PoreSize::minFeret) @@ -22,6 +24,7 @@ PYBIND11_MODULE(PyScaffolder, m) { .def(py::init<>()) .def_readwrite("is_build_inverse", &PyScaffolder::Parameter::is_build_inverse) .def_readwrite("is_intersect", &PyScaffolder::Parameter::is_intersect) + .def_readwrite("verbose", &PyScaffolder::Parameter::verbose) .def_readwrite("coff", &PyScaffolder::Parameter::coff) .def_readwrite("grid_offset", &PyScaffolder::Parameter::grid_offset) .def_readwrite("grid_size", &PyScaffolder::Parameter::grid_size) @@ -34,7 +37,7 @@ PYBIND11_MODULE(PyScaffolder, m) { .def_readwrite("fix_self_intersect", &PyScaffolder::Parameter::fix_self_intersect) .def_readwrite("surface_name", &PyScaffolder::Parameter::surface_name); - m.def("slice_test", &PyScaffolder::slice_test, "A function to slice input mesh into pore sizes", + m.def("slice_test", &PyScaffolder::slice_test, py::call_guard(), "A function to slice input mesh into pore sizes", py::arg("vertices"), py::arg("faces"), py::arg("k_slice") = 100, @@ -46,10 +49,19 @@ PYBIND11_MODULE(PyScaffolder, m) { PyScaffolder::Parameter default_parameters; - m.def("generate_mesh", &PyScaffolder::generate_mesh, "A function to generate isosurface from input mesh and parameters", + m.def("generate_scaffold", &PyScaffolder::generate_scaffold, py::call_guard(), "A function to generate isosurface from input mesh and parameters", py::arg("vertices"), py::arg("faces"), py::arg("params") = default_parameters, py::arg("callback") = py::none() ); + + m.def("marching_cubes", &PyScaffolder::marching_cubes, py::call_guard(), "A function to generate a triangular mesh (v, f) from isovalues (f)", + py::arg("f"), + py::arg("grid_size") = std::tuple(100, 100, 100), + py::arg("v_min") = std::tuple(0, 0, 0), + py::arg("delta") = 0.01, + py::arg("clean") = false, + py::arg("callback") = py::none() + ); } \ No newline at end of file diff --git a/src/python/PyScaffolder.hpp b/src/python/PyScaffolder.hpp index 6c17d21..8c3061e 100644 --- a/src/python/PyScaffolder.hpp +++ b/src/python/PyScaffolder.hpp @@ -1,5 +1,3 @@ -#pragma once - #ifndef PYSCAFFOLDER_INCLUDED_H #define PYSCAFFOLDER_INCLUDED_H @@ -10,6 +8,9 @@ #include "pybind11/stl_bind.h" #include "pybind11/eigen.h" #include "pybind11/functional.h" +#include "Scaffolder.h" + +namespace py = pybind11; namespace PyScaffolder { @@ -30,9 +31,10 @@ namespace PyScaffolder { struct Parameter { bool is_build_inverse = false; bool is_intersect = true; + bool verbose = false; uint16_t shell = 0; - uint16_t grid_offset = 5; - uint16_t smooth_step = 5; + uint16_t grid_offset = 3; + uint16_t smooth_step = 0; uint16_t k_slice = 100; uint16_t k_polygon = 4; uint16_t fix_self_intersect = 0; @@ -53,11 +55,20 @@ namespace PyScaffolder { const std::function& callback = NULL ); - MeshInfo generate_mesh( + MeshInfo generate_scaffold( Eigen::MatrixXd v, Eigen::MatrixXi f, Parameter params, const std::function& callback = NULL ); + + std::tuple marching_cubes( + Eigen::VectorXd& Fxyz, + py::object& grid_size, + py::object& Vmin, + py::object& delta, + bool clean = false, + const std::function& callback = NULL + ); } #endif \ No newline at end of file diff --git a/src/python/implements.cpp b/src/python/implements.cpp index 3ab4444..5320d84 100644 --- a/src/python/implements.cpp +++ b/src/python/implements.cpp @@ -1,15 +1,36 @@ #include "PyScaffolder.hpp" -#include "Scaffolder_2.h" -#undef LOG +#include "OptimalSlice.hpp" +#include "MeshOperation.h" +#include "implicit_function.h" +#include "QuadricSimplification.h" +#include "dualmc/dualmc.h" + #define LOG std::cout +#define IS_INSTANCE_ARRAYLIST(x) (py::isinstance(x) || py::isinstance(x) || py::isinstance(x)) using namespace PyScaffolder; +using namespace optimal_slice; +namespace py = pybind11; // format Eigen object to readable text Eigen::IOFormat CleanFmt(4, Eigen::DontAlignCols, ", ", "\n", "[", "]"); Eigen::IOFormat CSVFmt(-1, Eigen::DontAlignCols, ", ", ", "); -MeshInfo PyScaffolder::generate_mesh( +void clean_mesh(TMesh &mesh) { + // It's a good practice to clean STL mesh firstly (VCGLib) + mesh.face.EnableFFAdjacency(); + vcg::tri::Clean::RemoveDuplicateFace(mesh); + vcg::tri::Clean::RemoveDuplicateVertex(mesh); + vcg::tri::Clean::RemoveUnreferencedVertex(mesh); + vcg::tri::UpdateTopology::FaceFace(mesh); + vcg::tri::Clean::RemoveZeroAreaFace(mesh); + vcg::tri::Clean::RemoveNonManifoldFace(mesh); + vcg::tri::Clean::RemoveUnreferencedVertex(mesh); + vcg::tri::UpdateTopology::FaceFace(mesh); + mesh.face.DisableFFAdjacency(); +} + +MeshInfo PyScaffolder::generate_scaffold( Eigen::MatrixXd V1, Eigen::MatrixXi F1, Parameter params, @@ -17,8 +38,8 @@ MeshInfo PyScaffolder::generate_mesh( ) { MeshInfo mesh_info; TMesh mesh; - double volume1 = eps, volume2 = eps; - double area1 = eps, area2 = eps; + double volume1 = EPSIL, volume2 = EPSIL; + double area1 = EPSIL, area2 = EPSIL; { Eigen::MatrixXd V, Fxyz, IFxyz; Eigen::MatrixXi F; @@ -69,15 +90,16 @@ MeshInfo PyScaffolder::generate_mesh( // Create border offset from the original box V1min1 = V1min - params.grid_offset * grid_delta * Eigen::RowVector3d::Ones(); grid_size = (L / grid_delta).cast() + 2 * params.grid_offset * Eigen::RowVector3i::Ones(); - LOG << "-- Bounding Box: " - << V1min.format(CleanFmt) << ' ' << V1max.format(CleanFmt) << std::endl - << "-- Length: " << L.format(CleanFmt) << std::endl - << "-- Grid delta: " << grid_delta << std::endl; + if (params.verbose) + LOG << "-- Bounding Box: " + << V1min.format(CleanFmt) << ' ' << V1max.format(CleanFmt) << std::endl + << "-- Length: " << L.format(CleanFmt) << std::endl + << "-- Grid delta: " << grid_delta << std::endl; // Convert input STL from VCG::TMesh to Eigen matrixXd V1 mesh_to_eigen_vector(stl, V1, F1); } // Generate Grid index (x,y,z) - LOG << "[Generating grid] "; + if (params.verbose) LOG << "[Generating grid] "; Eigen::MatrixXd GV(grid_size.prod(), 3); for (size_t k = 0; k < grid_size(2); k++) { const double z = V1min1(2) + k * grid_delta; @@ -90,8 +112,10 @@ MeshInfo PyScaffolder::generate_mesh( } } } - LOG << "OK" << std::endl - << "-- Grid size: " << grid_size.prod() << " " << grid_size.format(CleanFmt) << std::endl; + if (params.verbose) + LOG << "OK" << std::endl + << "-- Grid size: " << grid_size.prod() << " " << grid_size.format(CleanFmt) << std::endl; + if (callback != NULL) callback(1); Eigen::VectorXd W, D, Signed_Distance; @@ -104,9 +128,10 @@ MeshInfo PyScaffolder::generate_mesh( igl::signed_distance_fast_winding_number(GV, V1, F1, tree, bvh, Signed_Distance); Signed_Distance *= -1; } - LOG << "OK" << std::endl - << "-- Sign Distance: [ " << Signed_Distance.minCoeff() << ", " << Signed_Distance.maxCoeff() << "] " + if (params.verbose) + LOG << "-- Sign Distance: [ " << Signed_Distance.minCoeff() << ", " << Signed_Distance.maxCoeff() << "] " << "Wind: [ " << W.minCoeff() << ", " << W.maxCoeff() << "]" << std::endl; + if (callback != NULL) callback(10); sol::state lua; @@ -238,7 +263,8 @@ MeshInfo PyScaffolder::generate_mesh( is_manifold = is_manifold && fix_non_manifold_edges(mesh, params.minimum_diameter, 5); } if (!is_manifold) { - std::cout << "[Warning] Mesh is not two manifold" << std::endl; + if (params.verbose) + LOG << "[Warning] Mesh is not two manifold" << std::endl; } else { area2 = vcg::tri::Stat::ComputeMeshArea(mesh); @@ -248,24 +274,94 @@ MeshInfo PyScaffolder::generate_mesh( mesh_info.surface_area_ratio = abs(area2 / area1); } } - // It's a good practice to clean STL mesh firstly (VCGLib) - mesh.face.EnableFFAdjacency(); - vcg::tri::Clean::RemoveDuplicateFace(mesh); - vcg::tri::Clean::RemoveDuplicateVertex(mesh); - vcg::tri::Clean::RemoveUnreferencedVertex(mesh); - vcg::tri::UpdateTopology::FaceFace(mesh); - vcg::tri::Clean::RemoveZeroAreaFace(mesh); - vcg::tri::Clean::RemoveNonManifoldFace(mesh); - vcg::tri::Clean::RemoveUnreferencedVertex(mesh); - vcg::tri::UpdateTopology::FaceFace(mesh); - mesh.face.DisableFFAdjacency(); + clean_mesh(mesh); mesh_to_eigen_vector(mesh, mesh_info.v, mesh_info.f); if (callback != NULL) callback(100); return mesh_info; } -using namespace optimal_slice; +std::tuple PyScaffolder::marching_cubes( + Eigen::VectorXd& Fxyz, + py::object& grid_size, + py::object& _Vmin, + py::object& _delta, + bool clean, + const std::function& callback +) { + dualmc::DualMC builder; + std::vector mc_vertices; + std::vector mc_quads; + if (callback != NULL) builder.callback = callback; + + // Type conversion + std::array g; + if (IS_INSTANCE_ARRAYLIST(grid_size)) { + g = py::cast< std::array >(grid_size); + } + else { + int32_t gs = py::cast(grid_size); + g[0] = gs; + g[1] = gs; + g[2] = gs; + } + std::array delta; + if (IS_INSTANCE_ARRAYLIST(_delta)) { + delta = py::cast< std::array >(_delta); + } + else { + double d = py::cast(_delta); + delta[0] = d; + delta[1] = d; + delta[2] = d; + } + std::array Vmin = py::cast>(_Vmin); + + // Dual-Marching cubes + builder.build((double const*)Fxyz.data(), g[0], g[1], g[2], 0, true, true, mc_vertices, mc_quads); + Eigen::MatrixXd v; + Eigen::MatrixXi f; + if (clean) { + TMesh mesh; + TMesh::VertexIterator vi = vcg::tri::Allocator::AddVertices(mesh, mc_vertices.size()); + TMesh::FaceIterator fi = vcg::tri::Allocator::AddFaces(mesh, mc_quads.size() * 2); + std::vector vp(mc_vertices.size()); + for (size_t i = 0, len = mc_vertices.size(); i < len; i++, ++vi) { + vp[i] = &(*vi); + vi->P() = TMesh::CoordType( + Vmin[0] + mc_vertices[i].x * delta[0], + Vmin[1] + mc_vertices[i].y * delta[1], + Vmin[2] + mc_vertices[i].z * delta[2] + ); + } + for (size_t i = 0, len = mc_quads.size(); i < len; i++, ++fi) { + fi->V(0) = vp[mc_quads[i].i0]; + fi->V(1) = vp[mc_quads[i].i1]; + fi->V(2) = vp[mc_quads[i].i2]; + ++fi; + fi->V(0) = vp[mc_quads[i].i2]; + fi->V(1) = vp[mc_quads[i].i3]; + fi->V(2) = vp[mc_quads[i].i0]; + } + clean_mesh(mesh); + mesh_to_eigen_vector(mesh, v, f); + } + else { + v.resize(mc_vertices.size(), 3); + for (size_t i = 0, len = mc_vertices.size(); i < len; i++) { + v.row(i) << Vmin[0] + mc_vertices[i].x * delta[0], Vmin[1] + mc_vertices[i].y * delta[1], Vmin[2] + mc_vertices[i].z * delta[2]; + } + f.resize(mc_quads.size()*2, 3); + for (size_t i = 0, j = 0, len = mc_quads.size(); i < len; i++) { + f.row(j) << mc_quads[i].i0, mc_quads[i].i1, mc_quads[i].i2; + j++; + f.row(j) << mc_quads[i].i2, mc_quads[i].i3, mc_quads[i].i0; + j++; + } + } + + return make_tuple(v, f); +} PoreSize PyScaffolder::slice_test(Eigen::MatrixXd v, Eigen::MatrixXi f, size_t k_slice, size_t k_polygon, int direction, const std::function& callback) { PoreSize pore_size; diff --git a/src/python/test.py b/src/python/test.py index e684c64..56316bf 100644 --- a/src/python/test.py +++ b/src/python/test.py @@ -1,42 +1,106 @@ import PyScaffolder +import unittest -# Load STL mesh with igl -v = [ - [0.0, 0.0, 0.0], - [20.0, 0.0, 0.0], - [0.0, -20.0, 0.0], - [20.0, -20.0, 0.0], - [0.0, 0.0, 20.0], - [20.0, 0.0, 20.0], - [0.0, -20.0, 20.0], - [20.0, -20.0, 20.0] -] - -f = [ - [1, 2, 0], - [2, 1, 3], - [7, 2, 3], - [2, 7, 6], - [1, 7, 3], - [7, 1, 5], - [7, 4, 6], - [4, 7, 5], - [4, 2, 6], - [2, 4, 0], - [4, 1, 0], - [1, 4, 5] -] - -# Test Slice test with all direction -a = PyScaffolder.slice_test(v, f, direction=3) -print("Vertices: ", len(v)) -print("Faces: ", len(f)) -print(len(a.minFeret)) - -# Test generate surface with default parameter -params = PyScaffolder.Parameter() -params.coff = 1.0 -a = PyScaffolder.generate_mesh(v, f, params) -print(a.porosity) -print(a.surface_area_ratio) -print(len(a.v), len(a.f)) \ No newline at end of file +from math import sqrt + +class TestPyScaffolder(unittest.TestCase): + + v = [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, -1.0, 0.0], + [1.0, -1.0, 0.0], + [0.0, 0.0, 1.0], + [1.0, 0.0, 1.0], + [0.0, -1.0, 1.0], + [1.0, -1.0, 1.0] + ] + + f = [ + [1, 2, 0], + [2, 1, 3], + [7, 2, 3], + [2, 7, 6], + [1, 7, 3], + [7, 1, 5], + [7, 4, 6], + [4, 7, 5], + [4, 2, 6], + [2, 4, 0], + [4, 1, 0], + [1, 4, 5] + ] + + counter = 0 + + def _callback(self, p): + self.counter = p + + def test_version(self): + self.assertIsNotNone(PyScaffolder.__version__) + + def test_slicer(self): + a = PyScaffolder.slice_test(self.v, self.f, direction=3) + self.assertEqual(len(self.v), 8) + self.assertEqual(len(self.f), 12) + self.assertEqual(len(a.minFeret), 0) + self.assertEqual(len(a.maxFeret), 0) + + def test_slicer_with_callback(self): + a = PyScaffolder.slice_test(self.v, self.f, direction=3, callback=self._callback) + self.assertEqual(self.counter, 100) + + def test_scaffolder(self): + # Test generate surface with default parameter + params = PyScaffolder.Parameter() + params.coff = 12.0 + params.verbose = False + a = PyScaffolder.generate_scaffold(self.v, self.f, params) + self.assertAlmostEqual(a.porosity, 0.453, places=2) + self.assertAlmostEqual(a.surface_area_ratio, 0.912, places=2) + self.assertGreaterEqual(len(a.v), 7.2e4) + self.assertEqual(len(a.v[0]), 3) + self.assertGreaterEqual(len(a.f), 1.4e5) + self.assertEqual(len(a.f[0]), 3) + + def test_scaffolder_with_callback(self): + params = PyScaffolder.Parameter() + params.coff = 12.0 + params.verbose = False + self.progress = 0 + def callback(p): + self.progress = p + PyScaffolder.generate_scaffold(self.v, self.f, params, callback=callback) + self.assertEqual(self.progress, 100) + + + def test_marching_cubes(self): + Fxyz = [] + for i in range(100): + for j in range(100): + for k in range(100): + Fxyz.append(sqrt((i-50)**2+(j-50)**2+(k-50)**2) - 2**2) + self.assertEqual(len(Fxyz), 1e6) + import numpy as np + (v, f) = PyScaffolder.marching_cubes(Fxyz, grid_size=(100, 100, 100), delta=0.02, v_min=(-1, -1, -1), clean=False) + self.assertGreaterEqual(len(v), 1080) + self.assertEqual(len(v[0]), 3) + self.assertGreaterEqual(len(f), 540) + self.assertEqual(len(f[0]), 3) + + def test_marching_cubes_with_callback(self): + Fxyz = [ + 1,1,1,1, + 1,-1,-1,1, + 1,-1,-1,1, + 1,1,1,1 + ]*4 + self.counter = 0 + (v, f) = PyScaffolder.marching_cubes(Fxyz, grid_size=4, delta=0.25, v_min=(-.5, -.5, -.5), clean=True, callback=self._callback) + self.assertEqual(len(v), 6) + self.assertEqual(len(f), 4) + self.assertEqual(self.counter, 100) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file