diff --git a/.gitignore b/.gitignore index 7e4bae9..c7c4f05 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,8 @@ __pycache__/ # Distribution / packaging dist/ +*.egg-info +build/ # coverage file .coverage @@ -26,7 +28,10 @@ env # Executables *.out *.exe +*.o # Shared object files *.so *.pyd +*.dylib +*.dll diff --git a/.travis.yml b/.travis.yml index e32e1b8..49cd6fe 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,24 +2,31 @@ notifications: email: false install: - - pip install . - python setup.py build_ext --inplace + - pip3 install . script: - - pip install pytest~=3.6.1 - - pytest + - pip3 install pytest pytest-cov + - pytest --cov=gdist jobs: include: - name: "Python 3.8 on Xenial Linux" language: python python: 3.8 + before_install: pip3 install codecov + after_success: codecov + - name: "Python 3.7.4 on macOS" + os: osx + osx_image: xcode11.2 + language: shell - name: "Python 3.8.0 on Windows" os: windows language: shell before_install: - choco install python --version 3.8.0 - - python -m pip install --upgrade pip + - python -m pip install -U pip - pip install --user wheel + - ./create_dll.bat env: PATH=/c/Python38:/c/Python38/Scripts:$PATH - stage: lint language: python diff --git a/create_dll.bat b/create_dll.bat new file mode 100644 index 0000000..75a8fd3 --- /dev/null +++ b/create_dll.bat @@ -0,0 +1,47 @@ +echo on + +if NOT DEFINED VCINSTALLDIR ( + if exist "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" ( + call "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" x86_amd64 + echo "USING VISUAL STUDIO 17" + ) +) + +if NOT DEFINED VCINSTALLDIR ( + if exist "C:\Program Files (x86)\Microsoft Visual Studio 15.0\VC\vcvarsall.bat" ( + call "C:\Program Files (x86)\Microsoft Visual Studio 15.0\VC\vcvarsall.bat" amd64 + echo "USING VISUAL STUDIO 15" + ) +) + +if NOT DEFINED VCINSTALLDIR ( + if exist "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" ( + call "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" amd64 + echo "USING VISUAL STUDIO 14" + ) +) + +if NOT DEFINED VCINSTALLDIR ( + if exist "C:\Program Files (x86)\Microsoft Visual Studio 13.0\VC\vcvarsall.bat" ( + call "C:\Program Files (x86)\Microsoft Visual Studio 13.0\VC\vcvarsall.bat" amd64 + echo "USING VISUAL STUDIO 13" + ) +) + +if NOT DEFINED VCINSTALLDIR ( + if exist "C:\Program Files (x86)\Microsoft Visual Studio 12.0\VC\vcvarsall.bat" ( + call "C:\Program Files (x86)\Microsoft Visual Studio 12.0\VC\vcvarsall.bat" amd64 + echo "USING VISUAL STUDIO 12" + ) +) + +if NOT DEFINED VCINSTALLDIR ( + echo "No compatible visual studio found! run vcvarsall.bat first!" +) + +mkdir build\lib.win32 + +cd build\lib.win32 + +cl.exe /LD /DDLL_EXPORTS /DNDEBUG ..\..\gdist_c_api.cpp +ls diff --git a/data/flat_triangular_mesh_1_indexed.txt b/data/flat_triangular_mesh_1_indexed.txt new file mode 100644 index 0000000..ffc98f9 --- /dev/null +++ b/data/flat_triangular_mesh_1_indexed.txt @@ -0,0 +1,322 @@ +121 200 +1.500000 -1.000000 0.000000 +1.500000 -0.800000 0.000000 +1.500000 -0.600000 0.000000 +1.500000 -0.400000 0.000000 +1.500000 -0.200000 0.000000 +1.500000 0.000000 0.000000 +1.500000 0.200000 0.000000 +1.500000 0.400000 0.000000 +1.500000 0.600000 0.000000 +1.500000 0.800000 0.000000 +1.500000 1.000000 0.000000 +1.700000 -1.000000 0.000000 +1.700000 -0.800000 0.000000 +1.700000 -0.600000 0.000000 +1.700000 -0.400000 0.000000 +1.700000 -0.200000 0.000000 +1.700000 0.000000 0.000000 +1.700000 0.200000 0.000000 +1.700000 0.400000 0.000000 +1.700000 0.600000 0.000000 +1.700000 0.800000 0.000000 +1.700000 1.000000 0.000000 +1.900000 -1.000000 0.000000 +1.900000 -0.800000 0.000000 +1.900000 -0.600000 0.000000 +1.900000 -0.400000 0.000000 +1.900000 -0.200000 0.000000 +1.900000 0.000000 0.000000 +1.900000 0.200000 0.000000 +1.900000 0.400000 0.000000 +1.900000 0.600000 0.000000 +1.900000 0.800000 0.000000 +1.900000 1.000000 0.000000 +2.100000 -1.000000 0.000000 +2.100000 -0.800000 0.000000 +2.100000 -0.600000 0.000000 +2.100000 -0.400000 0.000000 +2.100000 -0.200000 0.000000 +2.100000 0.000000 0.000000 +2.100000 0.200000 0.000000 +2.100000 0.400000 0.000000 +2.100000 0.600000 0.000000 +2.100000 0.800000 0.000000 +2.100000 1.000000 0.000000 +2.300000 -1.000000 0.000000 +2.300000 -0.800000 0.000000 +2.300000 -0.600000 0.000000 +2.300000 -0.400000 0.000000 +2.300000 -0.200000 0.000000 +2.300000 0.000000 0.000000 +2.300000 0.200000 0.000000 +2.300000 0.400000 0.000000 +2.300000 0.600000 0.000000 +2.300000 0.800000 0.000000 +2.300000 1.000000 0.000000 +2.500000 -1.000000 0.000000 +2.500000 -0.800000 0.000000 +2.500000 -0.600000 0.000000 +2.500000 -0.400000 0.000000 +2.500000 -0.200000 0.000000 +2.500000 0.000000 0.000000 +2.500000 0.200000 0.000000 +2.500000 0.400000 0.000000 +2.500000 0.600000 0.000000 +2.500000 0.800000 0.000000 +2.500000 1.000000 0.000000 +2.700000 -1.000000 0.000000 +2.700000 -0.800000 0.000000 +2.700000 -0.600000 0.000000 +2.700000 -0.400000 0.000000 +2.700000 -0.200000 0.000000 +2.700000 0.000000 0.000000 +2.700000 0.200000 0.000000 +2.700000 0.400000 0.000000 +2.700000 0.600000 0.000000 +2.700000 0.800000 0.000000 +2.700000 1.000000 0.000000 +2.900000 -1.000000 0.000000 +2.900000 -0.800000 0.000000 +2.900000 -0.600000 0.000000 +2.900000 -0.400000 0.000000 +2.900000 -0.200000 0.000000 +2.900000 0.000000 0.000000 +2.900000 0.200000 0.000000 +2.900000 0.400000 0.000000 +2.900000 0.600000 0.000000 +2.900000 0.800000 0.000000 +2.900000 1.000000 0.000000 +3.100000 -1.000000 0.000000 +3.100000 -0.800000 0.000000 +3.100000 -0.600000 0.000000 +3.100000 -0.400000 0.000000 +3.100000 -0.200000 0.000000 +3.100000 0.000000 0.000000 +3.100000 0.200000 0.000000 +3.100000 0.400000 0.000000 +3.100000 0.600000 0.000000 +3.100000 0.800000 0.000000 +3.100000 1.000000 0.000000 +3.300000 -1.000000 0.000000 +3.300000 -0.800000 0.000000 +3.300000 -0.600000 0.000000 +3.300000 -0.400000 0.000000 +3.300000 -0.200000 0.000000 +3.300000 0.000000 0.000000 +3.300000 0.200000 0.000000 +3.300000 0.400000 0.000000 +3.300000 0.600000 0.000000 +3.300000 0.800000 0.000000 +3.300000 1.000000 0.000000 +3.500000 -1.000000 0.000000 +3.500000 -0.800000 0.000000 +3.500000 -0.600000 0.000000 +3.500000 -0.400000 0.000000 +3.500000 -0.200000 0.000000 +3.500000 0.000000 0.000000 +3.500000 0.200000 0.000000 +3.500000 0.400000 0.000000 +3.500000 0.600000 0.000000 +3.500000 0.800000 0.000000 +3.500000 1.000000 0.000000 +1 2 12 +2 12 13 +2 3 13 +3 13 14 +3 4 14 +4 14 15 +4 5 15 +5 15 16 +5 6 16 +6 16 17 +6 7 17 +7 17 18 +7 8 18 +8 18 19 +8 9 19 +9 19 20 +9 10 20 +10 20 21 +10 11 21 +11 21 22 +12 13 23 +13 23 24 +13 14 24 +14 24 25 +14 15 25 +15 25 26 +15 16 26 +16 26 27 +16 17 27 +17 27 28 +17 18 28 +18 28 29 +18 19 29 +19 29 30 +19 20 30 +20 30 31 +20 21 31 +21 31 32 +21 22 32 +22 32 33 +23 24 34 +24 34 35 +24 25 35 +25 35 36 +25 26 36 +26 36 37 +26 27 37 +27 37 38 +27 28 38 +28 38 39 +28 29 39 +29 39 40 +29 30 40 +30 40 41 +30 31 41 +31 41 42 +31 32 42 +32 42 43 +32 33 43 +33 43 44 +34 35 45 +35 45 46 +35 36 46 +36 46 47 +36 37 47 +37 47 48 +37 38 48 +38 48 49 +38 39 49 +39 49 50 +39 40 50 +40 50 51 +40 41 51 +41 51 52 +41 42 52 +42 52 53 +42 43 53 +43 53 54 +43 44 54 +44 54 55 +45 46 56 +46 56 57 +46 47 57 +47 57 58 +47 48 58 +48 58 59 +48 49 59 +49 59 60 +49 50 60 +50 60 61 +50 51 61 +51 61 62 +51 52 62 +52 62 63 +52 53 63 +53 63 64 +53 54 64 +54 64 65 +54 55 65 +55 65 66 +56 57 67 +57 67 68 +57 58 68 +58 68 69 +58 59 69 +59 69 70 +59 60 70 +60 70 71 +60 61 71 +61 71 72 +61 62 72 +62 72 73 +62 63 73 +63 73 74 +63 64 74 +64 74 75 +64 65 75 +65 75 76 +65 66 76 +66 76 77 +67 68 78 +68 78 79 +68 69 79 +69 79 80 +69 70 80 +70 80 81 +70 71 81 +71 81 82 +71 72 82 +72 82 83 +72 73 83 +73 83 84 +73 74 84 +74 84 85 +74 75 85 +75 85 86 +75 76 86 +76 86 87 +76 77 87 +77 87 88 +78 79 89 +79 89 90 +79 80 90 +80 90 91 +80 81 91 +81 91 92 +81 82 92 +82 92 93 +82 83 93 +83 93 94 +83 84 94 +84 94 95 +84 85 95 +85 95 96 +85 86 96 +86 96 97 +86 87 97 +87 97 98 +87 88 98 +88 98 99 +89 90 100 +90 100 101 +90 91 101 +91 101 102 +91 92 102 +92 102 103 +92 93 103 +93 103 104 +93 94 104 +94 104 105 +94 95 105 +95 105 106 +95 96 106 +96 106 107 +96 97 107 +97 107 108 +97 98 108 +98 108 109 +98 99 109 +99 109 110 +100 101 111 +101 111 112 +101 102 112 +102 112 113 +102 103 113 +103 113 114 +103 104 114 +104 114 115 +104 105 115 +105 115 116 +105 106 116 +106 116 117 +106 107 117 +107 117 118 +107 108 118 +108 118 119 +108 109 119 +109 119 120 +109 110 120 +110 120 121 diff --git a/gdist.py b/gdist.py new file mode 100644 index 0000000..b20f3fa --- /dev/null +++ b/gdist.py @@ -0,0 +1,169 @@ +import ctypes +import glob +import os +import sys + +import numpy as np +import scipy.sparse + + +if sys.platform == 'win32': + libfile = os.path.abspath('build/lib.win32/gdist_c_api.dll') + lib = ctypes.windll.LoadLibrary(libfile) +elif sys.platform == 'darwin': + try: + libfile = glob.glob('build/*/gdist*.so')[0] + except IndexError: + libfile = glob.glob('build/*/gdist*.dylib')[0] + lib = ctypes.cdll.LoadLibrary(libfile) +else: + libfile = glob.glob('build/*/gdist*.so')[0] + lib = ctypes.cdll.LoadLibrary(libfile) + +lib.compute_gdist.argtypes = [ + ctypes.c_uint, + ctypes.c_uint, + np.ctypeslib.ndpointer(dtype=np.float64), + np.ctypeslib.ndpointer(dtype=np.int32), + ctypes.c_uint, + ctypes.c_uint, + np.ctypeslib.ndpointer(dtype=np.int32), + np.ctypeslib.ndpointer(dtype=np.int32), + np.ctypeslib.ndpointer(dtype=np.float64), + ctypes.c_double, + ctypes.c_uint, +] +lib.compute_gdist.restype = None + +lib.local_gdist_matrix.argtypes = [ + ctypes.c_uint, + ctypes.c_uint, + np.ctypeslib.ndpointer(dtype=np.float64), + np.ctypeslib.ndpointer(dtype=np.int32), + ctypes.POINTER(ctypes.c_uint), + ctypes.c_double, + ctypes.c_uint, +] +lib.local_gdist_matrix.restype = ctypes.POINTER(ctypes.c_double) + +lib.free_memory.argtypes = [ + ctypes.POINTER(ctypes.c_double), +] +lib.free_memory.restype = None + + +class Gdist(object): + def compute_gdist( + self, + number_of_vertices, + number_of_triangles, + vertices, + triangles, + number_of_source_indices, + number_of_target_indices, + source_indices_array, + target_indices_array, + distance_limit, + is_one_indexed, + ): + target_indices_size = target_indices_array.size + distance = np.empty(target_indices_size, dtype=np.float64) + lib.compute_gdist( + 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, + is_one_indexed, + ) + return distance + + def local_gdist_matrix( + self, + number_of_vertices, + number_of_triangles, + vertices, + triangles, + max_distance, + is_one_indexed, + ): + sparse_matrix_size = ctypes.c_uint(0) + data = lib.local_gdist_matrix( + number_of_vertices, + number_of_triangles, + vertices, + triangles, + ctypes.byref(sparse_matrix_size), + max_distance, + is_one_indexed, + ) + + np_data = np.fromiter( + data, + dtype=np.float64, + count=3 * sparse_matrix_size.value, + ) + lib.free_memory(data) + return np_data + + +def compute_gdist( + vertices, + triangles, + source_indices=None, + target_indices=None, + max_distance=1e100, + is_one_indexed=False, +): + vertices = vertices.ravel() + triangles = triangles.ravel() + source_indices = source_indices.ravel() + target_indices = target_indices.ravel() + + g = Gdist() + distance = g.compute_gdist( + number_of_vertices=vertices.size, + number_of_triangles=triangles.size, + vertices=vertices, + triangles=triangles, + number_of_source_indices=source_indices.size, + number_of_target_indices=target_indices.size, + source_indices_array=source_indices, + target_indices_array=target_indices, + distance_limit=max_distance, + is_one_indexed=int(is_one_indexed), + ) + return np.fromiter(distance, dtype=np.float64, count=target_indices.size) + + +def local_gdist_matrix( + vertices, + triangles, + max_distance=1e100, + is_one_indexed=False, +): + vertices = vertices.ravel() + triangles = triangles.ravel() + + g = Gdist() + data = g.local_gdist_matrix( + number_of_vertices=vertices.size, + number_of_triangles=triangles.size, + vertices=vertices, + triangles=triangles, + max_distance=max_distance, + is_one_indexed=is_one_indexed, + ) + sizes = data.size // 3 + rows = data[:sizes] + columns = data[sizes: 2 * sizes] + data = data[2 * sizes:] + + return scipy.sparse.csc_matrix( + (data, (rows, columns)), shape=(vertices.size // 3, vertices.size // 3) + ) diff --git a/gdist_c_api.cpp b/gdist_c_api.cpp new file mode 100644 index 0000000..a0f0939 --- /dev/null +++ b/gdist_c_api.cpp @@ -0,0 +1,154 @@ +#include "gdist_c_api.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, + unsigned is_one_indexed +) { + + 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, is_one_indexed); // 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, + unsigned is_one_indexed +) { + 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, is_one_indexed); // 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; +} + + +void free_memory_impl(double *ptr) { + delete[] ptr; +} + + +extern "C" { + void 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, + unsigned is_one_indexed + ) { + 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, + is_one_indexed + ); + } + + double* local_gdist_matrix( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + unsigned *triangles, + unsigned *sparse_matrix_size, + double max_distance, + unsigned is_one_indexed + ) { + return local_gdist_matrix_impl( + number_of_vertices, + number_of_triangles, + vertices, + triangles, + sparse_matrix_size, + max_distance, + is_one_indexed + ); + } + + void free_memory(double *ptr) { + free_memory_impl(ptr); + } +}; diff --git a/gdist_c_api.h b/gdist_c_api.h new file mode 100644 index 0000000..717d24e --- /dev/null +++ b/gdist_c_api.h @@ -0,0 +1,71 @@ +#include +#include +#include + +#include "geodesic_library/geodesic_algorithm_exact.h" + + +#if defined(_WIN32) +# if defined(DLL_EXPORTS) +# define DLL_EXPORT_API __declspec(dllexport) +# else +# define DLL_EXPORT_API __declspec(dllimport) +# endif +#else +# define DLL_EXPORT_API +#endif + + +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, + unsigned is_one_indexed +); + +double* local_gdist_matrix_impl( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + unsigned *triangles, + unsigned *sparse_matrix_size, + double max_distance, + unsigned is_one_indexed +); + +void free_memory_impl(double *ptr); + +extern "C" { + DLL_EXPORT_API void 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, + unsigned is_one_indexed + ); + + DLL_EXPORT_API double* local_gdist_matrix( + unsigned number_of_vertices, + unsigned number_of_triangles, + double *vertices, + unsigned *triangles, + unsigned *sparse_matrix_size, + double max_distance, + unsigned is_one_indexed + ); + + DLL_EXPORT_API void free_memory(double *ptr); +}; diff --git a/geodesic_library/geodesic_mesh.h b/geodesic_library/geodesic_mesh.h index a1a7a60..bf895b2 100644 --- a/geodesic_library/geodesic_mesh.h +++ b/geodesic_library/geodesic_mesh.h @@ -34,10 +34,11 @@ class Mesh void initialize_mesh_data(unsigned num_vertices, Points& p, unsigned num_faces, - Faces& tri); //build mesh from regular point-triangle representation + Faces& tri, + bool is_one_indexed=false); //build mesh from regular point-triangle representation template - void initialize_mesh_data(Points& p, Faces& tri); //build mesh from regular point-triangle representation + void initialize_mesh_data(Points& p, Faces& tri, bool is_one_indexed=false); //build mesh from regular point-triangle representation std::vector& vertices(){return m_vertices;}; std::vector& edges(){return m_edges;}; @@ -111,21 +112,22 @@ inline unsigned Mesh::closest_vertices(SurfacePoint* p, } template -void Mesh::initialize_mesh_data(Points& p, Faces& tri) //build mesh from regular point-triangle representation +void Mesh::initialize_mesh_data(Points& p, Faces& tri, bool is_one_indexed) //build mesh from regular point-triangle representation { assert(p.size() % 3 == 0); unsigned const num_vertices = p.size() / 3; assert(tri.size() % 3 == 0); unsigned const num_faces = tri.size() / 3; - initialize_mesh_data(num_vertices, p, num_faces, tri); + initialize_mesh_data(num_vertices, p, num_faces, tri, is_one_indexed); } template void Mesh::initialize_mesh_data(unsigned num_vertices, Points& p, unsigned num_faces, - Faces& tri) + Faces& tri, + bool is_one_indexed) { unsigned const approximate_number_of_internal_pointers = (num_vertices + num_faces)*4; unsigned const max_number_of_pointer_blocks = 100; @@ -154,7 +156,7 @@ void Mesh::initialize_mesh_data(unsigned num_vertices, unsigned shift = 3*i; for(unsigned j=0; j<3; ++j) { - unsigned vertex_index = tri[shift + j]; + unsigned vertex_index = tri[shift + j] - (unsigned)is_one_indexed; assert(vertex_index < num_vertices); f.adjacent_vertices()[j] = &m_vertices[vertex_index]; } diff --git a/pyproject.toml b/pyproject.toml index 83de03a..8b88f61 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,2 +1,2 @@ [build-system] -requires = ["setuptools", "wheel", "Cython", "numpy"] +requires = ["setuptools", "wheel", "numpy"] diff --git a/requirements.txt b/requirements.txt index 701f255..6bad103 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,2 @@ -cython numpy scipy diff --git a/setup.py b/setup.py index a7f5f3b..93547a2 100644 --- a/setup.py +++ b/setup.py @@ -41,33 +41,31 @@ """ import os -import numpy -import shutil +import sys import setuptools -from Cython.Distutils import build_ext GEODESIC_NAME = "gdist" -GEODESIC_MODULE = [ - setuptools.Extension( - 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=[('NDEBUG', 1)], - extra_compile_args=['--std=c++14'], - extra_link_args=['--std=c++14'], - ) -] +GEODESIC_MODULE = [] + +if sys.platform == 'darwin' or sys.platform == 'linux': + GEODESIC_MODULE = [ + setuptools.Extension( + name=GEODESIC_NAME, # Name of extension + sources=["gdist_c_api.cpp"], + language="c++", + extra_compile_args=['--std=c++11'], + extra_link_args=['--std=c++11'], + ) + ] INCLUDE_DIRS = [ - numpy.get_include(), # NumPy dtypes "geodesic_library", # geodesic distance, C++ library. ] TEAM = "Danil Kirsanov, Gaurav Malhotra and Stuart Knock" -INSTALL_REQUIREMENTS = ['numpy', 'scipy', 'cython'] +INSTALL_REQUIREMENTS = ['numpy', 'scipy'] with open(os.path.join(os.path.dirname(__file__), 'README.rst')) as fd: DESCRIPTION = fd.read() @@ -75,9 +73,10 @@ setuptools.setup( name="tvb-" + GEODESIC_NAME, version='2.0.2', + scripts=['gdist.py'], + py_modules=['gdist'], ext_modules=GEODESIC_MODULE, include_dirs=INCLUDE_DIRS, - cmdclass={'build_ext': build_ext}, install_requires=INSTALL_REQUIREMENTS, description="Compute geodesic distances", long_description=DESCRIPTION, @@ -87,8 +86,3 @@ url='https://github.com/the-virtual-brain/tvb-gdist', keywords="gdist geodesic distance geo tvb" ) - -shutil.rmtree('tvb_gdist.egg-info', True) -if os.path.exists(GEODESIC_NAME + '.cpp'): - os.remove(GEODESIC_NAME + '.cpp') -shutil.rmtree('build', True) diff --git a/tests/test_gdist.py b/tests/test_gdist.py index d00a9f9..bc26d94 100644 --- a/tests/test_gdist.py +++ b/tests/test_gdist.py @@ -1,67 +1,107 @@ -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 +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, + target + ) + np.testing.assert_array_almost_equal(distance, [0.2]) + + def test_flat_triangular_mesh_1_indexed(self): + data = np.loadtxt( + "data/flat_triangular_mesh_1_indexed.txt", + skiprows=1, + ) + vertices = data[0:121].astype(np.float64) + triangles = data[121:].astype(np.int32) + source = np.array([2], dtype=np.int32) + target = np.array([3], dtype=np.int32) + distance = gdist.compute_gdist( + vertices, + triangles, + source, + target, + is_one_indexed=True, + ) + 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, + 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_flat_triangular_mesh_1_indexed(self): + data = np.loadtxt( + "data/flat_triangular_mesh_1_indexed.txt", + skiprows=1, + ) + vertices = data[0:121].astype(np.float64) + triangles = data[121:].astype(np.int32) + distances = gdist.local_gdist_matrix( + vertices, + triangles, + is_one_indexed=True + ) + 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, True) + # 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