SM7: Read lon file (#1)

* 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)

Author: vinush-vignes

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/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
 

openep/io/readers.py

@@ -53,6 +53,7 @@
 """
 
 import os
+import re
 import scipy.io
 
 import numpy as np
@@ -63,9 +64,10 @@
 from ..data_structures.surface import extract_surface_data, Fields
 from ..data_structures.electric import extract_electric_data, Electric
 from ..data_structures.ablation import extract_ablation_data, Ablation
+from ..data_structures.arrows import Arrows
 from ..data_structures.case import Case
 
-__all__ = ["load_openep_mat", "_load_mat", "load_opencarp", "load_circle_cvi", "load_vtk"]
+__all__ = ["load_openep_mat", "_load_mat", "load_opencarp", "load_circle_cvi", "load_vtk", "load_igb"]
 
 
 def _check_mat_version_73(filename):
@@ -158,28 +160,43 @@ def load_opencarp(
 
     points_data = np.loadtxt(points, skiprows=1)
     points_data *= scale_points
-    indices_data = np.loadtxt(indices, skiprows=1, usecols=[1, 2, 3], dtype=int)
-    cell_region = np.loadtxt(indices, skiprows=1, usecols=4, dtype=int)
-
-    longitudinal_fibres = None
-    transverse_fibres = None
-    if fibres is not None:
-        fibres_data = np.loadtxt(fibres, skiprows=1, dtype=float)
-        longitudinal_fibres = fibres_data[:, :3]
-        if fibres_data.shape[1] == 6:
-            transverse_fibres = fibres_data[:, 3:]
-
-    # Create empty data structures for pass to Case
+    fibres_data = None if fibres is None else np.loadtxt(fibres)
+
+    indices_data, cell_region_data = [], []
+    linear_connection_data, linear_connection_regions = [], []
+
+    with open(indices) as elem_file:
+        data = elem_file.readlines()
+        for elem in data:
+            parts = elem.strip().split()
+            if parts[0] == 'Tr':
+                indices_data.append(list(map(int, parts[1:4])))
+                cell_region_data.append(int(parts[4]))
+            elif parts[0] == 'Ln':
+                linear_connection_data.append(list(map(int, parts[1:3])))
+                linear_connection_regions.append(int(parts[3]))
+
+    indices_data = np.array(indices_data)
+    cell_region = np.array(cell_region_data)
+    linear_connection_data = np.array(linear_connection_data)
+    linear_connection_regions = np.array(linear_connection_regions)
+
+    arrows = Arrows(
+        fibres=fibres_data,
+        linear_connections=linear_connection_data,
+        linear_connection_regions=linear_connection_regions
+    )
+
     fields = Fields(
         cell_region=cell_region,
-        longitudinal_fibres=longitudinal_fibres,
-        transverse_fibres=transverse_fibres,
+        longitudinal_fibres=None,
+        transverse_fibres=None,
     )
     electric = Electric()
     ablation = Ablation()
     notes = np.asarray([], dtype=object)
 
-    return Case(name, points_data, indices_data, fields, electric, ablation, notes)
+    return Case(name, points_data, indices_data, fields, electric, ablation, notes, arrows)
 
 
 def load_vtk(filename, name=None):
@@ -281,3 +298,43 @@ def load_circle_cvi(filename, dicoms_directory, extract_epi=True, extract_endo=T
         return epi_mesh
     elif extract_endo:
         return endo_mesh
+
+
+def load_igb(igb_filepath):
+    """
+    Reads an .igb file, returning the data and header information.
+
+    Args:
+        igb_filepath (str): Path to the .igb file.
+
+    Returns:
+        tuple:
+            - numpy.ndarray: 2D array of the file's data.
+            - dict: Contents of the header including 't' value (time steps) and other parameters.
+    """
+    with open(igb_filepath, 'rb') as file:
+        header = file.read(1024).decode('utf-8')
+        header = header.replace('\r', ' ').replace('\n', ' ').replace('\0', ' ')
+        hdr_content = {}
+
+        # Parse the header to dict format
+        for part in header.split():
+            key, value = part.split(':')
+            if key in ['x', 'y', 'z', 't', 'bin', 'num', 'lut', 'comp']:
+                hdr_content[key] = int(value)
+            elif key in ['facteur','zero','epais'] or key.startswith('org_') or key.startswith('dim_') or key.startswith('inc_'):
+                hdr_content[key] = float(value)
+            else:
+                hdr_content[key] = value
+
+        # Process file data
+        words = header.split()
+        word = [int(re.split(r"(\d+)", w)[1]) for w in words[:4]]
+        nnode = word[0] * word[1] * word[2]
+        size = os.path.getsize(igb_filepath) // 4 // nnode
+
+        file.seek(1024)
+        data = np.fromfile(file, dtype=np.float32, count=size * nnode)
+        data = data.reshape((size, nnode)).transpose()
+
+    return data, hdr_content