OE37 - vtk loading error due to CV and Divergence (#227)

* vtk load restored by stopping CV calc with vtk

* OE38-BUG: signalMaps index error (#228)

* added index error to signal maps import

* added index error to signal maps import

* handle empty ablation (#229)

* SM7: Read lon file (#1) (#231)

* created arrows class to store lines + vector visuals

* read lon to obtain fibres data

* write lon and elem files

* check linear connections and arrows exist before exporting

* create empty arrows class if arrows not provided

* export vtx files per pacing site

* Feature/SM-31: IGB data load (#2)

* load igb data

* added load igb to imports

* fibres filled with dummy values

* fibres filled with dummy values

* fibres filled with dummy values

* Feature/SM-45: Write to CSV (#3)

* export CSV writing to file using Panda

* Bugfix/SM-63 Empty signalmaps (#4)

* load empty signalMaps

* Feature/SM-18: Linear connections writer (#5)

* openCARP reader saves linear connection regions as separate data structure

* faster write function for openCARP data

* Feature/SM-45: Export CSV cell and histogram regions (#6)

* include histogram and cell region in csv export

* Bugfix/SM-86: Free boundaries as one actor (#8)

* combine free boundaries into a single actor

* not add rf to export openep if ablation == None (#9)

openep/__init__.py

@@ -18,7 +18,7 @@
 
 __all__ = ['case', 'mesh', 'draw']
 
-from .io.readers import load_openep_mat, load_opencarp, load_circle_cvi, load_vtk
+from .io.readers import load_openep_mat, load_opencarp, load_circle_cvi, load_vtk, load_igb
 from .io.writers import export_openCARP, export_openep_mat, export_vtk
 from .converters.pyvista_converters import from_pyvista, to_pyvista
 from . import case, mesh, draw

openep/analysis/analyse.py

@@ -17,7 +17,6 @@
 # with this program (LICENSE.txt).  If not, see <http://www.gnu.org/licenses/>
 
 """Module containing analysis classes"""
-
 from ._conduction_velocity import *
 from ..case.case_routines import interpolate_general_cloud_points_onto_surface
 
@@ -59,22 +58,8 @@ class ConductionVelocity:
     """
     def __init__(self, case):
         self._case = case
-        self._values = None
-        self._centers = None
-
-    @property
-    def values(self):
-        if self._values is None:
-            raise ValueError('Before accessing ``conduction_velocity.values`` '
-                             'run ``divergence.calculate_divergence()``')
-        return self._values
-
-    @property
-    def centers(self):
-        if self._centers is None:
-            raise ValueError('Before accessing ``conduction_velocity.centers`` '
-                             'run ``divergence.calculate_divergence()``')
-        return self._centers
+        self.values = None
+        self.centers = None
 
     def calculate_cv(
             self,
@@ -134,12 +119,19 @@ def calculate_cv(
         if method.lower() not in supported_cv_methods:
             raise ValueError(f"`method` must be one of {supported_cv_methods.keys()}.")
 
+        if include is None and self._case.electric.include is None:
+            raise TypeError(f"include object is of 'NoneType'")
+        else:
+            include = self._case.electric.include.astype(bool) if include is None else include
+
         interpolation_kws = dict() if interpolation_kws is None else interpolation_kws
-        include = self._case.electric.include.astype(bool) if include is None else include
-        lat, bipolar_egm_pts = preprocess_lat_egm(self._case, include)
+
+        bipolar_egm_pts = self._case.electric.bipolar_egm.points[include]
+        lat = (self._case.electric.annotations.local_activation_time[include]
+               - self._case.electric.annotations.reference_activation_time[include])
 
         cv_method = supported_cv_methods[method]
-        self._values, self._centers = cv_method(bipolar_egm_pts, lat, **method_kwargs)
+        self.values, self.centers = cv_method(bipolar_egm_pts, lat, **method_kwargs)
 
         if apply_scalar_field:
             self._case.fields.conduction_velocity = interpolate_general_cloud_points_onto_surface(
@@ -169,20 +161,8 @@ class Divergence:
     """
     def __init__(self, case):
         self._case = case
-        self._direction = None
-        self._values = None
-
-    @property
-    def values(self):
-        if self._values is None:
-            raise ValueError('Before accessing ``..divergence.values`` run ``..divergence.calculate_divergence()``')
-        return self._values
-
-    @property
-    def direction(self):
-        if self._direction is None:
-            raise ValueError('Before accessing ``divergence.direction`` run ``divergence.calculate_divergence()``')
-        return self._direction
+        self.direction = None
+        self.values = None
 
     def calculate_divergence(
         self,
@@ -224,10 +204,16 @@ def calculate_divergence(
             direction, values = cv.calculate_divergence(output_binary_field=True, apply_scalar_field=True)
         """
 
-        include = self._case.electric.include.astype(bool) if include is None else include
-        lat, bipolar_egm_pts = preprocess_lat_egm(self._case, include, lat_threshold=None)
+        if include is None and self._case.electric.include is None:
+            raise TypeError(f"include object is of 'NoneType'")
+        else:
+            include = self._case.electric.include.astype(bool) if include is None else include
+
+        bipolar_egm_pts = self._case.electric.bipolar_egm.points[include]
+        lat = (self._case.electric.annotations.local_activation_time[include]
+               - self._case.electric.annotations.reference_activation_time[include])
 
-        self._direction, self._values = divergence(
+        self.direction, self.values = divergence(
             case=self._case,
             bipolar_egm_pts=bipolar_egm_pts,
             local_activation_time=lat,

openep/data_structures/ablation.py

@@ -114,7 +114,13 @@ def extract_ablation_data(ablation_data):
             as well as the force applied.
     """
 
-    if isinstance(ablation_data, np.ndarray) or ablation_data['originaldata']['ablparams']['time'].size == 0:
+    if isinstance(ablation_data, np.ndarray):
+        return Ablation()
+
+    try:
+        if not ablation_data['originaldata']['ablparams']['time'].size:
+            return Ablation()
+    except KeyError as e:
         return Ablation()
 
     times = ablation_data['originaldata']['ablparams']['time'].astype(float)

openep/data_structures/arrows.py

@@ -0,0 +1,62 @@
+from attr import attrs
+import numpy as np
+
+__all__ = []
+
+
+@attrs(auto_attribs=True, auto_detect=True)
+class Arrows:
+    """
+    Class for storing information about arrows and lines on surface
+
+    Args:
+        fibres (np.ndarray): array of shape N_cells x 3
+        divergence (np.ndarray): array of shape N_cells x 3
+        linear_connections (np.ndarray): array of shape M x 3 (represents the linear connections between endo and epi)
+        linear_connection_regions (np.ndarray): array of shape N_cells
+    """
+
+    # TODO: move divergence arrows into Arrows class
+    # TODO: remove longitudinal and transversal arrows from Fields class
+    fibres: np.ndarray = None
+    divergence: np.ndarray = None
+    linear_connections: np.ndarray = None
+    linear_connection_regions: np.ndarray = None
+
+    def __repr__(self):
+        return f"arrows: {tuple(self.__dict__.keys())}"
+
+    def __getitem__(self, arrow):
+        try:
+            return self.__dict__[arrow]
+        except KeyError:
+            raise ValueError(f"There is no arrow '{arrow}'.")
+
+    def __setitem__(self, arrow, value):
+        if arrow not in self.__dict__.keys():
+            raise ValueError(f"'{arrow}' is not a valid arrow name.")
+        self.__dict__[arrow] = value
+
+    def __iter__(self):
+        return iter(self.__dict__.keys())
+
+    def __contains__(self, arrow):
+        return arrow in self.__dict__.keys()
+
+    @property
+    def linear_connection_regions_names(self):
+        if self.linear_connection_regions is None:
+            return []
+        regions = np.unique(self.linear_connection_regions).astype(str)
+        return regions.tolist()
+
+    def copy(self):
+        """Create a deep copy of Arrows"""
+
+        arrows = Arrows()
+        for arrow in self:
+            if self[arrow] is None:
+                continue
+            arrows[arrow] = np.array(self[arrow])
+
+        return arrows

openep/data_structures/case.py

@@ -83,6 +83,7 @@
 import pyvista
 
 from .surface import Fields
+from .arrows import Arrows
 from .electric import Electric, Electrogram, Annotations, ElectricSurface
 from .ablation import Ablation
 from ..analysis.analyse import Analyse
@@ -122,6 +123,7 @@ def __init__(
         electric: Electric,
         ablation: Optional[Ablation] = None,
         notes: Optional[List] = None,
+        arrows: Optional[Arrows] = None,
     ):
 
         self.name = name
@@ -131,6 +133,7 @@ def __init__(
         self.ablation = ablation
         self.electric = electric
         self.notes = notes
+        self.arrows = Arrows() if arrows is None else arrows
         self.analyse = Analyse(case=self)
 
     def __repr__(self):

openep/data_structures/surface.py

@@ -36,7 +36,7 @@ class Fields:
         local_activation_time (np.ndarray): array of shape N_points
         impedance (np.ndarray): array of shape N_points
         force (np.ndarray): array of shape N_points
-        region (np.ndarray): array of shape N_cells
+        cell_region (np.ndarray): array of shape N_cells
         longitudinal_fibres (np.ndarray): array of shape N_cells x 3
         transverse_fibres (np.ndarray): array of shape N_cells x 3
         pacing_site (np.ndarray): array of shape N_points
@@ -54,6 +54,7 @@ class Fields:
     pacing_site: np.ndarray = None
     conduction_velocity: np.ndarray = None
     cv_divergence: np.ndarray = None
+    histogram: np.ndarray = None
 
     def __repr__(self):
         return f"fields: {tuple(self.__dict__.keys())}"
@@ -187,18 +188,22 @@ def extract_surface_data(surface_data):
     if isinstance(pacing_site, np.ndarray):
         pacing_site = None if pacing_site.size == 0 else pacing_site.astype(int)
 
-    try:
-        conduction_velocity = surface_data['signalMaps']['conduction_velocity_field'].get('value', None)
-        if isinstance(conduction_velocity, np.ndarray):
-            conduction_velocity = None if conduction_velocity.size == 0 else conduction_velocity.astype(float)
-    except KeyError:
-        conduction_velocity = None
+    if surface_data.get('signalMaps'):
+        try:
+            conduction_velocity = surface_data['signalMaps']['conduction_velocity_field'].get('value', None)
+            if isinstance(conduction_velocity, np.ndarray):
+                conduction_velocity = None if conduction_velocity.size == 0 else conduction_velocity.astype(float)
+        except KeyError:
+            conduction_velocity = None
 
-    try:
-        cv_divergence = surface_data['signalMaps']['divergence_field'].get('value', None)
-        if isinstance(cv_divergence, np.ndarray):
-            cv_divergence = None if cv_divergence.size == 0 else cv_divergence.astype(float)
-    except KeyError:
+        try:
+            cv_divergence = surface_data['signalMaps']['divergence_field'].get('value', None)
+            if isinstance(cv_divergence, np.ndarray):
+                cv_divergence = None if cv_divergence.size == 0 else cv_divergence.astype(float)
+        except KeyError:
+            cv_divergence = None
+    else:
+        conduction_velocity = None
         cv_divergence = None
 
     fields = Fields(

openep/draw/draw_routines.py

@@ -63,9 +63,10 @@ def draw_free_boundaries(
     width: int = 5,
     plotter: pyvista.Plotter = None,
     names: List[str] = None,
+    combine: bool = False
 ):
     """
-    Draw the freeboundaries of a mesh.
+    Draw the free boundaries of a mesh.
 
     Args:
         free_boundaries (FreeBoundary): `FreeBoundary` object. Can be generated using
@@ -76,28 +77,47 @@ def draw_free_boundaries(
             If None, a new plotting object will be created.
         names (List(str)): List of names to associated with the actors. Default is None, in which
             case actors will be called 'free_boundary_n', where n is the index of the boundary.
+        combine (bool): Combines all free boundaries into one Actor (faster load time).
 
     Returns:
         plotter (pyvista.Plotter): Plotting object with the free boundaries added.
-
     """
-
+    combined_lines = pyvista.PolyData() if combine else None
     plotter = pyvista.Plotter() if plotter is None else plotter
     colours = [colour] * free_boundaries.n_boundaries if isinstance(colour, str) else colour
+
     if names is None:
         names = [f"free_boundary_{boundary_index:d}" for boundary_index in range(free_boundaries.n_boundaries)]
 
     for boundary_index, boundary in enumerate(free_boundaries.separate_boundaries()):
 
         points = free_boundaries.points[boundary[:, 0]]
         points = np.vstack([points, points[:1]])  # we need to close the loop
-        plotter.add_lines(
-            points,
-            color=colours[boundary_index],
-            width=width,
-            name=names[boundary_index],
-            connected=True
-        )
+
+        # store the lines to be added in later
+        if combine:
+            lines = pyvista.lines_from_points(points)
+            combined_lines = combined_lines.merge(lines)
+
+        # add each line one-by-one
+        else:
+            plotter.add_lines(
+                points,
+                color=colours[boundary_index],
+                width=width,
+                name=names[boundary_index],
+                connected=True
+            )
+
+    if combine:
+        # add the combined lines as a single Actor manually (modified code of add_lines)
+        actor = pyvista.Actor(mapper=pyvista.DataSetMapper(combined_lines))
+        actor.prop.line_width = width
+        actor.prop.show_edges = True
+        actor.prop.edge_color = colour
+        actor.prop.color = colour
+        actor.prop.lighting = False
+        plotter.add_actor(actor, reset_camera=False, name=names[0], pickable=False)
 
     return plotter