Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

TVB-2719 Improve documentation #50

Merged
merged 1 commit into from
Jul 24, 2020
Merged

TVB-2719 Improve documentation #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
127 changes: 118 additions & 9 deletions README.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,10 @@
=================
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

The **tvb-gdist** module is a Cython interface to a C++ library
(https://code.google.com/archive/p/geodesic/) for computing
geodesic distance which is the length of shortest line between two
Expand All @@ -11,25 +15,130 @@ The algorithm is due Mitchell, Mount and Papadimitriou, 1987; the implementation
is due to Danil Kirsanov and the Cython interface to Gaurav Malhotra and
Stuart Knock (TVB Team).


Original library (published under MIT license):
https://code.google.com/archive/p/geodesic/

We added a Python wrapped and made small fixes to the original library, to make it compatible with Cython.
We added a Python wrapped and made small fixes to the original library, to make
it compatible with Cython.

To install this, either run ``pip install tvb-gdist`` or download
sources from GitHub and run ``python setup.py install`` in current folder.

To install this, either run **pip install tvb-gdist** or download
sources from Github and run **python setup.py install** in current folder.
You can also use pip to directly install from GitHub:
``pip install git+https://github.com/the-virtual-brain/tvb-gdist``.

Basic test could be::

python
import gdist


Python 3, Cython, and a C++ compiler are required unless the Pypi whl files are compatible with your system.
Python 3, Cython, and a C++ compiler are required unless the Pypi whl files are
compatible with your system.

APIs
====

Current Build Status
=====================
.. image:: https://travis-ci.com/the-virtual-brain/tvb-gdist.svg?branch=trunk
:target: https://travis-ci.com/the-virtual-brain/tvb-gdist
The module exposes 2 APIs.

**gdist.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,
bool is_one_indexed = False)**

This is the wrapper function for computing geodesic distance between a
set of sources and targets on a mesh surface.

Args:
vertices (numpy.ndarray[numpy.float64_t, ndim=2]): Defines x,y,z
coordinates of the mesh's vertices.
triangles (numpy.ndarray[numpy.float64_t, ndim=2]): Defines faces of
the mesh as index triplets into vertices.
source_indices (numpy.ndarray[numpy.int32_t, ndim=1]): Index of the
source on the mesh.
target_indices (numpy.ndarray[numpy.int32_t, ndim=1]): Index of the
targets on the mesh.
max_distance (double): Propagation algorithm stops after reaching the
certain distance from the source.
is_one_indexed (bool): defines if the index of the triangles data is
one-indexed.

Returns:
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])


**gdist.local_gdist_matrix(numpy.ndarray[numpy.float64_t, ndim=2] vertices,
numpy.ndarray[numpy.int32_t, ndim=2] triangles,
double max_distance = GEODESIC_INF,
bool is_one_indexed = False)**

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.

Args:
vertices (numpy.ndarray[numpy.float64_t, ndim=2]): Defines x,y,z
coordinates of the mesh's vertices.
triangles (numpy.ndarray[numpy.float64_t, ndim=2]): Defines faces of
the mesh as index triplets into vertices.
max_distance (double): Propagation algorithm stops after reaching the
certain distance from the source.
is_one_indexed (bool): defines if the index of the triangles data is
one-indexed.

Returns:
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 '<type 'numpy.float64'>'
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.


Notes
=====

* The obtained matrix will be almost symmetrical due to floating point
imprecision.

* In order for the algorithm to work the mesh must not be numbered incorrectly
or disconnected or of somehow degenerate.
84 changes: 49 additions & 35 deletions gdist.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,12 +100,12 @@ cdef get_mesh(
# 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():
Expand All @@ -120,23 +120,30 @@ def compute_gdist(numpy.ndarray[numpy.float64_t, ndim=2] vertices,
numpy.ndarray[numpy.int32_t, ndim=1] target_indices = None,
double max_distance = GEODESIC_INF,
bool is_one_indexed = False):
"""
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``:
``is_one_indexed``: defines if the index of the triangles data is one-
indexed
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.
"""This is the wrapper function for computing geodesic distance between a
set of sources and targets on a mesh surface.

Args:
vertices (numpy.ndarray[numpy.float64_t, ndim=2]): Defines x,y,z
coordinates of the mesh's vertices.
triangles (numpy.ndarray[numpy.float64_t, ndim=2]): Defines faces of
the mesh as index triplets into vertices.
source_indices (numpy.ndarray[numpy.int32_t, ndim=1]): Index of the
source on the mesh.
target_indices (numpy.ndarray[numpy.int32_t, ndim=1]): Index of the
targets on the mesh.
max_distance (double): Propagation algorithm stops after reaching the
certain distance from the source.
is_one_indexed (bool): defines if the index of the triangles data is
one-indexed.

Returns:
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
Expand All @@ -151,9 +158,8 @@ def compute_gdist(numpy.ndarray[numpy.float64_t, ndim=2] vertices,
>>> 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])

>>> gdist.compute_gdist(vertices, triangles, source_indices=src, target_indices=trg)
array([0.2])
"""

cdef Mesh amesh
Expand Down Expand Up @@ -209,18 +215,24 @@ 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,
bool is_one_indexed = False):
"""
This is the wrapper function for computing geodesic distance from every
"""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``:
``is_one_indexed``: defines if the index of the triangles data is one-
indexed
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.
them.

Args:
vertices (numpy.ndarray[numpy.float64_t, ndim=2]): Defines x,y,z
coordinates of the mesh's vertices.
triangles (numpy.ndarray[numpy.float64_t, ndim=2]): Defines faces of
the mesh as index triplets into vertices.
max_distance (double): Propagation algorithm stops after reaching the
certain distance from the source.
is_one_indexed (bool): defines if the index of the triangles data is
one-indexed.

Returns:
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
Expand All @@ -239,8 +251,10 @@ def local_gdist_matrix(numpy.ndarray[numpy.float64_t, ndim=2] vertices,
[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

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
Expand Down