ComputeAdjacencyList Assignment Open3D

cpp/open3d/t/geometry/TriangleMesh.cpp

@@ -18,6 +18,7 @@
 #include <Eigen/Core>
 #include <string>
 #include <unordered_map>
+#include <set>
 
 #include "open3d/core/CUDAUtils.h"
 #include "open3d/core/Device.h"
@@ -65,6 +66,54 @@ TriangleMesh::TriangleMesh(const core::Tensor &vertex_positions,
     SetVertexPositions(vertex_positions);
     SetTriangleIndices(triangle_indices);
 }
+
+std::pair<core::Tensor, core::Tensor> TriangleMesh::ComputeAdjacencyList() {
+
+    if (IsEmpty()) {
+        utility::LogWarning("TriangleMesh is empty. No attributes computed!");
+        return std::make_pair(core::Tensor(), core::Tensor());
+    }
+
+    if(!HasTriangleIndices()) {
+        utility::LogWarning("TriangleMesh has no Indices. No attributes computed!");
+        return std::make_pair(core::Tensor(), core::Tensor());
+    }
+
+    core::Tensor tris_cpu =
+            GetTriangleIndices().To(core::Device()).Contiguous();
+
+    std::unordered_map< size_t, std::set<size_t> > adjacencyList;
+    auto insertEdge = [&adjacencyList](size_t s, size_t t){
+        adjacencyList[s].insert(t);
+    };
+
+    for(int idx = 0; idx < tris_cpu.GetLength(); idx++){
+        auto triangle_tensor = tris_cpu[idx];
+        auto *triangle = triangle_tensor.GetDataPtr<int>();
+        insertEdge(triangle[0], triangle[1]);
+        insertEdge(triangle[1], triangle[2]);
+        insertEdge(triangle[2], triangle[0]);
+    }
+
+    int num_vertices = GetVertexPositions().GetLength();
+    core::Tensor adjst = core::Tensor::Zeros({num_vertices+1}, core::Dtype::Int64);
+
+    int num_edges = tris_cpu.GetLength() * 3;
+    core::Tensor adjv = core::Tensor::Zeros({num_edges}, core::Dtype::Int64);
+
+    long prev_nnz = 0;
+    for(int idx = 1; idx <= num_vertices; idx++){
+        adjst[idx] = adjst[idx-1] + static_cast<long>(adjacencyList[idx-1].size());
+
+        int i = 0;
+        for(auto x : adjacencyList[idx-1]){
+            adjv[ prev_nnz + i] = x;
+            i++;
+        }
+        prev_nnz += static_cast<long>(adjacencyList[idx-1].size());
+    }
+    return std::make_pair(adjv, adjst);
+}
 
 std::string TriangleMesh::ToString() const {
     size_t num_vertices = 0;

cpp/open3d/t/geometry/TriangleMesh.h

@@ -828,6 +828,8 @@ class TriangleMesh : public Geometry, public DrawableGeometry {
     TriangleMesh BooleanDifference(const TriangleMesh &mesh,
                                    double tolerance = 1e-6) const;
 
+    std::pair<core::Tensor, core::Tensor> ComputeAdjacencyList();
+
     /// Create an axis-aligned bounding box from vertex attribute "positions".
     AxisAlignedBoundingBox GetAxisAlignedBoundingBox() const;
 

cpp/pybind/t/geometry/trianglemesh.cpp

@@ -734,11 +734,14 @@ This function always uses the CPU device.
 
         o3d.visualization.draw([{'name': 'difference', 'geometry': ans}])
 )");
-
+    triangle_mesh.def("compute_adjacency_list", &TriangleMesh::ComputeAdjacencyList,
+    "Return Mesh Adjacency Matrix in CSR format using Triangle indices attribute.");
+
     triangle_mesh.def("get_axis_aligned_bounding_box",
                       &TriangleMesh::GetAxisAlignedBoundingBox,
                       "Create an axis-aligned bounding box from vertex "
                       "attribute 'positions'.");
+
     triangle_mesh.def("get_oriented_bounding_box",
                       &TriangleMesh::GetOrientedBoundingBox,
                       "Create an oriented bounding box from vertex attribute "

cpp/tests/t/geometry/TriangleMesh.cpp

@@ -25,6 +25,70 @@ INSTANTIATE_TEST_SUITE_P(TriangleMesh,
                          TriangleMeshPermuteDevices,
                          testing::ValuesIn(PermuteDevices::TestCases()));
 
+TEST_P(TriangleMeshPermuteDevices, ComputeAdjacencyList_emptyMesh) {
+//Test the interface and the case when mesh is empty
+
+    t::geometry::TriangleMesh empty_mesh;
+
+    auto listCSR = empty_mesh.ComputeAdjacencyList();
+
+    core::Tensor adjacent_vertex = listCSR.first;
+    core::Tensor adjacent_index_start = listCSR.second;
+    EXPECT_TRUE(adjacent_vertex.GetLength() == 0);
+    EXPECT_TRUE(adjacent_index_start.GetLength() == 0);
+}
+
+TEST_P(TriangleMeshPermuteDevices, ComputeAdjacencyList_matchValues) {
+//Test the actual values computed in the function
+
+    core::Device device = GetParam();
+    core::Dtype float_dtype_custom = core::Float64;
+    core::Dtype int_dtype_custom = core::Int32;
+
+    t::geometry::TriangleMesh mesh =
+            t::geometry::TriangleMesh::CreateTetrahedron(
+                    2, float_dtype_custom, int_dtype_custom, device);
+
+    auto listCSR = mesh.ComputeAdjacencyList();
+    core::Tensor adjv = listCSR.first;
+    core::Tensor adjst = listCSR.second;
+
+    EXPECT_TRUE( adjv.GetLength() > 0);
+    EXPECT_TRUE( adjst.GetLength() > 0);
+
+    core::Tensor csr_col = core::Tensor::Init<int64_t>(
+            {1, 2, 3, 0, 2, 3, 0, 1, 3, 0, 1, 2}, device);
+
+    core::Tensor csr_row_idx = core::Tensor::Init<int64_t>(
+            {0, 3, 6, 9, 12}, device);
+
+    EXPECT_EQ(adjv.GetLength(), csr_col.GetLength());
+    EXPECT_EQ(adjst.GetLength(), csr_row_idx.GetLength());
+    EXPECT_TRUE(adjv.AllEqual(csr_col));
+    EXPECT_TRUE(adjst.AllEqual(csr_row_idx));
+}
+
+TEST_P(TriangleMeshPermuteDevices, ComputeAdjacencyList_expectedBehaviour) {
+//On a larger mesh, test the interface and few expected properties
+
+    core::Device device = GetParam();
+    core::Dtype float_dtype_custom = core::Float64;
+    core::Dtype int_dtype_custom = core::Int32;
+
+    t::geometry::TriangleMesh mesh =
+            t::geometry::TriangleMesh::CreateIcosahedron(
+                    2, float_dtype_custom, int_dtype_custom, device);
+
+
+    auto listCSR = mesh.ComputeAdjacencyList();
+    core::Tensor adjv = listCSR.first;
+    core::Tensor adjst = listCSR.second;
+
+    EXPECT_TRUE( adjv.GetLength() > 0);
+    EXPECT_TRUE( adjst.GetLength() > 0);
+    EXPECT_EQ(adjst.ToFlatVector<int64_t>()[adjst.GetLength()-1], adjv.GetLength());
+}
+
 TEST_P(TriangleMeshPermuteDevices, DefaultConstructor) {
     t::geometry::TriangleMesh mesh;
 

python/test/t/geometry/test_trianglemesh.py

@@ -241,6 +241,66 @@ def test_create_octahedron(device):
     assert octahedron_custom.vertex.positions.allclose(vertex_positions_custom)
     assert octahedron_custom.triangle.indices.allclose(triangle_indices_custom)
 
+@pytest.mark.parametrize("device", list_devices())
+def test_compute_adjacency_list(device):
+    # Test with custom parameters.
+    mesh = o3d.t.geometry.TriangleMesh.create_icosahedron(
+        2, o3c.float64, o3c.int32, device)
+    adjv, adjst = mesh.compute_adjacency_list()
+
+    # Check the number of edges at the end of row_index as per CRS format
+    assert adjst[-1] == len(adjv)
+
+    num_vertices = len(adjst)-1
+    # Create a adjacency matrix 
+    adjacencyMatrix = np.zeros((num_vertices, num_vertices))
+
+    for s in range(num_vertices):
+        start = adjst[s].item()
+        end = adjst[s+1].item()
+        for v in range(start, end):
+            t = adjv[v].item()
+            adjacencyMatrix[s,t] = 1
+
+    # Number of edges
+    assert len(adjv) == adjacencyMatrix.sum()
+
+    # Adjacency Matrix should be symmetric
+    assert np.array_equal(adjacencyMatrix, adjacencyMatrix.T)
+
+    #Triangle faces from computed adjacency matrix should match
+    actual_indices = [
+            [0, 1, 4],
+            [0, 1, 5],
+            [0, 4, 8],
+            [0, 5, 11],
+            [0, 8, 11],
+            [1, 4, 9],
+            [1, 5, 10],
+            [1, 9, 10],
+            [2, 3, 6],
+            [2, 3, 7],
+            [2, 6, 10],
+            [2, 7, 9],
+            [2, 9, 10],
+            [3, 6, 11],
+            [3, 7, 8],
+            [3, 8, 11],
+            [4, 7, 8],
+            [4, 7, 9],
+            [5, 6, 10],
+            [5, 6, 11]
+        ]
+
+    computed_triangles = []
+    for i in range(num_vertices):
+        for j in range(i+1, num_vertices):
+            for k in range(j+1, num_vertices):
+                if (adjacencyMatrix[i,j] + adjacencyMatrix[j,k] + adjacencyMatrix[k,i] == 3):
+                    computed_triangles.append([i,j,k])
+
+    assert len(computed_triangles) == len(actual_indices)
+    assert np.array_equal(np.array(actual_indices,int), np.array(computed_triangles,int))
 
 @pytest.mark.parametrize("device", list_devices())
 def test_create_icosahedron(device):