diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..1377554ebea6f98a2c748183bc5a96852af12ac2
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+*.swp
diff --git a/Examples/CreateEnclosingMesh/CreateEnclosingMesh.cpp b/Examples/CreateEnclosingMesh/CreateEnclosingMesh.cpp
index 9115c3ae57233230129c838c38bf4f9cecd4a4de..251f6b9169450c7deebd4df8458ad56669bd64b9 100644
--- a/Examples/CreateEnclosingMesh/CreateEnclosingMesh.cpp
+++ b/Examples/CreateEnclosingMesh/CreateEnclosingMesh.cpp
@@ -25,6 +25,7 @@
#include "imstkTetrahedralMesh.h"
#include "imstkMeshIO.h"
#include "imstkVTKMeshIO.h"
+#include "imstkGeometryUtilities.h"
using namespace imstk;
@@ -56,7 +57,7 @@ main()
// add the scene object to the scene
scene->addSceneObject(surfaceObject);
- auto tetMesh = TetrahedralMesh::createTetrahedralMeshCover(surfMesh, nx, ny, nz);
+ auto tetMesh = GeometryUtils::createTetrahedralMeshCover(*surfMesh, nx, ny, nz);
// add scene object for surface object
auto volObject = std::make_shared<VisualObject>("VolObj");
@@ -74,4 +75,4 @@ main()
simManager->setActiveScene(scene);
simManager->getViewer()->setBackgroundColors(Vec3d(0.3285, 0.3285, 0.6525), Vec3d(0.13836, 0.13836, 0.2748), true);
simManager->start();
-}
\ No newline at end of file
+}
diff --git a/Examples/RCM/CMakeLists.txt b/Examples/RCM/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..eac125b3b7c5653d399ad604dc950250c36ff47e
--- /dev/null
+++ b/Examples/RCM/CMakeLists.txt
@@ -0,0 +1,43 @@
+###########################################################################
+#
+# Copyright (c) Kitware, Inc.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0.txt
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+###########################################################################
+
+project(Example-RCMDragon)
+
+#-----------------------------------------------------------------------------
+# Create executable
+#-----------------------------------------------------------------------------
+imstk_add_executable(${PROJECT_NAME} RCM.cpp)
+
+#-----------------------------------------------------------------------------
+# Add the target to Examples folder
+#-----------------------------------------------------------------------------
+SET_TARGET_PROPERTIES (${PROJECT_NAME} PROPERTIES FOLDER Examples)
+
+#-----------------------------------------------------------------------------
+# Link libraries to executable
+#-----------------------------------------------------------------------------
+target_link_libraries(${PROJECT_NAME} apiUtilities)
+
+#-----------------------------------------------------------------------------
+# Set MSVC working directory to the install/bin directory
+#-----------------------------------------------------------------------------
+if(MSVC) # Configure running executable out of MSVC
+ set_property(TARGET ${PROJECT_NAME} PROPERTY VS_DEBUGGER_WORKING_DIRECTORY "${iMSTK_INSTALL_BIN_DIR}")
+endif()
+
+
diff --git a/Examples/RCM/RCM.cpp b/Examples/RCM/RCM.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1079abacae6e6db2378de054bbfc5b60273f8fa8
--- /dev/null
+++ b/Examples/RCM/RCM.cpp
@@ -0,0 +1,141 @@
+/*=========================================================================
+
+ Library: iMSTK
+
+ Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+ & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+ Licensed under the Apache License, Version 2.0 (the "License");
+ you may not use this file except in compliance with the License.
+ You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0.txt
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+
+=========================================================================*/
+
+#include "imstkMeshIO.h"
+#include "imstkTetrahedralMesh.h"
+#include "imstkGeometryUtilities.h"
+#include "imstkLogUtility.h"
+#include "bandwidth.h"
+#include <chrono>
+#include <thread>
+
+using namespace imstk;
+
+using QuadConn = std::array<size_t, 4>;
+
+///
+/// \brief create a quad mesh
+/// \retval pair of connectivity and num of vertices
+///
+std::pair<std::vector<QuadConn>, size_t> createConn();
+
+template<typename ElemConn>
+void testRCM(const std::vector<ElemConn>& conn, const size_t numVerts);
+
+int
+main(int argc, char** argv)
+{
+ auto logUtil = std::make_shared<LogUtility>();
+ logUtil->createLogger("simulation", "./");
+
+ // a 2D Cartesian mesh
+ {
+ auto p = createConn();
+ testRCM(p.first, p.second);
+ }
+
+ // 3D mesh
+ {
+ auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(MeshIO::read(iMSTK_DATA_ROOT "/asianDragon/asianDragon.veg"));
+ const auto numVerts = tetMesh->getNumVertices();
+ std::cout << "Number of vertices = " << numVerts << std::endl;
+ testRCM(tetMesh->getTetrahedraVertices(), numVerts);
+ }
+
+ // a surface mesh cover
+ {
+ auto surfMesh = std::dynamic_pointer_cast<SurfaceMesh>(MeshIO::read(iMSTK_DATA_ROOT "/asianDragon/asianDragon.obj"));
+ auto tetMesh = GeometryUtils::createTetrahedralMeshCover(*surfMesh, 80, 40, 60);
+ auto conn = tetMesh->getTetrahedraVertices();
+ auto numVerts = tetMesh->getNumVertices();
+ std::cout << "Number of vertices = " << numVerts << std::endl;
+ testRCM(conn, numVerts);
+ }
+ using namespace std::chrono_literals;
+ std::this_thread::sleep_for(5s);
+ return 0;
+}
+
+template<typename ElemConn>
+void
+testRCM(const std::vector<ElemConn>& conn, const size_t numVerts)
+{
+ std::cout << "Old bandwidth = " << bandwidth(conn, numVerts) << std::endl;
+
+ // new-to-old permutation
+ auto perm = GeometryUtils::reorderConnectivity(conn, numVerts);
+
+ // old-to-new permutation
+ std::vector<size_t> invPerm(perm.size());
+ for (size_t i = 0; i < perm.size(); ++i)
+ {
+ CHECK(perm[i] < numVerts) << "new vertex index should not be greater than number of vertices";
+ invPerm[perm[i]] = i;
+ }
+
+ auto newConn = conn;
+
+ for (auto& vertices : newConn)
+ {
+ for (auto& vid : vertices)
+ {
+ CHECK(vid < numVerts) << "Vertex id invalid since its greater than the number of vertices";
+ vid = invPerm[vid];
+ }
+ }
+
+ std::cout << "New bandwidth = " << bandwidth(newConn, numVerts) << "\n" << std::endl;
+
+ return;
+}
+
+std::pair<std::vector<QuadConn>, size_t>
+createConn()
+{
+ /**
+ 6-------9-------7-------8
+ | | | |
+ | 6 | 7 | 8 |
+ | | | |
+ 4------11-------5-------10
+ | | | |
+ | 3 | 4 | 5 |
+ | | | |
+ 2------13-------3-------12
+ | | | |
+ | 0 | 1 | 2 |
+ | | | |
+ 0------15-------1-------14
+ **/
+
+ std::vector<QuadConn> conn(9);
+ conn[0] = { 0, 15, 13, 2 };
+ conn[1] = { 15, 1, 3, 13 };
+ conn[2] = { 1, 14, 12, 3 };
+ conn[3] = { 2, 13, 11, 4 };
+ conn[4] = { 13, 3, 5, 11 };
+ conn[5] = { 3, 12, 10, 5 };
+ conn[6] = { 4, 11, 9, 6 };
+ conn[7] = { 11, 5, 7, 9 };
+ conn[8] = { 5, 10, 8, 7 };
+
+ return std::make_pair(conn, 16);
+}
diff --git a/Examples/RCM/bandwidth.h b/Examples/RCM/bandwidth.h
new file mode 100644
index 0000000000000000000000000000000000000000..3ddf6f1ed76733c427550bf7c88be31f745e5d99
--- /dev/null
+++ b/Examples/RCM/bandwidth.h
@@ -0,0 +1,107 @@
+///
+/// \brief Build the vertex-to-vertex connectivity of a map
+//
+/// \param conn element-to-vertex connectivity of the map
+/// \param numVerts number of vertice in the map
+/// \retval vertToVert vertex-to-vertex connectivity
+///
+template<typename ElemConn>
+static void
+buildVertToVert(const std::vector<ElemConn>& conn,
+ const size_t numVerts,
+ std::vector<std::unordered_set<size_t>>& vertToVert)
+{
+ // constexpr size_t numVertPerElem = ElemConn::size();
+ std::vector<size_t> vertToElemPtr(numVerts + 1, 0);
+ std::vector<size_t> vertToElem;
+
+ // find the number of adjacent elements for each vertex
+ for (const auto& vertices : conn)
+ {
+ for (auto vid : vertices)
+ {
+ vertToElemPtr[vid + 1] += 1;
+ }
+ }
+
+ // accumulate pointer
+ for (size_t i = 0; i < numVerts; ++i)
+ {
+ vertToElemPtr[i + 1] += vertToElemPtr[i];
+ }
+
+ // track the number
+ auto pos = vertToElemPtr;
+ size_t totNum = vertToElemPtr.back();
+
+ vertToElem.resize(totNum);
+
+ for (size_t eid = 0; eid < conn.size(); ++eid)
+ {
+ for (auto vid : conn[eid])
+ {
+ vertToElem[pos[vid]] = eid;
+ ++pos[vid];
+ }
+ }
+
+ // connectivity of vertex-to-vertex
+ vertToVert.resize(numVerts);
+ auto getVertexNbrs = [&vertToElem, &vertToElemPtr, &conn, &vertToVert](const size_t i) {
+ const auto ptr0 = vertToElemPtr[i];
+ const auto ptr1 = vertToElemPtr[i + 1];
+ size_t eid;
+
+ for (auto ptr = ptr0; ptr < ptr1; ++ptr)
+ {
+ eid = vertToElem[ptr];
+ for (auto vid : conn[eid])
+ {
+ // vertex-i itself is also included.
+ vertToVert[i].insert(vid);
+ }
+ }
+ };
+
+ for (size_t i = 0; i < numVerts; ++i)
+ {
+ getVertexNbrs(i);
+ }
+}
+
+///
+/// \brief Returns the bandwidth of a map
+///
+/// \param neighbors array of neighbors of each vertex; eg, neighbors[i] is a object containing
+///
+template<typename NBR>
+size_t
+bandwidth(const std::vector<NBR>& neighbors)
+{
+ size_t d = 0;
+ size_t dCur = 0;
+ for (size_t i = 0; i < neighbors.size(); ++i)
+ {
+ for (const auto& j : neighbors[i])
+ {
+ dCur = (i > j) ? (i - j) : (j - i);
+ d = std::max(d, dCur);
+ }
+ }
+ return d;
+}
+
+///
+/// \brief Returns the bandwidth of a map
+///
+/// \param conn element-to-vertex connectivity of the map
+/// \param numVerts number of vertices in the map
+///
+template<typename ElemConn>
+size_t
+bandwidth(const std::vector<ElemConn>& conn, const size_t numVerts)
+{
+ std::vector<std::unordered_set<size_t>> vertToVert;
+ buildVertToVert(conn, numVerts, vertToVert);
+ return bandwidth(vertToVert);
+}
diff --git a/Source/Geometry/Analytic/imstkSphere.cpp b/Source/Geometry/Analytic/imstkSphere.cpp
index 65af6c65b4aa4811870f89d2197f8ec6f6c2448a..2005977042b6a0f290e24e219b37d2c5ac66dd85 100644
--- a/Source/Geometry/Analytic/imstkSphere.cpp
+++ b/Source/Geometry/Analytic/imstkSphere.cpp
@@ -67,7 +67,7 @@ Sphere::setRadius(const double r)
}
void
-Sphere::computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent)
+Sphere::computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent) const
{
updatePostTransformData();
const Vec3d span = Vec3d(1, 1, 1) * m_radiusPostTransform;
diff --git a/Source/Geometry/Analytic/imstkSphere.h b/Source/Geometry/Analytic/imstkSphere.h
index aab7b3da6be7ebb8a859e807d827a1fdd4ff6ea8..9b5b848d8f21f9d3c8f055672d763d2a615cc922 100644
--- a/Source/Geometry/Analytic/imstkSphere.h
+++ b/Source/Geometry/Analytic/imstkSphere.h
@@ -62,7 +62,7 @@ public:
///
/// \brief Compute the bounding box for the geometry
///
- virtual void computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent = 0.0) override;
+ virtual void computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent = 0.0) const override;
protected:
friend class VTKSphereRenderDelegate;
diff --git a/Source/Geometry/Map/imstkTetraTriangleMap.cpp b/Source/Geometry/Map/imstkTetraTriangleMap.cpp
index 70d05df9d33783e50d37821a0adc90520fa3877a..84122c1f223700fe56f9bd728c2b8882ad94b506 100644
--- a/Source/Geometry/Map/imstkTetraTriangleMap.cpp
+++ b/Source/Geometry/Map/imstkTetraTriangleMap.cpp
@@ -246,10 +246,9 @@ TetraTriangleMap::updateBoundingBox()
m_bBoxMin.resize(tetMesh->getNumTetrahedra());
m_bBoxMax.resize(tetMesh->getNumTetrahedra());
- for (size_t idx = 0; idx < tetMesh->getNumTetrahedra(); ++idx)
- {
- tetMesh->computeTetrahedronBoundingBox(idx, m_bBoxMin[idx], m_bBoxMax[idx]);
- }
+ ParallelUtils::parallelFor(tetMesh->getNumTetrahedra(), [&](const size_t tid) {
+ tetMesh->computeTetrahedronBoundingBox(tid, m_bBoxMin[tid], m_bBoxMax[tid]);
+ });
m_boundingBoxAvailable = true;
}
diff --git a/Source/Geometry/Mesh/imstkPointSet.cpp b/Source/Geometry/Mesh/imstkPointSet.cpp
index b33d9af57cbcb7a54328e4f193d0293aeef13e44..c97e326d5038e6f7f4c73d6bc847e5e385d5f305 100644
--- a/Source/Geometry/Mesh/imstkPointSet.cpp
+++ b/Source/Geometry/Mesh/imstkPointSet.cpp
@@ -53,7 +53,7 @@ PointSet::print() const
}
void
-PointSet::computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent)
+PointSet::computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent) const
{
updatePostTransformData();
ParallelUtils::findAABB(m_vertexPositions, lowerCorner, upperCorner);
diff --git a/Source/Geometry/Mesh/imstkPointSet.h b/Source/Geometry/Mesh/imstkPointSet.h
index f47b909e989811956a17158c2334f4fe0beb1376..0521489d0ca5f1659716a0fd9f7304794723e2b6 100644
--- a/Source/Geometry/Mesh/imstkPointSet.h
+++ b/Source/Geometry/Mesh/imstkPointSet.h
@@ -64,7 +64,7 @@ public:
///
/// \brief Compute the bounding box for the entire mesh
///
- virtual void computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent = 0.0) override;
+ virtual void computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent = 0.0) const override;
// Accessors
diff --git a/Source/Geometry/Mesh/imstkSurfaceMesh.cpp b/Source/Geometry/Mesh/imstkSurfaceMesh.cpp
index 4cf59846585bd02e82179b9ee3e23b75f571bd84..0760cdf30773cb0376de4b7557b3ffcfcb09e1b7 100644
--- a/Source/Geometry/Mesh/imstkSurfaceMesh.cpp
+++ b/Source/Geometry/Mesh/imstkSurfaceMesh.cpp
@@ -609,220 +609,4 @@ SurfaceMesh::getMaxNumTriangles()
{
return m_maxNumTriangles;
}
-
-std::vector<bool>
-SurfaceMesh::markPointsInsideAndOut(const StdVectorOfVec3d& coords)
-{
- std::vector<bool> isInside;
- isInside.resize(coords.size(), false);
-
- Vec3d aabbMin, aabbMax;
- this->computeBoundingBox(aabbMin, aabbMax, 1.);
-
- auto genRandomDirection = [](Vec3d& dir)
- {
- for (int i = 0; i < 3; ++i)
- {
- dir[i] = rand();
- }
- double mag = dir.norm();
-
- for (int i = 0; i < 3; ++i)
- {
- dir[i] /= mag;
- }
- return;
- };
-
- auto intersectTriangle = [](const Vec3d& xyz, const Vec3d& xyz0, const Vec3d& xyz1, const Vec3d& xyz2, const Vec3d& dir)
- {
- // const double eps = 1e-15;
- constexpr const double eps = std::numeric_limits<double>::epsilon();
- Vec3d edge0 = xyz1 - xyz0;
- Vec3d edge1 = xyz2 - xyz0;
- Vec3d pvec = dir.cross(edge1);
- double det = edge0.dot(pvec);
-
- if (det > -eps && det < eps)
- {
- return false;
- }
- double inv_det = 1.0 / det;
- Vec3d tvec = xyz - xyz0;
- double u = tvec.dot(pvec) * inv_det;
- if (u < 0.0 || u > 1.0)
- {
- return false;
- }
- Vec3d qvec = tvec.cross(edge0);
- double v = dir.dot(qvec) * inv_det;
- if (v < 0.0 || u + v > 1.0)
- {
- return false;
- }
-
- double t = edge1.dot(qvec) * inv_det;
- if (t > 0.0)
- {
- return true;
- }
- else
- {
- return false;
- }
- };
-
-#if 1
- std::vector<Vec3d> bBoxMin;
- std::vector<Vec3d> bBoxMax;
-
- bBoxMin.resize(this->getNumTriangles());
- bBoxMax.resize(this->getNumTriangles());
-
- for (size_t idx = 0; idx < this->getNumTriangles(); ++idx)
- {
- const auto& verts = m_trianglesVertices.at(idx);
- const auto& xyz0 = m_vertexPositions[verts[0]];
- const auto& xyz1 = m_vertexPositions[verts[1]];
- const auto& xyz2 = m_vertexPositions[verts[2]];
-
- bBoxMin[idx][0] = xyz0[0];
- bBoxMin[idx][1] = xyz0[1];
- bBoxMin[idx][2] = xyz0[2];
- bBoxMax[idx][0] = xyz0[0];
- bBoxMax[idx][1] = xyz0[1];
- bBoxMax[idx][2] = xyz0[2];
-
- bBoxMin[idx][0] = std::min(bBoxMin[idx][0], xyz1[0]);
- bBoxMin[idx][1] = std::min(bBoxMin[idx][1], xyz1[1]);
- bBoxMin[idx][2] = std::min(bBoxMin[idx][2], xyz1[2]);
- bBoxMin[idx][0] = std::min(bBoxMin[idx][0], xyz2[0]);
- bBoxMin[idx][1] = std::min(bBoxMin[idx][1], xyz2[1]);
- bBoxMin[idx][2] = std::min(bBoxMin[idx][2], xyz2[2]);
-
- bBoxMax[idx][0] = std::max(bBoxMax[idx][0], xyz1[0]);
- bBoxMax[idx][1] = std::max(bBoxMax[idx][1], xyz1[1]);
- bBoxMax[idx][2] = std::max(bBoxMax[idx][2], xyz1[2]);
- bBoxMax[idx][0] = std::max(bBoxMax[idx][0], xyz2[0]);
- bBoxMax[idx][1] = std::max(bBoxMax[idx][1], xyz2[1]);
- bBoxMax[idx][2] = std::max(bBoxMax[idx][2], xyz2[2]);
- }
-
- auto rayTracingFunc = [&coords, &aabbMin, &aabbMax, &bBoxMin, &bBoxMax, &isInside, &intersectTriangle, &genRandomDirection, this](const size_t i) {
- bool outBox = coords[i][0] < aabbMin[0] || coords[i][0] > aabbMax[0]
- || coords[i][1] < aabbMin[1] || coords[i][1] > aabbMax[1]
- || coords[i][2] < aabbMin[2] || coords[i][2] > aabbMax[2];
- if (outBox)
- {
- return;
- }
-
- // TODO: generate a random direction?
- const Vec3d direction = { 0.0, 0.0, 1.0 };
- // Vec3d direction;
- // genRandomDirection(direction);
- // std::cout << direction << std::endl;
- int numIntersections = 0;
- const auto& xyz = m_vertexPositions;
-
- for (size_t j = 0; j < this->getNumTriangles(); ++j)
- {
- const auto& verts = m_trianglesVertices[j];
-
- // consider directed ray
- if (coords[i][2] > bBoxMax[j][2])
- {
- continue;
- }
-
- if (coords[i][0] > bBoxMax[j][0])
- {
- continue;
- }
- if (coords[i][0] < bBoxMin[j][0])
- {
- continue;
- }
- if (coords[i][1] > bBoxMax[j][1])
- {
- continue;
- }
- if (coords[i][1] < bBoxMin[j][1])
- {
- continue;
- }
-
- auto intersected = intersectTriangle(coords[i],
- xyz[verts[0]],
- xyz[verts[1]],
- xyz[verts[2]],
- direction);
- if (intersected)
- {
- ++numIntersections;
- }
- }
-
- if (numIntersections % 2 == 1)
- {
- isInside[i] = true;
- }
-
- return;
- };
-
-#else
- auto rayTracingFunc = [&coords, &aabbMin, &aabbMax, &isInside, &intersectTriangle, &genRandomDirection, this](const size_t i) {
- bool outBox = coords[i][0] < aabbMin[0] || coords[i][0] > aabbMax[0]
- || coords[i][1] < aabbMin[1] || coords[i][1] > aabbMax[1]
- || coords[i][2] < aabbMin[2] || coords[i][2] > aabbMax[2];
- if (outBox)
- {
- return;
- }
-
- // TODO: generate a random direction?
- // const Vec3d direction = {1.0, 0.0, 0.0};
- Vec3d direction;
- genRandomDirection(direction);
- // std::cout << direction << std::endl;
- int numIntersections = 0;
-
- for (size_t j = 0; j < this->getNumTriangles(); ++j)
- {
- const auto& verts = m_trianglesVertices[j];
- auto intersected = intersectTriangle(coords[i],
- m_vertexPositions[verts[0]],
- m_vertexPositions[verts[1]],
- m_vertexPositions[verts[2]],
- direction);
- if (intersected)
- {
- ++numIntersections;
- }
- }
-
- if (numIntersections % 2 == 1)
- {
- isInside[i] = true;
- }
-
- return;
- };
-
-#endif
- // for (size_t i = 0; i < coords.size(); ++i)
- // {
- // rayTracingFunc(i);
- // }
-
- // CpuTimer timer;
- StopWatch timer;
- timer.start();
- ParallelUtils::parallelFor(coords.size(), rayTracingFunc);
- timer.stop();
- timer.printTimeElapsed();
-
- return isInside;
-}
} // namespace imstk
diff --git a/Source/Geometry/Mesh/imstkSurfaceMesh.h b/Source/Geometry/Mesh/imstkSurfaceMesh.h
index fc091f322c8224f9117fe76c621f5ecc34380e30..6a53c5c7a9bc3d664d6db972a1cda35e9e2ff385 100644
--- a/Source/Geometry/Mesh/imstkSurfaceMesh.h
+++ b/Source/Geometry/Mesh/imstkSurfaceMesh.h
@@ -212,12 +212,6 @@ public:
///
size_t getMaxNumTriangles();
- ///
- /// \brief Given a set of points mark them as inside (true) and outside
- /// <Explain how this works; limitations etc.>
- ///
- /// \note this function cannot be const because PointSet::computeBoundingBox, called inside, is not.
- std::vector<bool> markPointsInsideAndOut(const StdVectorOfVec3d& coords);
protected:
friend class VTKSurfaceMeshRenderDelegate;
diff --git a/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp b/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp
index 815099c9abe03c0c7f8942535365fb4ee65c8317..2de16dbebeb1603a4629316dca330ba9322cb28e 100644
--- a/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp
+++ b/Source/Geometry/Mesh/imstkTetrahedralMesh.cpp
@@ -295,278 +295,4 @@ TetrahedralMesh::getNumTetrahedra() const
{
return m_tetrahedraVertices.size();
}
-
-// std::unique_ptr<TetrahedralMesh>
-std::shared_ptr<TetrahedralMesh>
-TetrahedralMesh::createUniformMesh(const Vec3d& aabbMin, const Vec3d& aabbMax, const size_t nx,
- const size_t ny, const size_t nz)
-{
- Vec3d h = { (aabbMax[0] - aabbMin[0]) / nx,
- (aabbMax[1] - aabbMin[1]) / ny,
- (aabbMax[2] - aabbMin[2]) / nz };
- LOG_IF(FATAL, (h[0] <= 0.0 || h[1] <= 0.0 || h[2] <= 0.0)) << "Invalid bounding box";
-
- size_t numVertices = (nx + 1) * (ny + 1) * (nz + 1);
-
- // std::vector<Vec3d> coords;
- StdVectorOfVec3d coords;
- coords.resize(numVertices);
- size_t cnt = 0;
-
- for (size_t k = 0; k < nz + 1; ++k)
- {
- double z = aabbMin[2] + k * h[2];
- for (size_t j = 0; j < ny + 1; ++j)
- {
- double y = aabbMin[1] + j * h[1];
- for (size_t i = 0; i < nx + 1; ++i)
- {
- double x = aabbMin[0] + i * h[0];
- coords[cnt] = { x, y, z };
- ++cnt;
- }
- }
- }
-
- const size_t numDiv = 6;
- size_t numTets = numDiv * nx * ny * nz;
-
- std::vector<TetrahedralMesh::TetraArray> vertices;
- vertices.resize(numTets);
- cnt = 0;
- std::vector<size_t> indx(8);
-
- for (size_t k = 0; k < nz; ++k)
- {
- for (size_t j = 0; j < ny; ++j)
- {
- for (size_t i = 0; i < nx; ++i)
- {
- indx[3] = i + j * (nx + 1) + k * (nx + 1) * (ny + 1);
- indx[2] = indx[3] + 1;
- indx[0] = indx[3] + nx + 1;
- indx[1] = indx[0] + 1;
- indx[4] = indx[0] + (nx + 1) * (ny + 1);
- indx[5] = indx[1] + (nx + 1) * (ny + 1);
- indx[6] = indx[2] + (nx + 1) * (ny + 1);
- indx[7] = indx[3] + (nx + 1) * (ny + 1);
- vertices[cnt + 0] = { indx[0], indx[2], indx[3], indx[6] };
- vertices[cnt + 1] = { indx[0], indx[3], indx[7], indx[6] };
- vertices[cnt + 2] = { indx[0], indx[7], indx[4], indx[6] };
- vertices[cnt + 3] = { indx[0], indx[5], indx[6], indx[4] };
- vertices[cnt + 4] = { indx[1], indx[5], indx[6], indx[0] };
- vertices[cnt + 5] = { indx[1], indx[6], indx[2], indx[0] };
- cnt += numDiv;
- }
- }
- }
-
- auto mesh = std::make_shared<TetrahedralMesh>();
- // auto mesh = std::unique_ptr<TetrahedralMesh>(new TetrahedralMesh());
-
- mesh->initialize(coords, vertices);
- return mesh;
-}
-
-#if 0
-std::shared_ptr<TetrahedralMesh>
-TetrahedralMesh::createEnclosingMesh(const SurfaceMesh& surfMesh, const size_t nx, const size_t ny,
- const size_t nz)
-{
- // Given the index of a tet, return the corresponding index of the hex
- auto tetIdToHexId = [](const size_t tetId) { return tetId / 5; }
-
- Vec3d aabbMin, aabbMax;
- const double paddingPerc = 10.0;
- PointSet.computeBoundingBox(aabbMin, aabbMax, paddingPerc);
- const Vec3d h = { (aabbMax[0] - aabbMin[0]) / nx,
- (aabbMax[1] - aabbMin[1]) / ny,
- (aabbMax[2] - aabbMin[2]) / nz };
-
- auto findHexId = [&aabbMin, &aabbMax, &h](const Vec3d& xyz) {
- size_t idX = (xyz[0] - aabbMin[0]) / h[0];
- size_t idY = (xyz[1] - aabbMin[1]) / h[1];
- size_t idZ = (xyz[2] - aabbMin[2]) / h[2];
- return { idX, idY, idZ };
- };
-
- auto uniformMesh = createUniformMesh(aabbMin, aabbMax, nx, ny, nz);
-
- std::vector<bool> enclosePoint(uniformMesh.getNumHexahedra(), false);
-
- // ParallelUtils::parallelFor(surfMesh.getNumVertices(), [&](const size_t vid) {
- // auto xyz = surfMesh.getVertexPosition(vid);
- // size_t idX = (xyz[0] - aabbMin[0]) / h[0];
- // size_t idY = (xyz[1] - aabbMin[1]) / h[1];
- // size_t idZ = (xyz[2] - aabbMin[2]) / h[2];
- // size_t hexId = idX + idY *nx + idZ * nx*ny;
- // size_t tetId0 = 5*hexId;
- // size_t tetId1 = tetId0 + 5;
- //
- //
- // static ParallelUtils::SpinLock lock;
- // lock.lock();
- // for (size_t id = tetId0; id<tetId1; ++id)
- // {
- // enclosePoint[id] = true;
- // }
- // lock.unlock();
- // });
-
- for (size_t vid = 0; vid < surfMesh.getNumVertices(); ++vid)
- {
- auto xyz = surfMesh.getVertexPosition(vid);
- size_t idX = (xyz[0] - aabbMin[0]) / h[0];
- size_t idY = (xyz[1] - aabbMin[1]) / h[1];
- size_t idZ = (xyz[2] - aabbMin[2]) / h[2];
- size_t hexId = idX + idY * nx + idZ * nx * ny;
- size_t tetId0 = 5 * hexId;
- size_t tetId1 = tetId0 + 5;
-
- for (size_t i = tetId0; i < tetId1; ++i)
- {
- enclosePoint[i] = true;
- }
- }
-}
-
-#endif
-
-std::shared_ptr<TetrahedralMesh>
-TetrahedralMesh::createTetrahedralMeshCover(std::shared_ptr<SurfaceMesh> surfMesh, const size_t nx, const size_t ny,
- const size_t nz)
-{
- Vec3d aabbMin, aabbMax;
-
- // create a background mesh
- surfMesh->computeBoundingBox(aabbMin, aabbMax, 1. /*percentage padding*/);
- auto uniformMesh = createUniformMesh(aabbMin, aabbMax, nx, ny, nz);
-
- // ray-tracing
- const auto& coords = uniformMesh->getVertexPositions();
- auto insideSurfMesh = surfMesh->markPointsInsideAndOut(coords);
-
- // label elements
- std::vector<bool> validTet(uniformMesh->getNumTetrahedra(), false);
- std::vector<bool> validVtx(uniformMesh->getNumVertices(), false);
-
- // TetrahedralMesh::WeightsArray weights = {0.0, 0.0, 0.0, 0.0};
- const Vec3d h = { (aabbMax[0] - aabbMin[0]) / nx,
- (aabbMax[1] - aabbMin[1]) / ny,
- (aabbMax[2] - aabbMin[2]) / nz };
-
- // TODO: can be parallelized by make NUM_THREADS copies of validTet, or use atomic op on validTet
- auto labelEnclosingTet = [&surfMesh, &uniformMesh, &aabbMin, &h, nx, ny, nz, &validTet](const size_t i)
- {
- const auto& xyz = surfMesh->getVertexPosition(i);
- // find the hex that encloses the point;
- size_t idX = (size_t)((xyz[0] - aabbMin[0]) / h[0]);
- size_t idY = (size_t)((xyz[1] - aabbMin[1]) / h[1]);
- size_t idZ = (size_t)((xyz[2] - aabbMin[2]) / h[2]);
- size_t hexId = idX + idY * nx + idZ * nx * ny;
-
- // the index range of tets inside the enclosing hex
- const int numDiv = 6;
- size_t tetId0 = numDiv * hexId;
- size_t tetId1 = tetId0 + numDiv;
-
- // loop over the tets to find the enclosing tets
- bool found = false;
- TetrahedralMesh::WeightsArray weights = { 0.0, 0.0, 0.0, 0.0 };
- for (size_t id = tetId0; id < tetId1; ++id)
- {
- uniformMesh->computeBarycentricWeights(id, xyz, weights);
-
- if ((weights[0] >= 0) && (weights[1] >= 0) && (weights[2] >= 0) && (weights[3] >= 0))
- {
- validTet[id] = true;
- found = true;
- break;
- }
- }
-
- // TODO: not sure how to do thread-safe logging
- CHECK(found) << "Failed to find the enclosing tetrahedron";
- };
-
- // a customized approach to find the enclosing tet for each surface points
- for (size_t i = 0; i < surfMesh->getNumVertices(); ++i)
- {
- labelEnclosingTet(i);
- }
-
- for (size_t i = 0; i < validTet.size(); ++i)
- {
- const auto& verts = uniformMesh->getTetrahedronVertices(i);
- if (insideSurfMesh[verts[0]]
- || insideSurfMesh[verts[1]]
- || insideSurfMesh[verts[2]]
- || insideSurfMesh[verts[3]])
- {
- validTet[i] = true;
- }
- }
-
- size_t numElems = 0;
- for (size_t i = 0; i < validTet.size(); ++i)
- {
- const auto& verts = uniformMesh->getTetrahedronVertices(i);
- if (validTet[i])
- {
- validVtx[verts[0]] = true;
- validVtx[verts[1]] = true;
- validVtx[verts[2]] = true;
- validVtx[verts[3]] = true;
- ++numElems;
- }
- }
-
- // discard useless vertices, and build old-to-new index map
- size_t numVerts = 0;
- for (auto b : validVtx)
- {
- if (b)
- {
- ++numVerts;
- }
- }
-
- StdVectorOfVec3d newCoords(numVerts);
- std::vector<size_t> oldToNew(coords.size(), std::numeric_limits<size_t>::max());
- size_t cnt = 0;
-
- for (size_t i = 0; i < validVtx.size(); ++i)
- {
- if (validVtx[i])
- {
- newCoords[cnt] = coords[i];
- oldToNew[i] = cnt;
- ++cnt;
- }
- }
-
- // update tet-to-vert connectivity
- std::vector<TetraArray> newIndices(numElems);
- cnt = 0;
- for (size_t i = 0; i < uniformMesh->getNumTetrahedra(); ++i)
- {
- if (validTet[i])
- {
- const auto oldIndices = uniformMesh->getTetrahedronVertices(i);
-
- newIndices[cnt][0] = oldToNew[oldIndices[0]];
- newIndices[cnt][1] = oldToNew[oldIndices[1]];
- newIndices[cnt][2] = oldToNew[oldIndices[2]];
- newIndices[cnt][3] = oldToNew[oldIndices[3]];
-
- ++cnt;
- }
- }
-
- // ready to create the final mesh
- auto mesh = std::make_shared<TetrahedralMesh>();
- mesh->initialize(newCoords, newIndices);
-
- return mesh;
-}
} // namespace imstk
diff --git a/Source/Geometry/Mesh/imstkTetrahedralMesh.h b/Source/Geometry/Mesh/imstkTetrahedralMesh.h
index b667fd2c0fddc045d19f44b9b8d39594feee1f2c..0aeb2c944596d030acb5fef413e7077a7bb0c33b 100644
--- a/Source/Geometry/Mesh/imstkTetrahedralMesh.h
+++ b/Source/Geometry/Mesh/imstkTetrahedralMesh.h
@@ -123,23 +123,6 @@ public:
void setTetrahedraAsRemoved(const unsigned int tetId) { m_removedMeshElems[tetId] = true; }
const std::vector<bool>& getRemovedTetrahedra() const { return m_removedMeshElems; }
- ///
- /// \brief Create a tetrahedral mesh based on a uniform Cartesian mesh
- /// \param aabbMin the small conner of a box
- /// \param aabbMax the large conner of a box
- /// \param nx number of elements in the x-direction
- /// \param ny number of elements in the y-direction
- /// \param nz number of elements in the z-direction
- ///
- /// \note Refer: Dompierre, Julien & Labbé, Paul & Vallet, Marie-Gabrielle & Camarero, Ricardo. (1999).
- /// How to Subdivide Pyramids, Prisms, and Hexahedra into Tetrahedra.. 195-204.
- static std::shared_ptr<TetrahedralMesh> createUniformMesh(const Vec3d& aabbMin, const Vec3d& aabbMax, const size_t nx, const size_t ny, const size_t nz);
-
- ///
- /// \brief Create a tetrahedral mesh cover
- /// \note surfMesh can't be const since the non-const member function rayTracing is called inside.
- ///
- static std::shared_ptr<TetrahedralMesh> createTetrahedralMeshCover(std::shared_ptr<SurfaceMesh> surfMesh, const size_t nx, const size_t ny, size_t nz);
protected:
friend class VTKTetrahedralMeshRenderDelegate;
diff --git a/Source/Geometry/imstkGeometry.cpp b/Source/Geometry/imstkGeometry.cpp
index b03022d1697b224a03cdc00f659482f5153d41bd..48d7cd0432220a2b6ac60f519e88afa43d85adfc 100644
--- a/Source/Geometry/imstkGeometry.cpp
+++ b/Source/Geometry/imstkGeometry.cpp
@@ -63,7 +63,7 @@ Geometry::print() const
}
void
-Geometry::computeBoundingBox(Vec3d&, Vec3d&, const double)
+Geometry::computeBoundingBox(Vec3d&, Vec3d&, const double) const
{
LOG(FATAL) << "computeBoundingBox() must be called from an instance of a specific geometry class";
}
diff --git a/Source/Geometry/imstkGeometry.h b/Source/Geometry/imstkGeometry.h
index fc9c0993c012b9a769dd7481b71b1e56559f84cf..901165bf6e863a596d39cfbffa850bf5487f1b47 100644
--- a/Source/Geometry/imstkGeometry.h
+++ b/Source/Geometry/imstkGeometry.h
@@ -103,7 +103,7 @@ public:
///
/// \brief Compute the bounding box for the geometry
///
- virtual void computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent = 0.0);
+ virtual void computeBoundingBox(Vec3d& lowerCorner, Vec3d& upperCorner, const double paddingPercent = 0.0) const;
///
/// \brief Translate the geometry in Cartesian space
diff --git a/Source/Geometry/imstkGeometryUtilities.cpp b/Source/Geometry/imstkGeometryUtilities.cpp
index 0000f7dc6063760a3ebdc106e862146ff49f0f14..02d47efe403e77ec3e17b3299b21087555f8672f 100644
--- a/Source/Geometry/imstkGeometryUtilities.cpp
+++ b/Source/Geometry/imstkGeometryUtilities.cpp
@@ -20,6 +20,7 @@
=========================================================================*/
#include "imstkGeometryUtilities.h"
+#include "Mesh/imstkSurfaceMesh.h"
#include "imstkHexahedralMesh.h"
#include "imstkLineMesh.h"
#include "imstkTetrahedralMesh.h"
@@ -358,4 +359,876 @@ GeometryUtils::loopSubdivideSurfaceMesh(const SurfaceMesh& surfaceMesh, const in
return convertVtkPolyDataToSurfaceMesh(filter->GetOutput());
}
+
+namespace // anonymous namespace
+{
+///
+/// \brief Build the vertex-to-vertex connectivity
+///
+/// \param[in] conn element-to-vertex connectivity
+/// \param[in] numVerts number of vertices
+/// \param[out] vertToVert the vertex-to-vertex connectivity on return
+///
+template<typename ElemConn>
+static void
+buildVertexToVertexConnectivity(const std::vector<ElemConn>& conn,
+ const size_t numVerts,
+ std::vector<std::unordered_set<size_t>>& vertToVert)
+{
+ // constexpr size_t numVertPerElem = ElemConn::size();
+ std::vector<size_t> vertToElemPtr(numVerts + 1, 0);
+ std::vector<size_t> vertToElem;
+
+ // find the number of adjacent elements for each vertex
+ for (const auto& vertices : conn)
+ {
+ for (auto vid : vertices)
+ {
+ vertToElemPtr[vid + 1] += 1;
+ }
+ }
+
+ // accumulate pointer
+ for (size_t i = 0; i < numVerts; ++i)
+ {
+ vertToElemPtr[i + 1] += vertToElemPtr[i];
+ }
+
+ // track the front position for each vertex in \p vertToElem
+ auto pos = vertToElemPtr;
+ // the total number of element neighbors of all vertices, ie, size of \p vertToElem
+ const size_t totNum = vertToElemPtr.back();
+
+ vertToElem.resize(totNum);
+
+ for (size_t eid = 0; eid < conn.size(); ++eid)
+ {
+ for (auto vid : conn[eid])
+ {
+ vertToElem[pos[vid]] = eid;
+ ++pos[vid];
+ }
+ }
+
+ // connectivity of vertex-to-vertex
+ vertToVert.resize(numVerts);
+ auto getVertexNbrs = [&vertToElem, &vertToElemPtr, &conn, &vertToVert](const size_t i)
+ {
+ const auto ptr0 = vertToElemPtr[i];
+ const auto ptr1 = vertToElemPtr[i + 1];
+ size_t eid;
+
+ for (auto ptr = ptr0; ptr < ptr1; ++ptr)
+ {
+ eid = vertToElem[ptr];
+ for (auto vid : conn[eid])
+ {
+ // vertex-i itself is also included.
+ vertToVert[i].insert(vid);
+ }
+ }
+ };
+
+ // for (size_t i = 0; i < numVerts; ++i)
+ // {
+ // getVertexNbrs(i);
+ // }
+ ParallelUtils::parallelFor(numVerts, getVertexNbrs);
}
+
+///
+/// \brief Reverse Cuthill-Mckee (RCM) based reording to reduce bandwidth
+///
+/// \param[in] neighbors array of neighbors of each vertex; eg, neighbors[i] is an object containing all neighbors of vertex-i
+/// \return the permutation vector that map from new indices to old indices
+///
+/// \cite https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm
+///
+template<typename NeighborContainer>
+static std::vector<size_t>
+RCM(const std::vector<NeighborContainer>& neighbors)
+{
+ const size_t invalid = std::numeric_limits<size_t>::max();
+
+ // is greater in terms of degrees
+ auto isGreater = [&neighbors](const size_t i, const size_t j) {
+ return neighbors[i].size() > neighbors[j].size() ? true : false;
+ };
+
+ const size_t numVerts = neighbors.size();
+
+ std::vector<size_t> P(numVerts);
+ for (size_t i = 0; i < numVerts; ++i)
+ {
+ P[i] = i;
+ }
+ std::sort(P.begin(), P.end(), isGreater);
+
+ // \todo an alternative is to use std::set for P
+ // std::set<size_t, isGreater> P;
+ // for (size_t i=0; i<numVerts; ++i)
+ // {
+ // P.insert(i);
+ // }
+
+ std::vector<size_t> R(numVerts, invalid); // permutation
+ std::queue<size_t> Q; // queue
+ std::vector<bool> isInP(numVerts, true); // if a vertex is still in P
+ size_t pos = 0; // track how many vertices are put into R
+
+ ///
+ /// \brief Move a vertex into R, and enque its neighbors into Q in ascending order.
+ /// \param vid index of vertex to be moved into R
+ ///
+ auto moveVertexIntoRAndItsNeighborsIntoQ = [&neighbors, &isInP, &pos, &R, &Q](const size_t vid)
+ {
+ R[pos] = vid;
+ ++pos;
+ isInP[vid] = false;
+
+ // Put the unordered neighbors into Q, in ascending order.
+ // first find unordered neighbors
+ std::set<size_t> unorderedNbrs;
+ for (auto nbr : neighbors[vid])
+ {
+ if (isInP[nbr])
+ {
+ unorderedNbrs.insert(nbr);
+ }
+ }
+
+ for (auto nbr : unorderedNbrs)
+ {
+ Q.push(nbr);
+ isInP[nbr] = false;
+ }
+ };
+
+ size_t pCur = 0;
+
+ // main loop
+ while (true)
+ {
+ size_t parent = invalid;
+ for (size_t vid = pCur; vid < isInP.size(); ++vid)
+ {
+ if (isInP[vid])
+ {
+ isInP[vid] = false;
+ pCur = vid;
+ parent = vid;
+ break;
+ }
+ }
+ if (parent == invalid)
+ {
+ break;
+ }
+
+ moveVertexIntoRAndItsNeighborsIntoQ(parent);
+
+ while (!Q.empty())
+ {
+ parent = Q.front();
+ Q.pop();
+ moveVertexIntoRAndItsNeighborsIntoQ(parent);
+ }
+
+ // here we have empty Q
+ }
+
+ CHECK(pos == numVerts) << "RCM ordering: we should never get here";
+
+ std::reverse(R.begin(), R.end());
+ return R;
+}
+
+///
+/// \brief Reorder using Reverse Cuthill-Mckee algorithm
+//
+/// \param conn element-to-vertex connectivity
+/// \param numVerts number of vertices
+///
+/// \return the permutation vector that maps from new indices to old indices
+///
+template<typename ElemConn>
+static std::vector<size_t>
+RCM(const std::vector<ElemConn>& conn, const size_t numVerts)
+{
+ // connectivity of vertex-to-vertex
+ std::vector<std::unordered_set<size_t>> vertToVert;
+ buildVertexToVertexConnectivity(conn, numVerts, vertToVert);
+ return RCM(vertToVert);
+}
+
+///
+/// \brief Given a set of points mark them as inside (true) and outside
+/// \param surfaceMesh a \ref SurfaceMesh
+/// \param coords a set of points to be tested
+///
+/// \note this function cannot be const because PointSet::computeBoundingBox, called inside, is not.
+///
+std::vector<bool>
+markPointsInsideAndOut(const SurfaceMesh& surfaceMesh, const StdVectorOfVec3d& coords)
+{
+ std::vector<bool> isInside;
+ isInside.resize(coords.size(), false);
+
+ Vec3d aabbMin, aabbMax;
+ surfaceMesh.computeBoundingBox(aabbMin, aabbMax, 1.);
+
+ auto genRandomDirection = [](Vec3d& direction)
+ {
+ for (int i = 0; i < 3; ++i)
+ {
+ direction[i] = rand();
+ }
+ double mag = direction.norm();
+
+ for (int i = 0; i < 3; ++i)
+ {
+ direction[i] /= mag;
+ }
+ return;
+ };
+
+ auto triangleRayIntersection = [](const Vec3d& xyz,
+ const Vec3d& triVert0,
+ const Vec3d& triVert1,
+ const Vec3d& triVert2,
+ const Vec3d& direction)
+ {
+ // const double eps = 1e-15;
+ constexpr const double eps = std::numeric_limits<double>::epsilon();
+ Vec3d edge0 = triVert1 - triVert0;
+ Vec3d edge1 = triVert2 - triVert0;
+ Vec3d pvec = direction.cross(edge1);
+ double det = edge0.dot(pvec);
+
+ if (det > -eps && det < eps)
+ {
+ return false;
+ }
+ double inv_det = 1.0 / det;
+ Vec3d tvec = xyz - triVert0;
+ double u = tvec.dot(pvec) * inv_det;
+ if (u < 0.0 || u > 1.0)
+ {
+ return false;
+ }
+ Vec3d qvec = tvec.cross(edge0);
+ double v = direction.dot(qvec) * inv_det;
+ if (v < 0.0 || u + v > 1.0)
+ {
+ return false;
+ }
+
+ double t = edge1.dot(qvec) * inv_det;
+ if (t > 0.0)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ };
+
+ std::vector<Vec3d> bBoxMin;
+ std::vector<Vec3d> bBoxMax;
+
+ bBoxMin.resize(surfaceMesh.getNumTriangles());
+ bBoxMax.resize(surfaceMesh.getNumTriangles());
+
+ for (size_t idx = 0; idx < surfaceMesh.getNumTriangles(); ++idx)
+ {
+ const auto& verts = surfaceMesh.getTrianglesVertices().at(idx);
+ const auto& xyz0 = surfaceMesh.getVertexPosition(verts[0]);
+ const auto& xyz1 = surfaceMesh.getVertexPosition(verts[1]);
+ const auto& xyz2 = surfaceMesh.getVertexPosition(verts[2]);
+ // const auto& xyz0 = m_vertexPositions[verts[0]];
+ // const auto& xyz1 = m_vertexPositions[verts[1]];
+ // const auto& xyz2 = m_vertexPositions[verts[2]];
+
+ bBoxMin[idx][0] = xyz0[0];
+ bBoxMin[idx][1] = xyz0[1];
+ bBoxMin[idx][2] = xyz0[2];
+ bBoxMax[idx][0] = xyz0[0];
+ bBoxMax[idx][1] = xyz0[1];
+ bBoxMax[idx][2] = xyz0[2];
+
+ bBoxMin[idx][0] = std::min(bBoxMin[idx][0], xyz1[0]);
+ bBoxMin[idx][1] = std::min(bBoxMin[idx][1], xyz1[1]);
+ bBoxMin[idx][2] = std::min(bBoxMin[idx][2], xyz1[2]);
+ bBoxMin[idx][0] = std::min(bBoxMin[idx][0], xyz2[0]);
+ bBoxMin[idx][1] = std::min(bBoxMin[idx][1], xyz2[1]);
+ bBoxMin[idx][2] = std::min(bBoxMin[idx][2], xyz2[2]);
+
+ bBoxMax[idx][0] = std::max(bBoxMax[idx][0], xyz1[0]);
+ bBoxMax[idx][1] = std::max(bBoxMax[idx][1], xyz1[1]);
+ bBoxMax[idx][2] = std::max(bBoxMax[idx][2], xyz1[2]);
+ bBoxMax[idx][0] = std::max(bBoxMax[idx][0], xyz2[0]);
+ bBoxMax[idx][1] = std::max(bBoxMax[idx][1], xyz2[1]);
+ bBoxMax[idx][2] = std::max(bBoxMax[idx][2], xyz2[2]);
+ }
+
+ auto rayTracingFunc = [&coords,
+ &aabbMin,
+ &aabbMax,
+ &bBoxMin,
+ &bBoxMax,
+ &isInside,
+ &triangleRayIntersection,
+ &genRandomDirection,
+ &surfaceMesh](const size_t i)
+ {
+ bool outBox = coords[i][0] < aabbMin[0] || coords[i][0] > aabbMax[0]
+ || coords[i][1] < aabbMin[1] || coords[i][1] > aabbMax[1]
+ || coords[i][2] < aabbMin[2] || coords[i][2] > aabbMax[2];
+ if (outBox)
+ {
+ return;
+ }
+
+ /// \todo generate a random direction?
+ const Vec3d direction = { 0.0, 0.0, 1.0 };
+ // Vec3d direction;
+ // genRandomDirection(direction);
+ // std::cout << direction << std::endl;
+ int numIntersections = 0;
+ const auto& xyz = surfaceMesh.getVertexPositions();
+ const auto& triVerts = surfaceMesh.getTrianglesVertices();
+
+ for (size_t j = 0; j < surfaceMesh.getNumTriangles(); ++j)
+ {
+ const auto& verts = triVerts[j];
+
+ // consider directed ray
+ if (coords[i][2] > bBoxMax[j][2])
+ {
+ continue;
+ }
+
+ if (coords[i][0] > bBoxMax[j][0])
+ {
+ continue;
+ }
+ if (coords[i][0] < bBoxMin[j][0])
+ {
+ continue;
+ }
+ if (coords[i][1] > bBoxMax[j][1])
+ {
+ continue;
+ }
+ if (coords[i][1] < bBoxMin[j][1])
+ {
+ continue;
+ }
+
+ auto intersected = triangleRayIntersection(coords[i],
+ xyz[verts[0]],
+ xyz[verts[1]],
+ xyz[verts[2]],
+ direction);
+ if (intersected)
+ {
+ ++numIntersections;
+ }
+ }
+
+ if (numIntersections % 2 == 1)
+ {
+ isInside[i] = true;
+ }
+
+ return;
+ };
+
+ // for (size_t i = 0; i < coords.size(); ++i)
+ // {
+ // rayTracingFunc(i);
+ // }
+
+ ParallelUtils::parallelFor(coords.size(), rayTracingFunc);
+
+ return isInside;
+}
+
+///
+/// \brief Given a set of uniformly spaced points, mark them as inside (true) and outside.
+/// It makes uses of ray-tracing but skips points based on the nearest distance between current point and the surface.
+///
+/// \param surfaceMesh a \ref SurfaceMesh
+/// \param coords a set of points to be tested
+/// \param nx number of points in x-direction
+/// \param ny number of points in y-direction
+/// \param nz number of points in z-direction
+///
+/// \note this function cannot be const because PointSet::computeBoundingBox, called inside, is not.
+///
+std::vector<bool>
+markPointsInsideAndOut(const SurfaceMesh& surfaceMesh, const StdVectorOfVec3d& coords, const size_t nx, const size_t ny, const size_t nz)
+{
+ std::vector<bool> isInside;
+ isInside.resize(coords.size(), false);
+
+ Vec3d aabbMin, aabbMax;
+ surfaceMesh.computeBoundingBox(aabbMin, aabbMax, 1.);
+ // space between two adjacent points
+ const Vec3d h = { coords[1][0] - coords[0][0], coords[nx][1] - coords[0][1], coords[nx * ny][2] - coords[0][2] };
+
+ /// \brief if a ray intersects with a triangle.
+ /// \param xyz starting point of ray
+ /// \param triVert0 triVert0, triVert1, triVert2 are 3 vertices of the triangle
+ /// \param direction direction of the ray
+ /// \param[out] distance on return it is the distance of the point and the triangle
+ auto triangleRayIntersection = [](const Vec3d& xyz, const Vec3d& triVert0, const Vec3d& triVert1, const Vec3d& triVert2, const Vec3d& direction, double& distance)
+ {
+ const double eps = std::numeric_limits<double>::epsilon();
+ const Vec3d edge0 = triVert1 - triVert0;
+ const Vec3d edge1 = triVert2 - triVert0;
+ const Vec3d pvec = direction.cross(edge1);
+ const double det = edge0.dot(pvec);
+
+ if (det > -eps && det < eps)
+ {
+ return false;
+ }
+ const double inv_det = 1.0 / det;
+ const Vec3d tvec = xyz - triVert0;
+ const double u = tvec.dot(pvec) * inv_det;
+ if (u < 0.0 || u > 1.0)
+ {
+ return false;
+ }
+ const Vec3d qvec = tvec.cross(edge0);
+ const double v = direction.dot(qvec) * inv_det;
+ if (v < 0.0 || u + v > 1.0)
+ {
+ return false;
+ }
+
+ distance = edge1.dot(qvec) * inv_det;
+ if (distance > 0.0)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ };
+
+ std::vector<Vec3d> bBoxMin;
+ std::vector<Vec3d> bBoxMax;
+
+ bBoxMin.resize(surfaceMesh.getNumTriangles());
+ bBoxMax.resize(surfaceMesh.getNumTriangles());
+
+ /// \brief find the bounding boxes of each surface triangle
+ auto findBoundingBox = [&](const size_t idx)
+ {
+ const auto& verts = surfaceMesh.getTrianglesVertices().at(idx);
+ const auto& vertXyz = surfaceMesh.getVertexPositions();
+ const auto& xyz0 = vertXyz[verts[0]];
+ const auto& xyz1 = vertXyz[verts[1]];
+ const auto& xyz2 = vertXyz[verts[2]];
+
+ bBoxMin[idx][0] = xyz0[0];
+ bBoxMin[idx][1] = xyz0[1];
+ bBoxMin[idx][2] = xyz0[2];
+ bBoxMax[idx][0] = xyz0[0];
+ bBoxMax[idx][1] = xyz0[1];
+ bBoxMax[idx][2] = xyz0[2];
+
+ bBoxMin[idx][0] = std::min(bBoxMin[idx][0], xyz1[0]);
+ bBoxMin[idx][1] = std::min(bBoxMin[idx][1], xyz1[1]);
+ bBoxMin[idx][2] = std::min(bBoxMin[idx][2], xyz1[2]);
+ bBoxMin[idx][0] = std::min(bBoxMin[idx][0], xyz2[0]);
+ bBoxMin[idx][1] = std::min(bBoxMin[idx][1], xyz2[1]);
+ bBoxMin[idx][2] = std::min(bBoxMin[idx][2], xyz2[2]);
+
+ bBoxMax[idx][0] = std::max(bBoxMax[idx][0], xyz1[0]);
+ bBoxMax[idx][1] = std::max(bBoxMax[idx][1], xyz1[1]);
+ bBoxMax[idx][2] = std::max(bBoxMax[idx][2], xyz1[2]);
+ bBoxMax[idx][0] = std::max(bBoxMax[idx][0], xyz2[0]);
+ bBoxMax[idx][1] = std::max(bBoxMax[idx][1], xyz2[1]);
+ bBoxMax[idx][2] = std::max(bBoxMax[idx][2], xyz2[2]);
+ };
+
+ ParallelUtils::parallelFor(surfaceMesh.getNumTriangles(), findBoundingBox);
+
+ // ray tracing for all points in the x-axis. These points are those start with indices (0,j,k)
+ // and jk = j + k*ny
+ auto rayTracingLine = [&](const size_t jk)
+ {
+ size_t idx0 = jk * nx;
+ bool outBox = coords[idx0][0] < aabbMin[0] || coords[idx0][0] > aabbMax[0]
+ || coords[idx0][1] < aabbMin[1] || coords[idx0][1] > aabbMax[1]
+ || coords[idx0][2] < aabbMin[2] || coords[idx0][2] > aabbMax[2];
+ if (outBox)
+ {
+ return;
+ }
+
+ const Vec3d direction = { 1.0, 0.0, 0.0 };
+ const auto& xyz = surfaceMesh.getVertexPositions();
+ const auto& triVerts = surfaceMesh.getTrianglesVertices();
+
+ size_t i = 0;
+ while (i < nx)
+ {
+ size_t idx = idx0 + i;
+ int numIntersections = 0;
+ double dist = 0.0;
+ double distMin = h[0] * (nz + 1);
+
+ for (size_t j = 0; j < surfaceMesh.getNumTriangles(); ++j)
+ {
+ const auto& verts = triVerts[j];
+
+ // consider directed ray
+ if (coords[idx][0] > bBoxMax[j][0])
+ {
+ continue;
+ }
+
+ if (coords[idx][1] > bBoxMax[j][1])
+ {
+ continue;
+ }
+ if (coords[idx][1] < bBoxMin[j][1])
+ {
+ continue;
+ }
+ if (coords[idx][2] > bBoxMax[j][2])
+ {
+ continue;
+ }
+ if (coords[idx][2] < bBoxMin[j][2])
+ {
+ continue;
+ }
+
+ auto intersected = triangleRayIntersection(coords[idx],
+ xyz[verts[0]],
+ xyz[verts[1]],
+ xyz[verts[2]],
+ direction,
+ dist);
+ if (intersected)
+ {
+ ++numIntersections;
+ distMin = std::min(dist, distMin);
+ }
+ }
+
+ // core of the algorithm: points between the current one and iEnd share the same label so we can skip them.
+ size_t iEnd = i + static_cast<size_t>(distMin / h[0]) + 1;
+ iEnd = std::min(iEnd, nx);
+
+ if (numIntersections % 2 == 1)
+ {
+ for (size_t ii = idx; ii < idx0 + iEnd; ++ii)
+ {
+ isInside[ii] = true;
+ }
+ }
+
+ i = iEnd;
+ }
+ };
+
+ ParallelUtils::parallelFor(ny * nz, rayTracingLine);
+
+ return isInside;
+}
+} // anonymous namespace
+
+std::shared_ptr<TetrahedralMesh>
+GeometryUtils::createUniformMesh(const Vec3d& aabbMin,
+ const Vec3d& aabbMax,
+ const size_t nx,
+ const size_t ny,
+ const size_t nz)
+{
+ const Vec3d h = { (aabbMax[0] - aabbMin[0]) / nx,
+ (aabbMax[1] - aabbMin[1]) / ny,
+ (aabbMax[2] - aabbMin[2]) / nz };
+ LOG_IF(FATAL, (h[0] <= 0.0 || h[1] <= 0.0 || h[2] <= 0.0)) << "Invalid bounding box";
+
+ const size_t numVertices = (nx + 1) * (ny + 1) * (nz + 1);
+
+ // std::vector<Vec3d> coords;
+ StdVectorOfVec3d coords;
+ coords.resize(numVertices);
+ size_t cnt = 0;
+
+ for (size_t k = 0; k < nz + 1; ++k)
+ {
+ double z = aabbMin[2] + k * h[2];
+ for (size_t j = 0; j < ny + 1; ++j)
+ {
+ double y = aabbMin[1] + j * h[1];
+ for (size_t i = 0; i < nx + 1; ++i)
+ {
+ double x = aabbMin[0] + i * h[0];
+ coords[cnt] = { x, y, z };
+ ++cnt;
+ }
+ }
+ }
+
+ const size_t numDiv = 6;
+ const size_t numTets = numDiv * nx * ny * nz;
+
+ std::vector<TetrahedralMesh::TetraArray> vertices;
+ vertices.resize(numTets);
+ cnt = 0;
+ std::vector<size_t> indx(8);
+
+ for (size_t k = 0; k < nz; ++k)
+ {
+ for (size_t j = 0; j < ny; ++j)
+ {
+ for (size_t i = 0; i < nx; ++i)
+ {
+ indx[3] = i + j * (nx + 1) + k * (nx + 1) * (ny + 1);
+ indx[2] = indx[3] + 1;
+ indx[0] = indx[3] + nx + 1;
+ indx[1] = indx[0] + 1;
+ indx[4] = indx[0] + (nx + 1) * (ny + 1);
+ indx[5] = indx[1] + (nx + 1) * (ny + 1);
+ indx[6] = indx[2] + (nx + 1) * (ny + 1);
+ indx[7] = indx[3] + (nx + 1) * (ny + 1);
+ vertices[cnt + 0] = { indx[0], indx[2], indx[3], indx[6] };
+ vertices[cnt + 1] = { indx[0], indx[3], indx[7], indx[6] };
+ vertices[cnt + 2] = { indx[0], indx[7], indx[4], indx[6] };
+ vertices[cnt + 3] = { indx[0], indx[5], indx[6], indx[4] };
+ vertices[cnt + 4] = { indx[1], indx[5], indx[6], indx[0] };
+ vertices[cnt + 5] = { indx[1], indx[6], indx[2], indx[0] };
+ cnt += numDiv;
+ }
+ }
+ }
+
+ auto mesh = std::make_shared<TetrahedralMesh>();
+ // auto mesh = std::unique_ptr<TetrahedralMesh>(new TetrahedralMesh());
+
+ mesh->initialize(coords, vertices);
+ return mesh;
+}
+
+std::shared_ptr<TetrahedralMesh>
+GeometryUtils::createTetrahedralMeshCover(const SurfaceMesh& surfMesh,
+ const size_t nx,
+ const size_t ny,
+ const size_t nz)
+{
+ Vec3d aabbMin, aabbMax;
+
+ // create a background mesh
+ surfMesh.computeBoundingBox(aabbMin, aabbMax, 1. /*percentage padding*/);
+ auto uniformMesh = createUniformMesh(aabbMin, aabbMax, nx, ny, nz);
+
+ // ray-tracing
+ const auto& coords = uniformMesh->getVertexPositions();
+ // auto insideSurfMesh = markPointsInsideAndOut(surfMesh, coords);
+ const auto insideSurfMesh = markPointsInsideAndOut(surfMesh, coords, nx + 1, ny + 1, nz + 1);
+
+ // label elements
+ std::vector<bool> validTet(uniformMesh->getNumTetrahedra(), false);
+ std::vector<bool> validVtx(uniformMesh->getNumVertices(), false);
+
+ // TetrahedralMesh::WeightsArray weights = {0.0, 0.0, 0.0, 0.0};
+ const Vec3d h = { (aabbMax[0] - aabbMin[0]) / nx,
+ (aabbMax[1] - aabbMin[1]) / ny,
+ (aabbMax[2] - aabbMin[2]) / nz };
+
+ // a customized approach to find the enclosing tet for each surface points
+ /// \todo can be parallelized by make NUM_THREADS copies of validTet, or use atomic op on validTet
+ auto labelEnclosingTet = [&aabbMin, &h, nx, ny, nz, &uniformMesh, &validTet](const Vec3d& xyz)
+ {
+ const size_t idX = (xyz[0] - aabbMin[0]) / h[0];
+ const size_t idY = (xyz[1] - aabbMin[1]) / h[1];
+ const size_t idZ = (xyz[2] - aabbMin[2]) / h[2];
+ const size_t hexId = idX + idY * nx + idZ * nx * ny;
+
+ // the index range of tets inside the enclosing hex
+ const int numDiv = 6;
+ const size_t tetId0 = numDiv * hexId;
+ const size_t tetId1 = tetId0 + numDiv;
+
+ static TetrahedralMesh::WeightsArray weights = { 0.0, 0.0, 0.0, 0.0 };
+
+ // loop over the tets to find the enclosing tets
+ for (size_t id = tetId0; id < tetId1; ++id)
+ {
+ if (validTet[id])
+ {
+ continue;
+ }
+ uniformMesh->computeBarycentricWeights(id, xyz, weights);
+
+ if ((weights[0] >= 0) && (weights[1] >= 0) && (weights[2] >= 0) && (weights[3] >= 0))
+ {
+ validTet[id] = true;
+ break;
+ }
+ }
+ };
+
+ auto labelEnclosingTetOfVertices = [&surfMesh, &uniformMesh, &aabbMin, &h, nx, ny, nz, &labelEnclosingTet, &validTet](const size_t i)
+ {
+ const auto& xyz = surfMesh.getVertexPosition(i);
+ labelEnclosingTet(xyz);
+ };
+
+ for (size_t i = 0; i < validTet.size(); ++i)
+ {
+ const auto& verts = uniformMesh->getTetrahedronVertices(i);
+ if (insideSurfMesh[verts[0]]
+ || insideSurfMesh[verts[1]]
+ || insideSurfMesh[verts[2]]
+ || insideSurfMesh[verts[3]])
+ {
+ validTet[i] = true;
+ }
+ }
+
+ // find the enclosing tets of a group of points on a surface triangle
+ auto labelEnclosingTetOfInteriorPnt = [&surfMesh, &labelEnclosingTet](const size_t fid)
+ {
+ auto verts = surfMesh.getTrianglesVertices()[fid];
+ const auto& vtx0 = surfMesh.getVertexPosition(verts[0]);
+ const auto& vtx1 = surfMesh.getVertexPosition(verts[1]);
+ const auto& vtx2 = surfMesh.getVertexPosition(verts[2]);
+ std::vector<Vec3d> pnts(12);
+
+ pnts[0] = 0.75 * vtx0 + 0.25 * vtx1;
+ pnts[1] = 0.50 * vtx0 + 0.50 * vtx1;
+ pnts[2] = 0.25 * vtx0 + 0.75 * vtx1;
+ pnts[3] = 0.75 * vtx1 + 0.25 * vtx2;
+ pnts[4] = 0.50 * vtx1 + 0.50 * vtx2;
+ pnts[5] = 0.25 * vtx1 + 0.75 * vtx2;
+ pnts[6] = 0.75 * vtx2 + 0.25 * vtx0;
+ pnts[7] = 0.50 * vtx2 + 0.50 * vtx0;
+ pnts[8] = 0.25 * vtx2 + 0.75 * vtx0;
+ pnts[9] = 2.0 / 3.0 * pnts[0] + 1.0 / 3.0 * pnts[5];
+ pnts[10] = 0.5 * (pnts[1] + pnts[4]);
+ pnts[11] = 0.5 * (pnts[4] + pnts[7]);
+
+ for (size_t i = 0; i < pnts.size(); ++i)
+ {
+ labelEnclosingTet(pnts[i]);
+ }
+ };
+
+ // enclose all vertices
+ for (size_t i = 0; i < surfMesh.getNumVertices(); ++i)
+ {
+ labelEnclosingTetOfVertices(i);
+ }
+
+ // enclose some interior points on triangles
+ for (size_t i = 0; i < surfMesh.getNumTriangles(); ++i)
+ {
+ labelEnclosingTetOfInteriorPnt(i);
+ }
+
+ size_t numElems = 0;
+ for (size_t i = 0; i < validTet.size(); ++i)
+ {
+ const auto& verts = uniformMesh->getTetrahedronVertices(i);
+ if (validTet[i])
+ {
+ validVtx[verts[0]] = true;
+ validVtx[verts[1]] = true;
+ validVtx[verts[2]] = true;
+ validVtx[verts[3]] = true;
+ ++numElems;
+ }
+ }
+
+ // discard useless vertices, and build old-to-new index map
+ size_t numVerts = 0;
+ for (auto b : validVtx)
+ {
+ if (b)
+ {
+ ++numVerts;
+ }
+ }
+
+ StdVectorOfVec3d newCoords(numVerts);
+ std::vector<size_t> oldToNew(coords.size(), std::numeric_limits<size_t>::max());
+
+ size_t cnt = 0;
+ for (size_t i = 0; i < validVtx.size(); ++i)
+ {
+ if (validVtx[i])
+ {
+ newCoords[cnt] = coords[i];
+ oldToNew[i] = cnt;
+ ++cnt;
+ }
+ }
+
+ // update tet-to-vert connectivity
+ std::vector<TetrahedralMesh::TetraArray> newIndices(numElems);
+ cnt = 0;
+ for (size_t i = 0; i < uniformMesh->getNumTetrahedra(); ++i)
+ {
+ if (validTet[i])
+ {
+ const auto oldIndices = uniformMesh->getTetrahedronVertices(i);
+
+ newIndices[cnt][0] = oldToNew[oldIndices[0]];
+ newIndices[cnt][1] = oldToNew[oldIndices[1]];
+ newIndices[cnt][2] = oldToNew[oldIndices[2]];
+ newIndices[cnt][3] = oldToNew[oldIndices[3]];
+
+ ++cnt;
+ }
+ }
+
+ // ready to create the final mesh
+ auto mesh = std::make_shared<TetrahedralMesh>();
+ mesh->initialize(newCoords, newIndices);
+
+ return mesh;
+}
+
+template<typename NeighborContainer>
+std::vector<size_t>
+GeometryUtils::reorderConnectivity(const std::vector<NeighborContainer>& neighbors, const MeshNodeRenumberingStrategy& method)
+{
+ switch (method)
+ {
+ case (MeshNodeRenumberingStrategy::ReverseCuthillMckee):
+ return RCM(neighbors);
+ default:
+ LOG(WARNING) << "Unrecognized reorder method; using RCM instead";
+ return RCM(neighbors);
+ }
+}
+
+template<typename ElemConn>
+std::vector<size_t>
+GeometryUtils::reorderConnectivity(const std::vector<ElemConn>& conn, const size_t numVerts, const MeshNodeRenumberingStrategy& method)
+{
+ switch (method)
+ {
+ case (MeshNodeRenumberingStrategy::ReverseCuthillMckee):
+ return RCM(conn, numVerts);
+ default:
+ LOG(WARNING) << "Unrecognized reorder method; using Reverse Cuthill-Mckee strategy instead";
+ return RCM(conn, numVerts);
+ }
+}
+} // namespace imstk
+
+template std::vector<size_t> imstk::GeometryUtils::reorderConnectivity<std::set<size_t>>(const std::vector<std::set<size_t>>&, const GeometryUtils::MeshNodeRenumberingStrategy&);
+
+template std::vector<size_t> imstk::GeometryUtils::reorderConnectivity<std::array<size_t, 4>>(const std::vector<std::array<size_t, 4>>&, const size_t, const MeshNodeRenumberingStrategy&);
diff --git a/Source/Geometry/imstkGeometryUtilities.h b/Source/Geometry/imstkGeometryUtilities.h
index 55db04a8b6cd86129df266de546442a0ecfb2832..717577c33fddcf9a37f792712ab03b89d8280136 100644
--- a/Source/Geometry/imstkGeometryUtilities.h
+++ b/Source/Geometry/imstkGeometryUtilities.h
@@ -21,9 +21,16 @@
#pragma once
+#include "g3log/g3log.hpp"
#include "imstkMath.h"
-#include <memory>
#include <vtkSmartPointer.h>
+#include "imstkParallelUtils.h"
+
+#include <memory>
+#include <numeric>
+#include <queue>
+#include <unordered_set>
+#include <set>
class vtkCellArray;
class vtkPolyData;
@@ -153,5 +160,53 @@ std::unique_ptr<SurfaceMesh> linearSubdivideSurfaceMesh(const SurfaceMesh& surfa
/// for more details
///
std::unique_ptr<SurfaceMesh> loopSubdivideSurfaceMesh(const SurfaceMesh& surfaceMesh, const int numSubdivisions = 1);
+
+///
+/// \brief Create a tetrahedral mesh based on a uniform Cartesian mesh
+/// \param aabbMin the small conner of a box
+/// \param aabbMax the large conner of a box
+/// \param nx number of elements in the x-direction
+/// \param ny number of elements in the y-direction
+/// \param nz number of elements in the z-direction
+///
+/// \note Refer: Dompierre, Julien & Labbé, Paul & Vallet, Marie-Gabrielle & Camarero, Ricardo. (1999).
+/// How to Subdivide Pyramids, Prisms, and Hexahedra into Tetrahedra.. 195-204.
+std::shared_ptr<TetrahedralMesh> createUniformMesh(const Vec3d& aabbMin, const Vec3d& aabbMax, const size_t nx, const size_t ny, const size_t nz);
+
+///
+/// \brief Create a tetrahedral mesh cover
+///
+std::shared_ptr<TetrahedralMesh> createTetrahedralMeshCover(const SurfaceMesh& surfMesh, const size_t nx, const size_t ny, size_t nz);
+
+///
+/// \brief Enumeration for reordering method
+///
+enum class MeshNodeRenumberingStrategy
+{
+ ReverseCuthillMckee // Reverse Cuthill-Mckee
};
-}
+
+///
+/// \brief Reorder indices in a connectivity to reduce bandwidth
+///
+/// \param[in] neighbors array of neighbors of each vertex; eg, neighbors[i] is an object containing all neighbors of vertex-i
+/// \param[i] method reordering method; see \ref ReorderMethod
+///
+/// \return the permutation vector that map from new indices to old indices
+///
+template<typename NeighborContainer>
+std::vector<size_t> reorderConnectivity(const std::vector<NeighborContainer>& neighbors, const MeshNodeRenumberingStrategy& method = MeshNodeRenumberingStrategy::ReverseCuthillMckee);
+
+///
+/// \brief Reorder using Reverse Cuthill-Mckee
+///
+/// \param[in] conn element-to-vertex connectivity
+/// \param[in] numVerts number of vertices
+/// \param[in] method reordering method; see \ref ReorderMethod
+///
+/// \return the permutation vector that maps from new indices to old indices
+///
+template<typename ElemConn>
+std::vector<size_t> reorderConnectivity(const std::vector<ElemConn>& conn, const size_t numVerts, const MeshNodeRenumberingStrategy& method = MeshNodeRenumberingStrategy::ReverseCuthillMckee);
+} // namespace GeometryUtils
+} // namespace imstk