|
| 1 | +# Code written by Jerin Rajan on 18th Mar 2021 |
| 2 | +# Functions to support OpenEp visualisation |
| 3 | + |
| 4 | +import scipy.io as sio |
| 5 | +from scipy.spatial import ConvexHull, Delaunay |
| 6 | +from scipy import ndimage as ndi |
| 7 | +import numpy as np |
| 8 | +import matplotlib.pyplot as plt |
| 9 | +import matplotlib as mp |
| 10 | +from matplotlib import cm |
| 11 | +from matplotlib.colors import ListedColormap, LinearSegmentedColormap |
| 12 | +from skimage import feature |
| 13 | +import trimesh as tm |
| 14 | +# from matplotlib.tri import Triangulation |
| 15 | + |
| 16 | + |
| 17 | +# Declarations |
| 18 | +file_path = '../openep/openep-examples-main/openep_dataset_1.mat' |
| 19 | +data_tri_X_lst = [] |
| 20 | +data_tri_T_lst = [] |
| 21 | +x_values = [] |
| 22 | +y_values = list() |
| 23 | +z_values = list() |
| 24 | +# voltage threshold to split the colorbar by solid and JET |
| 25 | +volt_threshold = 2 |
| 26 | + |
| 27 | + |
| 28 | + |
| 29 | +# Load the file |
| 30 | +# Main MAT File - userdata |
| 31 | +main_file = sio.loadmat(file_path) |
| 32 | + |
| 33 | +# Parsing TriRep Anatomy Data |
| 34 | +# .X & .Triangulation |
| 35 | + |
| 36 | +# Vertices/Points = .X |
| 37 | +data_tri_X = main_file['userdata']['surface_triRep_X'][0][0] |
| 38 | +x = data_tri_X[:,0] |
| 39 | +y = data_tri_X[:,1] |
| 40 | +z = data_tri_X[:,2] |
| 41 | + |
| 42 | +# Triangulation/faces - .T |
| 43 | +data_tri_T = main_file['userdata']['surface_triRep_Triangulation'][0][0] |
| 44 | +# Resolving indexing issue - matlab to python |
| 45 | +t = data_tri_T-1 |
| 46 | + |
| 47 | +data_act_bip = main_file['userdata']['surface'][0][0]['act_bip'][0][0] |
| 48 | +voltage_data = data_act_bip[:,1] |
| 49 | + |
| 50 | +# Checking for NaN values and replacing with 0 in voltage data |
| 51 | +voltage_new = np.asarray(list(map(lambda x:0 if np.isnan(x) else x, voltage_data))) |
| 52 | + |
| 53 | + |
| 54 | + |
| 55 | +# ColorShell |
| 56 | +# coloring the shell by the voltage value rather than z-axis height |
| 57 | +min_volt = min(voltage_new) |
| 58 | +max_volt = max(voltage_new) |
| 59 | + |
| 60 | +# normalisation - normalising the voltage - values |
| 61 | +norm = mp.colors.Normalize(vmin=min_volt, vmax=max_volt) |
| 62 | + |
| 63 | +# Plot the Mesh |
| 64 | +# Create a new fig window |
| 65 | +fig = plt.figure() |
| 66 | +# get the current 3d axis on the current fig |
| 67 | +ax1 = fig.gca(projection='3d') |
| 68 | + |
| 69 | + |
| 70 | +# defining a custom color scalar field |
| 71 | +# if only nodal values are known, then the signal need to be averaged by the triangles |
| 72 | +color = np.mean(voltage_new[t], axis=1) |
| 73 | + |
| 74 | +# Voltage Threshold - to split the colorbar axis |
| 75 | +# determining the size of the voltage values |
| 76 | +diff_volt = max_volt - min_volt |
| 77 | + |
| 78 | +dist_below_threshold = volt_threshold - min_volt |
| 79 | +dist_above_threshold = max_volt - volt_threshold |
| 80 | + |
| 81 | +# splitting the colorbar |
| 82 | +# color below threshold - size |
| 83 | +size_below_threshold = round((dist_below_threshold/diff_volt)*256) |
| 84 | +# color above threshold - size |
| 85 | +size_above_threshold = 256 - size_below_threshold |
| 86 | + |
| 87 | +# Define Colormap - Reverse JET |
| 88 | +cm_jet = plt.cm.get_cmap('jet_r',size_below_threshold) |
| 89 | +newcolors = cm_jet(np.linspace(0,1,256)) |
| 90 | +# Magenta - RGBA array |
| 91 | +magenta = np.array([255/255, 0/255, 255/255, 1]) |
| 92 | +# Gray - RGBA array |
| 93 | +gray = np.array([128/255, 128/255, 128/255, 1]) |
| 94 | + |
| 95 | +# Creating new colormap columns |
| 96 | +# Col 1 - Jet |
| 97 | +col1 = cm_jet(np.linspace(0,1,42)) |
| 98 | +# Col2 - Magenta |
| 99 | +col2 = np.zeros((214,4)) + magenta |
| 100 | +# Combining Col1 + Col2 |
| 101 | +final_col = np.concatenate((col1,col2)) |
| 102 | +# NewColormap |
| 103 | +newcmp = ListedColormap(final_col) |
| 104 | + |
| 105 | +# Points - vertices |
| 106 | +points = data_tri_X |
| 107 | +print('x-y-z points\n',points.shape) |
| 108 | + |
| 109 | + |
| 110 | +# Trimesh object |
| 111 | +mesh = tm.Trimesh(vertices=points,faces=t) |
| 112 | + |
| 113 | +# All edges |
| 114 | +print('total number of edges:\n',mesh.edges.shape) |
| 115 | + |
| 116 | + |
| 117 | + |
| 118 | +# Unique edges |
| 119 | +unique_edges = mesh.edges_unique |
| 120 | +print('unique edges of the mesh\n',unique_edges.shape) |
| 121 | + |
| 122 | +# Identify the unshared edges |
| 123 | +unshared_edges = mesh.face_adjacency_unshared |
| 124 | +print('unshared\n',unshared_edges.shape) |
| 125 | + |
| 126 | +unique_edges_in_unshared = [] |
| 127 | +for item in unshared_edges: |
| 128 | + if item in unique_edges: |
| 129 | + unique_edges_in_unshared.append(item) |
| 130 | + else: |
| 131 | + pass |
| 132 | +unique_edges_in_unshared = np.array(unique_edges_in_unshared) |
| 133 | +print('unique edges in unshared:\n',unique_edges_in_unshared) |
| 134 | + |
| 135 | +# Identify the unique unshared edges |
| 136 | +uniques = {} |
| 137 | +count = 0 |
| 138 | +for a,b in unique_edges_in_unshared: |
| 139 | + if (a,b) not in uniques and (b,a) not in uniques: |
| 140 | + uniques[(a,b)]=(a,b) |
| 141 | + else: |
| 142 | + if (a,b) in uniques: |
| 143 | + del uniques[(a,b)] |
| 144 | + else: |
| 145 | + del uniques[(b,a)] |
| 146 | +# print('uniques\n',uniques) |
| 147 | + |
| 148 | +unique_unshared_edges = np.array(list(uniques.values())) |
| 149 | +print('unique_unshared_edges:\n',unique_unshared_edges.shape) |
| 150 | + |
| 151 | +# Get the vertex points for unique unshared edges |
| 152 | + |
| 153 | + |
| 154 | +# extract the ponts/vertices of the unique unshared edges |
| 155 | +for item in unique_unshared_edges: |
| 156 | + for i in item: |
| 157 | + x_values.append(points[i][0]) |
| 158 | + y_values.append(points[i][1]) |
| 159 | + z_values.append(points[i][2]) |
| 160 | + # print('unshared_edge_coordinates(x);\n',x_values.append(points[i][0])) |
| 161 | + # print('unshared_edge_coordinates(y);\n',y_values.append(points[i][1])) |
| 162 | + # print('unshared_edge_coordinates(z);\n',z_values.append(points[i][2])) |
| 163 | + |
| 164 | +# Get the unique edge points |
| 165 | +x_values_array = np.array(x_values).reshape(len(x_values),1) |
| 166 | +y_values_array = np.array(y_values).reshape(len(y_values),1) |
| 167 | +z_values_array = np.array(z_values).reshape(len(z_values),1) |
| 168 | + |
| 169 | + |
| 170 | +xx = np.array(x_values).reshape(1,len(x_values)) |
| 171 | +yy = np.array(y_values).reshape(1,len(y_values)) |
| 172 | +zz = np.array(z_values).reshape(1,len(z_values)) |
| 173 | +print('xx-values\n',xx) |
| 174 | +print('yy-values\n',yy.shape) |
| 175 | +print('zz-values\n',zz.shape) |
| 176 | + |
| 177 | + |
| 178 | +print('x_values\n',x_values_array.shape) |
| 179 | +print('y_values\n',y_values_array.shape) |
| 180 | +print('z_values\n',z_values_array.shape) |
| 181 | + |
| 182 | +print('x-values\n',x) |
| 183 | +print('y-values\n',y.shape) |
| 184 | +print('z-values\n',z.shape) |
| 185 | + |
| 186 | +# # Join the unique edges together |
| 187 | +unique_unshared_edge_points = np.concatenate((x_values_array, |
| 188 | + y_values_array, |
| 189 | + z_values_array), |
| 190 | + axis=1) |
| 191 | +print('border_points:\n',unique_unshared_edge_points.shape) |
| 192 | + |
| 193 | +# tri = Triangulation(x_values_array,y_values_array,z_values_array) |
| 194 | + |
| 195 | +surf = ax1.plot_trisurf(xx[0], |
| 196 | + yy[0], |
| 197 | + zz[0], |
| 198 | + linewidth=0.1) |
| 199 | + |
| 200 | +# mesh1 = tm.Trimesh(vertices=unique_unshared_edge_points,faces=t) |
| 201 | +# # |
| 202 | +# mesh1.show() |
| 203 | + |
| 204 | +# Plotting the border to the axis |
| 205 | +# ax1.plot_trisurf(xx[0], |
| 206 | +# yy[0], |
| 207 | +# zz[0], |
| 208 | +# triangles=t, |
| 209 | +# linewidth=0.5) |
| 210 | + |
| 211 | +surf1 = ax1.plot_trisurf(x, |
| 212 | + y, |
| 213 | + z, |
| 214 | + triangles=t, |
| 215 | + linewidth=0.5, |
| 216 | + antialiased=True, |
| 217 | + cmap=newcmp) |
| 218 | + |
| 219 | + |
| 220 | +# print('dir:\n',dir(surf1)) |
| 221 | +# print('typeofsurf',type(surf1)) |
| 222 | +# surf1.set_edgecolor([0,0,0,1]) |
| 223 | +surf1.set_array(color) |
| 224 | +# surf1.set_markeredgecolor([0,0,0,1]) |
| 225 | +# surf1.set_fc([0.0,1,0]) |
| 226 | +ax1.set_title('OpenEP TriRep Anatomy Data') |
| 227 | +plt.axis('off') |
| 228 | + |
| 229 | +# Colour Pallette, Position Left |
| 230 | +cb = plt.colorbar(mp.cm.ScalarMappable(norm=norm, cmap=newcmp), |
| 231 | + ax=[ax1], |
| 232 | + location='left', |
| 233 | + label='Voltage (mV)') |
| 234 | + |
| 235 | +# Show plot |
| 236 | +plt.show() |
| 237 | + |
| 238 | + |
| 239 | + |
| 240 | +# DRAWMAP plots an OpenEP map |
| 241 | +# userdata is a Carto data structure |
| 242 | +def draw_map( |
| 243 | + userdata, *arg): |
| 244 | + pass |
0 commit comments