ENH: added Reverse Cuthill-Mckee

Examples/RCM/CMakeLists.txt

+###########################################################################
+#
+# 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()
+
+

Examples/RCM/RCM.cpp

+/*=========================================================================
+
+   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 "imstkGeometry.h"
+
+#include "bandwidth.h"
+
+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) {
+     
+    // a 2D Cartesian mesh
+    {
+        auto p = createConn();
+        testRCM(p.first, p.second);
+    }
+
+    // dragon mesh
+    {
+        auto tetMesh = std::dynamic_pointer_cast<TetrahedralMesh>(MeshIO::read(iMSTK_DATA_ROOT "/asianDragon/asianDragon.veg"));
+        auto conn = tetMesh->getTetrahedraVertices(); 
+        auto numVerts = tetMesh->getNumVertices();
+        testRCM(conn, numVerts);
+    }
+
+    // a mesh cover
+    {
+        const size_t nx = 80;
+        const size_t ny = 40;
+        const size_t nz = 60;
+        auto surfMesh = std::dynamic_pointer_cast<SurfaceMesh>(MeshIO::read(iMSTK_DATA_ROOT "/asianDragon/asianDragon.obj"));
+        auto tetMesh = TetrahedralMesh::createTetrahedralMeshCover(surfMesh, nx, ny, nz);
+        auto conn = tetMesh->getTetrahedraVertices(); 
+        auto numVerts = tetMesh->getNumVertices();
+        testRCM(conn, numVerts);
+    }
+
+    return 0;
+}
+
+template <typename ElemConn>
+void testRCM(const std::vector<ElemConn>& conn, const size_t numVerts) 
+{
+    std::cout << "bandwidth_old = " << bandwidth(conn, numVerts) << std::endl;
+
+    // new-to-old permutation
+    auto perm = GeometricUtils::RCM(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);
+        invPerm[perm[i]] = i;
+    }
+
+    auto newConn = conn;
+
+    for (auto& vertices : newConn)
+    {
+        for (auto& vid : vertices)
+        {
+            CHECK(vid < numVerts);
+            vid = invPerm[vid];
+        }
+    }
+
+    std::cout << "bandwidth_new = " << 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);
+}

Examples/RCM/bandwidth.h

+///
+/// \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);
+}

Source/Geometry/imstkGeometry.h

@@ -27,6 +27,11 @@
 #include "imstkParallelUtils.h"
 
 #include <string>
+#include <memory>
+#include <numeric>
+#include <unordered_set>
+#include <set>
+
 
 namespace imstk
 {
@@ -218,4 +223,207 @@ protected:
     RigidTransform3d m_transform = RigidTransform3d::Identity(); ///> Transformation matrix
     double m_scaling = 1.0;
 };
+
+class GeometricUtils {
+public:
+    ///
+    /// \brief Reverse Cuthill-Mckee; return the permutation vector that maps from new indices to old indices
+    //
+    /// \param conn element-to-vertex connectivity
+    /// \param numVerts number of vertices
+    ///
+    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;
+        buildVertToVert(conn, numVerts, vertToVert);
+        return RCM(vertToVert);
+    }
+
+    ///
+    /// \brief Reverse Cuthill-Mckee; returns the permutation vector that map from new indices to old indices
+    ///
+    /// \param neighbors array of neighbors of each vertex; eg, neighbors[i] is an object containing all neighbors of vertex-i
+    ///
+    template <typename NBR>
+    static std::vector<size_t>
+    RCM(const std::vector<NBR>& 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);
+
+        // 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 moveIntoRAndNbrIntoQ = [&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;
+
+        ///
+        /// \brief pop a vertex that is not ordered from \p P
+        ///
+        auto popFromP = [&pCur, &isInP]() {
+            for (size_t vid = pCur; vid < isInP.size(); ++vid)
+            {
+                if (isInP[vid])
+                {
+                    isInP[vid] = false;
+                    pCur       = vid;
+                    return vid;
+                }
+            }
+            return INVALID;
+        };
+
+        // main loop
+        while (true)
+        {
+            std::size_t parent = popFromP();
+            if (parent == INVALID)
+            {
+                break;
+            }
+
+            moveIntoRAndNbrIntoQ(parent);
+
+            while (!Q.empty())
+            {
+                parent = Q.front();
+                Q.pop();
+                moveIntoRAndNbrIntoQ(parent);
+            }
+
+            // here we have empty Q
+        }
+
+        CHECK(pos == numVerts);
+
+        std::reverse(R.begin(), R.end());
+        return R;
+    }
+
+private:
+    ///
+    /// \brief Build the vertex-to-vertex connectivity
+    /// 
+    /// \param conn element-to-vertex connectivity
+    /// \param numVerts number of vertices
+    /// \param vertToVert the vertex-to-vertex connectivity on return
+    ///
+    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);
+        }
+    }
+
+};  // class GeometricUtils 
+
+
 } //imstk