diff --git a/.travis.yml b/.travis.yml index b30dad7f..45df2c7c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,10 @@ notifications: email: false +branches: + only: + - trunk + install: - pip3 install . script: @@ -25,7 +29,7 @@ jobs: - ubuntu-toolchain-r-test packages: - g++-6 - before_install: pip3 install pytest-cov codecov + before_install: pip3 install pytest-cov before_script: python3 setup.py clean --all build_ext --force --inplace after_script: - sudo chmod +x cppci.sh @@ -56,6 +60,12 @@ jobs: os: osx osx_image: xcode11.2 language: shell + before_install: + - brew pin gdal; brew pin postgis; brew pin cairo; brew pin glib; brew pin gnutls; brew pin p11-kit; brew pin gnupg; brew pin poppler; + - brew install llvm + env: + - CC="/usr/local/opt/llvm/bin/clang -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" + - CXX="/usr/local/opt/llvm/bin/clang++ -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" - stage: lint language: python install: pip install flake8 diff --git a/README.rst b/README.rst index 573b0288..d8044227 100644 --- a/README.rst +++ b/README.rst @@ -4,6 +4,8 @@ Geodesic Library .. image:: https://travis-ci.com/the-virtual-brain/tvb-gdist.svg?branch=trunk :target: https://travis-ci.com/the-virtual-brain/tvb-gdist +.. image:: https://codecov.io/gh/the-virtual-brain/tvb-gdist/branch/trunk/graph/badge.svg + :target: https://codecov.io/gh/the-virtual-brain/tvb-gdist The **tvb-gdist** module is a Cython interface to a C++ library (https://code.google.com/archive/p/geodesic/) for computing @@ -182,3 +184,13 @@ Notes * In order for the algorithm to work the mesh must not be numbered incorrectly or disconnected or of somehow degenerate. + +* The c++ compiler must have OpenMP installed. This is generally not an issue + on g++ or Microsoft Visual Studio. However, in macOS one may need to install + llvm from brew in order to use OpenMP: ``brew install llvm``. Then, change the + ``CC`` and ``CXX`` environment variables to the following: + + .. code-block:: txt + + CC="/usr/local/opt/llvm/bin/clang -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" + CXX="/usr/local/opt/llvm/bin/clang++ -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" diff --git a/gdist.pyx b/gdist.pyx index 704b7583..8233143f 100644 --- a/gdist.pyx +++ b/gdist.pyx @@ -61,61 +61,19 @@ from libcpp.vector cimport vector ################################################################################ ############ Defininitions to access the C++ geodesic distance library ######### ################################################################################ -cdef extern from "geodesic_mesh_elements.h" namespace "geodesic": - cdef cppclass Vertex: - Vertex() - -cdef extern from "geodesic_mesh_elements.h" namespace "geodesic": - cdef cppclass SurfacePoint: - SurfacePoint() - SurfacePoint(Vertex*) - double& x() - double& y() - double& z() - -cdef extern from "geodesic_mesh.h" namespace "geodesic": - cdef cppclass Mesh: - Mesh() - void initialize_mesh_data(vector[double]&, vector[unsigned]&, bool) - vector[Vertex]& vertices() - cdef extern from "geodesic_utils.h": + cdef cppclass SparseMatrix: + vector[unsigned] rows + vector[unsigned] columns + vector[double] data vector[double] compute_gdist_impl(vector[double], vector[unsigned], vector[unsigned], vector[unsigned], double, bool, bool) - -cdef extern from "geodesic_algorithm_exact.h" namespace "geodesic": - cdef cppclass GeodesicAlgorithmExact: - GeodesicAlgorithmExact(Mesh*) - void propagate(vector[SurfacePoint]&, double, vector[SurfacePoint]*) - unsigned best_source(SurfacePoint&, double&) + SparseMatrix local_gdist_matrix_impl(vector[double], vector[unsigned], double, bool) cdef extern from "geodesic_constants_and_simple_functions.h" namespace "geodesic": double GEODESIC_INF ################################################################################ -cdef get_mesh( - numpy.ndarray[numpy.float64_t, ndim=2] vertices, - numpy.ndarray[numpy.int32_t, ndim=2] triangles, - Mesh &amesh, - bool is_one_indexed -): - # Define C++ vectors to contain the mesh surface components. - cdef vector[double] points - cdef vector[unsigned] faces - - # Map numpy array of mesh "vertices" to C++ vector of mesh "points" - cdef numpy.float64_t coord - for coord in vertices.flatten(): - points.push_back(coord) - - # Map numpy array of mesh "triangles" to C++ vector of mesh "faces" - cdef numpy.int32_t indx - for indx in triangles.flatten(): - faces.push_back(indx) - - amesh.initialize_mesh_data(points, faces, is_one_indexed) - - def compute_gdist(numpy.ndarray[numpy.float64_t, ndim=2] vertices, numpy.ndarray[numpy.int32_t, ndim=2] triangles, numpy.ndarray[numpy.int32_t, ndim=1] source_indices=None, @@ -245,43 +203,26 @@ def local_gdist_matrix(numpy.ndarray[numpy.float64_t, ndim=2] vertices, from the propgate step... """ - cdef Mesh amesh - get_mesh(vertices, triangles, amesh, is_one_indexed) + cdef Py_ssize_t N = vertices.shape[0] - # Define and create object for exact algorithm on that mesh: - cdef GeodesicAlgorithmExact *algorithm = new GeodesicAlgorithmExact(&amesh) + cdef vector[double] points + cdef vector[unsigned] faces - cdef vector[SurfacePoint] source, targets - cdef Py_ssize_t N = vertices.shape[0] - cdef Py_ssize_t k - cdef Py_ssize_t kk - cdef numpy.float64_t distance = 0 - - # Add all vertices as targets - for k in range(N): - targets.push_back(SurfacePoint(&amesh.vertices()[k])) - - rows = [] - columns = [] - data = [] - for k in range(N): - source.push_back(SurfacePoint(&amesh.vertices()[k])) - algorithm.propagate(source, max_distance, NULL) - source.pop_back() - - for kk in range(N): # TODO: Reduce to vertices reached during propagate. - algorithm.best_source(targets[kk], distance) - - if ( - distance is not GEODESIC_INF - and distance is not 0 - and distance <= max_distance - ): - rows.append(k) - columns.append(kk) - data.append(distance) - - return scipy.sparse.csc_matrix((data, (rows, columns)), shape=(N, N)) + for k in vertices.flatten(): + points.push_back(k) + for k in triangles.flatten(): + faces.push_back(k) + + cdef SparseMatrix distances = local_gdist_matrix_impl( + points, + faces, + max_distance, + is_one_indexed, + ) + + return scipy.sparse.csc_matrix( + (distances.data, (distances.rows, distances.columns)), shape=(N, N) + ) def distance_matrix_of_selected_points( diff --git a/geodesic_library/geodesic_utils.h b/geodesic_library/geodesic_utils.h index c2dedc6a..67c21ed9 100644 --- a/geodesic_library/geodesic_utils.h +++ b/geodesic_library/geodesic_utils.h @@ -1,7 +1,15 @@ #include +#include #include "geodesic_algorithm_exact.h" +class SparseMatrix { +public: + std::vector rows, columns; + std::vector data; +}; + + std::vector compute_gdist_impl( std::vector points, std::vector faces, @@ -35,3 +43,48 @@ std::vector compute_gdist_impl( } return distances; } + +SparseMatrix local_gdist_matrix_impl( + std::vector points, + std::vector faces, + double max_distance, + bool is_one_indexed +) { + geodesic::Mesh mesh; + mesh.initialize_mesh_data(points, faces, is_one_indexed=is_one_indexed); // create internal mesh data structure including edges + + SparseMatrix distances; + + std::vector rows, columns; + std::vector data; + #pragma omp parallel + { + geodesic::GeodesicAlgorithmExact algorithm(&mesh); + std::vector rows_private, columns_private; + std::vector data_private; + #pragma omp for nowait + for (int i = 0; i < (int)mesh.vertices().size(); ++i) { + std::vector sources {&mesh.vertices()[i]}; + algorithm.propagate(sources, geodesic::GEODESIC_INF, NULL); // cover the whole mesh + for(int j = 0; j < (int)mesh.vertices().size(); ++j) { + geodesic::SurfacePoint p(&mesh.vertices()[j]); + + double distance; + unsigned best_source = algorithm.best_source(p, distance); // for a given surface point, find closest source and distance to this source + + if (distance != 0 && distance <= geodesic::GEODESIC_INF && distance <= max_distance) { + rows_private.push_back(i); + columns_private.push_back(j); + data_private.push_back(distance); + } + } + } + rows.insert(rows.end(), rows_private.begin(), rows_private.end()); + columns.insert(columns.end(), columns_private.begin(), columns_private.end()); + data.insert(data.end(), data_private.begin(), data_private.end()); + } + distances.rows = rows; + distances.columns = columns; + distances.data = data; + return distances; +} diff --git a/setup.py b/setup.py index 98981880..bef6e252 100644 --- a/setup.py +++ b/setup.py @@ -41,6 +41,7 @@ """ import os +import sys import setuptools import numpy @@ -66,9 +67,16 @@ name=GEODESIC_NAME, # Name of extension sources=["gdist.pyx"], # Filename of Cython source language="c++", # Cython create C++ source + # Disable assertions; one is failing geodesic_mesh.h:405 define_macros=define_macros, - extra_compile_args=["--std=c++14"], - extra_link_args=["--std=c++14"], + extra_compile_args=[ + '--std=c++14', + '/openmp' if sys.platform == 'win32' else '-fopenmp', + ], + extra_link_args=[ + '--std=c++14', + '/openmp' if sys.platform == 'win32' else '-fopenmp', + ], include_dirs=[numpy.get_include(), "geodesic_library"], ) ] diff --git a/tests/test_geodesic_utils.cpp b/tests/test_geodesic_utils.cpp index 4b118db1..5aae584b 100644 --- a/tests/test_geodesic_utils.cpp +++ b/tests/test_geodesic_utils.cpp @@ -5,6 +5,13 @@ #include "../geodesic_library/geodesic_utils.h" +std::vector> sparse_to_matrix(SparseMatrix spareseMatrix, unsigned size) { + std::vector> matrix(size, std::vector(size)); + for (unsigned i = 0; i < spareseMatrix.rows.size(); ++i) { + matrix[spareseMatrix.rows[i]][spareseMatrix.columns[i]] = spareseMatrix.data[i]; + } + return matrix; +} TEST(compute_gdist_impl, flat_traingular_mesh_test) { std::vector points; @@ -27,6 +34,20 @@ TEST(compute_gdist_impl, flat_traingular_mesh_test) { } } +TEST(local_gdist_matrix_impl, flat_triangular_mesh_test) { + std::vector points; + std::vector faces; + geodesic::read_mesh_from_file("../data/flat_triangular_mesh.txt", points, faces); + SparseMatrix gdist_matrix = local_gdist_matrix_impl( + points, + faces, + geodesic::GEODESIC_INF, + false + ); + std::vector> matrix = sparse_to_matrix(gdist_matrix, points.size() / 3); + EXPECT_NEAR(matrix[1][0], 0.2, 1e-6); +} + int main(int argc, char **argv) { testing::InitGoogleTest(&argc, argv); return RUN_ALL_TESTS();