diff --git a/.travis.yml b/.travis.yml index ec4aed0..7c9e8dd 100644 --- a/.travis.yml +++ b/.travis.yml @@ -46,6 +46,9 @@ 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 - stage: lint language: python install: pip install flake8 diff --git a/gdist.pyx b/gdist.pyx index 704b758..e9e0d9c 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, @@ -244,44 +202,27 @@ def local_gdist_matrix(numpy.ndarray[numpy.float64_t, ndim=2] vertices, this will first require the efficient extraction of this information from the propgate step... """ + + cdef Py_ssize_t N = vertices.shape[0] - cdef Mesh amesh - get_mesh(vertices, triangles, amesh, is_one_indexed) + cdef vector[double] points + cdef vector[unsigned] faces - # Define and create object for exact algorithm on that mesh: - cdef GeodesicAlgorithmExact *algorithm = new GeodesicAlgorithmExact(&amesh) + for k in vertices.flatten(): + points.push_back(k) + for k in triangles.flatten(): + faces.push_back(k) - 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)) + 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 c2dedc6..67c21ed 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 9898188..a9bea6e 100644 --- a/setup.py +++ b/setup.py @@ -41,6 +41,7 @@ """ import os +import sys import setuptools import numpy @@ -59,6 +60,11 @@ compiler_directives["linetrace"] = True define_macros.append(("CYTHON_TRACE_NOGIL", "1")) + +if sys.platform == "darwin": + os.environ["CC"] = "/usr/local/opt/llvm/bin/clang -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" + os.environ["CXX"] = "/usr/local/opt/llvm/bin/clang++ -L/usr/local/opt/llvm/lib -Wl,-rpath,/usr/local/opt/llvm/lib" + GEODESIC_NAME = "gdist" GEODESIC_MODULE = [ @@ -66,9 +72,16 @@ name=GEODESIC_NAME, # Name of extension sources=["gdist.pyx"], # Filename of Cython source language="c++", # Cython create C++ source - define_macros=define_macros, - extra_compile_args=["--std=c++14"], - extra_link_args=["--std=c++14"], + # Disable assertions; one is failing geodesic_mesh.h:405 + define_macros=[('NDEBUG', 1)], + 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"], ) ]