First draft of the brainprintpython package

.gitignore

@@ -0,0 +1,5 @@
+__pycache__
+*pyc
+deprecated
+tests
+package

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2019 Image Analysis, DZNE e.V.
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,39 @@
+# BrainPrint-python
+
+This is the `brainprintpython` package, an experimental derivative of the 
+original BrainPrint scripts, with the following goals and changes:
+
+## Goals
+
+- provide a Python-only version of the BrainPrint scripts (except some 
+  Freesurfer dependencies)
+- remove dependencies on third-party software (shapeDNA binaries, gmsh, meshfix)
+- provide a light-weight version of the original scripts that contains only the
+  most frequently used submodules
+- integrate the post-processing module (for computing distances etc.)
+- create a fully modularized package whose functions can be called from other
+  python scripts without the need of spawning subprocesses, while still
+  maintaining the command-line interface of the scripts
+- provide additional files (setup.py, LICENSE, README) so that it can be 
+  packaged and distributed as a stand-alone Python package
+- revision, and, where possible, simplification and reduction of the original 
+  code base for easier maintainability
+- allow for future integration of code from the `lapy` package
+
+## Changes
+
+- no more support for analyses of cortical parcellation or label files
+- no more Python 2.x compatibility
+- currently no more support for tetrahedral meshes
+
+## Installation
+
+Use the following code to download, build and install a package from this 
+repository into your local Python package directory:
+
+`pip3 install --user git+https://github.com/reuter-lab/BrainPrint-python.git@master#egg=brainprintpython`
+
+Use the following code to install the package in editable mode to a location of
+your choice:
+
+`pip3 install --user --src /my/preferred/location --editable git+https://github.com/reuter-lab/BrainPrint-python.git@master#egg=brainprintpython`

brainprintpython/__init__.py

@@ -0,0 +1,2 @@
+name = "brainprintpython"
+

brainprintpython/computeABtria.py

@@ -0,0 +1,82 @@
+def computeABtria(v, t, lump = False):
+    """
+    computeABtria(v,t) computes the two sparse symmetric matrices representing
+           the Laplace Beltrami Operator for a given triangle mesh using
+           the linear finite element method (assuming a closed mesh or 
+           the Neumann boundary condition).
+
+    Inputs:   v - vertices : list of lists of 3 floats
+              t - triangles: list of lists of 3 int of indices (>=0) into v array
+
+    Outputs:  A - sparse sym. (n x n) positive semi definite numpy matrix 
+              B - sparse sym. (n x n) positive definite numpy matrix (inner product)
+
+    Can be used to solve sparse generalized Eigenvalue problem: A x = lambda B x
+    or to solve Poisson equation: A x = B f (where f is function on mesh vertices)
+    or to solve Laplace equation: A x = 0
+    or to model the operator's action on a vector x:   y = B\(Ax) 
+    """
+
+    import sys
+    import numpy as np
+    from scipy import sparse
+
+    v = np.array(v)
+    t = np.array(t)
+
+    # Compute vertex coordinates and a difference vector for each triangle:
+    t1 = t[:, 0]
+    t2 = t[:, 1]
+    t3 = t[:, 2]
+    v1 = v[t1, :]
+    v2 = v[t2, :]
+    v3 = v[t3, :]
+    v2mv1 = v2 - v1
+    v3mv2 = v3 - v2
+    v1mv3 = v1 - v3
+
+    # Compute cross product and 4*vol for each triangle:
+    cr  = np.cross(v3mv2,v1mv3)
+    vol = 2 * np.sqrt(np.sum(cr*cr, axis=1))
+    # zero vol will cause division by zero below, so set to small value:
+    vol_mean = 0.0001*np.mean(vol)
+    vol[vol<sys.float_info.epsilon] = vol_mean 
+
+    # compute cotangents for A
+    # using that v2mv1 = - (v3mv2 + v1mv3) this can also be seen by 
+    # summing the local matrix entries in the old algorithm
+    A12 = np.sum(v3mv2*v1mv3,axis=1)/vol
+    A23 = np.sum(v1mv3*v2mv1,axis=1)/vol
+    A31 = np.sum(v2mv1*v3mv2,axis=1)/vol
+    # compute diagonals (from row sum = 0)
+    A11 = -A12-A31
+    A22 = -A12-A23
+    A33 = -A31-A23
+    # stack columns to assemble data
+    localA = np.column_stack((A12, A12, A23, A23, A31, A31, A11, A22, A33))
+    I = np.column_stack((t1, t2, t2, t3, t3, t1, t1, t2, t3))
+    J = np.column_stack((t2, t1, t3, t2, t1, t3, t1, t2, t3))
+    # Flatten arrays:
+    I = I.flatten()
+    J = J.flatten()
+    localA = localA.flatten()
+    # Construct sparse matrix:
+    A = sparse.csc_matrix((localA, (I, J)))
+
+    if not lump:
+        # create b matrix data (account for that vol is 4 times area)
+        Bii = vol / 24
+        Bij = vol / 48
+        localB = np.column_stack((Bij, Bij, Bij, Bij, Bij, Bij, Bii, Bii, Bii))
+        localB = localB.flatten()
+        B = sparse.csc_matrix((localB, (I, J)))
+    else:
+        # when lumping put all onto diagonal  (area/3 for each vertex)
+        Bii = vol / 12
+        localB = np.column_stack((Bii, Bii, Bii))
+        I = np.column_stack((t1, t2, t3))
+        I = I.flatten()
+        localB = localB.flatten()
+        B = sparse.csc_matrix((localB, (I, I)))
+
+    return A, B

brainprintpython/laplaceTria.py

@@ -0,0 +1,34 @@
+def laplaceTria(v, t, k=10, lump=False):
+    """
+    Compute linear finite-element method Laplace-Beltrami spectrum
+    """
+
+    from scipy.sparse.linalg import LinearOperator, eigsh, splu
+    from brainprintpython.computeABtria import computeABtria
+
+    useCholmod = True
+    try:
+        from sksparse.cholmod import cholesky
+    except ImportError:
+        useCholmod = False
+    if useCholmod:
+        print("Solver: cholesky decomp - performance optimal ...")
+    else:
+        print("Package scikit-sparse not found (Cholesky decomp)")
+        print("Solver: spsolve (LU decomp) - performance not optimal ...")
+
+    A, M = computeABtria(v,t,lump=lump)
+
+    # turns out it is much faster to use cholesky and pass operator
+    sigma=-0.01
+    if useCholmod:
+        chol = cholesky(A-sigma*M)
+        OPinv = LinearOperator(matvec=chol,shape=A.shape, dtype=A.dtype)
+    else:
+        lu = splu(A-sigma*M)
+        OPinv = LinearOperator(matvec=lu.solve,shape=A.shape, dtype=A.dtype)
+
+    eigenvalues, eigenvectors = eigsh(A, k, M, sigma=sigma,OPinv=OPinv)
+
+    return eigenvalues, eigenvectors
+