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 4fe9937..5e035c0 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -3,28 +3,36 @@ notifications:
install:
- pip3 install .
- # - python setup.py build_ext --inplace
script:
- - pip3 install pytest~=3.6.1
- - pytest
+ - pip3 install nose2 nose2-cov
+ - nose2 --with-coverage
jobs:
include:
- name: "Python 3.8 on Xenial Linux"
language: python
python: 3.8
+ before_install:
+ - pip3 install codecov
+ - sudo chmod +x linux_so.sh
+ - ./linux_so.sh
+ after_success: codecov
+ - name: "Python 3.7.4 on macOS"
+ os: osx
+ osx_image: xcode11.2
+ language: shell
+ before_install:
+ - chmod 755 macos_dylib.sh
+ - ./macos_dylib.sh
- 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
- - name: "Python 3.7.4 on macOS"
- os: osx
- osx_image: xcode11.2
- language: shell
- stage: lint
language: python
python: 3.8
diff --git a/alternative_geodesic.py b/alternative_geodesic.py
deleted file mode 100644
index 4ae05aa..0000000
--- a/alternative_geodesic.py
+++ /dev/null
@@ -1,224 +0,0 @@
-# -*- coding: utf-8 -*-
-#
-#
-# TheVirtualBrain-Framework Package. This package holds all Data Management, and
-# Web-UI helpful to run brain-simulations. To use it, you also need do download
-# TheVirtualBrain-Scientific Package (for simulators). See content of the
-# documentation-folder for more details. See also http://www.thevirtualbrain.org
-#
-# (c) 2012-2020, Baycrest Centre for Geriatric Care ("Baycrest") and others
-#
-# This program is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software Foundation,
-# either version 3 of the License, or (at your option) any later version.
-# This program is distributed in the hope that it will be useful, but WITHOUT ANY
-# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
-# PARTICULAR PURPOSE. See the GNU General Public License for more details.
-# You should have received a copy of the GNU General Public License along with this
-# program. If not, see .
-#
-#
-# CITATION:
-# When using The Virtual Brain for scientific publications, please cite it as follows:
-#
-# Paula Sanz Leon, Stuart A. Knock, M. Marmaduke Woodman, Lia Domide,
-# Jochen Mersmann, Anthony R. McIntosh, Viktor Jirsa (2013)
-# The Virtual Brain: a simulator of primate brain network dynamics.
-# Frontiers in Neuroinformatics (7:10. doi: 10.3389/fninf.2013.00010)
-#
-#
-
-"""
-Translation of the C++ geodesic library to Python
-
-mw 18/10/2013
-
-"""
-
-import time
-import numpy
-
-# geodesic_constants_and_simple_functions.h
-GEODESIC_INF = 1e100
-SMALLEST_INTERVAL_RATIO = 1e-6
-
-
-def cos_from_edges(a, b, c):
- assert all(a > 1e-50) and all(b > 1e-50) and all(c > 1e-50)
- return numpy.clip((b * b + c * c - a * a) / (2.0 * b * c), -1.0, 1.0)
-
-
-def angle_from_edges(a, b, c):
- return numpy.arccos(cos_from_edges(a, b, c))
-
-
-def read_mesh_from(filename):
- raise NotImplemented
-
-
-# geodesic_mesh_elements.h
-class MeshElement(object):
- def __init__(self):
- self.adjacent_vertices = []
- self.adjacent_faces = []
- self.adjacent_edges = []
-
-
-class Point3D(object):
- x, y, z = 0., 0., 0.
-
- def distance(self, v):
- x, y, z = v
- dx, dy, dz = self.x - x, self.y - y, self.z - z
- return numpy.sqrt(dx * dx + dy * dy + dz * dz)
-
- def set(self, x, y, z):
- self.x, self.y, self.z = 0., 0., 0.
-
- def iadd(self, v):
- self.x += v[0]
- self.y += v[1]
- self.z += v[2]
-
- def imul(self, v):
- self.x *= v
- self.y *= v
- self.z *= v
-
-
-class Vertex(MeshElement, Point3D):
- saddle_or_boundary = None
-
-
-class Face(MeshElement):
- corner_angles = [None] * 3
-
- def opposite_edge(vertex):
- raise NotImplemented
-
- def opposite_vertex(edge):
- raise NotImplemented
-
- def next_edge(edge, vertex):
- raise NotImplemented
-
- def vertex_angle(self, vertex):
- for v, a in zip(self.adjacent_vertices, self.corner_angles):
- if v == vertex:
- return a
-
-
-class Edge(MeshElement):
- length = 0.0
-
- def opposite_face(self, face):
- raise NotImplemented
-
- def opposite_vertex(self, vertex):
- raise NotImplemented
-
- def belongs(self, vertex):
- raise NotImplemented
-
- def is_boundary(self):
- return 1 == len(self.adjacent_faces)
-
- def local_coordinates(point, x, y):
- raise NotImplemented
-
-
-class SurfacePoint(Point3D):
- """
- A point lying anywhere on the surface of the mesh
-
- """
-
- def __init__(self, p3, a=0.5):
- if isinstance(p3, Vertex):
- self.p = p3
- elif isinstance(p3, Face):
- self.set(0., 0., 0.)
- [self.iadd(vi) for vi in p3.adjacent_vertices]
- self.imul(1. / 3)
- elif isinstance(p3, Edge):
- b = 1 - a
- v0 = p3.adjacent_vertices[0]
- v1 = p3.adjacent_vertices[1]
- self.x = b * v0.x + a * v1.x
- self.y = b * v0.y + a * v1.y
- self.z = b * v0.z + a * v1.z
-
-
-class HalfEdge(object):
- # ints in C++
- face_id, vertex_0, vertex_1 = None, None, None
-
- def __lt__(x, y):
- if (x.vertex_0 == y.vertex_0):
- return x.vertex_1 < y.vertex_1
- else:
- return x.vertex_0 < y.vertex_0
-
- def __ne__(x, y):
- return x.vertex_0 != y.vertex_0 or x.vertex_1 != y.vertex_1
-
- def __eq__(x, y):
- return x.vertex_0 == y.vertex_0 and x.vertex_1 == y.vertex_1
-
-
-class SurfacePath(object):
- # std::vector& path
- path = []
-
- def length(self):
- raise NotImplemented
-
- def print_info_about_path(self):
- raise NotImplemented
-
-
-# geodesic_algorithm_base.h
-
-
-class Base(object):
- """
- Base algorithm, from geodesic_algorithm_base.h
-
- """
-
- def __init__(self):
- self.max_propagation_distance = 1e100
- self.mesh = None
-
- def propagate(self, sources, max_propagation_distance, stop_points):
- raise NotImplemented
-
- def trace_back(self, destination, path):
- raise NotImplemented
-
- def geodesic(self, source, destination, path):
- raise NotImplemented
-
- def best_source(self, point, distance):
- raise NotImplemented
-
- def print_statistics(self):
- raise NotImplemented
-
- def set_stop_conditions(self, stop_points, stop_distance):
- raise NotImplemented
-
- def stop_distance(self):
- return self.max_propagation_distance
-
-
-class Dijkstra(Base):
- pass
-
-
-class Subdivision(Base):
- pass
-
-
-class Exact(Base):
- pass
diff --git a/create_dll.bat b/create_dll.bat
new file mode 100644
index 0000000..7230c8c
--- /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 /EHsc ..\..\geodesic_library\gdist_c_api.cpp
+ls
diff --git a/gdist.py b/gdist.py
new file mode 100644
index 0000000..bc11a3f
--- /dev/null
+++ b/gdist.py
@@ -0,0 +1,161 @@
+import ctypes
+import glob
+import os
+import sys
+
+import numpy as np
+import scipy.sparse
+
+
+if sys.platform == 'win32':
+ libfile = glob.glob('build/*/gdist_c_api.dll')[0]
+ libfile = os.path.abspath(libfile)
+ lib = ctypes.CDLL(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(libfile)
+else:
+ libfile = glob.glob('build/*/gdist*.so')[0]
+ lib = ctypes.CDLL(libfile)
+
+lib.compute_gdist.argtypes = [
+ ctypes.c_uint,
+ ctypes.c_uint,
+ np.ctypeslib.ndpointer(dtype=np.float64),
+ np.ctypeslib.ndpointer(dtype=np.uint32),
+ ctypes.c_uint,
+ ctypes.c_uint,
+ np.ctypeslib.ndpointer(dtype=np.uint32),
+ np.ctypeslib.ndpointer(dtype=np.uint32),
+ np.ctypeslib.ndpointer(dtype=np.float64),
+ ctypes.c_double,
+]
+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.uint32),
+ ctypes.POINTER(ctypes.c_uint),
+ ctypes.c_double,
+]
+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,
+ ):
+ 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,
+ )
+ return distance
+
+ def local_gdist_matrix(
+ self,
+ number_of_vertices,
+ number_of_triangles,
+ vertices,
+ triangles,
+ max_distance,
+ ):
+ 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,
+ )
+
+ 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,
+):
+ 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,
+ )
+ return np.fromiter(distance, dtype=np.float64, count=target_indices.size)
+
+
+def local_gdist_matrix(
+ vertices,
+ triangles,
+ max_distance=1e100,
+):
+ vertices = vertices.ravel()
+ triangles = triangles.ravel()
+
+ g = Gdist()
+ data = g.local_gdist_matrix(
+ vertices.size,
+ triangles.size,
+ vertices,
+ triangles,
+ max_distance,
+ )
+ assert data.size % 3 == 0
+ 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.pyx b/gdist.pyx
deleted file mode 100644
index a2898c0..0000000
--- a/gdist.pyx
+++ /dev/null
@@ -1,278 +0,0 @@
-# -*- coding: utf-8 -*-
-# cython: language_level=3
-#
-#
-# TheVirtualBrain-Framework Package. This package holds all Data Management, and
-# Web-UI helpful to run brain-simulations. To use it, you also need do download
-# TheVirtualBrain-Scientific Package (for simulators). See content of the
-# documentation-folder for more details. See also http://www.thevirtualbrain.org
-#
-# (c) 2012-2020, Baycrest Centre for Geriatric Care ("Baycrest") and others
-#
-# This program is free software: you can redistribute it and/or modify it under the
-# terms of the GNU General Public License as published by the Free Software Foundation,
-# either version 3 of the License, or (at your option) any later version.
-# This program is distributed in the hope that it will be useful, but WITHOUT ANY
-# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
-# PARTICULAR PURPOSE. See the GNU General Public License for more details.
-# You should have received a copy of the GNU General Public License along with this
-# program. If not, see .
-#
-#
-# CITATION:
-# When using The Virtual Brain for scientific publications, please cite it as follows:
-#
-# Paula Sanz Leon, Stuart A. Knock, M. Marmaduke Woodman, Lia Domide,
-# Jochen Mersmann, Anthony R. McIntosh, Viktor Jirsa (2013)
-# The Virtual Brain: a simulator of primate brain network dynamics.
-# Frontiers in Neuroinformatics (7:10. doi: 10.3389/fninf.2013.00010)
-#
-#
-
-"""
-This module defines a Cython wrapper for the geodesic distance C++ library.
-The algorithm (implemented in C++) extends Mitchell, Mount and Papadimitriou
-(1987) and was written by Danil Kirsanov
-(http://code.google.com/archive/p/geodesic/).
-
-The interface definitions and wrapper functions are written in Cython syntax
-(http://cython.org/) and provide an API for some of the classes, functions and
-constants useful for calculating the geodesic distance.
-
-To compile, and build gdist.so using Cython::
-
- python setup.py build_ext --inplace
-
-.. moduleauthor:: Gaurav Malhotra
-.. moduleauthor:: Stuart A. Knock
-
-"""
-
-#For compatible datatypes
-import numpy
-cimport numpy
-
-# For csc_matrix returned by local_gdist_matrix()
-import scipy.sparse
-
-# Pre-cdef'd containers from the C++ standard library
-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]&)
- vector[Vertex]& vertices()
-
-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&)
-
-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
-):
- # 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)
-
-
-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,
- numpy.ndarray[numpy.int32_t, ndim=1] target_indices = None,
- double max_distance = GEODESIC_INF):
- """
- This is the wrapper function for computing geodesic distance between a set
- of sources and targets on a mesh surface. This function accepts five
- arguments:
- ``vertices``: defines x,y,z coordinates of the mesh's vertices.
- ``triangles``: defines faces of the mesh as index triplets into vertices.
- ``source_indices``: Index of the source on the mesh.
- ``target_indices``: Index of the targets on the mesh.
- ``max_distance``:
- and returns a numpy.ndarray((len(target_indices), ), dtype=numpy.float64)
- specifying the shortest distance to the target vertices from the nearest
- source vertex on the mesh. If no target_indices are provided, all vertices
- of the mesh are considered as targets, however, in this case, specifying
- max_distance will limit the targets to those vertices within max_distance of
- a source.
-
- NOTE: This is the function to use when specifying localised stimuli and
- parameter variations. For efficiently using the whole mesh as sources, such
- as is required to represent local connectivity within the cortex, see the
- local_gdist_matrix() function.
-
- Basic usage then looks like::
- >>> import numpy
- >>> temp = numpy.loadtxt("data/flat_triangular_mesh.txt", skiprows=1)
- >>> vertices = temp[0:121].astype(numpy.float64)
- >>> triangles = temp[121:321].astype(numpy.int32)
- >>> src = numpy.array([1], dtype=numpy.int32)
- >>> trg = numpy.array([2], dtype=numpy.int32)
- >>> import gdist
- >>> gdist.compute_gdist(vertices, triangles, source_indices = src, target_indices = trg)
- array([ 0.2])
-
- """
-
- cdef Mesh amesh
- get_mesh(vertices, triangles, amesh)
-
- # Define and create object for exact algorithm on that mesh:
- cdef GeodesicAlgorithmExact *algorithm = new GeodesicAlgorithmExact(&amesh)
-
- # Define a C++ vector for the source vertices
- cdef vector[SurfacePoint] all_sources
-
- # Define a C++ vector for the target vertices
- cdef vector[SurfacePoint] stop_points
-
- # Define a NumPy array to hold the results
- cdef numpy.ndarray[numpy.float64_t, ndim=1] distances
-
- cdef Py_ssize_t k
-
- if source_indices is None: #Default to 0
- source_indices = numpy.arange(0, dtype=numpy.int32)
- # Map the NumPy sources of targets to a C++ vector
- for k in source_indices:
- all_sources.push_back(SurfacePoint(&amesh.vertices()[k]))
-
- if target_indices is None:
- # Propagate purely on max_distance
- algorithm.propagate(all_sources, max_distance, NULL)
- # Make all vertices targets
- # NOTE: this is inefficient in the best_source step below, but not sure
- # how to avoid.
- target_indices = numpy.arange(vertices.shape[0], dtype=numpy.int32)
- # Map the NumPy array of targets to a C++ vector
- for k in target_indices:
- stop_points.push_back(SurfacePoint(&amesh.vertices()[k]))
- else:
- # Map the NumPy array of targets to a C++ vector
- for k in target_indices:
- stop_points.push_back(SurfacePoint(&amesh.vertices()[k]))
- # Propagate to the specified targets
- algorithm.propagate(all_sources, max_distance, &stop_points)
-
- distances = numpy.zeros((len(target_indices), ), dtype=numpy.float64)
- for k in range(stop_points.size()):
- algorithm.best_source(stop_points[k], distances[k])
-
- distances[distances==GEODESIC_INF] = numpy.inf
-
- return distances
-
-
-def local_gdist_matrix(numpy.ndarray[numpy.float64_t, ndim=2] vertices,
- numpy.ndarray[numpy.int32_t, ndim=2] triangles,
- double max_distance = GEODESIC_INF):
- """
- This is the wrapper function for computing geodesic distance from every
- vertex on the surface to all those within a distance ``max_distance`` of
- them. The function accepts three arguments:
- ``vertices``: defines x,y,z coordinates of the mesh's vertices
- ``triangles``: defines faces of the mesh as index triplets into vertices.
- ``max_distance``:
- and returns a scipy.sparse.csc_matrix((N, N), dtype=numpy.float64), where N
- is the number of vertices, specifying the shortest distance from all
- vertices to all the vertices within max_distance.
-
- Basic usage then looks like::
- >>> import numpy
- >>> temp = numpy.loadtxt("data/flat_triangular_mesh.txt", skiprows=1)
- >>> import gdist
- >>> vertices = temp[0:121].astype(numpy.float64)
- >>> triangles = temp[121:321].astype(numpy.int32)
- >>> gdist.local_gdist_matrix(vertices, triangles, max_distance=1.0)
- <121x121 sparse matrix of type ''
- with 6232 stored elements in Compressed Sparse Column format>
-
- Runtime and guesstimated memory usage as a function of max_distance for the
- reg_13 cortical surface mesh, ie containing 2**13 vertices per hemisphere.
- ::
- [[10, 20, 30, 40, 50, 60, 70, 80, 90, 100], # mm
- [19, 28, 49, 81, 125, 181, 248, 331, 422, 522], # s
- [ 3, 13, 30, 56, 89, 129, 177, 232, 292, 358]] # MB]
-
- where memory is a min-guestimate given by: mem_req = nnz * 8 / 1024 / 1024
-
- NOTE: The best_source loop could be sped up considerably by restricting
- targets to those vertices within max_distance of the source, however,
- this will first require the efficient extraction of this information
- from the propgate step...
- """
-
- cdef Mesh amesh
- get_mesh(vertices, triangles, amesh)
-
- # Define and create object for exact algorithm on that mesh:
- cdef GeodesicAlgorithmExact *algorithm = new GeodesicAlgorithmExact(&amesh)
-
- 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))
diff --git a/geodesic.py b/geodesic.py
deleted file mode 100644
index c3f452b..0000000
--- a/geodesic.py
+++ /dev/null
@@ -1,196 +0,0 @@
-
-"""
-Translation of the C++ geodesic library to Python
-
-mw 18/10/2013
-
-"""
-
-import time
-import numpy
-
-# geodesic_constants_and_simple_functions.h
-GEODESIC_INF = 1e100
-SMALLEST_INTERVAL_RATIO = 1e-6
-
-
-def cos_from_edges(a, b, c):
- assert all(a > 1e-50) and all(b > 1e-50) and all(c > 1e-50)
- return numpy.clip((b*b + c*c - a*a) / (2.0 * b * c), -1.0, 1.0)
-
-
-def angle_from_edges(a, b, c):
- return numpy.arccos(cos_from_edges(a, b, c))
-
-
-def read_mesh_from(filename):
- raise NotImplemented
-
-# geodesic_mesh_elements.h
-
-
-class MeshElement(object):
- def __init__(self):
- self.adjacent_vertices = []
- self.adjacent_faces = []
- self.adjacent_edges = []
-
-
-class Point3D(object):
- x, y, z = 0., 0., 0.
-
- def distance(self, v):
- x, y, z = v
- dx, dy, dz = self.x - x, self.y - y, self.z - z
- return numpy.sqrt(dx*dx + dy*dy + dz*dz)
-
- def set(self, x, y, z):
- self.x, self.y, self.z = 0., 0., 0.
-
- def iadd(self, v):
- self.x += v[0]
- self.y += v[1]
- self.z += v[2]
-
- def imul(self, v):
- self.x *= v
- self.y *= v
- self.z *= v
-
-
-class Vertex(MeshElement, Point3D):
- saddle_or_boundary = None
-
-
-class Face(MeshElement):
- corner_angles = [None]*3
-
- def opposite_edge(vertex):
- raise NotImplemented
-
- def opposite_vertex(edge):
- raise NotImplemented
-
- def next_edge(edge, vertex):
- raise NotImplemented
-
- def vertex_angle(self, vertex):
- for v, a in zip(self.adjacent_vertices, self.corner_angles):
- if v == vertex:
- return a
-
-
-class Edge(MeshElement):
- length = 0.0
-
- def opposite_face(self, face):
- raise NotImplemented
-
- def opposite_vertex(self, vertex):
- raise NotImplemented
-
- def belongs(self, vertex):
- raise NotImplemented
-
- def is_boundary(self):
- return 1 == len(self.adjacent_faces)
-
- def local_coordinates(point, x, y):
- raise NotImplemented
-
-
-class SurfacePoint(Point3D):
- """
- A point lying anywhere on the surface of the mesh
-
- """
-
- def __init__(self, p3, a=0.5):
- if isinstance(p3, Vertex):
- self.p = p3
- elif isinstance(p3, Face):
- self.set(0., 0., 0.)
- [self.iadd(vi) for vi in p3.adjacent_vertices]
- self.imul(1./3)
- elif isinstance(p3, Edge):
- b = 1 - a
- v0 = p3.adjacent_vertices[0]
- v1 = p3.adjacent_vertices[1]
- self.x = b*v0.x + a*v1.x
- self.y = b*v0.y + a*v1.y
- self.z = b*v0.z + a*v1.z
-
-
-class HalfEdge(object):
- # ints in C++
- face_id, vertex_0, vertex_1 = None, None, None
-
- def __lt__(x, y):
- if (x.vertex_0 == y.vertex_0):
- return x.vertex_1 < y.vertex_1
- else:
- return x.vertex_0 < y.vertex_0
-
- def __ne__(x, y):
- return x.vertex_0 != y.vertex_0 or x.vertex_1 != y.vertex_1
-
- def __eq__(x, y):
- return x.vertex_0 == y.vertex_0 and x.vertex_1 == y.vertex_1
-
-
-class SurfacePath(object):
- # std::vector& path
- path = []
-
- def length(self):
- raise NotImplemented
-
- def print_info_about_path(self):
- raise NotImplemented
-
-
-# geodesic_algorithm_base.h
-
-
-class Base(object):
- """
- Base algorithm, from geodesic_algorithm_base.h
-
- """
-
- def __init__(self):
- self.max_propagation_distance = 1e100
- self.mesh = None
-
- def propagate(self, sources, max_propagation_distance, stop_points):
- raise NotImplemented
-
- def trace_back(self, destination, path):
- raise NotImplemented
-
- def geodesic(self, source, destination, path):
- raise NotImplemented
-
- def best_source(self, point, distance):
- raise NotImplemented
-
- def print_statistics(self):
- raise NotImplemented
-
- def set_stop_conditions(self, stop_points, stop_distance):
- raise NotImplemented
-
- def stop_distance(self):
- return self.max_propagation_distance
-
-
-class Dijkstra(Base):
- pass
-
-
-class Subdivision(Base):
- pass
-
-
-class Exact(Base):
- pass
diff --git a/geodesic_library/gdist_c_api.cpp b/geodesic_library/gdist_c_api.cpp
new file mode 100644
index 0000000..2bbb192
--- /dev/null
+++ b/geodesic_library/gdist_c_api.cpp
@@ -0,0 +1,148 @@
+#include "gdist_c_api.h"
+
+
+void compute_gdist_impl(
+ unsigned number_of_vertices,
+ unsigned number_of_triangles,
+ double *vertices,
+ unsigned *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()];
+ assert (data != NULL); // memory allocation should not fail
+ *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,
+ unsigned *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
+ );
+ }
+
+ void free_memory(double *ptr) {
+ free_memory_impl(ptr);
+ }
+};
diff --git a/geodesic_library/gdist_c_api.h b/geodesic_library/gdist_c_api.h
new file mode 100644
index 0000000..c31f106
--- /dev/null
+++ b/geodesic_library/gdist_c_api.h
@@ -0,0 +1,67 @@
+#include
+#include
+#include
+
+#include "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,
+ unsigned *triangles,
+ unsigned number_of_source_indices,
+ unsigned number_of_target_indices,
+ unsigned *source_indices_array,
+ unsigned *target_indices_array,
+ double *distance,
+ double distance_limit
+);
+
+double* local_gdist_matrix_impl(
+ unsigned number_of_vertices,
+ unsigned number_of_triangles,
+ double *vertices,
+ unsigned *triangles,
+ unsigned *sparse_matrix_size,
+ double max_distance
+);
+
+void free_memory_impl(double *ptr);
+
+extern "C" {
+ DLL_EXPORT_API void compute_gdist(
+ unsigned number_of_vertices,
+ unsigned number_of_triangles,
+ double *vertices,
+ unsigned *triangles,
+ unsigned number_of_source_indices,
+ unsigned number_of_target_indices,
+ unsigned *source_indices_array,
+ unsigned *target_indices_array,
+ double *distance,
+ double distance_limit
+ );
+
+ 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
+ );
+
+ DLL_EXPORT_API void free_memory(double *ptr);
+};
diff --git a/geodesic_library/geodesic_algorithm_exact.h b/geodesic_library/geodesic_algorithm_exact.h
index fad8c41..0428bf7 100644
--- a/geodesic_library/geodesic_algorithm_exact.h
+++ b/geodesic_library/geodesic_algorithm_exact.h
@@ -1140,7 +1140,7 @@ inline void GeodesicAlgorithmExact::construct_propagated_intervals(bool invert,
p->min() = 0.0; //it will be changed later on
- assert(p->start() < p->stop());
+ // assert(p->start() < p->stop());
}
}
else //now we have to invert the intervals
@@ -1163,7 +1163,7 @@ inline void GeodesicAlgorithmExact::construct_propagated_intervals(bool invert,
p->min() = 0;
- assert(p->start() < p->stop());
+ // assert(p->start() < p->stop());
assert(p->start() >= 0.0);
assert(p->stop() <= edge->length());
}
diff --git a/linux_so.sh b/linux_so.sh
new file mode 100644
index 0000000..be95d99
--- /dev/null
+++ b/linux_so.sh
@@ -0,0 +1,5 @@
+mkdir build
+cd build
+mkdir lib.linux
+cd lib.linux
+g++ -std=c++11 -shared -fPIC ../../geodesic_library/gdist_c_api.cpp -o gdist_c_api.so
diff --git a/macos_dylib.sh b/macos_dylib.sh
new file mode 100644
index 0000000..1f84f4f
--- /dev/null
+++ b/macos_dylib.sh
@@ -0,0 +1,5 @@
+mkdir build
+cd build
+mkdir lib.macos
+cd lib.macos
+clang++ -std=c++11 -shared -fPIC ../../geodesic_library/gdist_c_api.cpp -o gdist_c_api.dylib
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 9826b81..1a46cbf 100644
--- a/setup.py
+++ b/setup.py
@@ -41,34 +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'],
- include_dirs=[numpy.get_include(), "geodesic_library"],
- )
-]
+GEODESIC_MODULE = []
+
+if sys.platform == 'darwin' or sys.platform == 'linux':
+ GEODESIC_MODULE = [
+ setuptools.Extension(
+ name=GEODESIC_NAME, # Name of extension
+ sources=['geodesic_library/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()
@@ -76,9 +73,10 @@
setuptools.setup(
name="tvb-" + GEODESIC_NAME,
version='2.0.2',
- ext_modules=GEODESIC_MODULE,
+ 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,
@@ -88,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_equality_with_stable.py b/tests/test_equality_with_stable.py
index 594a6fc..f828563 100644
--- a/tests/test_equality_with_stable.py
+++ b/tests/test_equality_with_stable.py
@@ -14,7 +14,7 @@ def test_equality_with_stable():
)
triangles = np.loadtxt(
f'data/surface_data/{surface_data}/triangles.txt',
- dtype=np.int32,
+ dtype=np.uint32,
)
actual = gdist.local_gdist_matrix(
vertices=vertices,
diff --git a/tests/test_gdist.py b/tests/test_gdist.py
index d00a9f9..f32e733 100644
--- a/tests/test_gdist.py
+++ b/tests/test_gdist.py
@@ -1,67 +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)
- # 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.uint32)
+ source = np.array([1], dtype=np.uint32)
+ target = np.array([2], dtype=np.uint32)
+ distance = gdist.compute_gdist(
+ vertices,
+ triangles,
+ source,
+ 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.uint32)
+ source = np.array([0], dtype=np.uint32)
+ target = np.array([1], dtype=np.uint32)
+ 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.uint32)
+ 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.uint32)
+ 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