🔧 redo install to follow PEP517 standard

Author: vanandrew

README.md

@@ -13,14 +13,17 @@ Smith SM. Fast robust automated brain extraction. Hum Brain Mapp.
 To install, simply use `pip` to install this repo:
 
 ```
+# install from pypispheresphere
+pip install brainextractor
+
 # install repo with pip
 pip install git+https://github.com/vanandrew/brainextractor@main
 
 # install from local copy
 pip install /path/to/local/repo
 ```
 
-Note that is reccomended to use `brainextractor` on python 3.7+
+Note that it is recommended to use `brainextractor` on python 3.7+
 
 ## Usage
 

brainextractor/__init__.py

@@ -1 +1 @@
-from .main import BrainExtractor
+from .main import BrainExtractor

brainextractor/helpers.py

@@ -5,9 +5,10 @@
 import trimesh
 from numba import jit
 
+
 def sphere(shape: list, radius: float, position: list):
     """
-        Creates a binary sphere
+    Creates a binary sphere
     """
     # assume shape and position are both a 3-tuple of int or float
     # the units are pixels / voxels (px for short)
@@ -30,90 +31,22 @@ def sphere(shape: list, radius: float, position: list):
     # the inner part of the sphere will have distance below 1
     return arr <= 1.0
 
-def cartesian(arrays, out=None):
-    """
-        Generate a cartesian product of input arrays
-    """
-    arrays = [np.asarray(x) for x in arrays]
-    dtype = arrays[0].dtype
-
-    n = np.prod([x.size for x in arrays])
-    if out is None:
-        out = np.zeros([n, len(arrays)], dtype=dtype)
-
-    m = int(n / arrays[0].size)
-    out[:,0] = np.repeat(arrays[0], m)
-    if arrays[1:]:
-        cartesian(arrays[1:], out=out[0:m, 1:])
-        for j in range(1, arrays[0].size):
-            out[j*m:(j+1)*m, 1:] = out[0:m, 1:]
-    return out
-
-def find_enclosure(surface: trimesh.Trimesh, data_shape: tuple):
-    """
-        Finds all voxels inside of a surface
-
-        This function stores all the surface vertices in a k-d tree
-        and uses it to quickly look up the closest vertex to each
-        volume voxel in the image.
-
-        Once the closest vertex is found, a vector is created between
-        the voxel location and the vertex. The resulting vector is dot
-        product with the corresponding vertex normal. Negative values
-        indicate that the voxel lies exterior to the surface (since it
-        is anti-parallel to the vertex normal), while positive values
-        indicate that they are interior to the surface (parallel to
-        the vertex normal).
-    """
-    # get vertex normals for each vertex on the surface
-    normals = surface.vertex_normals
-
-    # get KDTree over surface vertices
-    searcher = surface.kdtree
-
-    # get bounding box around surface
-    max_loc = np.ceil(np.max(surface.vertices, axis=0)).astype(np.int64)
-    min_loc = np.floor(np.min(surface.vertices, axis=0)).astype(np.int64)
-
-    # build a list of locations representing the volume grid
-    # within the bounding box
-    locs = cartesian([
-        np.arange(min_loc[0], max_loc[0]),
-        np.arange(min_loc[1], max_loc[1]),
-        np.arange(min_loc[2], max_loc[2])])
-
-    # find the nearest vertex to each voxel
-    # searcher.query returns a list of vertices corresponding
-    # to the closest vertex to the given voxel location
-    _, nearest_idx = searcher.query(locs, n_jobs=6)
-    nearest_vertices = surface.vertices[nearest_idx]
-
-    # get the directional vector from each voxel location to it's nearest vertex
-    direction_vectors = nearest_vertices - locs
-
-    # find it's direction by taking the dot product with vertex normal
-    # this is done row-by-row between directional vectors and the vertex normals
-    dot_products = np.einsum('ij,ij->i', direction_vectors, normals[nearest_idx])
-
-    # get the interior (where dot product is > 0)
-    interior = (dot_products > 0).reshape((max_loc - min_loc).astype(np.int64))
-
-    # create mask
-    mask = np.zeros(data_shape)
-    mask[min_loc[0]:max_loc[0],min_loc[1]:max_loc[1],min_loc[2]:max_loc[2]] = interior
-
-    # return the mask
-    return mask
 
 @jit(nopython=True, cache=True)
 def closest_integer_point(vertex: np.ndarray):
     """
-        Gives the closest integer point based on euclidean distance
+    Gives the closest integer point based on euclidean distance
     """
     # get neighboring grid points to search
-    x = vertex[0]; y = vertex[1]; z = vertex[2]
-    x0 = np.floor(x); y0 = np.floor(y); z0 = np.floor(z)
-    x1 = x0 + 1; y1 = y0 + 1; z1 = z0 + 1
+    x = vertex[0]
+    y = vertex[1]
+    z = vertex[2]
+    x0 = np.floor(x)
+    y0 = np.floor(y)
+    z0 = np.floor(z)
+    x1 = x0 + 1
+    y1 = y0 + 1
+    z1 = z0 + 1
 
     # initialize min euclidean distance
     min_euclid = 99
@@ -132,12 +65,13 @@ def closest_integer_point(vertex: np.ndarray):
     # return the final coords
     return final_coords.astype(np.int64)
 
+
 @jit(nopython=True, cache=True)
 def bresenham3d(v0: np.ndarray, v1: np.ndarray):
     """
-        Bresenham's algorithm for a 3-D line
+    Bresenham's algorithm for a 3-D line
 
-        https://www.geeksforgeeks.org/bresenhams-algorithm-for-3-d-line-drawing/
+    https://www.geeksforgeeks.org/bresenhams-algorithm-for-3-d-line-drawing/
     """
     # initialize axis differences
 
@@ -150,32 +84,50 @@ def bresenham3d(v0: np.ndarray, v1: np.ndarray):
 
     # determine the driving axis
     if dx >= dy and dx >= dz:
-        d0 = dx; d1 = dy; d2 = dz
-        s0 = xs; s1 = ys; s2 = zs
-        a0 = 0; a1 = 1; a2 = 2
+        d0 = dx
+        d1 = dy
+        d2 = dz
+        s0 = xs
+        s1 = ys
+        s2 = zs
+        a0 = 0
+        a1 = 1
+        a2 = 2
     elif dy >= dx and dy >= dz:
-        d0 = dy; d1 = dx; d2 = dz
-        s0 = ys; s1 = xs; s2 = zs
-        a0 = 1; a1 = 0; a2 = 2
+        d0 = dy
+        d1 = dx
+        d2 = dz
+        s0 = ys
+        s1 = xs
+        s2 = zs
+        a0 = 1
+        a1 = 0
+        a2 = 2
     elif dz >= dx and dz >= dy:
-        d0 = dz; d1 = dx; d2 = dy
-        s0 = zs; s1 = xs; s2 = ys
-        a0 = 2; a1 = 0; a2 = 1
+        d0 = dz
+        d1 = dx
+        d2 = dy
+        s0 = zs
+        s1 = xs
+        s2 = ys
+        a0 = 2
+        a1 = 0
+        a2 = 1
 
     # create line array
     line = np.zeros((d0 + 1, 3), dtype=np.int64)
     line[0] = v0
 
     # get points
-    p1 = 2*d1 - d0
-    p2 = 2*d2 - d0
+    p1 = 2 * d1 - d0
+    p2 = 2 * d2 - d0
     for i in range(d0):
         c = line[i].copy()
         c[a0] += s0
-        if (p1 >= 0):
+        if p1 >= 0:
             c[a1] += s1
             p1 -= 2 * d0
-        if (p2 >= 0):
+        if p2 >= 0:
             c[a2] += s2
             p2 -= 2 * d0
         p1 += 2 * d1
@@ -185,27 +137,30 @@ def bresenham3d(v0: np.ndarray, v1: np.ndarray):
     # return list
     return line
 
+
 @jit(nopython=True, cache=True)
 def l2norm(vec: np.ndarray):
     """
-        Computes the l2 norm for 3d vector
+    Computes the l2 norm for 3d vector
     """
-    return np.sqrt(vec[0]**2 + vec[1]**2 + vec[2]**2)
+    return np.sqrt(vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2)
+
 
 @jit(nopython=True, cache=True)
 def l2normarray(array: np.ndarray):
     """
-        Computes the l2 norm for several 3d vectors
+    Computes the l2 norm for several 3d vectors
     """
-    return np.sqrt(array[:, 0]**2 + array[:, 1]**2 + array[:, 2]**2)
+    return np.sqrt(array[:, 0] ** 2 + array[:, 1] ** 2 + array[:, 2] ** 2)
+
 
 def diagonal_dot(a: np.ndarray, b: np.ndarray):
     """
-        Dot product by row of a and b.
-        There are a lot of ways to do this though
-        performance varies very widely. This method
-        uses a dot product to sum the row and avoids
-        function calls if at all possible.
+    Dot product by row of a and b.
+    There are a lot of ways to do this though
+    performance varies very widely. This method
+    uses a dot product to sum the row and avoids
+    function calls if at all possible.
     """
     a = np.asanyarray(a)
-    return np.dot(a * b, [1.0] * a.shape[1])
+    return np.dot(a * b, [1.0] * a.shape[1])