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/data_structures/arrows.py

@@ -13,13 +13,15 @@ class Arrows:
         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())}"
@@ -41,6 +43,13 @@ def __iter__(self):
     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"""
 

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

@@ -162,7 +162,9 @@ def load_opencarp(
     points_data *= scale_points
     fibres_data = None if fibres is None else np.loadtxt(fibres)
 
-    indices_data, cell_region_data, linear_connection_data = [], [], []
+    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:
@@ -171,14 +173,18 @@ def load_opencarp(
                 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:4])))
+                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_connections=linear_connection_data,
+        linear_connection_regions=linear_connection_regions
     )
 
     fields = Fields(

openep/io/writers.py

@@ -56,7 +56,6 @@
 import scipy.io
 import pathlib
 
-
 from openep.data_structures.ablation import Ablation
 from openep.data_structures.case import Case
 from openep.data_structures.surface import Fields
@@ -109,22 +108,41 @@ def export_openCARP(
         comments='',
     )
 
-    # Save elements info
-    n_lines = case.indices.shape[0]
+    # Total number of lines
+    n_triangles = case.indices.shape[0]
+    n_lin_conns = 0 if case.arrows.linear_connections is None else case.arrows.linear_connections.shape[0]
+    n_lines = n_triangles + n_lin_conns
+
+    # Save triangulation elements info
+    cell_type = np.full(n_triangles, fill_value="Tr")
     cell_region = case.fields.cell_region if case.fields.cell_region is not None else np.zeros(n_triangles, dtype=int)
 
-    elem_data = ""
-    for indices, region in zip(case.indices, cell_region):
-        elem_data += 'Tr ' + ' '.join(map(str, indices)) + ' ' + str(region) + '\n'
+    combined_cell_data = [cell_type[:, np.newaxis], case.indices, cell_region[:, np.newaxis]]
+    combined_cell_data = np.concatenate(combined_cell_data, axis=1, dtype=object)
 
+    np.savetxt(
+        output_path.with_suffix(".elem"),
+        combined_cell_data,
+        fmt="%s %d %d %d %d",
+        header=str(n_lines),
+        comments='',
+    )
+
+    # Save linear connections info
     if case.arrows.linear_connections is not None:
-        n_lines += case.arrows.linear_connections.shape[0]
-        for indices in case.arrows.linear_connections:
-            elem_data += 'Ln ' + ' '.join(map(str, indices)) + '\n'
+        conn_type = np.full(n_lin_conns, fill_value="Ln")
+        ln_region = case.arrows.linear_connection_regions if case.arrows.linear_connection_regions is not None else np.zeros(n_lin_conns, dtype=int)
+
+        combined_ln_conn_data = [conn_type[:, np.newaxis], case.arrows.linear_connections, ln_region[:, np.newaxis]]
+        combined_ln_conn_data = np.concatenate(combined_ln_conn_data, axis=1, dtype=object)
 
-    with open(output_path.with_suffix('.elem'), 'w') as file:
-        file.write(f'{n_lines}\n')
-        file.write(elem_data)
+        # append to the elem file
+        with open(output_path.with_suffix(".elem"), 'a') as elem_file:
+            np.savetxt(
+                elem_file,
+                combined_ln_conn_data,
+                fmt="%s %d %d %d",
+            )
 
     # Save fibres
     if case.arrows.fibres is not None:
@@ -141,7 +159,7 @@ def export_openCARP(
 
     # Ignore all -1 values, as these points do not belong to any pacing site
     for site_index in np.unique(case.fields.pacing_site)[1:]:
-        
+
         pacing_site_points = np.nonzero(case.fields.pacing_site == site_index)[0]
         n_points = pacing_site_points.size
 
@@ -240,7 +258,9 @@ def export_openep_mat(
             divergence=case.analyse.divergence
         )
 
-    userdata['rf'] = _export_ablation_data(ablation=case.ablation)
+    if case.ablation is not None:
+        userdata['rf'] = _export_ablation_data(ablation=case.ablation)
+
     scipy.io.savemat(
         file_name=filename,
         mdict={'userdata': userdata},
@@ -276,7 +296,7 @@ def _add_surface_maps(surface_data, **kwargs):
 
 
 def export_csv(
-        case: Case,
+        system,
         filename: str,
         selections: dict,
 ):
@@ -286,10 +306,17 @@ def export_csv(
     Saves selected field data.
 
     Args:
-        case (Case): dataset to be exported
+        system (model.system_manager.System): System with dataset to be exported
         filename (str): name of file to be written
         selections (dict): the data field name and flag whether to export or not
     """
+    case = system.case
+    mesh = system.mesh
+
+    point_region = _convert_cell_to_point(
+        cell_data=case.fields.cell_region,
+        mesh=mesh
+    )
 
     available_exports = {
         'Bipolar voltage': case.fields.bipolar_voltage,
@@ -298,10 +325,11 @@ def export_csv(
         'Force': case.fields.force,
         'Impedance': case.fields.impedance,
         'Thickness': case.fields.thickness,
-        'Cell region': case.fields.cell_region,
+        'Cell region': point_region,
         'Pacing site': case.fields.pacing_site,
         'Conduction velocity': case.fields.conduction_velocity,
         'Divergence': case.fields.cv_divergence,
+        'Histogram': case.fields.histogram,
     }
 
     df = pd.DataFrame()
@@ -555,4 +583,19 @@ def _export_ablation_data(ablation: Ablation):
     ablation_data['originaldata']['force']['lateralangle'] = ablation.force.lateral_angle if ablation.force.lateral_angle is not None else empty_float_array
     ablation_data['originaldata']['force']['position'] = ablation.force.points if ablation.force.points is not None else empty_float_array
 
-    return ablation_data
+    return ablation_data
+
+
+def _convert_cell_to_point(cell_data, mesh):
+    """Convert cell data to point data by applying the cell value to its vertices"""
+    point_data = np.zeros(mesh.n_points, dtype=int)
+    for cell_i in range(mesh.n_cells):
+
+        cell_value = cell_data[cell_i]
+        point_ids = mesh.get_cell(cell_i).point_ids
+
+        # Add cell value to point ids
+        for pid in point_ids:
+            point_data[pid] = cell_value
+
+    return point_data