Source code for navis.core.volumes

#    This script is part of navis (http://www.github.com/navis-org/navis).
#    Copyright (C) 2018 Philipp Schlegel
#
#    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.

import csv
import functools
import json
import math
import numbers
import os
import uuid

import numpy as np
import scipy.spatial
import trimesh

from typing import Union, Optional, Sequence, List, Dict, Any
from typing_extensions import Literal

from .. import utils, config

# Set up logging
logger = config.get_logger(__name__)


[docs] class Volume(trimesh.Trimesh): """Mesh consisting of vertices and faces. Subclass of ``trimesh.Trimesh`` with a few additional methods. Parameters ---------- vertices : list | array | mesh-like `(N, 3)` vertices coordinates or an object that has ``.vertices`` and ``.faces`` attributes in which case ``faces`` parameter will be ignored. faces : list | array `(M, 3)` array of indexed triangle faces. name : str, optional A name for the volume. color : tuple, optional RGB(A) color. id : int, optional If not provided, neuron will be assigned a random UUID as ``.id``. **kwargs Keyword arguments passed through to ``trimesh.Trimesh`` See Also -------- :func:`~navis.example_volume` Loads example volume(s). """
[docs] def __init__(self, vertices: Union[list, np.ndarray], faces: Union[list, np.ndarray] = None, name: Optional[str] = None, color: Union[str, Sequence[Union[int, float]]] = (.85, .85, .85, .2), id: Optional[int] = None, **kwargs): if hasattr(vertices, 'vertices') and hasattr(vertices, 'faces'): vertices, faces = vertices.vertices, vertices.faces super().__init__(vertices=vertices, faces=faces, **kwargs) self.name: Optional[str] = name self.color: Union[str, Sequence[Union[int, float]]] = color self.id: Optional[int] = id if id else uuid.uuid4() # This is very hackish but we want to make sure that parent methods of # Trimesh return a navis.Volume instead of a trimesh.Trimesh for f in dir(trimesh.Trimesh): # Don't mess with magic/private methods if f.startswith('_'): continue # Skip properties if not callable(getattr(trimesh.Trimesh, f)): continue setattr(self, f, _force_volume(getattr(self, f)))
@property def name(self): """Name of this volume.""" return self.metadata.get('name') @name.setter def name(self, value): self.metadata['name'] = value @property def color(self): """Color used for plotting.""" return self.metadata.get('color') @color.setter def color(self, value): self.metadata['color'] = value @property def id(self): """ID of this volume.""" return self.metadata.get('id') @id.setter def id(self, value): self.metadata['id'] = value @classmethod def from_csv(cls, vertices: str, faces: str, name: Optional[str] = None, color: Union[str, Sequence[Union[int, float]]] = (.85, .85, .85, .2), volume_id: Optional[int] = None, **kwargs) -> 'Volume': """Load volume from csv files containing vertices and faces. Parameters ---------- vertices : filepath | file-like CSV file containing vertices. faces : filepath | file-like CSV file containing faces. **kwargs Keyword arguments passed to ``csv.reader``. Returns ------- navis.Volume """ if not os.path.isfile(vertices) or not os.path.isfile(faces): raise ValueError('File(s) not found.') with open(vertices, 'r') as f: reader = csv.reader(f, **kwargs) vertices = np.array([r for r in reader]).astype(float) with open(faces, 'r') as f: reader = csv.reader(f, **kwargs) faces = np.array([r for r in reader]).astype(int) return cls(faces=faces, vertices=vertices, name=name, color=color, volume_id=volume_id) def to_csv(self, filename: str, **kwargs) -> None: """Save volume as two separated csv files containing vertices and faces. Parameters ---------- filename : str Filename to use. Will get a ``_vertices.csv`` and ``_faces.csv`` suffix. **kwargs Keyword arguments passed to ``csv.reader``. """ for data, suffix in zip([self.faces, self.vertices], ['_faces.csv', '_vertices.csv']): with open(filename + suffix, 'w') as csvfile: writer = csv.writer(csvfile) writer.writerows(data) @classmethod def from_json(cls, filename: str, import_kwargs: Dict = {}, **init_kwargs) -> 'Volume': """Load volume from json file containing vertices and faces. Parameters ---------- filename import_kwargs Keyword arguments passed to ``json.load``. **init_kwargs Keyword arguments passed to navis.Volume upon initialization. Returns ------- navis.Volume """ if not os.path.isfile(filename): raise ValueError('File not found.') with open(filename, 'r') as f: data = json.load(f, **import_kwargs) return cls(faces=data['faces'], vertices=data['vertices'], **init_kwargs) @classmethod def from_object(cls, obj: Any, **init_kwargs) -> 'Volume': """Load volume from generic object that has ``.vertices`` and ``.faces`` attributes. Parameters ---------- obj **init_kwargs Keyword arguments passed to navis.Volume upon initialization. Returns ------- navis.Volume """ if not hasattr(obj, 'vertices') or not hasattr(obj, 'faces'): raise ValueError('Object must have faces and vertices attributes.') return cls(faces=obj.faces, vertices=obj.vertices, **init_kwargs) @classmethod def from_file(cls, filename: str, import_kwargs: Dict = {}, **init_kwargs) -> 'Volume': """Load volume from file. Parameters ---------- filename : str File to load from. import_kwargs Keyword arguments passed to importer: - ``json.load`` for JSON file - ``trimesh.load_mesh`` for OBJ and STL files **init_kwargs Keyword arguments passed to navis.Volume upon initialization. Returns ------- navis.Volume """ if not os.path.isfile(filename): raise ValueError('File not found.') f, ext = os.path.splitext(filename) if ext == '.json': return cls.from_json(filename=filename, import_kwargs=import_kwargs, **init_kwargs) try: import trimesh except ImportError: raise ImportError('Unable to import: trimesh missing - please ' 'install: "pip install trimesh"') except BaseException: raise tm = trimesh.load_mesh(filename, **import_kwargs) return cls.from_object(tm, **init_kwargs) def to_json(self, filename: str) -> None: """Save volume as json file. Parameters ---------- filename : str Filename to use. """ with open(filename, 'w') as f: json.dump({'vertices': self.vertices.tolist(), 'faces': self.faces.tolist()}, f)
[docs] @classmethod def combine(cls, x: Sequence['Volume'], name: str = 'comb_vol', color: Union[str, Sequence[Union[int, float]]] = (.85, .85, .85, .2) ) -> 'Volume': """Merge multiple volumes into a single object. Parameters ---------- x : list or dict of Volumes name : str, optional Name of the combined volume. color : tuple | str, optional Color of the combined volume. Returns ------- :class:`~navis.Volume` """ if isinstance(x, Volume): return x if isinstance(x, dict): x = list(x.values()) if not utils.is_iterable(x): x = [x] # type: ignore if False in [isinstance(v, Volume) for v in x]: raise TypeError('Input must be list of volumes') vertices: np.ndarray = np.empty((0, 3)) faces: List[List[int]] = [] # Reindex faces for vol in x: offs = len(vertices) vertices = np.append(vertices, vol.vertices, axis=0) faces += [[f[0] + offs, f[1] + offs, f[2] + offs] for f in vol.faces] return cls(vertices=vertices, faces=faces, name=name, color=color)
@property def bbox(self) -> np.ndarray: """Bounding box of this volume.""" return self.bounds @property def verts(self) -> np.ndarray: """Legacy access to ``.vertices``.""" return self.vertices @verts.setter def verts(self, v): self.vertices = v @property def center(self) -> np.ndarray: """Center of mass.""" return np.mean(self.vertices, axis=0) def __getstate__(self): """Get state (used e.g. for pickling).""" return {k: v for k, v in self.__dict__.items() if not callable(v)} def __setstate__(self, d): """Update state (used e.g. for pickling).""" self.__dict__.update(d) def __str__(self): return self.__repr__() def __repr__(self): """Return quick summary of the current geometry. Avoids computing properties. """ elements = [] if hasattr(self, 'name'): # for Trimesh elements.append(f'name={self.name}') if hasattr(self, 'id') and not isinstance(self.id, uuid.UUID): # for Trimesh elements.append(f'id={self.id}') if hasattr(self, 'color'): # for Trimesh elements.append(f'color={self.color}') if hasattr(self, 'vertices'): # for Trimesh and PointCloud elements.append(f'vertices.shape={self.vertices.shape}') if hasattr(self, 'faces'): # for Trimesh elements.append(f'faces.shape={self.faces.shape}') return f'<navis.Volume({", ".join(elements)})>' def __truediv__(self, other): """Implement division for vertices.""" if isinstance(other, numbers.Number) or utils.is_iterable(other): n = self.copy() _ = np.divide(n.vertices, other, out=n.vertices, casting='unsafe') return n return NotImplemented def __mul__(self, other): """Implement multiplication for vertices.""" if isinstance(other, numbers.Number) or utils.is_iterable(other): n = self.copy() _ = np.multiply(n.vertices, other, out=n.vertices, casting='unsafe') return n return NotImplemented
[docs] def resize(self, x: Union[float, int], method: Union[Literal['center'], Literal['centroid'], Literal['normals'], Literal['origin']] = 'center', inplace: bool = False) -> Optional['Volume']: """Resize volume. Parameters ---------- x : int | float Resizing factor. For methods "center", "centroid" and "origin" this is the fraction of original size (e.g. ``.5`` for half size). For method "normals", this is is the absolute displacement (e.g. ``-1000`` to shrink volume by that many units)! method : "center" | "centroid" | "normals" | "origin" Point in space to use for resizing. .. list-table:: :widths: 15 75 :header-rows: 1 * - method - explanation * - center - average of all vertices * - centroid - average of the triangle centroids weighted by the area of each triangle. * - origin - resizes relative to origin of coordinate system (0, 0, 0) * - normals - resize using face normals. Note that this method uses absolute displacement for parameter ``x``. inplace : bool, optional If False, will return resized copy. Returns ------- :class:`navis.Volume` Resized copy of original volume. Only if ``inplace=False``. None If ``inplace=True``. """ assert isinstance(method, str) method = method.lower() perm_methods = ['center', 'origin', 'normals', 'centroid'] if method not in perm_methods: raise ValueError(f'Unknown method "{method}". Allowed ' f'methods: {", ".join(perm_methods)}') if not inplace: v = self.copy() else: v = self if method == 'normals': v.vertices = v.vertices + (v.vertex_normals * x) else: # Get the center if method == 'center': cn = np.mean(v.vertices, axis=0) elif method == 'centroid': cn = v.centroid elif method == 'origin': cn = np.array([0, 0, 0]) # Get vector from center to each vertex vec = v.vertices - cn # Multiply vector by resize factor vec *= x # Recalculate vertex positions v.vertices = vec + cn # Make sure to reset any pyoctree data on this volume if hasattr(v, 'pyoctree'): delattr(v, 'pyoctree') if not inplace: return v
[docs] def plot3d(self, **kwargs): """Plot volume using :func:`navis.plot3d`. Parameters ---------- **kwargs Keyword arguments. Will be passed to :func:`navis.plot3d`. See ``help(navis.plot3d)`` for a list of keywords. See Also -------- :func:`navis.plot3d` Function called to generate 3d plot. Examples -------- >>> import navis >>> vol = navis.example_volume('LH') >>> v = vol.plot3d(color = (255, 0, 0)) """ from .. import plotting if 'color' in kwargs: self.color = kwargs['color'] return plotting.plot3d(self, **kwargs)
def show(self, **kwargs): """See ``.plot3d``.""" # This is mostly to override trimesh.Trimesh method return self.plot3d(**kwargs) def _outlines_3d(self, view='xy', **kwargs): """Generate 3d outlines along a given view (see ``.to_2d()``). Parameters ---------- **kwargs Keyword arguments passed to :func:`~navis.Volume.to_2d`. Returns ------- list Coordinates of 2d circumference. e.g. ``[(x1, y1, z1), (x2, y2, z2), (x3, y3, z3), ...]`` Third dimension is averaged. """ co2d = np.array(self.to_2d(view=view, **kwargs)) if view in ['xy', 'yx']: third = np.repeat(self.center[2], co2d.shape[0]) elif view in ['xz', 'zx']: third = np.repeat(self.center[1], co2d.shape[0]) elif view in ['yz', 'zy']: third = np.repeat(self.center[0], co2d.shape[0]) return np.append(co2d, third.reshape(co2d.shape[0], 1), axis=1) def to_2d(self, alpha: float = 0.00017, view: tuple = ('x', 'y'), invert_y: bool = False) -> Sequence[Union[float, int]]: """Compute the 2d alpha shape (concave hull) this volume. Uses Scipy Delaunay and shapely. Parameters ---------- alpha: float, optional Alpha value to influence the gooeyness of the border. Smaller numbers don't fall inward as much as larger numbers. Too large, and you lose everything! view : tuple Determines axis. Can be prefixed with a '-' to invert the axis. Returns ------- list Coordinates of 2d circumference e.g. ``[(x1, y1), (x2, y2), (x3, y3), ...]`` """ def add_edge(edges, edge_points, coords, i, j): """Add line between the i-th and j-th points.""" if (i, j) in edges or (j, i) in edges: # already added return edges.add((i, j)) edge_points.append(coords[[i, j]]) accepted_views = ['x', 'z', 'y', '-x', '-z', '-y'] for ax in view: if ax not in accepted_views: raise ValueError(f'Unable to parse "{ax}" view') try: from shapely.ops import unary_union, polygonize # type: ignore import shapely.geometry as geometry # type: ignore except ImportError: raise ImportError('This function needs the shapely>=1.8.0') coords: np.ndarray map = {'x': 0, 'y': 1, 'z': 2} x_ix = map[view[0].replace('-', '').replace('+', '')] y_ix = map[view[1].replace('-', '').replace('+', '')] xmod = -1 if '-' in view[0] else 1 ymod = -1 if '-' in view[1] else 1 coords = self.vertices[:, [x_ix, y_ix]] * np.array([xmod, ymod]) tri = scipy.spatial.Delaunay(coords) edges: set = set() edge_points: list = [] # loop over triangles: # ia, ib, ic = indices of corner points of the triangle # Note that "vertices" property was renamed to "simplices" for ia, ib, ic in getattr(tri, 'simplices', getattr(tri, 'vertices', [])): pa: np.ndarray = coords[ia] # type: ignore pb: np.ndarray = coords[ib] # type: ignore pc: np.ndarray = coords[ic] # type: ignore # Lengths of sides of triangle a = math.sqrt((pa[0] - pb[0])**2 + (pa[1] - pb[1])**2) # type: ignore b = math.sqrt((pb[0] - pc[0])**2 + (pb[1] - pc[1])**2) # type: ignore c = math.sqrt((pc[0] - pa[0])**2 + (pc[1] - pa[1])**2) # type: ignore # Semiperimeter of triangle s = (a + b + c) / 2.0 # Area of triangle by Heron's formula area = math.sqrt(s * (s - a) * (s - b) * (s - c)) circum_r = a * b * c / (4.0 * area) # Here's the radius filter. if circum_r < 1.0 / alpha: add_edge(edges, edge_points, coords, ia, ib) add_edge(edges, edge_points, coords, ib, ic) add_edge(edges, edge_points, coords, ic, ia) m = geometry.MultiLineString(edge_points) triangles = list(polygonize(m)) concave_hull = unary_union(triangles) # Try with current settings, if this is not successful, try again # with lower alpha try: return list(concave_hull.exterior.coords) except AttributeError: return self.to_2d(alpha=alpha / 10, view=view, invert_y=invert_y) except BaseException: raise
[docs] def validate(self): """Use trimesh to try and fix issues (holes/normals).""" if not self.is_volume: logger.info("Mesh not valid, attempting to fix") self.fill_holes() self.fix_normals() if not self.is_volume: raise utils.VolumeError("Mesh is not a volume " "(e.g. not watertight, incorrect " "winding) and could not be fixed.")
def _force_volume(f): """Convert result from ``trimesh.Trimesh`` to ``navis.Volume``.""" @functools.wraps(f) def wrapper(*args, **kwargs): res = f(*args, **kwargs) if isinstance(res, trimesh.Trimesh): res = Volume(res.vertices, res.faces, **res.metadata) return res return wrapper