From 51c88ea31dcde420332ef47bf04159806f252c1b Mon Sep 17 00:00:00 2001 From: root Date: Wed, 17 Jun 2020 11:50:19 +0530 Subject: [PATCH] Add an extra condition --- gdist_c_api.cpp | 284 ++++++++++++++++++++++---------------------- tests/test_gdist.py | 130 ++++++++++---------- 2 files changed, 209 insertions(+), 205 deletions(-) diff --git a/gdist_c_api.cpp b/gdist_c_api.cpp index 7c0e4e7..5d038e9 100644 --- a/gdist_c_api.cpp +++ b/gdist_c_api.cpp @@ -1,142 +1,142 @@ -#include -#include - -#include "geodesic_library/geodesic_algorithm_exact.h" - - -void compute_gdist_impl( - unsigned number_of_vertices, - unsigned number_of_triangles, - double *vertices, - int *triangles, - unsigned number_of_source_indices, - unsigned number_of_target_indices, - unsigned *source_indices_array, - unsigned *target_indices_array, - double *distance, - double distance_limit -) { - - std::vector points (vertices, vertices + number_of_vertices); - std::vector faces (triangles, triangles + number_of_triangles); - std::vector source_indices (source_indices_array, source_indices_array + number_of_source_indices); - std::vector target_indices (target_indices_array, target_indices_array + number_of_target_indices); - - geodesic::Mesh mesh; - mesh.initialize_mesh_data(points, faces); // create internal mesh data structure including edges - - geodesic::GeodesicAlgorithmExact algorithm(&mesh); // create exact algorithm for the mesh - - std::vector all_sources, stop_points; - - for (unsigned i = 0; i < number_of_source_indices; ++i) { - all_sources.push_back(geodesic::SurfacePoint(&mesh.vertices()[source_indices[i]])); - } - - for (unsigned i = 0; i < number_of_target_indices; ++i) { - stop_points.push_back(geodesic::SurfacePoint(&mesh.vertices()[target_indices[i]])); - } - - algorithm.propagate(all_sources, distance_limit, &stop_points); - - for (unsigned i = 0; i < stop_points.size(); ++i) { - algorithm.best_source(stop_points[i], distance[i]); - } -} - -double* local_gdist_matrix_impl( - unsigned number_of_vertices, - unsigned number_of_triangles, - double *vertices, - unsigned *triangles, - unsigned *sparse_matrix_size, - double max_distance -) { - std::vector points (vertices, vertices + number_of_vertices); - std::vector faces (triangles, triangles + number_of_triangles); - - geodesic::Mesh mesh; - mesh.initialize_mesh_data(points, faces); // create internal mesh data structure including edges - geodesic::GeodesicAlgorithmExact algorithm(&mesh); // create exact algorithm for the mesh - std::vector rows_vector, columns_vector; - std::vector data_vector; - - double distance = 0; - - std::vector targets(number_of_vertices), source; - - for (unsigned i = 0; i < number_of_vertices; ++i) { - targets[i] = geodesic::SurfacePoint(&mesh.vertices()[i]); - } - for (unsigned i = 0; i < number_of_vertices / 3; ++i) { - source.push_back(geodesic::SurfacePoint(&mesh.vertices()[i])); - algorithm.propagate(source, max_distance, NULL); - source.pop_back(); - for (unsigned j = 0; j < number_of_vertices / 3; ++j) { - algorithm.best_source(targets[j], distance); - if (distance != geodesic::GEODESIC_INF && distance != 0) { - rows_vector.push_back(i); - columns_vector.push_back(j); - data_vector.push_back(distance); - } - } - } - - double *data; - data = new double[3 * rows_vector.size()]; - - *sparse_matrix_size = rows_vector.size(); - - std::copy(rows_vector.begin(), rows_vector.end(), data); - std::copy(columns_vector.begin(), columns_vector.end(), data + data_vector.size()); - std::copy(data_vector.begin(), data_vector.end(), data + 2 * data_vector.size()); - - return data; -} - - -extern "C" { - double* compute_gdist( - unsigned number_of_vertices, - unsigned number_of_triangles, - double *vertices, - int *triangles, - unsigned number_of_source_indices, - unsigned number_of_target_indices, - unsigned *source_indices_array, - unsigned *target_indices_array, - double *distance, - double distance_limit - ) { - compute_gdist_impl( - number_of_vertices, - number_of_triangles, - vertices, - triangles, - number_of_source_indices, - number_of_target_indices, - source_indices_array, - target_indices_array, - distance, - distance_limit - ); - } - - double* local_gdist_matrix( - unsigned number_of_vertices, - unsigned number_of_triangles, - double *vertices, - unsigned *triangles, - unsigned *sparse_matrix_size, - double max_distance - ) { - return local_gdist_matrix_impl( - number_of_vertices, - number_of_triangles, - vertices, - triangles, - sparse_matrix_size, - max_distance - ); - } -}; +#include +#include + +#include "geodesic_library/geodesic_algorithm_exact.h" + + +void compute_gdist_impl( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + int *triangles, + unsigned number_of_source_indices, + unsigned number_of_target_indices, + unsigned *source_indices_array, + unsigned *target_indices_array, + double *distance, + double distance_limit +) { + + std::vector points (vertices, vertices + number_of_vertices); + std::vector faces (triangles, triangles + number_of_triangles); + std::vector source_indices (source_indices_array, source_indices_array + number_of_source_indices); + std::vector target_indices (target_indices_array, target_indices_array + number_of_target_indices); + + geodesic::Mesh mesh; + mesh.initialize_mesh_data(points, faces); // create internal mesh data structure including edges + + geodesic::GeodesicAlgorithmExact algorithm(&mesh); // create exact algorithm for the mesh + + std::vector all_sources, stop_points; + + for (unsigned i = 0; i < number_of_source_indices; ++i) { + all_sources.push_back(geodesic::SurfacePoint(&mesh.vertices()[source_indices[i]])); + } + + for (unsigned i = 0; i < number_of_target_indices; ++i) { + stop_points.push_back(geodesic::SurfacePoint(&mesh.vertices()[target_indices[i]])); + } + + algorithm.propagate(all_sources, distance_limit, &stop_points); + + for (unsigned i = 0; i < stop_points.size(); ++i) { + algorithm.best_source(stop_points[i], distance[i]); + } +} + +double* local_gdist_matrix_impl( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + unsigned *triangles, + unsigned *sparse_matrix_size, + double max_distance +) { + std::vector points (vertices, vertices + number_of_vertices); + std::vector faces (triangles, triangles + number_of_triangles); + + geodesic::Mesh mesh; + mesh.initialize_mesh_data(points, faces); // create internal mesh data structure including edges + geodesic::GeodesicAlgorithmExact algorithm(&mesh); // create exact algorithm for the mesh + std::vector rows_vector, columns_vector; + std::vector data_vector; + + double distance = 0; + + std::vector targets(number_of_vertices), source; + + for (unsigned i = 0; i < number_of_vertices; ++i) { + targets[i] = geodesic::SurfacePoint(&mesh.vertices()[i]); + } + for (unsigned i = 0; i < number_of_vertices / 3; ++i) { + source.push_back(geodesic::SurfacePoint(&mesh.vertices()[i])); + algorithm.propagate(source, max_distance, NULL); + source.pop_back(); + for (unsigned j = 0; j < number_of_vertices / 3; ++j) { + algorithm.best_source(targets[j], distance); + if (distance != geodesic::GEODESIC_INF && distance != 0 && distance <= max_distance) { + rows_vector.push_back(i); + columns_vector.push_back(j); + data_vector.push_back(distance); + } + } + } + + double *data; + data = new double[3 * rows_vector.size()]; + + *sparse_matrix_size = rows_vector.size(); + + std::copy(rows_vector.begin(), rows_vector.end(), data); + std::copy(columns_vector.begin(), columns_vector.end(), data + data_vector.size()); + std::copy(data_vector.begin(), data_vector.end(), data + 2 * data_vector.size()); + + return data; +} + + +extern "C" { + double* compute_gdist( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + int *triangles, + unsigned number_of_source_indices, + unsigned number_of_target_indices, + unsigned *source_indices_array, + unsigned *target_indices_array, + double *distance, + double distance_limit + ) { + compute_gdist_impl( + number_of_vertices, + number_of_triangles, + vertices, + triangles, + number_of_source_indices, + number_of_target_indices, + source_indices_array, + target_indices_array, + distance, + distance_limit + ); + } + + double* local_gdist_matrix( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + unsigned *triangles, + unsigned *sparse_matrix_size, + double max_distance + ) { + return local_gdist_matrix_impl( + number_of_vertices, + number_of_triangles, + vertices, + triangles, + sparse_matrix_size, + max_distance + ); + } +}; diff --git a/tests/test_gdist.py b/tests/test_gdist.py index 5c3feb3..3c0472c 100644 --- a/tests/test_gdist.py +++ b/tests/test_gdist.py @@ -1,63 +1,67 @@ -import numpy as np - -import gdist - - -class TestComputeGdist(): - def test_flat_triangular_mesh(self): - data = np.loadtxt("data/flat_triangular_mesh.txt", skiprows=1) - vertices = data[0:121].astype(np.float64) - triangles = data[121:].astype(np.int32) - source = np.array([1], dtype=np.int32) - target = np.array([2], dtype=np.int32) - distance = gdist.compute_gdist( - vertices, - triangles, - source_indices=source, - target_indices=target - ) - np.testing.assert_array_almost_equal(distance, [0.2]) - - def test_hedgehog_mesh(self): - data = np.loadtxt("data/hedgehog_mesh.txt", skiprows=1) - vertices = data[0:300].astype(np.float64) - triangles = data[300:].astype(np.int32) - source = np.array([0], dtype=np.int32) - target = np.array([1], dtype=np.int32) - distance = gdist.compute_gdist( - vertices, - triangles, - source_indices=source, - target_indices=target - ) - np.testing.assert_array_almost_equal(distance, [1.40522]) - - -class TestLocalGdistMatrix: - def test_flat_triangular_mesh(self): - data = np.loadtxt("data/flat_triangular_mesh.txt", skiprows=1) - vertices = data[0:121].astype(np.float64) - triangles = data[121:].astype(np.int32) - distances = gdist.local_gdist_matrix(vertices, triangles) - epsilon = 1e-6 # the default value used in `assert_array_almost_equal` - # test if the obtained matrix is symmetric - assert (abs(distances - distances.T) > epsilon).nnz == 0 - np.testing.assert_array_almost_equal(distances.toarray()[1][0], 0.2) - # set max distance as 0.3 - distances = gdist.local_gdist_matrix(vertices, triangles, 0.3) - assert np.max(distances) <= 0.3 - - def test_hedgehog_mesh(self): - data = np.loadtxt("data/hedgehog_mesh.txt", skiprows=1) - vertices = data[0:300].astype(np.float64) - triangles = data[300:].astype(np.int32) - distances = gdist.local_gdist_matrix(vertices, triangles) - epsilon = 1e-6 # the default value used in `assert_array_almost_equal` - # test if the obtained matrix is symmetric - assert (abs(distances - distances.T) > epsilon).nnz == 0 - np.testing.assert_array_almost_equal( - distances.toarray()[1][0], 1.40522 - ) - # set max distance as 1.45 - distances = gdist.local_gdist_matrix(vertices, triangles, 1.45) - assert np.max(distances) <= 1.45 +import numpy as np + +import gdist + + +class TestComputeGdist(): + def test_flat_triangular_mesh(self): + data = np.loadtxt("data/flat_triangular_mesh.txt", skiprows=1) + vertices = data[0:121].astype(np.float64) + triangles = data[121:].astype(np.int32) + source = np.array([1], dtype=np.int32) + target = np.array([2], dtype=np.int32) + distance = gdist.compute_gdist( + vertices, + triangles, + source_indices=source, + target_indices=target + ) + np.testing.assert_array_almost_equal(distance, [0.2]) + + def test_hedgehog_mesh(self): + data = np.loadtxt("data/hedgehog_mesh.txt", skiprows=1) + vertices = data[0:300].astype(np.float64) + triangles = data[300:].astype(np.int32) + source = np.array([0], dtype=np.int32) + target = np.array([1], dtype=np.int32) + distance = gdist.compute_gdist( + vertices, + triangles, + source_indices=source, + target_indices=target + ) + np.testing.assert_array_almost_equal(distance, [1.40522]) + + +class TestLocalGdistMatrix: + def test_flat_triangular_mesh(self): + data = np.loadtxt("data/flat_triangular_mesh.txt", skiprows=1) + vertices = data[0:121].astype(np.float64) + triangles = data[121:].astype(np.int32) + distances = gdist.local_gdist_matrix(vertices, triangles) + epsilon = 1e-6 # the default value used in `assert_array_almost_equal` + # test if the obtained matrix is symmetric + assert (abs(distances - distances.T) > epsilon).nnz == 0 + np.testing.assert_array_almost_equal(distances.toarray()[1][0], 0.2) + # set max distance as 0.3 + distances = gdist.local_gdist_matrix(vertices, triangles, 0.3) + # test if the obtained matrix is symmetric + assert (abs(distances - distances.T) > epsilon).nnz == 0 + assert np.max(distances) <= 0.3 + + def test_hedgehog_mesh(self): + data = np.loadtxt("data/hedgehog_mesh.txt", skiprows=1) + vertices = data[0:300].astype(np.float64) + triangles = data[300:].astype(np.int32) + distances = gdist.local_gdist_matrix(vertices, triangles) + epsilon = 1e-6 # the default value used in `assert_array_almost_equal` + # test if the obtained matrix is symmetric + assert (abs(distances - distances.T) > epsilon).nnz == 0 + np.testing.assert_array_almost_equal( + distances.toarray()[1][0], 1.40522 + ) + # set max distance as 1.45 + distances = gdist.local_gdist_matrix(vertices, triangles, 1.45) + # test if the obtained matrix is symmetric + assert (abs(distances - distances.T) > epsilon).nnz == 0 + assert np.max(distances) <= 1.45