From 06c6d4e25b49aefc3c8ef3dfe6b139481c5671ba Mon Sep 17 00:00:00 2001 From: gharaviri Date: Tue, 3 Oct 2023 01:51:35 +0100 Subject: [PATCH 01/31] Add CV filed to the userdata --- .DS_Store | Bin 0 -> 6148 bytes openep/Analyses/CV.py | 568 ++++++++++++++++++++++++++++++ openep/Analyses/run.py | 70 ++++ openep/data_structures/surface.py | 8 + openep/io/writers.py | 2 + 5 files changed, 648 insertions(+) create mode 100644 .DS_Store create mode 100644 openep/Analyses/CV.py create mode 100644 openep/Analyses/run.py diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..0036caa8f2050a2c94eb804ca92bbd022344a24f GIT binary patch literal 6148 zcmeHKyH3ME5S)bwMWRSb%KHJ3_ya3SYKr^-f`kak@5lzWAWJ+|oR z1O8y)V}0n_+U@)>%9F{R-0Zs zW=sWAfmGmI0ewFdx?&U9JKC#*!A1b$glRKg+blsWnjkiTy(2R;aVpWN5+jB Triangulation process') + surf = egmX_point_cloud.delaunay_3d(alpha=alpha) + print('---> Triangulation process') + print('---> CV calculation (Triangulation method started ...)') + for num_point in range(len(surf.cells_dict[5])): + + vtx_id = [] + lat = [] + + for i in range (3): + vtx_id.append(surf.cells_dict[5][num_point][i]) + lat.append(LAT[surf.cells_dict[5][num_point][i]]) + + id_lat_sorted = np.argsort(lat) + + O = [egmX[int(vtx_id[id_lat_sorted[0]])],lat[id_lat_sorted[0]]] + A = [egmX[int(vtx_id[id_lat_sorted[1]])],lat[id_lat_sorted[1]]] + B = [egmX[int(vtx_id[id_lat_sorted[2]])],lat[id_lat_sorted[2]]] + + OA = np.sqrt(sum(np.power(np.subtract(O[0],A[0]),2))) + OB = np.sqrt(sum(np.power(np.subtract(O[0],B[0]),2))) + AB = np.sqrt(sum(np.power(np.subtract(A[0],B[0]),2))) + + tOA = A[1] - O[1] + tOB = B[1] - O[1] + + theta = np.arccos((np.power(OA,2)+ np.power(OB,2) - np.power(AB,2))/(2 * OA * OB)) + # check if the conditions are meet to accept the triangle set as a viable acceptable one + if (math.degrees(theta) >= min_theta and OA >= min_elec_distance and OA <= max_elec_distance and + OB >= min_elec_distance and OB <= max_elec_distance and tOA >= min_lat_difference and tOB >= min_lat_difference): + + alpha = np.arctan((tOB * OA - tOA * OB * np.cos(theta)) / (tOA * OB * np.sin(theta))) + cv_temp = (OA/tOA) * np.cos(alpha) + cv.append(cv_temp) + cv_centers_x.append(O[0][0]) + cv_centers_y.append(O[0][1]) + cv_centers_z.append(O[0][2]) + cv_centers_id.append(vtx_id[id_lat_sorted[0]]) + #if cv_temp >= 0.2 and cv_temp <=2: + # cv.append(cv_temp) + # cv_centers_x.append(O[0][0]) + # cv_centers_y.append(O[0][1]) + # cv_centers_z.append(O[0][2]) + # cv_centers_id.append(vtx_id[id_lat_sorted[0]]) + #else: + # cv_temp = (OA/tOA) * np.cos(alpha) + # cv_centers_x.append(O[0][0]) + # cv_centers_y.append(O[0][1]) + # cv_centers_z.append(O[0][2]) + # cv_centers_id.append(vtx_id[id_lat_sorted[0]]) + # cv.append(cv_temp) + ####### Create a triangulation mesh from egmX pointcloud ##### + cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] + print('---> CV calculation ended') + return cv, cv_centers, cv_centers_id + + def plane_fitting(self,egmX, LAT): + + cv = [] + cv_centers_x = [] + cv_centers_y = [] + cv_centers_z = [] + cv_centers_id = [] + cent = np.random.random((len(egmX), 3)) + direction = np.random.random((len(egmX), 3)) + + tree = KDTree(egmX, leaf_size=5) + all_nn_indices = tree.query_radius(egmX, r=10) + all_nns = [[egmX[idx] for idx in nn_indices] for nn_indices in all_nn_indices] + + print(len(egmX)) + for k in range(len(egmX)): + print(k,end='\r') + if len(all_nns[k]) >= 5: + cv_data = np.zeros(shape= (len(all_nns[k]) + 1,4)) + cv_data[0][:3] = egmX[k] + cv_data[0][3] = LAT[k] + for i in range(len(all_nns[k])): + cv_data[i+1][:3] = all_nns[k][i] + cv_data[i+1][3] = LAT[all_nn_indices[k][i]] + + x1 = cv_data[:,0] + x2 = cv_data[:,1] + x3 = cv_data[:,2] + m = cv_data[:,3] + #mesh_original = pv.PolyData(Points,Face) + #p2.add_points(cv_data[:,:3], render_points_as_spheres=True, point_size=5.0, color = 'yellow') + #p2.add_mesh(mesh_original,opacity=0.4) + #p2.show() + + res = minimize(func, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], args=(x1, x2, x3, m)) + a, b, c, d,e, f, g, h, q, l= res.x + #print(a*x1**2+b*x2**2+c*x3**2+d*x1*x2+e*x1*x3+f*x2*x3+g*x1+h*x2+q*x3+l -m) + + + x = egmX[k,0] + y = egmX[k,1] + z = egmX[k,2] + + tx = 2*a*x + d*y + e*z + g + ty = 2*b*y + d*x + f*z + h + tz = 2*c*z + e*x + f*y + q + + vx = tx/(tx**2 + ty**2 + tz**2) + vy = ty/(tx**2 + ty**2 + tz**2) + vz = tz/(tx**2 + ty**2 + tz**2) + + + cent[k,0] = x + cent[k,1] = y + cent[k,2] = z + + direction[k,0] = vx + direction[k,1] = vy + direction[k,2] = vz + + cv_temp = np.sqrt(np.sum(np.power(vx,2)+np.power(vy,2)+np.power(vz,2))) + direction[k,0] = vx + direction[k,1] = vy + direction[k,2] = vz + + if cv_temp >=0.2 and cv_temp <=2: + cv.append(cv_temp) + cv_centers_x.append(x) + cv_centers_y.append(y) + cv_centers_z.append(z) + cv_centers_id.append(k) + else: + + cv.append(cv_temp) + cv_centers_x.append(x) + cv_centers_y.append(y) + cv_centers_z.append(z) + else: + x = egmX[k,0] + y = egmX[k,1] + z = egmX[k,2] + cv.append(np.nan) + cv_centers_x.append(x) + cv_centers_y.append(y) + cv_centers_z.append(z) + cv_centers_id.append(k) + cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] + + return cv, cv_centers, cv_centers_id + + def RBF_method(self, egmX, LAT, points, face,Fiber_data='',): + + #rbfi = Rbf(egmX[:,0],egmX[:,1],egmX[:,2],LAT,kernel='gaussian',epsilon=0.5) + rbfi = Rbf(egmX[:,0],egmX[:,1],egmX[:,2],LAT) + lat_i = [] + for i in range(len(points)): + lat_i.append(rbfi(points[i][0], points[i][1], points[i][2])) + + mesh_original_RBF = pv.PolyData(points,face) + mesh_original_RBF['values'] = lat_i + deriv = mesh_original_RBF.compute_derivative('values') + df = deriv['gradient'] + + cvx = deriv['gradient'][:,0]/np.sum(np.power( deriv['gradient'],2),axis=1) + cvy = deriv['gradient'][:,1]/np.sum(np.power( deriv['gradient'],2),axis=1) + cvz = deriv['gradient'][:,2]/np.sum(np.power( deriv['gradient'],2),axis=1) + cv = np.sqrt(np.power(cvx,2) + np.power(cvy,2) + np.power(cvz,2)) + + cv_data = [] + cv_x = [] + cv_y = [] + cv_z = [] + cv_center_id =[] + cv_centers_x = [] + cv_centers_y = [] + cv_centers_z = [] + tree = KDTree(points, leaf_size=5) + dist, ind = tree.query(egmX, k=1) + for k in range(len(ind)): + + if cv[ind[k][0]] > 0 and cv[ind[k][0]] <= 10: + cv_data.append(cv[ind[k][0]]) + cv_x.append(cvx[ind[k][0]]) + cv_y.append(cvy[ind[k][0]]) + cv_z.append(cvz[ind[k][0]]) + cv_center_id.append(ind[k][0]) + #cv_centers_x.append(points[ind[k][0]][0]) + #cv_centers_y.append(points[ind[k][0]][1]) + #cv_centers_z.append(points[ind[k][0]][2]) + cv_centers_x.append(egmX[k][0]) + cv_centers_y.append(egmX[k][1]) + cv_centers_z.append(egmX[k][2]) + + + else: + cv_data.append(np.nan) + cv_centers_x.append(points[ind[k][0]][0]) + cv_centers_y.append(points[ind[k][0]][1]) + cv_centers_z.append(points[ind[k][0]][2]) + if 0: + cv_data.append(cv[ind[k][0]]) + cv_center_id.append(ind[k][0]) + cv_centers_x.append(points[ind[k][0]][0]) + cv_centers_y.append(points[ind[k][0]][1]) + cv_centers_z.append(points[ind[k][0]][2]) + cv_x.append(cvx[ind[k][0]]) + cv_y.append(cvy[ind[k][0]]) + cv_z.append(cvz[ind[k][0]]) + + direction = [cv_x,cv_y,cv_z] + cv_centers = [cv_centers_x,cv_centers_y,cv_centers_z] + + + + + return cv_data, cv_centers, direction, cv_center_id + + def visualization(self,lat_data = '', + cv_data = '', + mesh = '', + egmX = '', + method = '', + file_name = '', + direction='', + direction_centers = '', + show_histogram= False, + show_cv_fibrosis = True, + interpolation_radius = 5): + + boring_cmap = plt.cm.get_cmap("jet", 256) + boring_cmap_r = plt.cm.get_cmap("jet_r", 256) + boring_cmap_2 = plt.cm.get_cmap("Paired", 4) + + boring_cmap_2.colors[0][0]= 34/255 #,237,222#1 + boring_cmap_2.colors[0][1]=94/255 + boring_cmap_2.colors[0][2]=168/255 + boring_cmap_2.colors[1][0]= 62/255 #,190,133#0 + boring_cmap_2.colors[1][1]=182/255 #1 + boring_cmap_2.colors[1][2]=196/255 #1 + boring_cmap_2.colors[2][0]=253/255 + boring_cmap_2.colors[2][1]=141/255 + boring_cmap_2.colors[2][2]= 60/255 + boring_cmap_2.colors[3][0]=217/255 + boring_cmap_2.colors[3][1]=71/255 + boring_cmap_2.colors[3][2]=1/255 + + + + + + pv.rcParams['transparent_background'] = True + + + + if method == 'LAT': + #pl = pv.Plotter(off_screen=True) + pl = pv.Plotter() + egmX_point_cloud = pv.PolyData(egmX) + egmX_point_cloud['LAT[ms]'] = lat_data + interpolated_mesh_lat = mesh.interpolate(egmX_point_cloud, radius=interpolation_radius) + pl.add_mesh(egmX_point_cloud,scalars="LAT[ms]",render_points_as_spheres=True, + point_size=20, cmap = boring_cmap_r) + pl.add_mesh(interpolated_mesh_lat,scalars="LAT[ms]", clim = [np.min(lat_data), np.max(lat_data)], + below_color='brown', above_color = 'magenta', cmap = boring_cmap_r) + if direction != '' and direction_centers != '': + pl.add_arrows(direction_centers, direction, mag=2, color=[90, 90, 90]) + #here + pl.show(screenshot=f'/Users/aligharaviri/Desktop/CV_Data/Figures_picsLAT_{file_name}.png') + if method == 'CV': + pl = pv.Plotter(shape=(1,2),off_screen=False) + sargs = dict( + title_font_size=32, + label_font_size=32, + shadow=True, + n_labels=11, + italic=True, + fmt="%.1f", + font_family="arial", + color='black' + ) + egmX_point_cloud = pv.PolyData(egmX) + egmX_point_cloud['CV[m/s]'] = cv_data + ######## Ver 1 + #interpolated_mesh_cv = mesh.interpolate(egmX_point_cloud,radius=5 ,n_points=2 ,sharpness=2) + ######## Ver 1 + interpolated_mesh_cv = mesh.interpolate(egmX_point_cloud,radius=10 ,n_points=None,null_value=np.nan) + + pl.subplot(0,0) + + pl.add_mesh(egmX_point_cloud,scalars='CV[m/s]',render_points_as_spheres=True, + point_size=20, cmap = boring_cmap,clim=[0,1.5],below_color='navy', above_color = 'darkred') + pl.add_mesh(interpolated_mesh_cv,scalars='CV[m/s]', clim=[0,1.5], + below_color='navy', above_color = 'darkred', cmap = boring_cmap,scalar_bar_args=sargs) + + if direction != '': + pl.add_arrows(egmX, np.transpose(direction), mag=4, color='black')#color=[90, 90, 90]) + pl.camera.roll += 20 + pl.camera.elevation -=10 + pl.camera.azimuth = 10 + pl.camera.zoom (1.2) + pl.subplot(0,1) + + pl.add_mesh(egmX_point_cloud,scalars='CV[m/s]',render_points_as_spheres=True, + point_size=20, cmap = boring_cmap,clim=[0,1.2],below_color='navy', above_color = 'darkred') + pl.add_mesh(interpolated_mesh_cv,scalars='CV[m/s]', clim=[0,1.5], + below_color='navy', above_color = 'darkred', cmap = boring_cmap,scalar_bar_args=sargs) + if direction != '': + pl.add_arrows(egmX, np.transpose(direction), mag=4, color='black')#color=[90, 90, 90]) + #here + pl.camera.azimuth += 190 + pl.camera.roll -=30 + pl.camera.elevation -=40 + pl.camera.zoom (1.2) + #pl.show(screenshot=f'/Users/aligharaviri/Desktop/CV_Data/Figures_pics/{file_name}.png') + pl.show() + #pl.close() + + if method== 'CV' and show_histogram: + plt.rcParams['font.size'] = '14' + plt.rcParams['font.family'] = 'arial' + bin_width = 0.01 + plt.hist(cv_data, bins=np.arange(np.nanmin(cv_data), np.nanmax(cv_data) + bin_width, bin_width)) + + plt.xlabel('CV (m/s)') + plt.ylabel('Number of Events') + + #plt.savefig(f'/Volumes/Simulations/CV_Data/Figures_pics/hist_{file_name}.pdf',bbox_inches='tight',transparent=True, dpi=400) + #here + plt.show() + #plt.close() + if method == 'CV' and show_cv_fibrosis: + + interpolated_mesh_cv_face = interpolated_mesh_cv.point_data_to_cell_data() + correlation_data = [] + + low_cv_mask = interpolated_mesh_cv_face['CV[m/s]'] <= 0.3 + high_cv_mask = interpolated_mesh_cv_face['CV[m/s]'] > 0.3 + fibrosis_mask = mesh['regions'] == 34 + normal_mask = mesh['regions'] != 34 + low_cv_fibrosis = np.multiply(low_cv_mask,fibrosis_mask) + high_cv_no_fibrosis = np.multiply(high_cv_mask,normal_mask) + low_cv_no_fibrosis = np.multiply(low_cv_mask,normal_mask) + high_cv_fibrosis = np.multiply(high_cv_mask,fibrosis_mask) + correlation_data = (np.multiply(low_cv_fibrosis,3) + np.multiply(high_cv_fibrosis,2) + np.multiply(low_cv_no_fibrosis,1) + np.multiply(high_cv_no_fibrosis,0)) + accuracy = np.sum(np.multiply(low_cv_fibrosis,1) + np.multiply(high_cv_no_fibrosis,1)) / len(interpolated_mesh_cv_face['CV[m/s]']) + print('Accuracy= ',accuracy) + mesh['correlation'] = correlation_data + pl = pv.Plotter(off_screen=True) + pl.add_mesh(mesh,scalars='correlation', cmap = boring_cmap_2) + pl.show(screenshot=f'/Users/aligharaviri/Desktop/CV_Data/Figures_pics/correlation_{file_name}.png') + pl.close() + cv_not_nan = interpolated_mesh_cv_face['CV[m/s]'][~np.isnan(interpolated_mesh_cv_face['CV[m/s]'])] + region_not_nan = fibrosis_mask[~np.isnan(interpolated_mesh_cv_face['CV[m/s]'])] + reg_id = np.multiply(region_not_nan,1) + elif method == 'CV' and not show_cv_fibrosis: + reg_id = [] + cv_not_nan = [] + #cv_corr = [] + #for i in range(len(mesh['regions'])): + # print(i,'end:\r') + # if mesh['regions'][i] == 34: + # if ~np.isnan(interpolated_mesh_cv_face['CV[m/s]'][i]): + # reg_id.append(1) + # cv_corr.append(interpolated_mesh_cv_face['CV[m/s]'][i]) + # else: + # if ~np.isnan(interpolated_mesh_cv_face['CV[m/s]'][i]): + # reg_id.append(0) + # cv_corr.append(interpolated_mesh_cv_face['CV[m/s]'][i]) + + + if method == 'CV': + return interpolated_mesh_cv['CV[m/s]'], reg_id, cv_not_nan + + def gradient_method(self, egmX, LAT, points, face, Fiber_data = ''): + '''#sampling_points = Points[sampling_points_ID]''' + + mesh_original = pv.PolyData(points,face) + mesh_original['LAT[ms]'] = LAT + deriv = mesh_original.compute_derivative('LAT[ms]') + df = deriv['gradient'] + d_lat_x = deriv['gradient'][:,0]#/np.sum(np.power( deriv['gradient'],2),axis=1) + d_lat_y = deriv['gradient'][:,1]#/np.sum(np.power( deriv['gradient'],2),axis=1) + d_lat_z = deriv['gradient'][:,2]#/np.sum(np.power( deriv['gradient'],2),axis=1) + #cv_data = 1 / (np.power(d_lat_x,2) + np.power(d_lat_y,2) + np.power(d_lat_z,2)) + cvx = d_lat_x/np.sum(np.power( deriv['gradient'],2),axis=1) + cvy = d_lat_y/np.sum(np.power( deriv['gradient'],2),axis=1) + cvz = d_lat_z/np.sum(np.power( deriv['gradient'],2),axis=1) + cv_data = np.sqrt(np.power(cvx,2) + np.power(cvy,2) + np.power(cvz,2))#/np.sum(np.power( deriv['gradient'],2),axis=1) + cv_data_sampled = cv_data[egmX] + + return cv_data, cv_data_sampled + + def calculate_divergence (self, egmX, LAT, points, face, sampling_data,collision_threshold=-1,focal_threshold=1,interpolation_radius=5): + boring_cmap = plt.cm.get_cmap("jet", 256) + boring_cmap_r = plt.cm.get_cmap("jet_r", 256) + egmX_point_cloud = pv.PolyData(egmX) + egmX_point_cloud['LAT[ms]'] = LAT + mesh_original = pv.PolyData(points,face) + interpolated = mesh_original.interpolate(egmX_point_cloud, radius=interpolation_radius,n_points=None) + pl = pv.Plotter() + pl.add_mesh(interpolated, cmap=boring_cmap_r) + #pl.show() + pl.close() + deriv = interpolated.compute_derivative('LAT[ms]') + + cv_x = deriv['gradient'][:,0]/np.sum(np.power( deriv['gradient'],2),axis=1) + cv_y = deriv['gradient'][:,1]/np.sum(np.power( deriv['gradient'],2),axis=1) + cv_z = deriv['gradient'][:,2]/np.sum(np.power( deriv['gradient'],2),axis=1) + + direction = np.ndarray(shape=(len(cv_x),3)) + cv_direction = np.ndarray(shape=(len(cv_x),3)) + for i in range(len(cv_x)): + direction[i][0]=cv_x[i] /np.sqrt(np.sum(np.power(cv_x[i],2) + np.power(cv_y[i],2) + np.power(cv_z[i],2))) + direction[i][1]=cv_y[i] /np.sqrt(np.sum(np.power(cv_x[i],2) + np.power(cv_y[i],2) + np.power(cv_z[i],2))) + direction[i][2]=cv_z[i] /np.sqrt(np.sum(np.power(cv_x[i],2) + np.power(cv_y[i],2) + np.power(cv_z[i],2))) + cv_direction[i][0]=cv_x[i] + cv_direction[i][1]=cv_y[i] + cv_direction[i][2]=cv_z[i] + + mesh_original['activation_direction'] = cv_direction + div = mesh_original.compute_derivative(scalars='activation_direction',divergence=True) + sargs = dict( + title_font_size=32, + label_font_size=32, + shadow=True, + n_labels=11, + italic=True, + fmt="%.1f", + font_family="arial", + color='black' + ) + pv.rcParams['transparent_background'] = True + pl = pv.Plotter(shape=(1,2)) + pl.subplot(0,0) + pl.add_mesh(div.cell_data_to_point_data(), scalars='divergence',cmap=plt.cm.get_cmap("PRGn", 9),clim=[-1.5,1.5],below_color=[63, 1, 44],above_color=[1, 50,32], + scalar_bar_args=sargs) + pl.add_arrows(points[sampling_data],direction[sampling_data],mag=2,color=[90, 90, 90]) + pl.set_background('white') + pl.camera.roll += 20 + pl.camera.elevation -=10 + pl.camera.azimuth = 10 + pl.camera.zoom (1.5) + + pl.subplot(0,1) + pl.add_mesh(div.cell_data_to_point_data(), scalars='divergence',cmap=plt.cm.get_cmap("PRGn", 9),clim=[-1.5,1.5],below_color=[63, 1, 44],above_color=[1, 50,32], + scalar_bar_args=sargs) + pl.add_arrows(points[sampling_data],direction[sampling_data],mag=2,color=[90, 90, 90]) + pl.set_background('white') + pl.camera.azimuth += 190 + pl.camera.roll -=30 + pl.camera.elevation -=40 + pl.camera.zoom (1.5) + #pl.show(screenshot='/Volumes/Simulations/CV_Data/Figures_pics/foo.png') + #pl.show() + pl.close() + + div_value = [] + for i in range(len(div['divergence'])): + if div['divergence'][i] < collision_threshold or div['divergence'][i] > focal_threshold : + div_value.append(1) + + else: + div_value.append(0) + mesh_divergence = pv.PolyData(points,face) + mesh_divergence ['div'] = div_value + divergence_per_point = mesh_divergence.cell_data_to_point_data() + pl = pv.Plotter() + pl.add_mesh(divergence_per_point,cmap=plt.cm.get_cmap("jet", 2)) + pl.add_arrows(points[sampling_data],direction[sampling_data],mag=2,color=[90, 90, 90]) + #here + #pl.show() + pl.close() + + + return direction, mesh_original, div_value + + def exclude_collision_points(self, conduction_velocity_centers = '', conduction_velocity= '', + divergence_value = '', upper_threshold = '',mesh_points = ''): + tree = KDTree(mesh_points, leaf_size=2) + ind = tree.query_radius(conduction_velocity_centers, r=0.5) + non_colision_points = [] + collision_detection = [] + cv_data = [] + for i in range(len(ind)): + print(i,end='\r') + for k in range(len(ind[i])): + if divergence_value[ind[i][k]] == 1: + collision_detection.append(1) + else: + collision_detection.append(0) + + if np.sum(collision_detection) or (conduction_velocity[i] > upper_threshold): + non_colision_points.append(False) + cv_data.append(np.nan) + else: + non_colision_points.append(True) + cv_data.append(conduction_velocity[i]) + collision_detection = [] + + #bin_width = 0.01 + #plt.hist(cv_data, bins=np.arange(np.nanmin(cv_data), np.nanmax(cv_data) + bin_width, bin_width)) + + #plt.xlabel('CV (m/s)') + #plt.ylabel('Number of Events') + #here + #plt.show() + return non_colision_points, cv_data + #pl = pv.Plotter() + #pl.add_mesh(interpolated,cmap = boring_cmap_r) + #pl.add_arrows(points,direction,mag=2,color=[90, 90, 90]) + #pl.show() + + + + + #mesh_original['divergence'] = div_value + #pl.add_mesh(mesh_original, scalars='divergence',cmap=plt.cm.get_cmap("jet", 2)) + + diff --git a/openep/Analyses/run.py b/openep/Analyses/run.py new file mode 100644 index 00000000..db795855 --- /dev/null +++ b/openep/Analyses/run.py @@ -0,0 +1,70 @@ +import openep +from CV import CV +import pyvista as pv +import numpy as np +def run_cv (openep_file_name,cv_method, visualsiation): + case = openep.load_openep_mat(openep_file_name) + egmX = case.electric.bipolar_egm.points + voltage = case.electric.bipolar_egm.voltage + print('---> read electrode data') + LAT = case.electric.annotations.local_activation_time - case.electric.annotations.reference_activation_time + print('---> read bi-polar LAT data') + surface = case.create_mesh() + + valid_id = np.where(LAT!=-10000) + invalid_id = np.where(LAT==-10000) + egmX_clean = np.delete(egmX, np.where(LAT == -10000),0) + LAT_clean = np.delete(LAT, np.where(LAT == -10000)) + Points = surface.points + faces = surface.faces + Face = faces.reshape(-1,4) + temp_face = Face[:,1:4] + mesh_original = pv.PolyData(Points,Face) + cal_cv = CV() + if cv_method == 'Plane': + cv_plane, cv_centers_plane, cv_centers_id_plane = cal_cv.plane_fitting(egmX_clean, LAT_clean) + cv_total = np.zeros((len(egmX),1)) + for i in range(len(cv_plane)): + cv_total[valid_id[0][i]] = cv_plane[i] + + for i in range(len(invalid_id)): + cv_total[invalid_id[i]] = np.nan + + if visulaisation: + cal_cv.visualization(lat_data = '', cv_data = cv_plane, mesh = mesh_original, egmX = egmX_clean, method = 'CV',file_name = '', + show_histogram= True, direction_centers = Points,show_cv_fibrosis=False) + cv = cv_plane + + if cv_method == 'RBF': + cv_rbf, cv_centers_rbf, direction_rbf, cv_centers_id_rbf = cal_cv.RBF_method(egmX_clean, LAT_clean,Points, Face) + if visulaisation: + cal_cv.visualization(lat_data = '', cv_data = cv_rbf, mesh = mesh_original, egmX = np.transpose(cv_centers_rbf), method = 'CV',file_name = '', + show_histogram= True, direction_centers = Points,show_cv_fibrosis=False) + cv = cv_rbf + + if cv_method == 'Tri': + cv_tri, cv_centers_tri, cv_centers_id_tri = cal_cv.Triangulation(egmX_clean,LAT_clean) + if visulaisation: + cal_cv.visualization(lat_data = '', cv_data = cv_tri, mesh = mesh_original, egmX = np.transpose(cv_centers_tri), method = 'CV',file_name = '', + show_histogram= True, direction_centers = Points,show_cv_fibrosis=False) + cv = cv_tri + ######## Save the CV data as openep ####### + case.fields['conduction_velocity'] = cv + return case + + + + + + + +if __name__=="__main__": + file_path = '/Users/aligharaviri/Downloads/LAWT_CartoAF_ EP_Studies/' + file_name = '517_cv_cv.mat' + openep_file_name = f"{file_path}{file_name}" + out_put_file_name = f"{file_name[:-4]}_cv.mat" + method = "RBF" ######## methods are Plane, RBF, and Tri + visulaisation = True + case = run_cv (openep_file_name,method, visulaisation) + openep.export_openep_mat(case,out_put_file_name) + diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index 543e4a35..4fb46a17 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -51,6 +51,7 @@ class Fields: longitudinal_fibres: np.ndarray = None transverse_fibres: np.ndarray = None pacing_site: np.ndarray = None + conduction_velocity: np.ndarray = None def __repr__(self): return f"fields: {tuple(self.__dict__.keys())}" @@ -190,6 +191,11 @@ def extract_surface_data(surface_data): if isinstance(pacing_site, np.ndarray) and pacing_site.size == 0: pacing_site = None + + try: + conduction_velocity = surface_data['conduction_velocity'].astype(float) + except KeyError as e: + conduction_velocity = None fields = Fields( bipolar_voltage=bipolar_voltage, @@ -202,6 +208,7 @@ def extract_surface_data(surface_data): longitudinal_fibres=longitudinal_fibres, transverse_fibres=transverse_fibres, pacing_site=pacing_site, + conduction_velocity = conduction_velocity, ) return points, indices, fields @@ -237,6 +244,7 @@ def empty_fields(n_points=0, n_cells=0): longitudinal_fibres, transverse_fibres, pacing_site, + conudction_vlocity, ) return fields diff --git a/openep/io/writers.py b/openep/io/writers.py index 06490d58..290583c2 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -300,6 +300,8 @@ def _extract_surface_data( surface_data['longitudinal'] = fields.longitudinal_fibres if fields.longitudinal_fibres is not None else empty_float_array surface_data['transverse'] = fields.transverse_fibres if fields.transverse_fibres is not None else empty_float_array surface_data['pacing_site'] = fields.pacing_site if fields.pacing_site is not None else empty_int_array + surface_data['conduction_velocity'] = fields.conduction_velocity if fields.conduction_velocity is not None else empty_int_array + # Remove arrays that are full of NaNs for field_name, field in surface_data.items(): From 256d8ab5e66eb56c4a3f25d095b2b0c322b8428e Mon Sep 17 00:00:00 2001 From: gharaviri Date: Tue, 3 Oct 2023 01:57:16 +0100 Subject: [PATCH 02/31] Minor change --- openep/Analyses/run.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openep/Analyses/run.py b/openep/Analyses/run.py index db795855..91800710 100644 --- a/openep/Analyses/run.py +++ b/openep/Analyses/run.py @@ -60,10 +60,10 @@ def run_cv (openep_file_name,cv_method, visualsiation): if __name__=="__main__": file_path = '/Users/aligharaviri/Downloads/LAWT_CartoAF_ EP_Studies/' - file_name = '517_cv_cv.mat' + file_name = '521.mat' openep_file_name = f"{file_path}{file_name}" out_put_file_name = f"{file_name[:-4]}_cv.mat" - method = "RBF" ######## methods are Plane, RBF, and Tri + method = "Plane" ######## methods are Plane, RBF, and Tri visulaisation = True case = run_cv (openep_file_name,method, visulaisation) openep.export_openep_mat(case,out_put_file_name) From 5bf02be36185cbcb507c014d2f46429dfe48ac15 Mon Sep 17 00:00:00 2001 From: gharaviri Date: Tue, 3 Oct 2023 12:59:29 +0100 Subject: [PATCH 03/31] Add normals to the userdata --- openep/Analyses/run.py | 4 +++- openep/data_structures/surface.py | 8 ++++++++ openep/io/writers.py | 2 +- 3 files changed, 12 insertions(+), 2 deletions(-) diff --git a/openep/Analyses/run.py b/openep/Analyses/run.py index 91800710..3e0f9a4f 100644 --- a/openep/Analyses/run.py +++ b/openep/Analyses/run.py @@ -50,6 +50,8 @@ def run_cv (openep_file_name,cv_method, visualsiation): cv = cv_tri ######## Save the CV data as openep ####### case.fields['conduction_velocity'] = cv + case.fields['mesh_normals'] = mesh_original.point_normals + return case @@ -60,7 +62,7 @@ def run_cv (openep_file_name,cv_method, visualsiation): if __name__=="__main__": file_path = '/Users/aligharaviri/Downloads/LAWT_CartoAF_ EP_Studies/' - file_name = '521.mat' + file_name = '521_cv_cv.mat' openep_file_name = f"{file_path}{file_name}" out_put_file_name = f"{file_name[:-4]}_cv.mat" method = "Plane" ######## methods are Plane, RBF, and Tri diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index 4fb46a17..3cd50afa 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -52,6 +52,7 @@ class Fields: transverse_fibres: np.ndarray = None pacing_site: np.ndarray = None conduction_velocity: np.ndarray = None + mesh_normals : np.ndarray = None def __repr__(self): return f"fields: {tuple(self.__dict__.keys())}" @@ -196,6 +197,11 @@ def extract_surface_data(surface_data): conduction_velocity = surface_data['conduction_velocity'].astype(float) except KeyError as e: conduction_velocity = None + + try: + mesh_normals = surface_data['mesh_normals'].astype(float) + except KeyError as e: + mesh_normals = None fields = Fields( bipolar_voltage=bipolar_voltage, @@ -209,6 +215,7 @@ def extract_surface_data(surface_data): transverse_fibres=transverse_fibres, pacing_site=pacing_site, conduction_velocity = conduction_velocity, + mesh_normals = mesh_normals ) return points, indices, fields @@ -245,6 +252,7 @@ def empty_fields(n_points=0, n_cells=0): transverse_fibres, pacing_site, conudction_vlocity, + mesh_normals, ) return fields diff --git a/openep/io/writers.py b/openep/io/writers.py index 290583c2..31852185 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -301,7 +301,7 @@ def _extract_surface_data( surface_data['transverse'] = fields.transverse_fibres if fields.transverse_fibres is not None else empty_float_array surface_data['pacing_site'] = fields.pacing_site if fields.pacing_site is not None else empty_int_array surface_data['conduction_velocity'] = fields.conduction_velocity if fields.conduction_velocity is not None else empty_int_array - + surface_data['mesh_normals'] = fields.mesh_normals if fields.mesh_normals is not None else empty_int_array # Remove arrays that are full of NaNs for field_name, field in surface_data.items(): From c810f355f8042689766062be84f57f44e3ecf3e4 Mon Sep 17 00:00:00 2001 From: gharaviri Date: Tue, 3 Oct 2023 16:42:25 +0100 Subject: [PATCH 04/31] Small changes in Tri method --- openep/Analyses/CV.py | 14 +++++++++----- openep/Analyses/run.py | 4 ++-- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/openep/Analyses/CV.py b/openep/Analyses/CV.py index b3f0d603..3ef3ad57 100644 --- a/openep/Analyses/CV.py +++ b/openep/Analyses/CV.py @@ -85,9 +85,9 @@ def Triangulation(self,egmX,LAT,data_type=''): alpha = np.arctan((tOB * OA - tOA * OB * np.cos(theta)) / (tOA * OB * np.sin(theta))) cv_temp = (OA/tOA) * np.cos(alpha) cv.append(cv_temp) - cv_centers_x.append(O[0][0]) - cv_centers_y.append(O[0][1]) - cv_centers_z.append(O[0][2]) + cv_centers_x.append((O[0][0] + A[0][0] + B[0][0])/3) + cv_centers_y.append((O[0][1] + A[0][1] + B[0][1])/3) + cv_centers_z.append((O[0][2] + A[0][2] + B[0][2])/3) cv_centers_id.append(vtx_id[id_lat_sorted[0]]) #if cv_temp >= 0.2 and cv_temp <=2: # cv.append(cv_temp) @@ -95,13 +95,17 @@ def Triangulation(self,egmX,LAT,data_type=''): # cv_centers_y.append(O[0][1]) # cv_centers_z.append(O[0][2]) # cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - #else: + else: # cv_temp = (OA/tOA) * np.cos(alpha) # cv_centers_x.append(O[0][0]) # cv_centers_y.append(O[0][1]) # cv_centers_z.append(O[0][2]) # cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - # cv.append(cv_temp) + cv.append(np.nan) + cv_centers_x.append((O[0][0] + A[0][0] + B[0][0])/3) + cv_centers_y.append((O[0][1] + A[0][1] + B[0][1])/3) + cv_centers_z.append((O[0][2] + A[0][2] + B[0][2])/3) + cv_centers_id.append(vtx_id[id_lat_sorted[0]]) ####### Create a triangulation mesh from egmX pointcloud ##### cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] print('---> CV calculation ended') diff --git a/openep/Analyses/run.py b/openep/Analyses/run.py index 3e0f9a4f..22953c23 100644 --- a/openep/Analyses/run.py +++ b/openep/Analyses/run.py @@ -62,10 +62,10 @@ def run_cv (openep_file_name,cv_method, visualsiation): if __name__=="__main__": file_path = '/Users/aligharaviri/Downloads/LAWT_CartoAF_ EP_Studies/' - file_name = '521_cv_cv.mat' + file_name = '521.mat' openep_file_name = f"{file_path}{file_name}" out_put_file_name = f"{file_name[:-4]}_cv.mat" - method = "Plane" ######## methods are Plane, RBF, and Tri + method = "Tri" ######## methods are Plane, RBF, and Tri visulaisation = True case = run_cv (openep_file_name,method, visulaisation) openep.export_openep_mat(case,out_put_file_name) From 1f0817c9c3a91d901d594dbc12ec1f7496e29b8b Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:42:00 +0000 Subject: [PATCH 05/31] analysis class with cv and div calculator --- openep/analysis/analyse.py | 8 + openep/analysis/conduction_velocity.py | 297 +++++++++++++++++++++++++ openep/analysis/divergence.py | 132 +++++++++++ 3 files changed, 437 insertions(+) create mode 100644 openep/analysis/analyse.py create mode 100644 openep/analysis/conduction_velocity.py create mode 100644 openep/analysis/divergence.py diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py new file mode 100644 index 00000000..16aa6989 --- /dev/null +++ b/openep/analysis/analyse.py @@ -0,0 +1,8 @@ +from .conduction_velocity import ConductionVelocity +from .divergence import Divergence + +class Analyse: + + def __init__(self, case): + self.conduction_velocity = ConductionVelocity(case) + self.divergence = Divergence(case) \ No newline at end of file diff --git a/openep/analysis/conduction_velocity.py b/openep/analysis/conduction_velocity.py new file mode 100644 index 00000000..f927bec5 --- /dev/null +++ b/openep/analysis/conduction_velocity.py @@ -0,0 +1,297 @@ +import openep +from openep.analysis.run import run_cv +import sys +import os +import math +import pyvista as pv +import numpy as np +import scipy.io +import scipy +from scipy.optimize import minimize +import matplotlib.pyplot as plt +from sklearn.neighbors import KDTree +import pdb +from scipy.interpolate import Rbf +from vedo import * +import pyvistaqt +from sklearn import metrics + +from ..case.case_routines import ( + interpolate_general_cloud_points_onto_surface +) + +class ConductionVelocity(): + + + def __init__(self, case): + self.case = case + self.prep_egmX, self.prep_LAT = self._preprocess_data() + self.include = None + self.values = None + self.points = None + + def calculate( + self, + method='triangulation', + method_kws=dict() + ): + + """Calculate conduction velocity and interpolates each point on a mesh, + stores in conduction_velocity fields. + + Args: + method (str): method to use for interpolation. This should be one of: + - 'triangualtion' + - 'plane_fitting' + - 'rbf' + """ + NAME_TO_FUNC = { + "triangulation": self.calculate_with_triangualtion, + "plane_fitting": self.calculate_with_plane_fitting, + "rbf": self.calculate_with_rbf, + } + + method = method.lower() + if method not in NAME_TO_FUNC: + raise ValueError(f"`method` must be one of {NAME_TO_FUNC.keys()}.") + + cv_calc_method = NAME_TO_FUNC[method] + out = cv_calc_method(method_kws) + + self.case.fields.conduction_velocity = interpolate_general_cloud_points_onto_surface( + case=self.case, + cloud_values=self.case.analyse.conduction_velocity.values, + cloud_points=self.case.analyse.conduction_velocity.points + ) + return out + + def _preprocess_data(self): + egmX = self.case.electric.bipolar_egm.points + voltage = self.case.electric.bipolar_egm.voltage + LAT = self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time + surface = self.case.create_mesh() + + valid_id, invalid_id = np.where(LAT != -10000), np.where(LAT == -10000) + egmX_clean = np.delete(egmX, np.where(LAT == -10000), 0) + LAT_clean = np.delete(LAT, np.where(LAT == -10000)) + + return egmX_clean, LAT_clean + + def method_trial(self, method): + filename = r'C:\Users\VINUS\Documents\45_OPEN_EP\openep-py\openep-py\openep\_datasets\OpenEP-MATLAB\openep_dataset_2.mat' + case, x_mesh, cv = run_cv(openep_file_name=filename, + cv_method=method, + visualisation=True) + return x_mesh.point_data.active_scalars + + + def calculate_with_triangualtion(self, method_kws=dict()): + + print('---> Starting triangulation') + egmX, LAT = self.prep_egmX, self.prep_LAT + cv = [] + min_theta = method_kws.get('min_theta', 30) + max_elec_distance = method_kws.get('max_elec_dist', 10) + min_elec_distance = method_kws.get('min_elec_dist', 1.5) + min_lat_difference = method_kws.get('min_lat_diff', 2) + cv_centers_x = [] + cv_centers_y = [] + cv_centers_z = [] + cv_centers_id = [] + alpha = method_kws.get('alpha', 5) + + egmX_point_cloud = pv.PolyData(egmX) + print('---> Triangulation process') + surf = egmX_point_cloud.delaunay_3d(alpha=alpha) + print('---> Triangulation process') + print('---> CV calculation (Triangulation method started ...)') + + for num_point in range(len(surf.cells_dict[5])): + + vtx_id = [] + lat = [] + + for i in range(3): + vtx_id.append(surf.cells_dict[5][num_point][i]) + lat.append(LAT[surf.cells_dict[5][num_point][i]]) + + id_lat_sorted = np.argsort(lat) + + O = [egmX[int(vtx_id[id_lat_sorted[0]])], lat[id_lat_sorted[0]]] + A = [egmX[int(vtx_id[id_lat_sorted[1]])], lat[id_lat_sorted[1]]] + B = [egmX[int(vtx_id[id_lat_sorted[2]])], lat[id_lat_sorted[2]]] + + OA = np.sqrt(sum(np.power(np.subtract(O[0], A[0]), 2))) + OB = np.sqrt(sum(np.power(np.subtract(O[0], B[0]), 2))) + AB = np.sqrt(sum(np.power(np.subtract(A[0], B[0]), 2))) + + tOA = A[1] - O[1] + tOB = B[1] - O[1] + + theta = np.arccos((np.power(OA, 2) + np.power(OB, 2) - np.power(AB, 2)) / (2 * OA * OB)) + # check if the conditions are meet to accept the triangle set as a viable acceptable one + if (math.degrees(theta) >= min_theta and OA >= min_elec_distance and OA <= max_elec_distance and + OB >= min_elec_distance and OB <= max_elec_distance and tOA >= min_lat_difference and tOB >= min_lat_difference): + + alpha = np.arctan((tOB * OA - tOA * OB * np.cos(theta)) / (tOA * OB * np.sin(theta))) + cv_temp = (OA / tOA) * np.cos(alpha) + cv.append(cv_temp) + cv_centers_x.append((O[0][0] + A[0][0] + B[0][0]) / 3) + cv_centers_y.append((O[0][1] + A[0][1] + B[0][1]) / 3) + cv_centers_z.append((O[0][2] + A[0][2] + B[0][2]) / 3) + cv_centers_id.append(vtx_id[id_lat_sorted[0]]) + else: + cv.append(np.nan) + cv_centers_x.append((O[0][0] + A[0][0] + B[0][0]) / 3) + cv_centers_y.append((O[0][1] + A[0][1] + B[0][1]) / 3) + cv_centers_z.append((O[0][2] + A[0][2] + B[0][2]) / 3) + cv_centers_id.append(vtx_id[id_lat_sorted[0]]) + ####### Create a triangulation mesh from egmX pointcloud ##### + cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] + print('---> CV calculation ended') + print(cv) + + egmX_point_cloud = pv.PolyData(np.transpose(cv_centers)) + egmX_point_cloud['CV[m/s]'] = cv + # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) + + + egmX_point_cloud = pv.PolyData(np.transpose(cv_centers)) + egmX_point_cloud['CV[m/s]'] = cv + + + # interpolated_cv_field = interpolate_conduction_velocity_onto_surface( + # self.case, + # np.array(cv), + # egmX_point_cloud + # ) + # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) + + self.values = np.array(cv) + self.points = egmX_point_cloud.points + + return np.array(cv), egmX_point_cloud.points + + def calculate_with_plane_fitting(self, method_kws=dict()): + egmX, LAT = self.prep_egmX, self.prep_LAT + + cv = [] + cv_centers_x = [] + cv_centers_y = [] + cv_centers_z = [] + cv_centers_id = [] + cent = np.random.random((len(egmX), 3)) + direction = np.random.random((len(egmX), 3)) + + leaf_size = method_kws.get('leaf_size', 5) + tree = KDTree(egmX, leaf_size=leaf_size) + all_nn_indices = tree.query_radius(egmX, r=10) # 10 closest neighnours + all_nns = [[egmX[idx] for idx in nn_indices] for nn_indices in all_nn_indices] + + min_n_nearest_neighbours = method_kws.get('min_n_nearest_neighbours', 5) + print(len(egmX)) + for k in range(len(egmX)): + print(k, end='\r') + if len(all_nns[k]) >= min_n_nearest_neighbours: + cv_data = np.zeros(shape=(len(all_nns[k]) + 1, 4)) + cv_data[0][:3] = egmX[k] + cv_data[0][3] = LAT[k] + for i in range(len(all_nns[k])): + cv_data[i + 1][:3] = all_nns[k][i] + cv_data[i + 1][3] = LAT[all_nn_indices[k][i]] + + x1 = cv_data[:, 0] + x2 = cv_data[:, 1] + x3 = cv_data[:, 2] + m = cv_data[:, 3] + + res = minimize(self.func, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], args=(x1, x2, x3, m)) + a, b, c, d, e, f, g, h, q, l = res.x + + x = egmX[k, 0] + y = egmX[k, 1] + z = egmX[k, 2] + + tx = 2 * a * x + d * y + e * z + g + ty = 2 * b * y + d * x + f * z + h + tz = 2 * c * z + e * x + f * y + q + + vx = tx / (tx ** 2 + ty ** 2 + tz ** 2) + vy = ty / (tx ** 2 + ty ** 2 + tz ** 2) + vz = tz / (tx ** 2 + ty ** 2 + tz ** 2) + + cent[k, 0] = x + cent[k, 1] = y + cent[k, 2] = z + + direction[k, 0] = vx + direction[k, 1] = vy + direction[k, 2] = vz + + cv_temp = np.sqrt(np.sum(np.power(vx, 2) + np.power(vy, 2) + np.power(vz, 2))) + direction[k, 0] = vx + direction[k, 1] = vy + direction[k, 2] = vz + + if cv_temp >= 0.2 and cv_temp <= 2: + cv.append(cv_temp) + cv_centers_x.append(x) + cv_centers_y.append(y) + cv_centers_z.append(z) + cv_centers_id.append(k) + else: + + cv.append(cv_temp) + cv_centers_x.append(x) + cv_centers_y.append(y) + cv_centers_z.append(z) + else: + x = egmX[k, 0] + y = egmX[k, 1] + z = egmX[k, 2] + cv.append(np.nan) + cv_centers_x.append(x) + cv_centers_y.append(y) + cv_centers_z.append(z) + cv_centers_id.append(k) + cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] + + egmX_point_cloud = pv.PolyData(np.transpose(cv_centers)) + egmX_point_cloud['CV[m/s]'] = cv + + # interpolated_cv_field = interpolate_conduction_velocity_onto_surface( + # self.case, + # np.array(cv), + # egmX_point_cloud + # ) + # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) + + self.values = np.array(cv) + self.points = egmX_point_cloud.points + + return np.array(cv), egmX_point_cloud.points + + + # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) + # return interpolated_mesh_cv.point_data.active_scalars + + def func(self, coef, x1, x2, x3, m): + ''' This the function we use to fit a surface to the activations in plane + fitting method''' + a = coef[0] + b = coef[1] + c = coef[2] + d = coef[3] + e = coef[4] + f = coef[5] + g = coef[6] + h = coef[7] + q = coef[8] + l = coef[9] + # print(np.linalg.norm(m-(a*x1**2 + b*x2**2 + c*x3**2 + d*x1*x2 + e*x1*x3 + f*x2*x3 + g*x1 +h*x2 + q*x3 + l))) + norm = np.linalg.norm(m - ( + a * x1 ** 2 + b * x2 ** 2 + c * x3 ** 2 + d * x1 * x2 + e * x1 * x3 + f * x2 * x3 + g * x1 + h * x2 + q * x3 + l)) + return norm + + def calculate_with_rbf(self): + pass \ No newline at end of file diff --git a/openep/analysis/divergence.py b/openep/analysis/divergence.py new file mode 100644 index 00000000..2a0e3fe6 --- /dev/null +++ b/openep/analysis/divergence.py @@ -0,0 +1,132 @@ +import openep +from openep.analysis.run import run_cv +import sys +import os +import math +import pyvista as pv +import numpy as np +import scipy.io +import scipy +from scipy.optimize import minimize +import matplotlib.pyplot as plt +from sklearn.neighbors import KDTree +import pdb +from scipy.interpolate import Rbf +from vedo import * +import pyvistaqt +from sklearn import metrics + +from ..case.case_routines import ( + interpolate_general_cloud_points_onto_surface, + interpolate_activation_time_onto_surface +) + +_AVAILABLE_CV_METHODS = { + "triangulation", + "plane_fitting", + "rbf", +} + +class Divergence(): + + + def __init__(self, case): + self.case = case #.copy() change to copy later + self.include = None, + self.direction = None, + self.mesh_original = None, + self.div_value = None + self.divergence_per_point=None + + def calculate( + self, + output_binary_field=False, + method_kws=dict() + ): + """Calculate conduction velocity and interpolates each point on a mesh, + stores in conduction_velocity fields. + + Args: + method (str): method to use for interpolation. This should be one of: + - 'triangualtion' + - 'plane_fitting' + - 'rbf' + """ + + out = self.calculate_divergence(method_kws) + + + def calculate_divergence(self, method_kws, output_binary_field=False): + + if output_binary_field: + collision_threshold = method_kws.get("collision_threshold",-1) + focal_threshold = method_kws.get("collision_threshold", 1) + + surface = self.case.create_mesh() + egmX=self.case.electric.bipolar_egm.points + LAT=self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time + points=surface.points + face=surface.faces + + egmX_point_cloud = pv.PolyData(egmX) + egmX_point_cloud['LAT[ms]'] = LAT + mesh_original = pv.PolyData(points, face) + + # interpolated_old = mesh_original.interpolate(egmX_point_cloud, radius=5, n_points=None) + interpolated_scalar=interpolate_general_cloud_points_onto_surface( + case=self.case, + cloud_values=LAT, + cloud_points=egmX, + ) + + mesh_original['LAT_scalar'] = interpolated_scalar + deriv = mesh_original.compute_derivative(scalars='LAT_scalar') + + cv_x = deriv['gradient'][:, 0] / np.sum(np.power(deriv['gradient'], 2), axis=1) + cv_y = deriv['gradient'][:, 1] / np.sum(np.power(deriv['gradient'], 2), axis=1) + cv_z = deriv['gradient'][:, 2] / np.sum(np.power(deriv['gradient'], 2), axis=1) + + direction = np.ndarray(shape=(len(cv_x), 3)) + cv_direction = np.ndarray(shape=(len(cv_x), 3)) + for i in range(len(cv_x)): + direction[i][0] = cv_x[i] / np.sqrt( + np.sum(np.power(cv_x[i], 2) + np.power(cv_y[i], 2) + np.power(cv_z[i], 2))) + direction[i][1] = cv_y[i] / np.sqrt( + np.sum(np.power(cv_x[i], 2) + np.power(cv_y[i], 2) + np.power(cv_z[i], 2))) + direction[i][2] = cv_z[i] / np.sqrt( + np.sum(np.power(cv_x[i], 2) + np.power(cv_y[i], 2) + np.power(cv_z[i], 2))) + cv_direction[i][0] = cv_x[i] + cv_direction[i][1] = cv_y[i] + cv_direction[i][2] = cv_z[i] + + mesh_original['activation_direction'] = cv_direction + div = mesh_original.compute_derivative(scalars='activation_direction', divergence=True) + sargs = dict( + title_font_size=32, + label_font_size=32, + shadow=True, + n_labels=11, + italic=True, + fmt="%.1f", + font_family="arial", + color='black' + ) + + if output_binary_field: + div_value = [] + for i in range(len(div['divergence'])): + if div['divergence'][i] < collision_threshold or div['divergence'][i] > focal_threshold: + div_value.append(1) + else: + div_value.append(0) + else: + div_value = div['divergence'] + + mesh_divergence = pv.PolyData(points, face) + mesh_divergence['div'] = div_value + divergence_per_point = mesh_divergence.cell_data_to_point_data() + + self.case.fields.divergence = div_value + self.divergence_per_point=divergence_per_point + + self.direction, self.mesh_original, self.div_value = direction, mesh_original, div_value From adf2e2d3f8181dc6117c4af7306eaeb4604037be Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:42:53 +0000 Subject: [PATCH 06/31] plot add line connected=True --- openep/draw/draw_routines.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openep/draw/draw_routines.py b/openep/draw/draw_routines.py index e586be3c..cb8d1065 100644 --- a/openep/draw/draw_routines.py +++ b/openep/draw/draw_routines.py @@ -96,6 +96,7 @@ def draw_free_boundaries( color=colours[boundary_index], width=width, name=names[boundary_index], + connected=True ) return plotter From 93d0ccd6791782c6bab68bc09f09bd34b31b4fa5 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:43:19 +0000 Subject: [PATCH 07/31] exporting cv and div data --- openep/io/writers.py | 55 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/openep/io/writers.py b/openep/io/writers.py index 31852185..1c21b079 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -189,8 +189,17 @@ def export_openep_mat( indices=case.indices, fields=case.fields, ) + userdata['surface'] = _add_surface_maps( + surface_data=userdata['surface'], + cv_field=case.fields.conduction_velocity, + divergence_field=case.fields.divergence + ) userdata['electric'] = _extract_electric_data(electric=case.electric) + userdata['electric'] = _add_electric_signal_props( + electric_data=userdata['electric'], + conduction_velocity=case.analyse.conduction_velocity + ) userdata['rf'] = _export_ablation_data(ablation=case.ablation) scipy.io.savemat( @@ -201,6 +210,52 @@ def export_openep_mat( oned_as='column', ) +def _add_surface_maps(surface_data, **kwargs): + cv_field = kwargs.get('cv_field') + div_field = kwargs.get('divergence_field') + + if not surface_data.get('signalMaps'): + surface_data['signalMaps'] = {} + + # TODO: connect propSetting from user setting + if cv_field is not None: + surface_data['signalMaps']['conduction_velocity_field'] = { + 'name' : 'Conduction Velocity Field', + 'value': cv_field, + # 'propSettings':None, + } + + if div_field is not None: + surface_data['signalMaps']['divergence_field'] = { + 'name' : 'Divergence Field', + 'value': div_field, + # 'propSettings':None, + } + + return surface_data + +def _add_electric_signal_props(electric_data, **kwargs): + conduction_velocity = kwargs.get('conduction_velocity') + + if not electric_data.get('annotations'): + electric_data['annotations'] = {} + + if not electric_data['annotations'].get('signalProps'): + electric_data['annotations']['signalProps'] = {} + + signal_props = electric_data['annotations']['signalProps'] + + if conduction_velocity: + #TODO: connect propSetting from user setting + signal_props['conduction_velocity'] = { + 'name' : 'Conduction Velocity Values', + 'value': conduction_velocity.values, + } + signal_props['cvX'] = { + 'name' : 'Conduction Velocity Coordinates', + 'value': conduction_velocity.points, + } + return electric_data def export_vtk( case: Case, From 05a718fc7c5af51fd6e7a2710858404dc5d1af83 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:43:53 +0000 Subject: [PATCH 08/31] adding cv field attributes --- openep/data_structures/surface.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index 3cd50afa..5b262c22 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -52,6 +52,7 @@ class Fields: transverse_fibres: np.ndarray = None pacing_site: np.ndarray = None conduction_velocity: np.ndarray = None + divergence: np.ndarray = None mesh_normals : np.ndarray = None def __repr__(self): @@ -251,7 +252,7 @@ def empty_fields(n_points=0, n_cells=0): longitudinal_fibres, transverse_fibres, pacing_site, - conudction_vlocity, + conduction_velocity, mesh_normals, ) From 80c2b169fec13c3eb2c39cdfbe9f7da4a00d1976 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:48:03 +0000 Subject: [PATCH 09/31] interpolate cloud to surface function --- openep/case/__init__.py | 1 + openep/case/case_routines.py | 50 ++++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/openep/case/__init__.py b/openep/case/__init__.py index 61a81044..e6e4d3b5 100644 --- a/openep/case/__init__.py +++ b/openep/case/__init__.py @@ -9,5 +9,6 @@ Interpolator, interpolate_activation_time_onto_surface, interpolate_voltage_onto_surface, + interpolate_general_cloud_points_onto_surface, bipolar_from_unipolar_surface_points, ) diff --git a/openep/case/case_routines.py b/openep/case/case_routines.py index d3bbf431..d924ac11 100644 --- a/openep/case/case_routines.py +++ b/openep/case/case_routines.py @@ -453,6 +453,56 @@ def __call__(self, surface_points, max_distance=None): def __repr__(self): return f"Interpolator: method={self.method}, kws={self.method_kws}" +def interpolate_general_cloud_points_onto_surface( + case, + cloud_values, + cloud_points, + method=scipy.interpolate.RBFInterpolator, + method_kws=None, + max_distance=None, + include=None, +): + """Interpolate general point data onto the points of a mesh. + + Args: + case (openep.case.Case): case from which the conduction velocity will be interpolated + method (callable): method to use for interpolation. The default is + scipy.interpolate.RBFInterpolator. + method_kws (dict): dictionary of keyword arguments to pass to `method` + when creating the interpolator. + max_distance (float, optional): If provided, any points on the surface of the mesh + further than this distance to all mapping coordinates will have their + interpolated activation times set NaN. The default it None, in which case + the distance from surface points to mapping points is not considered. + include (np.ndarray, optional): Flag for which mapping points to include when creating + the interpolator. If None, `case.electric.include` will be used. + + Returns: + interpolated_lat (ndarray): local activation times interpolated onto the surface of the mesh, + one value per point on the mesh. + """ + + surface_points = case.points + + include = [not np.isnan(x) for x in cloud_values] if include is None else include + cloud_points = cloud_points[include] + cloud_values = cloud_values[include] + + interpolator = Interpolator( + cloud_points, + cloud_values, + method=method, + method_kws=method_kws, + ) + + interpolated = interpolator(surface_points, max_distance=max_distance) + + # Any points that are not part of the mesh faces should have CV set to NaN + n_surface_points = surface_points.shape[0] + not_on_surface = ~np.in1d(np.arange(n_surface_points), case.indices) + interpolated[not_on_surface] = np.NaN + + return interpolated def interpolate_activation_time_onto_surface( case, From 9e9cc826c27ad0242177dd6b7a6a566ad5a2077a Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:49:19 +0000 Subject: [PATCH 10/31] adding analysis attribute to case --- openep/data_structures/case.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openep/data_structures/case.py b/openep/data_structures/case.py index f4ced3b6..57c64b01 100644 --- a/openep/data_structures/case.py +++ b/openep/data_structures/case.py @@ -85,6 +85,7 @@ from .surface import Fields from .electric import Electric, Electrogram, Annotations, ElectricSurface from .ablation import Ablation +from ..analysis.analyse import Analyse from ..case.case_routines import ( bipolar_from_unipolar_surface_points, calculate_distance, @@ -130,6 +131,7 @@ def __init__( self.ablation = ablation self.electric = electric self.notes = notes + self.analyse = Analyse(case=self) def __repr__(self): return f"{self.name}( nodes: {self.points.shape} indices: {self.indices.shape} {self.fields} )" From e05d2b11fdbd2b3b29f7f95f31e49e7f6bfd9659 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:49:49 +0000 Subject: [PATCH 11/31] cv and divergence for pyvista converters --- openep/converters/pyvista_converters.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openep/converters/pyvista_converters.py b/openep/converters/pyvista_converters.py index 0c602607..18f8501e 100644 --- a/openep/converters/pyvista_converters.py +++ b/openep/converters/pyvista_converters.py @@ -93,5 +93,7 @@ def to_pyvista( mesh.point_data["LAT"] = case.fields.local_activation_time mesh.point_data["Impedance"] = case.fields.impedance mesh.point_data["Force"] = case.fields.force + mesh.point_data["Conduction Velocity"] = case.fields.conduction_velocity + mesh.point_data["Divergence"] = case.fields.divergence return mesh From bfd17c181083117bcb3f0454079828025d319f98 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 09:50:32 +0000 Subject: [PATCH 12/31] analysis package --- openep/analysis/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 openep/analysis/__init__.py diff --git a/openep/analysis/__init__.py b/openep/analysis/__init__.py new file mode 100644 index 00000000..e69de29b From f5ce7f4b59e5747e286667ce381269ef003952ce Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 10:39:29 +0000 Subject: [PATCH 13/31] remove run_cv --- openep/analysis/conduction_velocity.py | 9 --------- openep/analysis/divergence.py | 1 - 2 files changed, 10 deletions(-) diff --git a/openep/analysis/conduction_velocity.py b/openep/analysis/conduction_velocity.py index f927bec5..26bb052c 100644 --- a/openep/analysis/conduction_velocity.py +++ b/openep/analysis/conduction_velocity.py @@ -1,5 +1,4 @@ import openep -from openep.analysis.run import run_cv import sys import os import math @@ -77,14 +76,6 @@ def _preprocess_data(self): return egmX_clean, LAT_clean - def method_trial(self, method): - filename = r'C:\Users\VINUS\Documents\45_OPEN_EP\openep-py\openep-py\openep\_datasets\OpenEP-MATLAB\openep_dataset_2.mat' - case, x_mesh, cv = run_cv(openep_file_name=filename, - cv_method=method, - visualisation=True) - return x_mesh.point_data.active_scalars - - def calculate_with_triangualtion(self, method_kws=dict()): print('---> Starting triangulation') diff --git a/openep/analysis/divergence.py b/openep/analysis/divergence.py index 2a0e3fe6..e3c3a8bd 100644 --- a/openep/analysis/divergence.py +++ b/openep/analysis/divergence.py @@ -1,5 +1,4 @@ import openep -from openep.analysis.run import run_cv import sys import os import math From ae13f554566fa08ddcb79a17b5ebbf4fe60500da Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Tue, 21 Nov 2023 11:10:46 +0000 Subject: [PATCH 14/31] added divergence values to export openep --- openep/analysis/divergence.py | 10 ++-------- openep/io/writers.py | 13 ++++++++++++- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/openep/analysis/divergence.py b/openep/analysis/divergence.py index e3c3a8bd..f67391f5 100644 --- a/openep/analysis/divergence.py +++ b/openep/analysis/divergence.py @@ -20,12 +20,6 @@ interpolate_activation_time_onto_surface ) -_AVAILABLE_CV_METHODS = { - "triangulation", - "plane_fitting", - "rbf", -} - class Divergence(): @@ -34,7 +28,7 @@ def __init__(self, case): self.include = None, self.direction = None, self.mesh_original = None, - self.div_value = None + self.values = None self.divergence_per_point=None def calculate( @@ -128,4 +122,4 @@ def calculate_divergence(self, method_kws, output_binary_field=False): self.case.fields.divergence = div_value self.divergence_per_point=divergence_per_point - self.direction, self.mesh_original, self.div_value = direction, mesh_original, div_value + self.direction, self.mesh_original, self.values = direction, mesh_original, div_value diff --git a/openep/io/writers.py b/openep/io/writers.py index 1c21b079..ca76ae71 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -198,7 +198,8 @@ def export_openep_mat( userdata['electric'] = _extract_electric_data(electric=case.electric) userdata['electric'] = _add_electric_signal_props( electric_data=userdata['electric'], - conduction_velocity=case.analyse.conduction_velocity + conduction_velocity=case.analyse.conduction_velocity, + divergence=case.analyse.divergence ) userdata['rf'] = _export_ablation_data(ablation=case.ablation) @@ -218,6 +219,7 @@ def _add_surface_maps(surface_data, **kwargs): surface_data['signalMaps'] = {} # TODO: connect propSetting from user setting + # TODO: handle with None values as matlab cannot handle if cv_field is not None: surface_data['signalMaps']['conduction_velocity_field'] = { 'name' : 'Conduction Velocity Field', @@ -236,6 +238,7 @@ def _add_surface_maps(surface_data, **kwargs): def _add_electric_signal_props(electric_data, **kwargs): conduction_velocity = kwargs.get('conduction_velocity') + divergence = kwargs.get('divergence') if not electric_data.get('annotations'): electric_data['annotations'] = {} @@ -255,6 +258,14 @@ def _add_electric_signal_props(electric_data, **kwargs): 'name' : 'Conduction Velocity Coordinates', 'value': conduction_velocity.points, } + + if divergence: + #TODO: connect propSetting from user setting + signal_props['divergence'] = { + 'name' : 'Divergence Values', + 'value': divergence.values, + } + return electric_data def export_vtk( From 1a57cfaa2651fef3556a3fe7790a01f911e77cd0 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Wed, 22 Nov 2023 07:52:43 +0000 Subject: [PATCH 15/31] quick fix to get remove_points working --- openep/analysis/conduction_velocity.py | 24 +++++++++++++++++++----- openep/analysis/divergence.py | 22 ++++++++++++++++------ 2 files changed, 35 insertions(+), 11 deletions(-) diff --git a/openep/analysis/conduction_velocity.py b/openep/analysis/conduction_velocity.py index 26bb052c..60e56556 100644 --- a/openep/analysis/conduction_velocity.py +++ b/openep/analysis/conduction_velocity.py @@ -24,15 +24,16 @@ class ConductionVelocity(): def __init__(self, case): self.case = case - self.prep_egmX, self.prep_LAT = self._preprocess_data() self.include = None self.values = None self.points = None + self.prep_egmX, self.prep_LAT = self._preprocess_data() def calculate( self, method='triangulation', - method_kws=dict() + method_kws=dict(), + include=None ): """Calculate conduction velocity and interpolates each point on a mesh, @@ -44,6 +45,10 @@ def calculate( - 'plane_fitting' - 'rbf' """ + if include is not None: + self.include=include + self.prep_egmX, self.prep_LAT = self._preprocess_data() + NAME_TO_FUNC = { "triangulation": self.calculate_with_triangualtion, "plane_fitting": self.calculate_with_plane_fitting, @@ -65,9 +70,18 @@ def calculate( return out def _preprocess_data(self): - egmX = self.case.electric.bipolar_egm.points - voltage = self.case.electric.bipolar_egm.voltage - LAT = self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time + if self.include is None: + egmX = self.case.electric.bipolar_egm.points + #TODO: voltage is not even used! + voltage = self.case.electric.bipolar_egm.voltage + LAT = self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time + else: + egmX = self.case.electric.bipolar_egm.points[self.include] + voltage = self.case.electric.bipolar_egm.voltage[self.include] + LAT = self.case.electric.annotations.local_activation_time[self.include] - self.case.electric.annotations.reference_activation_time[self.include] + + + print('len_cv:', len(egmX)) surface = self.case.create_mesh() valid_id, invalid_id = np.where(LAT != -10000), np.where(LAT == -10000) diff --git a/openep/analysis/divergence.py b/openep/analysis/divergence.py index f67391f5..5b36b577 100644 --- a/openep/analysis/divergence.py +++ b/openep/analysis/divergence.py @@ -25,16 +25,17 @@ class Divergence(): def __init__(self, case): self.case = case #.copy() change to copy later - self.include = None, - self.direction = None, - self.mesh_original = None, + self.include = None + self.direction = None + self.mesh_original = None self.values = None self.divergence_per_point=None def calculate( self, output_binary_field=False, - method_kws=dict() + method_kws=dict(), + include=None, ): """Calculate conduction velocity and interpolates each point on a mesh, stores in conduction_velocity fields. @@ -45,6 +46,8 @@ def calculate( - 'plane_fitting' - 'rbf' """ + if include is not None: + self.include=include out = self.calculate_divergence(method_kws) @@ -56,8 +59,15 @@ def calculate_divergence(self, method_kws, output_binary_field=False): focal_threshold = method_kws.get("collision_threshold", 1) surface = self.case.create_mesh() - egmX=self.case.electric.bipolar_egm.points - LAT=self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time + + if self.include is None: + egmX = self.case.electric.bipolar_egm.points + LAT = self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time + else: + egmX=self.case.electric.bipolar_egm.points[self.include] + LAT=self.case.electric.annotations.local_activation_time[self.include] - self.case.electric.annotations.reference_activation_time[self.include] + print("len_div:", len(egmX)) + points=surface.points face=surface.faces From 62585292afb9fa309f1a01360a217723042506bb Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 10:39:53 +0000 Subject: [PATCH 16/31] cleaned up general interpolator docstring --- openep/case/case_routines.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/openep/case/case_routines.py b/openep/case/case_routines.py index d924ac11..12fe6e68 100644 --- a/openep/case/case_routines.py +++ b/openep/case/case_routines.py @@ -453,6 +453,7 @@ def __call__(self, surface_points, max_distance=None): def __repr__(self): return f"Interpolator: method={self.method}, kws={self.method_kws}" + def interpolate_general_cloud_points_onto_surface( case, cloud_values, @@ -465,7 +466,9 @@ def interpolate_general_cloud_points_onto_surface( """Interpolate general point data onto the points of a mesh. Args: - case (openep.case.Case): case from which the conduction velocity will be interpolated + case (openep.case.Case): case from which the cloud values and points will be interpolated + cloud_values (ndarray): Array of values to be interpolated, corresponding to cloud_points. + cloud_points (ndarray): Array of points where cloud_values are defined. method (callable): method to use for interpolation. The default is scipy.interpolate.RBFInterpolator. method_kws (dict): dictionary of keyword arguments to pass to `method` @@ -497,7 +500,7 @@ def interpolate_general_cloud_points_onto_surface( interpolated = interpolator(surface_points, max_distance=max_distance) - # Any points that are not part of the mesh faces should have CV set to NaN + # Any points that are not part of the mesh faces should have its value set to NaN n_surface_points = surface_points.shape[0] not_on_surface = ~np.in1d(np.arange(n_surface_points), case.indices) interpolated[not_on_surface] = np.NaN From b683bee1b43d33b19f43eaad161b955177fa685d Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 10:40:37 +0000 Subject: [PATCH 17/31] added conduction velocity and divergence class to analyse.py file --- openep/analysis/analyse.py | 225 ++++++++++++++++++++++++++++++++++++- 1 file changed, 222 insertions(+), 3 deletions(-) diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py index 16aa6989..d4258515 100644 --- a/openep/analysis/analyse.py +++ b/openep/analysis/analyse.py @@ -1,8 +1,227 @@ -from .conduction_velocity import ConductionVelocity -from .divergence import Divergence +# OpenEP +# Copyright (c) 2021 OpenEP Collaborators +# +# This file is part of OpenEP. +# +# OpenEP 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. +# +# OpenEP 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. +# +# You should have received a copy of the GNU General Public License along +# with this program (LICENSE.txt). If not, see + +"""Module containing analysis classes""" + +from ._conduction_velocity import * +from ..case.case_routines import interpolate_general_cloud_points_onto_surface + class Analyse: + """ + Analyse class for managing analysis of case data. + + Args: + case (Case): The case data to be analysed. + + Attributes: + conduction_velocity (ConductionVelocity): An instance of ConductionVelocity for conducting + velocity analysis on the case. + + divergence (Divergence): An instance of Divergence for conducting divergence analysis on the case. + """ def __init__(self, case): self.conduction_velocity = ConductionVelocity(case) - self.divergence = Divergence(case) \ No newline at end of file + self.divergence = Divergence(case) + + +class ConductionVelocity: + """ + Class for calculating and analyzing conduction velocity. + + This class provides methods to calculate conduction velocity from case data and optionally + interpolate these values onto a mesh. It handles different methods of calculation and manages + the resulting data. + + Attributes: + _case (Case): Internal reference to the case data. + values (np.ndarray): The calculated conduction velocity values. + points (np.ndarray): The points corresponding to the conduction velocity values. + + Parameters: + case (Case): The case data from which conduction velocity is to be calculated. + """ + def __init__(self, case): + self._case = case + self._values = None + self._points = None + + @property + def values(self): + if not self._values: + raise ValueError('Before accessing ``conduction_velocity.values`` ' + 'run ``divergence.calculate_divergence()``') + return self._values + + @property + def points(self): + if not self._points: + raise ValueError('Before accessing ``conduction_velocity.points`` ' + 'run ``divergence.calculate_divergence()``') + return self._points + + def calculate_cv( + self, + method='triangulation', + include=None, + apply_scalar_field=True, + **kwargs + ): + """ + Calculate the conduction velocity and interpolate points onto a mesh. + + This method calculates conduction velocities using one of the specified methods and + optionally interpolates the calculated values onto a mesh, storing the results in the + 'conduction_velocity' field of the case object. + + Args: + method (str): The method to use for calculating conduction velocity. Options are: + - 'triangulation' + - 'plane_fitting' + - 'rbf' + Defaults to 'triangulation'. + + include (np.ndarray, optional): A boolean array specifying which data points to include + in the preprocessing for bipolar EGM data and LAT. + Defaults to "includes" specified in the imported data. + + apply_scalar_field (bool, optional): If True, applies the calculated conduction velocity + values as a scalar field on the case object's mesh. + Defaults to True. + **kwargs: Additional keyword arguments specific to the chosen calculation method. + + Returns: + tuple: A tuple containing two elements: + - values (np.ndarray): An array of calculated conduction velocity values. + - points (np.ndarray): An array of points corresponding to the calculated values. + + Raises: + ValueError: If the specified method is not among the available options. + + Example: + >> cv = ConductionVelocity(case) + >> values, points = cv.calculate(method='plane_fitting', apply_scalar_field=True) + """ + + supported_cv_methods = { + "triangulation": triangulation, + "plane_fitting": plane_fitting, + "rbf": radial_basis_function + } + + if method.lower() not in supported_cv_methods: + raise ValueError(f"`method` must be one of {supported_cv_methods.keys()}.") + + include = self._case.electric.include.astype(bool) if include is None else include + lat, bipolar_egm_pts = preprocess_lat_egm(self._case, include) + + cv_method = supported_cv_methods[method] + self._values, self._points = cv_method(bipolar_egm_pts, lat, **kwargs) + + if apply_scalar_field: + self._case.fields.conduction_velocity = interpolate_general_cloud_points_onto_surface( + case=self._case, + cloud_values=self.values, + cloud_points=self.points, + ) + + return self.values, self.points + + +class Divergence: + """ + Class for calculating and analyzing divergence of conduction velocity. + + This class provides methods to calculate the divergence of conduction velocity from case data. + It can output binary fields indicating regions of divergence and manage the resulting data. + + Attributes: + _case (Case): Internal reference to the case data. + direction (np.ndarray): The calculated divergence direction vectors. + values (np.ndarray): The calculated divergence values. + + Args: + case (Case): The case data from which divergence is to be calculated. + """ + def __init__(self, case): + self._case = case + self._direction = None + self._values = None + + @property + def values(self): + if not self._values: + raise ValueError('Before accessing ``..divergence.values`` run ``..divergence.calculate_divergence()``') + return self._values + + @property + def direction(self): + if not self._direction: + raise ValueError('Before accessing ``divergence.direction`` run ``divergence.calculate_divergence()``') + return self._direction + + def calculate_divergence( + self, + include=None, + output_binary_field=False, + apply_scalar_field=True + ): + """ + Calculate the divergence of conduction velocity and optionally apply the result as a scalar field. + + Parameters: + include (np.ndarray, optional): A boolean array specifying which data points to include in the + calculation. If None, the 'include' array from the case's electric + data is used. Defaults to None. + + output_binary_field (bool, optional): If True, the output will be a binary field where 1 indicates + regions of divergence and 0 otherwise. + Defaults to False. + + apply_scalar_field (bool, optional): If True, the calculated divergence values (or binary field, if + 'output_binary_field' is True) are applied as a scalar field to + the case object. + Defaults to True. + + Returns: + tuple of (np.ndarray, np.ndarray): A tuple containing the divergence direction and divergence values. + The direction is a 2D array of shape (N, 3), where N is the number + of points, and each row represents the (x, y, z) components of the + divergence direction. The values are a 1D array of divergence values. + + Example: + case = load_case_data('path/to/case_data') + cv = ConductionVelocity(case) + 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) + + self._direction, self._values = divergence( + case=self._case, + bipolar_egm_pts=bipolar_egm_pts, + local_activation_time=lat, + output_binary_field=output_binary_field + ) + + if apply_scalar_field: + self._case.fields.divergence = self.values + + return self.direction, self.values From 2990654b59b39665de69b5284bc535d1f57d84aa Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 10:40:59 +0000 Subject: [PATCH 18/31] removed rough code for CV + div development --- openep/Analyses/CV.py | 572 ------------------------- openep/Analyses/run.py | 72 ---- openep/analysis/conduction_velocity.py | 302 ------------- openep/analysis/divergence.py | 135 ------ 4 files changed, 1081 deletions(-) delete mode 100644 openep/Analyses/CV.py delete mode 100644 openep/Analyses/run.py delete mode 100644 openep/analysis/conduction_velocity.py delete mode 100644 openep/analysis/divergence.py diff --git a/openep/Analyses/CV.py b/openep/Analyses/CV.py deleted file mode 100644 index 3ef3ad57..00000000 --- a/openep/Analyses/CV.py +++ /dev/null @@ -1,572 +0,0 @@ -import sys -import os -import math -import pyvista as pv -import numpy as np -import scipy.io -import scipy -from scipy.optimize import minimize -import matplotlib.pyplot as plt -from sklearn.neighbors import KDTree -import pdb -from scipy.interpolate import Rbf -from vedo import * -import pyvistaqt -from sklearn import metrics - -def func(coef, x1, x2, x3, m): - ''' This the function we use to fit a surface to the activations in plane - fitting method''' - a = coef[0] - b = coef[1] - c = coef[2] - d = coef[3] - e = coef[4] - f = coef[5] - g = coef[6] - h = coef[7] - q = coef[8] - l = coef[9] - #print(np.linalg.norm(m-(a*x1**2 + b*x2**2 + c*x3**2 + d*x1*x2 + e*x1*x3 + f*x2*x3 + g*x1 +h*x2 + q*x3 + l))) - return np.linalg.norm(m-(a*x1**2 + b*x2**2 + c*x3**2 + d*x1*x2 + e*x1*x3 + f*x2*x3 + g*x1 +h*x2 + q*x3 + l)) -class CV: - - def Triangulation(self,egmX,LAT,data_type=''): - ''' This function calculates local atrial conduction velocity from recorded - electro-anatomical data, using triangulation method. (C.D. Cantwell et.al. 2015 - ''Techniques for automated local activation time annotation and conduction velocity estimation in cardiac mapping") - egmX: Recorded Electrode locations (x,y,z) - LAT : Recorded local activation local activation time''' - ##### initiation ###### - cv = [] - min_theta = 30 - max_elec_distance = 10 - min_elec_distance = 1.5 - min_lat_difference = 2 - cv_centers_x = [] - cv_centers_y = [] - cv_centers_z = [] - cv_centers_id = [] - alpha = 5 - - ######## Create a triangulation mesh from the recording electrodes - egmX_point_cloud = pv.PolyData(egmX) - print('---> Triangulation process') - surf = egmX_point_cloud.delaunay_3d(alpha=alpha) - print('---> Triangulation process') - print('---> CV calculation (Triangulation method started ...)') - for num_point in range(len(surf.cells_dict[5])): - - vtx_id = [] - lat = [] - - for i in range (3): - vtx_id.append(surf.cells_dict[5][num_point][i]) - lat.append(LAT[surf.cells_dict[5][num_point][i]]) - - id_lat_sorted = np.argsort(lat) - - O = [egmX[int(vtx_id[id_lat_sorted[0]])],lat[id_lat_sorted[0]]] - A = [egmX[int(vtx_id[id_lat_sorted[1]])],lat[id_lat_sorted[1]]] - B = [egmX[int(vtx_id[id_lat_sorted[2]])],lat[id_lat_sorted[2]]] - - OA = np.sqrt(sum(np.power(np.subtract(O[0],A[0]),2))) - OB = np.sqrt(sum(np.power(np.subtract(O[0],B[0]),2))) - AB = np.sqrt(sum(np.power(np.subtract(A[0],B[0]),2))) - - tOA = A[1] - O[1] - tOB = B[1] - O[1] - - theta = np.arccos((np.power(OA,2)+ np.power(OB,2) - np.power(AB,2))/(2 * OA * OB)) - # check if the conditions are meet to accept the triangle set as a viable acceptable one - if (math.degrees(theta) >= min_theta and OA >= min_elec_distance and OA <= max_elec_distance and - OB >= min_elec_distance and OB <= max_elec_distance and tOA >= min_lat_difference and tOB >= min_lat_difference): - - alpha = np.arctan((tOB * OA - tOA * OB * np.cos(theta)) / (tOA * OB * np.sin(theta))) - cv_temp = (OA/tOA) * np.cos(alpha) - cv.append(cv_temp) - cv_centers_x.append((O[0][0] + A[0][0] + B[0][0])/3) - cv_centers_y.append((O[0][1] + A[0][1] + B[0][1])/3) - cv_centers_z.append((O[0][2] + A[0][2] + B[0][2])/3) - cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - #if cv_temp >= 0.2 and cv_temp <=2: - # cv.append(cv_temp) - # cv_centers_x.append(O[0][0]) - # cv_centers_y.append(O[0][1]) - # cv_centers_z.append(O[0][2]) - # cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - else: - # cv_temp = (OA/tOA) * np.cos(alpha) - # cv_centers_x.append(O[0][0]) - # cv_centers_y.append(O[0][1]) - # cv_centers_z.append(O[0][2]) - # cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - cv.append(np.nan) - cv_centers_x.append((O[0][0] + A[0][0] + B[0][0])/3) - cv_centers_y.append((O[0][1] + A[0][1] + B[0][1])/3) - cv_centers_z.append((O[0][2] + A[0][2] + B[0][2])/3) - cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - ####### Create a triangulation mesh from egmX pointcloud ##### - cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] - print('---> CV calculation ended') - return cv, cv_centers, cv_centers_id - - def plane_fitting(self,egmX, LAT): - - cv = [] - cv_centers_x = [] - cv_centers_y = [] - cv_centers_z = [] - cv_centers_id = [] - cent = np.random.random((len(egmX), 3)) - direction = np.random.random((len(egmX), 3)) - - tree = KDTree(egmX, leaf_size=5) - all_nn_indices = tree.query_radius(egmX, r=10) - all_nns = [[egmX[idx] for idx in nn_indices] for nn_indices in all_nn_indices] - - print(len(egmX)) - for k in range(len(egmX)): - print(k,end='\r') - if len(all_nns[k]) >= 5: - cv_data = np.zeros(shape= (len(all_nns[k]) + 1,4)) - cv_data[0][:3] = egmX[k] - cv_data[0][3] = LAT[k] - for i in range(len(all_nns[k])): - cv_data[i+1][:3] = all_nns[k][i] - cv_data[i+1][3] = LAT[all_nn_indices[k][i]] - - x1 = cv_data[:,0] - x2 = cv_data[:,1] - x3 = cv_data[:,2] - m = cv_data[:,3] - #mesh_original = pv.PolyData(Points,Face) - #p2.add_points(cv_data[:,:3], render_points_as_spheres=True, point_size=5.0, color = 'yellow') - #p2.add_mesh(mesh_original,opacity=0.4) - #p2.show() - - res = minimize(func, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], args=(x1, x2, x3, m)) - a, b, c, d,e, f, g, h, q, l= res.x - #print(a*x1**2+b*x2**2+c*x3**2+d*x1*x2+e*x1*x3+f*x2*x3+g*x1+h*x2+q*x3+l -m) - - - x = egmX[k,0] - y = egmX[k,1] - z = egmX[k,2] - - tx = 2*a*x + d*y + e*z + g - ty = 2*b*y + d*x + f*z + h - tz = 2*c*z + e*x + f*y + q - - vx = tx/(tx**2 + ty**2 + tz**2) - vy = ty/(tx**2 + ty**2 + tz**2) - vz = tz/(tx**2 + ty**2 + tz**2) - - - cent[k,0] = x - cent[k,1] = y - cent[k,2] = z - - direction[k,0] = vx - direction[k,1] = vy - direction[k,2] = vz - - cv_temp = np.sqrt(np.sum(np.power(vx,2)+np.power(vy,2)+np.power(vz,2))) - direction[k,0] = vx - direction[k,1] = vy - direction[k,2] = vz - - if cv_temp >=0.2 and cv_temp <=2: - cv.append(cv_temp) - cv_centers_x.append(x) - cv_centers_y.append(y) - cv_centers_z.append(z) - cv_centers_id.append(k) - else: - - cv.append(cv_temp) - cv_centers_x.append(x) - cv_centers_y.append(y) - cv_centers_z.append(z) - else: - x = egmX[k,0] - y = egmX[k,1] - z = egmX[k,2] - cv.append(np.nan) - cv_centers_x.append(x) - cv_centers_y.append(y) - cv_centers_z.append(z) - cv_centers_id.append(k) - cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] - - return cv, cv_centers, cv_centers_id - - def RBF_method(self, egmX, LAT, points, face,Fiber_data='',): - - #rbfi = Rbf(egmX[:,0],egmX[:,1],egmX[:,2],LAT,kernel='gaussian',epsilon=0.5) - rbfi = Rbf(egmX[:,0],egmX[:,1],egmX[:,2],LAT) - lat_i = [] - for i in range(len(points)): - lat_i.append(rbfi(points[i][0], points[i][1], points[i][2])) - - mesh_original_RBF = pv.PolyData(points,face) - mesh_original_RBF['values'] = lat_i - deriv = mesh_original_RBF.compute_derivative('values') - df = deriv['gradient'] - - cvx = deriv['gradient'][:,0]/np.sum(np.power( deriv['gradient'],2),axis=1) - cvy = deriv['gradient'][:,1]/np.sum(np.power( deriv['gradient'],2),axis=1) - cvz = deriv['gradient'][:,2]/np.sum(np.power( deriv['gradient'],2),axis=1) - cv = np.sqrt(np.power(cvx,2) + np.power(cvy,2) + np.power(cvz,2)) - - cv_data = [] - cv_x = [] - cv_y = [] - cv_z = [] - cv_center_id =[] - cv_centers_x = [] - cv_centers_y = [] - cv_centers_z = [] - tree = KDTree(points, leaf_size=5) - dist, ind = tree.query(egmX, k=1) - for k in range(len(ind)): - - if cv[ind[k][0]] > 0 and cv[ind[k][0]] <= 10: - cv_data.append(cv[ind[k][0]]) - cv_x.append(cvx[ind[k][0]]) - cv_y.append(cvy[ind[k][0]]) - cv_z.append(cvz[ind[k][0]]) - cv_center_id.append(ind[k][0]) - #cv_centers_x.append(points[ind[k][0]][0]) - #cv_centers_y.append(points[ind[k][0]][1]) - #cv_centers_z.append(points[ind[k][0]][2]) - cv_centers_x.append(egmX[k][0]) - cv_centers_y.append(egmX[k][1]) - cv_centers_z.append(egmX[k][2]) - - - else: - cv_data.append(np.nan) - cv_centers_x.append(points[ind[k][0]][0]) - cv_centers_y.append(points[ind[k][0]][1]) - cv_centers_z.append(points[ind[k][0]][2]) - if 0: - cv_data.append(cv[ind[k][0]]) - cv_center_id.append(ind[k][0]) - cv_centers_x.append(points[ind[k][0]][0]) - cv_centers_y.append(points[ind[k][0]][1]) - cv_centers_z.append(points[ind[k][0]][2]) - cv_x.append(cvx[ind[k][0]]) - cv_y.append(cvy[ind[k][0]]) - cv_z.append(cvz[ind[k][0]]) - - direction = [cv_x,cv_y,cv_z] - cv_centers = [cv_centers_x,cv_centers_y,cv_centers_z] - - - - - return cv_data, cv_centers, direction, cv_center_id - - def visualization(self,lat_data = '', - cv_data = '', - mesh = '', - egmX = '', - method = '', - file_name = '', - direction='', - direction_centers = '', - show_histogram= False, - show_cv_fibrosis = True, - interpolation_radius = 5): - - boring_cmap = plt.cm.get_cmap("jet", 256) - boring_cmap_r = plt.cm.get_cmap("jet_r", 256) - boring_cmap_2 = plt.cm.get_cmap("Paired", 4) - - boring_cmap_2.colors[0][0]= 34/255 #,237,222#1 - boring_cmap_2.colors[0][1]=94/255 - boring_cmap_2.colors[0][2]=168/255 - boring_cmap_2.colors[1][0]= 62/255 #,190,133#0 - boring_cmap_2.colors[1][1]=182/255 #1 - boring_cmap_2.colors[1][2]=196/255 #1 - boring_cmap_2.colors[2][0]=253/255 - boring_cmap_2.colors[2][1]=141/255 - boring_cmap_2.colors[2][2]= 60/255 - boring_cmap_2.colors[3][0]=217/255 - boring_cmap_2.colors[3][1]=71/255 - boring_cmap_2.colors[3][2]=1/255 - - - - - - pv.rcParams['transparent_background'] = True - - - - if method == 'LAT': - #pl = pv.Plotter(off_screen=True) - pl = pv.Plotter() - egmX_point_cloud = pv.PolyData(egmX) - egmX_point_cloud['LAT[ms]'] = lat_data - interpolated_mesh_lat = mesh.interpolate(egmX_point_cloud, radius=interpolation_radius) - pl.add_mesh(egmX_point_cloud,scalars="LAT[ms]",render_points_as_spheres=True, - point_size=20, cmap = boring_cmap_r) - pl.add_mesh(interpolated_mesh_lat,scalars="LAT[ms]", clim = [np.min(lat_data), np.max(lat_data)], - below_color='brown', above_color = 'magenta', cmap = boring_cmap_r) - if direction != '' and direction_centers != '': - pl.add_arrows(direction_centers, direction, mag=2, color=[90, 90, 90]) - #here - pl.show(screenshot=f'/Users/aligharaviri/Desktop/CV_Data/Figures_picsLAT_{file_name}.png') - if method == 'CV': - pl = pv.Plotter(shape=(1,2),off_screen=False) - sargs = dict( - title_font_size=32, - label_font_size=32, - shadow=True, - n_labels=11, - italic=True, - fmt="%.1f", - font_family="arial", - color='black' - ) - egmX_point_cloud = pv.PolyData(egmX) - egmX_point_cloud['CV[m/s]'] = cv_data - ######## Ver 1 - #interpolated_mesh_cv = mesh.interpolate(egmX_point_cloud,radius=5 ,n_points=2 ,sharpness=2) - ######## Ver 1 - interpolated_mesh_cv = mesh.interpolate(egmX_point_cloud,radius=10 ,n_points=None,null_value=np.nan) - - pl.subplot(0,0) - - pl.add_mesh(egmX_point_cloud,scalars='CV[m/s]',render_points_as_spheres=True, - point_size=20, cmap = boring_cmap,clim=[0,1.5],below_color='navy', above_color = 'darkred') - pl.add_mesh(interpolated_mesh_cv,scalars='CV[m/s]', clim=[0,1.5], - below_color='navy', above_color = 'darkred', cmap = boring_cmap,scalar_bar_args=sargs) - - if direction != '': - pl.add_arrows(egmX, np.transpose(direction), mag=4, color='black')#color=[90, 90, 90]) - pl.camera.roll += 20 - pl.camera.elevation -=10 - pl.camera.azimuth = 10 - pl.camera.zoom (1.2) - pl.subplot(0,1) - - pl.add_mesh(egmX_point_cloud,scalars='CV[m/s]',render_points_as_spheres=True, - point_size=20, cmap = boring_cmap,clim=[0,1.2],below_color='navy', above_color = 'darkred') - pl.add_mesh(interpolated_mesh_cv,scalars='CV[m/s]', clim=[0,1.5], - below_color='navy', above_color = 'darkred', cmap = boring_cmap,scalar_bar_args=sargs) - if direction != '': - pl.add_arrows(egmX, np.transpose(direction), mag=4, color='black')#color=[90, 90, 90]) - #here - pl.camera.azimuth += 190 - pl.camera.roll -=30 - pl.camera.elevation -=40 - pl.camera.zoom (1.2) - #pl.show(screenshot=f'/Users/aligharaviri/Desktop/CV_Data/Figures_pics/{file_name}.png') - pl.show() - #pl.close() - - if method== 'CV' and show_histogram: - plt.rcParams['font.size'] = '14' - plt.rcParams['font.family'] = 'arial' - bin_width = 0.01 - plt.hist(cv_data, bins=np.arange(np.nanmin(cv_data), np.nanmax(cv_data) + bin_width, bin_width)) - - plt.xlabel('CV (m/s)') - plt.ylabel('Number of Events') - - #plt.savefig(f'/Volumes/Simulations/CV_Data/Figures_pics/hist_{file_name}.pdf',bbox_inches='tight',transparent=True, dpi=400) - #here - plt.show() - #plt.close() - if method == 'CV' and show_cv_fibrosis: - - interpolated_mesh_cv_face = interpolated_mesh_cv.point_data_to_cell_data() - correlation_data = [] - - low_cv_mask = interpolated_mesh_cv_face['CV[m/s]'] <= 0.3 - high_cv_mask = interpolated_mesh_cv_face['CV[m/s]'] > 0.3 - fibrosis_mask = mesh['regions'] == 34 - normal_mask = mesh['regions'] != 34 - low_cv_fibrosis = np.multiply(low_cv_mask,fibrosis_mask) - high_cv_no_fibrosis = np.multiply(high_cv_mask,normal_mask) - low_cv_no_fibrosis = np.multiply(low_cv_mask,normal_mask) - high_cv_fibrosis = np.multiply(high_cv_mask,fibrosis_mask) - correlation_data = (np.multiply(low_cv_fibrosis,3) + np.multiply(high_cv_fibrosis,2) + np.multiply(low_cv_no_fibrosis,1) + np.multiply(high_cv_no_fibrosis,0)) - accuracy = np.sum(np.multiply(low_cv_fibrosis,1) + np.multiply(high_cv_no_fibrosis,1)) / len(interpolated_mesh_cv_face['CV[m/s]']) - print('Accuracy= ',accuracy) - mesh['correlation'] = correlation_data - pl = pv.Plotter(off_screen=True) - pl.add_mesh(mesh,scalars='correlation', cmap = boring_cmap_2) - pl.show(screenshot=f'/Users/aligharaviri/Desktop/CV_Data/Figures_pics/correlation_{file_name}.png') - pl.close() - cv_not_nan = interpolated_mesh_cv_face['CV[m/s]'][~np.isnan(interpolated_mesh_cv_face['CV[m/s]'])] - region_not_nan = fibrosis_mask[~np.isnan(interpolated_mesh_cv_face['CV[m/s]'])] - reg_id = np.multiply(region_not_nan,1) - elif method == 'CV' and not show_cv_fibrosis: - reg_id = [] - cv_not_nan = [] - #cv_corr = [] - #for i in range(len(mesh['regions'])): - # print(i,'end:\r') - # if mesh['regions'][i] == 34: - # if ~np.isnan(interpolated_mesh_cv_face['CV[m/s]'][i]): - # reg_id.append(1) - # cv_corr.append(interpolated_mesh_cv_face['CV[m/s]'][i]) - # else: - # if ~np.isnan(interpolated_mesh_cv_face['CV[m/s]'][i]): - # reg_id.append(0) - # cv_corr.append(interpolated_mesh_cv_face['CV[m/s]'][i]) - - - if method == 'CV': - return interpolated_mesh_cv['CV[m/s]'], reg_id, cv_not_nan - - def gradient_method(self, egmX, LAT, points, face, Fiber_data = ''): - '''#sampling_points = Points[sampling_points_ID]''' - - mesh_original = pv.PolyData(points,face) - mesh_original['LAT[ms]'] = LAT - deriv = mesh_original.compute_derivative('LAT[ms]') - df = deriv['gradient'] - d_lat_x = deriv['gradient'][:,0]#/np.sum(np.power( deriv['gradient'],2),axis=1) - d_lat_y = deriv['gradient'][:,1]#/np.sum(np.power( deriv['gradient'],2),axis=1) - d_lat_z = deriv['gradient'][:,2]#/np.sum(np.power( deriv['gradient'],2),axis=1) - #cv_data = 1 / (np.power(d_lat_x,2) + np.power(d_lat_y,2) + np.power(d_lat_z,2)) - cvx = d_lat_x/np.sum(np.power( deriv['gradient'],2),axis=1) - cvy = d_lat_y/np.sum(np.power( deriv['gradient'],2),axis=1) - cvz = d_lat_z/np.sum(np.power( deriv['gradient'],2),axis=1) - cv_data = np.sqrt(np.power(cvx,2) + np.power(cvy,2) + np.power(cvz,2))#/np.sum(np.power( deriv['gradient'],2),axis=1) - cv_data_sampled = cv_data[egmX] - - return cv_data, cv_data_sampled - - def calculate_divergence (self, egmX, LAT, points, face, sampling_data,collision_threshold=-1,focal_threshold=1,interpolation_radius=5): - boring_cmap = plt.cm.get_cmap("jet", 256) - boring_cmap_r = plt.cm.get_cmap("jet_r", 256) - egmX_point_cloud = pv.PolyData(egmX) - egmX_point_cloud['LAT[ms]'] = LAT - mesh_original = pv.PolyData(points,face) - interpolated = mesh_original.interpolate(egmX_point_cloud, radius=interpolation_radius,n_points=None) - pl = pv.Plotter() - pl.add_mesh(interpolated, cmap=boring_cmap_r) - #pl.show() - pl.close() - deriv = interpolated.compute_derivative('LAT[ms]') - - cv_x = deriv['gradient'][:,0]/np.sum(np.power( deriv['gradient'],2),axis=1) - cv_y = deriv['gradient'][:,1]/np.sum(np.power( deriv['gradient'],2),axis=1) - cv_z = deriv['gradient'][:,2]/np.sum(np.power( deriv['gradient'],2),axis=1) - - direction = np.ndarray(shape=(len(cv_x),3)) - cv_direction = np.ndarray(shape=(len(cv_x),3)) - for i in range(len(cv_x)): - direction[i][0]=cv_x[i] /np.sqrt(np.sum(np.power(cv_x[i],2) + np.power(cv_y[i],2) + np.power(cv_z[i],2))) - direction[i][1]=cv_y[i] /np.sqrt(np.sum(np.power(cv_x[i],2) + np.power(cv_y[i],2) + np.power(cv_z[i],2))) - direction[i][2]=cv_z[i] /np.sqrt(np.sum(np.power(cv_x[i],2) + np.power(cv_y[i],2) + np.power(cv_z[i],2))) - cv_direction[i][0]=cv_x[i] - cv_direction[i][1]=cv_y[i] - cv_direction[i][2]=cv_z[i] - - mesh_original['activation_direction'] = cv_direction - div = mesh_original.compute_derivative(scalars='activation_direction',divergence=True) - sargs = dict( - title_font_size=32, - label_font_size=32, - shadow=True, - n_labels=11, - italic=True, - fmt="%.1f", - font_family="arial", - color='black' - ) - pv.rcParams['transparent_background'] = True - pl = pv.Plotter(shape=(1,2)) - pl.subplot(0,0) - pl.add_mesh(div.cell_data_to_point_data(), scalars='divergence',cmap=plt.cm.get_cmap("PRGn", 9),clim=[-1.5,1.5],below_color=[63, 1, 44],above_color=[1, 50,32], - scalar_bar_args=sargs) - pl.add_arrows(points[sampling_data],direction[sampling_data],mag=2,color=[90, 90, 90]) - pl.set_background('white') - pl.camera.roll += 20 - pl.camera.elevation -=10 - pl.camera.azimuth = 10 - pl.camera.zoom (1.5) - - pl.subplot(0,1) - pl.add_mesh(div.cell_data_to_point_data(), scalars='divergence',cmap=plt.cm.get_cmap("PRGn", 9),clim=[-1.5,1.5],below_color=[63, 1, 44],above_color=[1, 50,32], - scalar_bar_args=sargs) - pl.add_arrows(points[sampling_data],direction[sampling_data],mag=2,color=[90, 90, 90]) - pl.set_background('white') - pl.camera.azimuth += 190 - pl.camera.roll -=30 - pl.camera.elevation -=40 - pl.camera.zoom (1.5) - #pl.show(screenshot='/Volumes/Simulations/CV_Data/Figures_pics/foo.png') - #pl.show() - pl.close() - - div_value = [] - for i in range(len(div['divergence'])): - if div['divergence'][i] < collision_threshold or div['divergence'][i] > focal_threshold : - div_value.append(1) - - else: - div_value.append(0) - mesh_divergence = pv.PolyData(points,face) - mesh_divergence ['div'] = div_value - divergence_per_point = mesh_divergence.cell_data_to_point_data() - pl = pv.Plotter() - pl.add_mesh(divergence_per_point,cmap=plt.cm.get_cmap("jet", 2)) - pl.add_arrows(points[sampling_data],direction[sampling_data],mag=2,color=[90, 90, 90]) - #here - #pl.show() - pl.close() - - - return direction, mesh_original, div_value - - def exclude_collision_points(self, conduction_velocity_centers = '', conduction_velocity= '', - divergence_value = '', upper_threshold = '',mesh_points = ''): - tree = KDTree(mesh_points, leaf_size=2) - ind = tree.query_radius(conduction_velocity_centers, r=0.5) - non_colision_points = [] - collision_detection = [] - cv_data = [] - for i in range(len(ind)): - print(i,end='\r') - for k in range(len(ind[i])): - if divergence_value[ind[i][k]] == 1: - collision_detection.append(1) - else: - collision_detection.append(0) - - if np.sum(collision_detection) or (conduction_velocity[i] > upper_threshold): - non_colision_points.append(False) - cv_data.append(np.nan) - else: - non_colision_points.append(True) - cv_data.append(conduction_velocity[i]) - collision_detection = [] - - #bin_width = 0.01 - #plt.hist(cv_data, bins=np.arange(np.nanmin(cv_data), np.nanmax(cv_data) + bin_width, bin_width)) - - #plt.xlabel('CV (m/s)') - #plt.ylabel('Number of Events') - #here - #plt.show() - return non_colision_points, cv_data - #pl = pv.Plotter() - #pl.add_mesh(interpolated,cmap = boring_cmap_r) - #pl.add_arrows(points,direction,mag=2,color=[90, 90, 90]) - #pl.show() - - - - - #mesh_original['divergence'] = div_value - #pl.add_mesh(mesh_original, scalars='divergence',cmap=plt.cm.get_cmap("jet", 2)) - - diff --git a/openep/Analyses/run.py b/openep/Analyses/run.py deleted file mode 100644 index 22953c23..00000000 --- a/openep/Analyses/run.py +++ /dev/null @@ -1,72 +0,0 @@ -import openep -from CV import CV -import pyvista as pv -import numpy as np -def run_cv (openep_file_name,cv_method, visualsiation): - case = openep.load_openep_mat(openep_file_name) - egmX = case.electric.bipolar_egm.points - voltage = case.electric.bipolar_egm.voltage - print('---> read electrode data') - LAT = case.electric.annotations.local_activation_time - case.electric.annotations.reference_activation_time - print('---> read bi-polar LAT data') - surface = case.create_mesh() - - valid_id = np.where(LAT!=-10000) - invalid_id = np.where(LAT==-10000) - egmX_clean = np.delete(egmX, np.where(LAT == -10000),0) - LAT_clean = np.delete(LAT, np.where(LAT == -10000)) - Points = surface.points - faces = surface.faces - Face = faces.reshape(-1,4) - temp_face = Face[:,1:4] - mesh_original = pv.PolyData(Points,Face) - cal_cv = CV() - if cv_method == 'Plane': - cv_plane, cv_centers_plane, cv_centers_id_plane = cal_cv.plane_fitting(egmX_clean, LAT_clean) - cv_total = np.zeros((len(egmX),1)) - for i in range(len(cv_plane)): - cv_total[valid_id[0][i]] = cv_plane[i] - - for i in range(len(invalid_id)): - cv_total[invalid_id[i]] = np.nan - - if visulaisation: - cal_cv.visualization(lat_data = '', cv_data = cv_plane, mesh = mesh_original, egmX = egmX_clean, method = 'CV',file_name = '', - show_histogram= True, direction_centers = Points,show_cv_fibrosis=False) - cv = cv_plane - - if cv_method == 'RBF': - cv_rbf, cv_centers_rbf, direction_rbf, cv_centers_id_rbf = cal_cv.RBF_method(egmX_clean, LAT_clean,Points, Face) - if visulaisation: - cal_cv.visualization(lat_data = '', cv_data = cv_rbf, mesh = mesh_original, egmX = np.transpose(cv_centers_rbf), method = 'CV',file_name = '', - show_histogram= True, direction_centers = Points,show_cv_fibrosis=False) - cv = cv_rbf - - if cv_method == 'Tri': - cv_tri, cv_centers_tri, cv_centers_id_tri = cal_cv.Triangulation(egmX_clean,LAT_clean) - if visulaisation: - cal_cv.visualization(lat_data = '', cv_data = cv_tri, mesh = mesh_original, egmX = np.transpose(cv_centers_tri), method = 'CV',file_name = '', - show_histogram= True, direction_centers = Points,show_cv_fibrosis=False) - cv = cv_tri - ######## Save the CV data as openep ####### - case.fields['conduction_velocity'] = cv - case.fields['mesh_normals'] = mesh_original.point_normals - - return case - - - - - - - -if __name__=="__main__": - file_path = '/Users/aligharaviri/Downloads/LAWT_CartoAF_ EP_Studies/' - file_name = '521.mat' - openep_file_name = f"{file_path}{file_name}" - out_put_file_name = f"{file_name[:-4]}_cv.mat" - method = "Tri" ######## methods are Plane, RBF, and Tri - visulaisation = True - case = run_cv (openep_file_name,method, visulaisation) - openep.export_openep_mat(case,out_put_file_name) - diff --git a/openep/analysis/conduction_velocity.py b/openep/analysis/conduction_velocity.py deleted file mode 100644 index 60e56556..00000000 --- a/openep/analysis/conduction_velocity.py +++ /dev/null @@ -1,302 +0,0 @@ -import openep -import sys -import os -import math -import pyvista as pv -import numpy as np -import scipy.io -import scipy -from scipy.optimize import minimize -import matplotlib.pyplot as plt -from sklearn.neighbors import KDTree -import pdb -from scipy.interpolate import Rbf -from vedo import * -import pyvistaqt -from sklearn import metrics - -from ..case.case_routines import ( - interpolate_general_cloud_points_onto_surface -) - -class ConductionVelocity(): - - - def __init__(self, case): - self.case = case - self.include = None - self.values = None - self.points = None - self.prep_egmX, self.prep_LAT = self._preprocess_data() - - def calculate( - self, - method='triangulation', - method_kws=dict(), - include=None - ): - - """Calculate conduction velocity and interpolates each point on a mesh, - stores in conduction_velocity fields. - - Args: - method (str): method to use for interpolation. This should be one of: - - 'triangualtion' - - 'plane_fitting' - - 'rbf' - """ - if include is not None: - self.include=include - self.prep_egmX, self.prep_LAT = self._preprocess_data() - - NAME_TO_FUNC = { - "triangulation": self.calculate_with_triangualtion, - "plane_fitting": self.calculate_with_plane_fitting, - "rbf": self.calculate_with_rbf, - } - - method = method.lower() - if method not in NAME_TO_FUNC: - raise ValueError(f"`method` must be one of {NAME_TO_FUNC.keys()}.") - - cv_calc_method = NAME_TO_FUNC[method] - out = cv_calc_method(method_kws) - - self.case.fields.conduction_velocity = interpolate_general_cloud_points_onto_surface( - case=self.case, - cloud_values=self.case.analyse.conduction_velocity.values, - cloud_points=self.case.analyse.conduction_velocity.points - ) - return out - - def _preprocess_data(self): - if self.include is None: - egmX = self.case.electric.bipolar_egm.points - #TODO: voltage is not even used! - voltage = self.case.electric.bipolar_egm.voltage - LAT = self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time - else: - egmX = self.case.electric.bipolar_egm.points[self.include] - voltage = self.case.electric.bipolar_egm.voltage[self.include] - LAT = self.case.electric.annotations.local_activation_time[self.include] - self.case.electric.annotations.reference_activation_time[self.include] - - - print('len_cv:', len(egmX)) - surface = self.case.create_mesh() - - valid_id, invalid_id = np.where(LAT != -10000), np.where(LAT == -10000) - egmX_clean = np.delete(egmX, np.where(LAT == -10000), 0) - LAT_clean = np.delete(LAT, np.where(LAT == -10000)) - - return egmX_clean, LAT_clean - - def calculate_with_triangualtion(self, method_kws=dict()): - - print('---> Starting triangulation') - egmX, LAT = self.prep_egmX, self.prep_LAT - cv = [] - min_theta = method_kws.get('min_theta', 30) - max_elec_distance = method_kws.get('max_elec_dist', 10) - min_elec_distance = method_kws.get('min_elec_dist', 1.5) - min_lat_difference = method_kws.get('min_lat_diff', 2) - cv_centers_x = [] - cv_centers_y = [] - cv_centers_z = [] - cv_centers_id = [] - alpha = method_kws.get('alpha', 5) - - egmX_point_cloud = pv.PolyData(egmX) - print('---> Triangulation process') - surf = egmX_point_cloud.delaunay_3d(alpha=alpha) - print('---> Triangulation process') - print('---> CV calculation (Triangulation method started ...)') - - for num_point in range(len(surf.cells_dict[5])): - - vtx_id = [] - lat = [] - - for i in range(3): - vtx_id.append(surf.cells_dict[5][num_point][i]) - lat.append(LAT[surf.cells_dict[5][num_point][i]]) - - id_lat_sorted = np.argsort(lat) - - O = [egmX[int(vtx_id[id_lat_sorted[0]])], lat[id_lat_sorted[0]]] - A = [egmX[int(vtx_id[id_lat_sorted[1]])], lat[id_lat_sorted[1]]] - B = [egmX[int(vtx_id[id_lat_sorted[2]])], lat[id_lat_sorted[2]]] - - OA = np.sqrt(sum(np.power(np.subtract(O[0], A[0]), 2))) - OB = np.sqrt(sum(np.power(np.subtract(O[0], B[0]), 2))) - AB = np.sqrt(sum(np.power(np.subtract(A[0], B[0]), 2))) - - tOA = A[1] - O[1] - tOB = B[1] - O[1] - - theta = np.arccos((np.power(OA, 2) + np.power(OB, 2) - np.power(AB, 2)) / (2 * OA * OB)) - # check if the conditions are meet to accept the triangle set as a viable acceptable one - if (math.degrees(theta) >= min_theta and OA >= min_elec_distance and OA <= max_elec_distance and - OB >= min_elec_distance and OB <= max_elec_distance and tOA >= min_lat_difference and tOB >= min_lat_difference): - - alpha = np.arctan((tOB * OA - tOA * OB * np.cos(theta)) / (tOA * OB * np.sin(theta))) - cv_temp = (OA / tOA) * np.cos(alpha) - cv.append(cv_temp) - cv_centers_x.append((O[0][0] + A[0][0] + B[0][0]) / 3) - cv_centers_y.append((O[0][1] + A[0][1] + B[0][1]) / 3) - cv_centers_z.append((O[0][2] + A[0][2] + B[0][2]) / 3) - cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - else: - cv.append(np.nan) - cv_centers_x.append((O[0][0] + A[0][0] + B[0][0]) / 3) - cv_centers_y.append((O[0][1] + A[0][1] + B[0][1]) / 3) - cv_centers_z.append((O[0][2] + A[0][2] + B[0][2]) / 3) - cv_centers_id.append(vtx_id[id_lat_sorted[0]]) - ####### Create a triangulation mesh from egmX pointcloud ##### - cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] - print('---> CV calculation ended') - print(cv) - - egmX_point_cloud = pv.PolyData(np.transpose(cv_centers)) - egmX_point_cloud['CV[m/s]'] = cv - # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) - - - egmX_point_cloud = pv.PolyData(np.transpose(cv_centers)) - egmX_point_cloud['CV[m/s]'] = cv - - - # interpolated_cv_field = interpolate_conduction_velocity_onto_surface( - # self.case, - # np.array(cv), - # egmX_point_cloud - # ) - # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) - - self.values = np.array(cv) - self.points = egmX_point_cloud.points - - return np.array(cv), egmX_point_cloud.points - - def calculate_with_plane_fitting(self, method_kws=dict()): - egmX, LAT = self.prep_egmX, self.prep_LAT - - cv = [] - cv_centers_x = [] - cv_centers_y = [] - cv_centers_z = [] - cv_centers_id = [] - cent = np.random.random((len(egmX), 3)) - direction = np.random.random((len(egmX), 3)) - - leaf_size = method_kws.get('leaf_size', 5) - tree = KDTree(egmX, leaf_size=leaf_size) - all_nn_indices = tree.query_radius(egmX, r=10) # 10 closest neighnours - all_nns = [[egmX[idx] for idx in nn_indices] for nn_indices in all_nn_indices] - - min_n_nearest_neighbours = method_kws.get('min_n_nearest_neighbours', 5) - print(len(egmX)) - for k in range(len(egmX)): - print(k, end='\r') - if len(all_nns[k]) >= min_n_nearest_neighbours: - cv_data = np.zeros(shape=(len(all_nns[k]) + 1, 4)) - cv_data[0][:3] = egmX[k] - cv_data[0][3] = LAT[k] - for i in range(len(all_nns[k])): - cv_data[i + 1][:3] = all_nns[k][i] - cv_data[i + 1][3] = LAT[all_nn_indices[k][i]] - - x1 = cv_data[:, 0] - x2 = cv_data[:, 1] - x3 = cv_data[:, 2] - m = cv_data[:, 3] - - res = minimize(self.func, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], args=(x1, x2, x3, m)) - a, b, c, d, e, f, g, h, q, l = res.x - - x = egmX[k, 0] - y = egmX[k, 1] - z = egmX[k, 2] - - tx = 2 * a * x + d * y + e * z + g - ty = 2 * b * y + d * x + f * z + h - tz = 2 * c * z + e * x + f * y + q - - vx = tx / (tx ** 2 + ty ** 2 + tz ** 2) - vy = ty / (tx ** 2 + ty ** 2 + tz ** 2) - vz = tz / (tx ** 2 + ty ** 2 + tz ** 2) - - cent[k, 0] = x - cent[k, 1] = y - cent[k, 2] = z - - direction[k, 0] = vx - direction[k, 1] = vy - direction[k, 2] = vz - - cv_temp = np.sqrt(np.sum(np.power(vx, 2) + np.power(vy, 2) + np.power(vz, 2))) - direction[k, 0] = vx - direction[k, 1] = vy - direction[k, 2] = vz - - if cv_temp >= 0.2 and cv_temp <= 2: - cv.append(cv_temp) - cv_centers_x.append(x) - cv_centers_y.append(y) - cv_centers_z.append(z) - cv_centers_id.append(k) - else: - - cv.append(cv_temp) - cv_centers_x.append(x) - cv_centers_y.append(y) - cv_centers_z.append(z) - else: - x = egmX[k, 0] - y = egmX[k, 1] - z = egmX[k, 2] - cv.append(np.nan) - cv_centers_x.append(x) - cv_centers_y.append(y) - cv_centers_z.append(z) - cv_centers_id.append(k) - cv_centers = [cv_centers_x, cv_centers_y, cv_centers_z] - - egmX_point_cloud = pv.PolyData(np.transpose(cv_centers)) - egmX_point_cloud['CV[m/s]'] = cv - - # interpolated_cv_field = interpolate_conduction_velocity_onto_surface( - # self.case, - # np.array(cv), - # egmX_point_cloud - # ) - # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) - - self.values = np.array(cv) - self.points = egmX_point_cloud.points - - return np.array(cv), egmX_point_cloud.points - - - # interpolated_mesh_cv = self.mesh_original.interpolate(egmX_point_cloud, radius=10, n_points=None, null_value=np.nan) - # return interpolated_mesh_cv.point_data.active_scalars - - def func(self, coef, x1, x2, x3, m): - ''' This the function we use to fit a surface to the activations in plane - fitting method''' - a = coef[0] - b = coef[1] - c = coef[2] - d = coef[3] - e = coef[4] - f = coef[5] - g = coef[6] - h = coef[7] - q = coef[8] - l = coef[9] - # print(np.linalg.norm(m-(a*x1**2 + b*x2**2 + c*x3**2 + d*x1*x2 + e*x1*x3 + f*x2*x3 + g*x1 +h*x2 + q*x3 + l))) - norm = np.linalg.norm(m - ( - a * x1 ** 2 + b * x2 ** 2 + c * x3 ** 2 + d * x1 * x2 + e * x1 * x3 + f * x2 * x3 + g * x1 + h * x2 + q * x3 + l)) - return norm - - def calculate_with_rbf(self): - pass \ No newline at end of file diff --git a/openep/analysis/divergence.py b/openep/analysis/divergence.py deleted file mode 100644 index 5b36b577..00000000 --- a/openep/analysis/divergence.py +++ /dev/null @@ -1,135 +0,0 @@ -import openep -import sys -import os -import math -import pyvista as pv -import numpy as np -import scipy.io -import scipy -from scipy.optimize import minimize -import matplotlib.pyplot as plt -from sklearn.neighbors import KDTree -import pdb -from scipy.interpolate import Rbf -from vedo import * -import pyvistaqt -from sklearn import metrics - -from ..case.case_routines import ( - interpolate_general_cloud_points_onto_surface, - interpolate_activation_time_onto_surface -) - -class Divergence(): - - - def __init__(self, case): - self.case = case #.copy() change to copy later - self.include = None - self.direction = None - self.mesh_original = None - self.values = None - self.divergence_per_point=None - - def calculate( - self, - output_binary_field=False, - method_kws=dict(), - include=None, - ): - """Calculate conduction velocity and interpolates each point on a mesh, - stores in conduction_velocity fields. - - Args: - method (str): method to use for interpolation. This should be one of: - - 'triangualtion' - - 'plane_fitting' - - 'rbf' - """ - if include is not None: - self.include=include - - out = self.calculate_divergence(method_kws) - - - def calculate_divergence(self, method_kws, output_binary_field=False): - - if output_binary_field: - collision_threshold = method_kws.get("collision_threshold",-1) - focal_threshold = method_kws.get("collision_threshold", 1) - - surface = self.case.create_mesh() - - if self.include is None: - egmX = self.case.electric.bipolar_egm.points - LAT = self.case.electric.annotations.local_activation_time - self.case.electric.annotations.reference_activation_time - else: - egmX=self.case.electric.bipolar_egm.points[self.include] - LAT=self.case.electric.annotations.local_activation_time[self.include] - self.case.electric.annotations.reference_activation_time[self.include] - print("len_div:", len(egmX)) - - points=surface.points - face=surface.faces - - egmX_point_cloud = pv.PolyData(egmX) - egmX_point_cloud['LAT[ms]'] = LAT - mesh_original = pv.PolyData(points, face) - - # interpolated_old = mesh_original.interpolate(egmX_point_cloud, radius=5, n_points=None) - interpolated_scalar=interpolate_general_cloud_points_onto_surface( - case=self.case, - cloud_values=LAT, - cloud_points=egmX, - ) - - mesh_original['LAT_scalar'] = interpolated_scalar - deriv = mesh_original.compute_derivative(scalars='LAT_scalar') - - cv_x = deriv['gradient'][:, 0] / np.sum(np.power(deriv['gradient'], 2), axis=1) - cv_y = deriv['gradient'][:, 1] / np.sum(np.power(deriv['gradient'], 2), axis=1) - cv_z = deriv['gradient'][:, 2] / np.sum(np.power(deriv['gradient'], 2), axis=1) - - direction = np.ndarray(shape=(len(cv_x), 3)) - cv_direction = np.ndarray(shape=(len(cv_x), 3)) - for i in range(len(cv_x)): - direction[i][0] = cv_x[i] / np.sqrt( - np.sum(np.power(cv_x[i], 2) + np.power(cv_y[i], 2) + np.power(cv_z[i], 2))) - direction[i][1] = cv_y[i] / np.sqrt( - np.sum(np.power(cv_x[i], 2) + np.power(cv_y[i], 2) + np.power(cv_z[i], 2))) - direction[i][2] = cv_z[i] / np.sqrt( - np.sum(np.power(cv_x[i], 2) + np.power(cv_y[i], 2) + np.power(cv_z[i], 2))) - cv_direction[i][0] = cv_x[i] - cv_direction[i][1] = cv_y[i] - cv_direction[i][2] = cv_z[i] - - mesh_original['activation_direction'] = cv_direction - div = mesh_original.compute_derivative(scalars='activation_direction', divergence=True) - sargs = dict( - title_font_size=32, - label_font_size=32, - shadow=True, - n_labels=11, - italic=True, - fmt="%.1f", - font_family="arial", - color='black' - ) - - if output_binary_field: - div_value = [] - for i in range(len(div['divergence'])): - if div['divergence'][i] < collision_threshold or div['divergence'][i] > focal_threshold: - div_value.append(1) - else: - div_value.append(0) - else: - div_value = div['divergence'] - - mesh_divergence = pv.PolyData(points, face) - mesh_divergence['div'] = div_value - divergence_per_point = mesh_divergence.cell_data_to_point_data() - - self.case.fields.divergence = div_value - self.divergence_per_point=divergence_per_point - - self.direction, self.mesh_original, self.values = direction, mesh_original, div_value From f5b2f329615a1897f652a02e1605779223c26bd9 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 10:41:24 +0000 Subject: [PATCH 19/31] all CV and divergence methods in one python file --- openep/analysis/_conduction_velocity.py | 320 ++++++++++++++++++++++++ 1 file changed, 320 insertions(+) create mode 100644 openep/analysis/_conduction_velocity.py diff --git a/openep/analysis/_conduction_velocity.py b/openep/analysis/_conduction_velocity.py new file mode 100644 index 00000000..aca6654e --- /dev/null +++ b/openep/analysis/_conduction_velocity.py @@ -0,0 +1,320 @@ +# OpenEP +# Copyright (c) 2021 OpenEP Collaborators +# +# This file is part of OpenEP. +# +# OpenEP 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. +# +# OpenEP 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. +# +# You should have received a copy of the GNU General Public License along +# with this program (LICENSE.txt). If not, see + +"""Module containing conduction velocity (CV) calculation methods and CV divergence method""" + +import math +from vedo import * +import pyvista as pv +from scipy.optimize import minimize +from sklearn.neighbors import KDTree + +from ..case.case_routines import interpolate_general_cloud_points_onto_surface + + +__all__ = ['preprocess_lat_egm', 'plane_fitting', 'divergence', + 'triangulation', 'radial_basis_function'] + + +def preprocess_lat_egm( + case, + include=None, + lat_threshold=-10000 +): + """ + Preprocess case data for conduction velocity calculation by removing local activation times + and corresponding bipolar electrogram points based on threshold and include. + + Args: + case (Case): The case data containing electrogram and annotation information. + + include (np.ndarray, optional): A boolean array specifying which data points to include + in the preprocessing for bipolar EGM data and LAT. + Defaults to "includes" specified in the imported data. + + lat_threshold (int, optional): Threshold for local activation times. Points with local + activation times less than this value will be removed. + Defaults to -1000. + If None, no threshold is applied. + + Returns: + tuple: A tuple containing two numpy arrays: + - local_activation_time (np.ndarray): Array of cleaned local activation times. + - bipolar_egm_pts (np.ndarray): Array of bipolar electrogram points corresponding + to the cleaned local activation times. + + """ + bipolar_egm_pts = case.electric.bipolar_egm.points[include] + local_activation_time = (case.electric.annotations.local_activation_time[include] + - case.electric.annotations.reference_activation_time[include]) + + if lat_threshold: + bipolar_egm_pts = np.delete(bipolar_egm_pts, np.where(local_activation_time == lat_threshold), 0) + local_activation_time = np.delete(local_activation_time, np.where(local_activation_time == lat_threshold)) + + return local_activation_time, bipolar_egm_pts + + +def plane_fitting( + bipolar_egm_pts, + local_activation_time, + leaf_size=5, + min_n_nearest_neighbours=5, + tree_query_radius=10 +): + """ + Fit a plane to bipolar electrogram points and calculate conduction velocity values and centroids. + + Args: + bipolar_egm_pts (array): + Array of bipolar electrogram points. + + local_activation_time (array): + Array of local activation times. + + leaf_size (int, optional): + Leaf size for KDTree. Defaults to 5. + + min_n_nearest_neighbours (int, optional): + Minimum number of nearest neighbours. Defaults to 5. + + tree_query_radius (int, optional): + Radius for tree query. Defaults to 10. + + Returns: + tuple: Array of conduction velocities and array of triangle center points. + + Example 1: + import openep + + case = openep.load_openep_mat("path/to/openep.mat") + mesh = case.create_mesh() + case.analyse.conduction_velocity.calculate(method='plane_fitting') + + Example 2: + >> import numpy as np + >> bipolar_egm_pts = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], ...]) + >> local_activation_times = np.array([0.1, 0.2, 0.15, ...]) + >> cv_values, cv_centroids = plane_fitting(bipolar_egm_pts, local_activation_times) + + """ + cv_values = np.full(bipolar_egm_pts.shape[0], np.nan) + cv_centroids = np.empty((bipolar_egm_pts.shape)) + ones_arr = np.ones(10) + + tree = KDTree(bipolar_egm_pts, leaf_size=leaf_size) + all_nn_indices = tree.query_radius(bipolar_egm_pts, r=tree_query_radius) + all_nns = [[bipolar_egm_pts[idx] for idx in nn_indices] for nn_indices in all_nn_indices] + + for egm_i in range(bipolar_egm_pts.shape[0]): + + if len(all_nns[egm_i]) >= min_n_nearest_neighbours: + x1, x2, x3 = np.transpose(np.vstack((bipolar_egm_pts[egm_i], all_nns[egm_i]))) + m = np.append(local_activation_time[egm_i], local_activation_time[all_nn_indices[egm_i]]) + + res = minimize(_obj_func, ones_arr, args=(x1, x2, x3, m)) + a, b, c, d, e, f, g, h, q, l = res.x + + x, y, z = bipolar_egm_pts[egm_i] + + tx = 2 * a * x + d * y + e * z + g + ty = 2 * b * y + d * x + f * z + h + tz = 2 * c * z + e * x + f * y + q + + vx = tx / (tx ** 2 + ty ** 2 + tz ** 2) + vy = ty / (tx ** 2 + ty ** 2 + tz ** 2) + vz = tz / (tx ** 2 + ty ** 2 + tz ** 2) + + cv_values[egm_i] = np.sqrt(vx ** 2 + vy ** 2 + vz ** 2) + cv_centroids[egm_i] = [x, y, z] + else: + cv_centroids[egm_i] = bipolar_egm_pts[egm_i] + + return cv_values, cv_centroids + + +def triangulation( + bipolar_egm_pts, + local_activation_time, + min_theta=30, + electrode_distance_range=(1.5, 10), + min_lat_difference=2, + del_alpha=5 +): + """ + Perform triangulation on electrogram data to calculate conduction velocities. + + Args: + bipolar_egm_pts (array-like): + Array of bipolar electrogram points of size Nx3 + + local_activation_time (array-like): + Local activation times corresponding to bipolar_egm_pts. + + min_theta (float): + Minimum angle threshold for triangle cell acceptance (degrees). + + electrode_distance_range (tuple): + Range of acceptable electrode distances. + + min_lat_difference (float): + Minimum local activation time difference. + + del_alpha (float): + Alpha value for Delaunay triangulation. + + Returns: + tuple: Array of conduction velocities and array of triangle center points. + + Example: + import openep + + case = openep.load_openep_mat("path/to/openep.mat") + mesh = case.create_mesh() + case.analyse.conduction_velocity.calculate(method='triangulation') + + Reference: + This method is based on the paper: + Age-Related Changes in Human Left and Right Atrial Conduction by PIPIN KOJODJOJO M.R.C.P. + """ + + bi_egm_point_cloud = pv.PolyData(bipolar_egm_pts) + bi_egm_surface = bi_egm_point_cloud.delaunay_3d(alpha=del_alpha) + + # Only extract triangle cell types from mesh, see pyvista.CellType.TRIANGLE (indexed at 5 in cells_dict) + bi_egm_cells = bi_egm_surface.cells_dict[5] + + cv_values = np.full(bi_egm_cells.shape[0], np.nan) + cv_centroids = np.empty(bi_egm_cells.shape) + + for cell_i in range(bi_egm_cells.shape[0]): + + # Find vtx id of cell and sort by LAT value + vtx_id = bi_egm_cells[cell_i].astype(int) + lat = local_activation_time[vtx_id] + lat_sorted_id = np.argsort(lat) + + # Assign point O, A and B to the cell vtx + lat_O, lat_A, lat_B = lat[lat_sorted_id] + point_O, point_A, point_B = bipolar_egm_pts[vtx_id[lat_sorted_id]] + + # Calculate distance between OA, OB and AB + dist_OA = np.linalg.norm(np.subtract(point_O, point_A)) + dist_OB = np.linalg.norm(np.subtract(point_O, point_B)) + dist_AB = np.linalg.norm(np.subtract(point_A, point_B)) + + # Calculate difference in LAT between OA and OB + tOA = lat_A - lat_O + tOB = lat_B - lat_O + + # Calculate angles subtended by OA and OB + theta = np.arccos((dist_OA**2 + dist_OB**2 - dist_AB**2) / (2 * dist_OA * dist_OB)) + + # Conditions to check if triangle is acceptable + if (math.degrees(theta) >= min_theta + and electrode_distance_range[0] <= dist_OA <= electrode_distance_range[1] + and electrode_distance_range[0] <= dist_OB <= electrode_distance_range[1] + and tOA >= min_lat_difference + and tOB >= min_lat_difference + ): + # Calculate angle subtended by OA and the direction of wavefront propagation + alpha = np.arctan((tOB * dist_OA - tOA * dist_OB * np.cos(theta)) / (tOA * dist_OB * np.sin(theta))) + cv_temp = (dist_OA / tOA) * np.cos(alpha) + cv_values[cell_i] = cv_temp + else: + cv_values[cell_i] = np.nan + + centroid_i = np.mean([point_O, point_A, point_B], axis=0) + cv_centroids[cell_i] = centroid_i + + # return np.array(cv_values), np.array(cv_centroids) + return cv_values, cv_centroids + + +def _obj_func(coef, x1, x2, x3, m): + """ + Objective function for plane fitting method. + + This function is used to fit a surface to the activations in the plane fitting method. + It calculates the norm of the difference between the measured values 'm' and the + values predicted by the quadratic surface defined by the coefficients. + + Parameters: + coef (array): Coefficients of the quadratic surface (a, b, c, d, e, f, g, h, q, l). + x1 (array): X-coordinates of the data points. + x2 (array): Y-coordinates of the data points. + x3 (array): Z-coordinates of the data points. + m (array): Measured values at the data points. + + Returns: + float: The norm of the difference between measured and predicted values. + """ + a, b, c, d, e, f, g, h, q, l = coef + predicted = (a * x1 ** 2 + + b * x2 ** 2 + + c * x3 ** 2 + + d * x1 * x2 + + e * x1 * x3 + + f * x2 * x3 + + g * x1 + + h * x2 + + q * x3 + + l) + return np.linalg.norm(m - predicted) + + +def radial_basis_function( + bipolar_egm_pts, + local_activation_time +): + raise NotImplementedError("Method: radial_basis_function has not been implemented yet for estimating CV!") + + +def divergence( + case, + bipolar_egm_pts, + local_activation_time, + collision_threshold=-1, + focal_threshold=1, + output_binary_field=False +): + + case = case.copy() + original_mesh = case.create_mesh() + + interpolated_scalar = interpolate_general_cloud_points_onto_surface( + case=case, + cloud_values=local_activation_time, + cloud_points=bipolar_egm_pts, + ) + + original_mesh['LAT_scalar'] = interpolated_scalar + derivative = original_mesh.compute_derivative(scalars='LAT_scalar') + + cv_direction = derivative['gradient'] / np.sum(derivative['gradient'] ** 2, axis=1)[:, np.newaxis] + magnitude = np.sqrt(np.sum(cv_direction ** 2, axis=1)) + norm_cv_direction = cv_direction / magnitude[:, np.newaxis] + + original_mesh['activation_direction'] = cv_direction + div = original_mesh.compute_derivative(scalars='activation_direction', divergence=True) + divergence = div['divergence'] + + if output_binary_field: + divergence = np.where((divergence < collision_threshold) | (divergence > focal_threshold), 1, 0) + + return norm_cv_direction, divergence From a54b3f3cf72e21037754627cbe2126bd6e4b92b8 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 15:29:51 +0000 Subject: [PATCH 20/31] case copy error - changed to temp mesh --- openep/analysis/_conduction_velocity.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/openep/analysis/_conduction_velocity.py b/openep/analysis/_conduction_velocity.py index aca6654e..117ff16d 100644 --- a/openep/analysis/_conduction_velocity.py +++ b/openep/analysis/_conduction_velocity.py @@ -294,8 +294,8 @@ def divergence( output_binary_field=False ): - case = case.copy() - original_mesh = case.create_mesh() + temp_mesh = case.create_mesh() + basic_mesh = pv.PolyData(temp_mesh.points, temp_mesh.faces) interpolated_scalar = interpolate_general_cloud_points_onto_surface( case=case, @@ -303,15 +303,15 @@ def divergence( cloud_points=bipolar_egm_pts, ) - original_mesh['LAT_scalar'] = interpolated_scalar - derivative = original_mesh.compute_derivative(scalars='LAT_scalar') + basic_mesh['LAT_scalar'] = interpolated_scalar + derivative = basic_mesh.compute_derivative(scalars='LAT_scalar') cv_direction = derivative['gradient'] / np.sum(derivative['gradient'] ** 2, axis=1)[:, np.newaxis] magnitude = np.sqrt(np.sum(cv_direction ** 2, axis=1)) norm_cv_direction = cv_direction / magnitude[:, np.newaxis] - original_mesh['activation_direction'] = cv_direction - div = original_mesh.compute_derivative(scalars='activation_direction', divergence=True) + basic_mesh['activation_direction'] = cv_direction + div = basic_mesh.compute_derivative(scalars='activation_direction', divergence=True) divergence = div['divergence'] if output_binary_field: From 608df711d16463f305b186632e1740af8e8a95fe Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 15:30:36 +0000 Subject: [PATCH 21/31] array evaluation for None changed --- openep/analysis/analyse.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py index d4258515..6498ad5e 100644 --- a/openep/analysis/analyse.py +++ b/openep/analysis/analyse.py @@ -64,14 +64,14 @@ def __init__(self, case): @property def values(self): - if not self._values: + if self._values is None: raise ValueError('Before accessing ``conduction_velocity.values`` ' 'run ``divergence.calculate_divergence()``') return self._values @property def points(self): - if not self._points: + if self._points is None: raise ValueError('Before accessing ``conduction_velocity.points`` ' 'run ``divergence.calculate_divergence()``') return self._points @@ -166,13 +166,13 @@ def __init__(self, case): @property def values(self): - if not self._values: + if self._values is None: raise ValueError('Before accessing ``..divergence.values`` run ``..divergence.calculate_divergence()``') return self._values @property def direction(self): - if not self._direction: + if self._direction is None: raise ValueError('Before accessing ``divergence.direction`` run ``divergence.calculate_divergence()``') return self._direction From 072234435879de6467e2943e772d4f0f5d3600ad Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 15:32:28 +0000 Subject: [PATCH 22/31] divergence -> cv_divergence --- openep/analysis/analyse.py | 2 +- openep/data_structures/surface.py | 7 ++++--- openep/io/writers.py | 2 +- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py index 6498ad5e..0d3adb00 100644 --- a/openep/analysis/analyse.py +++ b/openep/analysis/analyse.py @@ -222,6 +222,6 @@ def calculate_divergence( ) if apply_scalar_field: - self._case.fields.divergence = self.values + self._case.fields.cv_divergence = self.values return self.direction, self.values diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index 5b262c22..8fe633a9 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -52,7 +52,7 @@ class Fields: transverse_fibres: np.ndarray = None pacing_site: np.ndarray = None conduction_velocity: np.ndarray = None - divergence: np.ndarray = None + cv_divergence: np.ndarray = None mesh_normals : np.ndarray = None def __repr__(self): @@ -215,8 +215,9 @@ def extract_surface_data(surface_data): longitudinal_fibres=longitudinal_fibres, transverse_fibres=transverse_fibres, pacing_site=pacing_site, - conduction_velocity = conduction_velocity, - mesh_normals = mesh_normals + conduction_velocity=conduction_velocity, + cv_divergence=cv_divergence, + mesh_normals=mesh_normals ) return points, indices, fields diff --git a/openep/io/writers.py b/openep/io/writers.py index ca76ae71..2fe36b5f 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -192,7 +192,7 @@ def export_openep_mat( userdata['surface'] = _add_surface_maps( surface_data=userdata['surface'], cv_field=case.fields.conduction_velocity, - divergence_field=case.fields.divergence + divergence_field=case.fields.cv_divergence ) userdata['electric'] = _extract_electric_data(electric=case.electric) From c04050ee2972d6e36c6af9386ee7403beffde72d Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 15:32:43 +0000 Subject: [PATCH 23/31] divergence -> cv_divergence --- openep/converters/pyvista_converters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openep/converters/pyvista_converters.py b/openep/converters/pyvista_converters.py index 18f8501e..e670b8aa 100644 --- a/openep/converters/pyvista_converters.py +++ b/openep/converters/pyvista_converters.py @@ -94,6 +94,6 @@ def to_pyvista( mesh.point_data["Impedance"] = case.fields.impedance mesh.point_data["Force"] = case.fields.force mesh.point_data["Conduction Velocity"] = case.fields.conduction_velocity - mesh.point_data["Divergence"] = case.fields.divergence + mesh.point_data["Divergence"] = case.fields.cv_divergence return mesh From 60c4b749bd9eb368aa8673b3808859b34960f0b2 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 15:33:24 +0000 Subject: [PATCH 24/31] cv + divergence import and export fixed --- openep/data_structures/surface.py | 11 +++++++---- openep/io/writers.py | 17 ++++++++--------- 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index 8fe633a9..bb9988d1 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -195,10 +195,15 @@ def extract_surface_data(surface_data): pacing_site = None try: - conduction_velocity = surface_data['conduction_velocity'].astype(float) + conduction_velocity = surface_data['signalMaps']['conduction_velocity_field']['value'].astype(float) except KeyError as e: conduction_velocity = None - + + try: + cv_divergence = surface_data['signalMaps']['divergence_field']['value'].astype(float) + except KeyError as e: + cv_divergence = None + try: mesh_normals = surface_data['mesh_normals'].astype(float) except KeyError as e: @@ -253,8 +258,6 @@ def empty_fields(n_points=0, n_cells=0): longitudinal_fibres, transverse_fibres, pacing_site, - conduction_velocity, - mesh_normals, ) return fields diff --git a/openep/io/writers.py b/openep/io/writers.py index 2fe36b5f..7b503371 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -211,6 +211,7 @@ def export_openep_mat( oned_as='column', ) + def _add_surface_maps(surface_data, **kwargs): cv_field = kwargs.get('cv_field') div_field = kwargs.get('divergence_field') @@ -218,24 +219,24 @@ def _add_surface_maps(surface_data, **kwargs): if not surface_data.get('signalMaps'): surface_data['signalMaps'] = {} - # TODO: connect propSetting from user setting - # TODO: handle with None values as matlab cannot handle + # TODO: connect propSetting from user setting in UI if cv_field is not None: surface_data['signalMaps']['conduction_velocity_field'] = { - 'name' : 'Conduction Velocity Field', + 'name': 'Conduction Velocity Field', 'value': cv_field, - # 'propSettings':None, + 'propSettings': {}, } if div_field is not None: surface_data['signalMaps']['divergence_field'] = { - 'name' : 'Divergence Field', + 'name': 'Divergence Field', 'value': div_field, - # 'propSettings':None, + 'propSettings': {}, } return surface_data + def _add_electric_signal_props(electric_data, **kwargs): conduction_velocity = kwargs.get('conduction_velocity') divergence = kwargs.get('divergence') @@ -249,7 +250,6 @@ def _add_electric_signal_props(electric_data, **kwargs): signal_props = electric_data['annotations']['signalProps'] if conduction_velocity: - #TODO: connect propSetting from user setting signal_props['conduction_velocity'] = { 'name' : 'Conduction Velocity Values', 'value': conduction_velocity.values, @@ -260,7 +260,6 @@ def _add_electric_signal_props(electric_data, **kwargs): } if divergence: - #TODO: connect propSetting from user setting signal_props['divergence'] = { 'name' : 'Divergence Values', 'value': divergence.values, @@ -268,6 +267,7 @@ def _add_electric_signal_props(electric_data, **kwargs): return electric_data + def export_vtk( case: Case, filename: str, @@ -366,7 +366,6 @@ def _extract_surface_data( surface_data['longitudinal'] = fields.longitudinal_fibres if fields.longitudinal_fibres is not None else empty_float_array surface_data['transverse'] = fields.transverse_fibres if fields.transverse_fibres is not None else empty_float_array surface_data['pacing_site'] = fields.pacing_site if fields.pacing_site is not None else empty_int_array - surface_data['conduction_velocity'] = fields.conduction_velocity if fields.conduction_velocity is not None else empty_int_array surface_data['mesh_normals'] = fields.mesh_normals if fields.mesh_normals is not None else empty_int_array # Remove arrays that are full of NaNs From 0c534da94a543b82f0dcb579b0f6e7ba265fd897 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Thu, 18 Jan 2024 17:09:30 +0000 Subject: [PATCH 25/31] removed mesh_normals var --- openep/data_structures/surface.py | 7 ------- openep/io/writers.py | 2 -- 2 files changed, 9 deletions(-) diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index bb9988d1..77a05053 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -53,7 +53,6 @@ class Fields: pacing_site: np.ndarray = None conduction_velocity: np.ndarray = None cv_divergence: np.ndarray = None - mesh_normals : np.ndarray = None def __repr__(self): return f"fields: {tuple(self.__dict__.keys())}" @@ -204,11 +203,6 @@ def extract_surface_data(surface_data): except KeyError as e: cv_divergence = None - try: - mesh_normals = surface_data['mesh_normals'].astype(float) - except KeyError as e: - mesh_normals = None - fields = Fields( bipolar_voltage=bipolar_voltage, unipolar_voltage=unipolar_voltage, @@ -222,7 +216,6 @@ def extract_surface_data(surface_data): pacing_site=pacing_site, conduction_velocity=conduction_velocity, cv_divergence=cv_divergence, - mesh_normals=mesh_normals ) return points, indices, fields diff --git a/openep/io/writers.py b/openep/io/writers.py index 7b503371..854f2ff6 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -172,7 +172,6 @@ def export_openCARP( def export_openep_mat( case: Case, filename: str, - separate_regions: bool = False, ): """Export data in OpenEP format. @@ -366,7 +365,6 @@ def _extract_surface_data( surface_data['longitudinal'] = fields.longitudinal_fibres if fields.longitudinal_fibres is not None else empty_float_array surface_data['transverse'] = fields.transverse_fibres if fields.transverse_fibres is not None else empty_float_array surface_data['pacing_site'] = fields.pacing_site if fields.pacing_site is not None else empty_int_array - surface_data['mesh_normals'] = fields.mesh_normals if fields.mesh_normals is not None else empty_int_array # Remove arrays that are full of NaNs for field_name, field in surface_data.items(): From d4fac8f284f6bf1e44cfe8623e167cd565afd475 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Wed, 24 Jan 2024 18:16:29 +0000 Subject: [PATCH 26/31] changed var name points to center --- openep/analysis/analyse.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py index 0d3adb00..4fb50af2 100644 --- a/openep/analysis/analyse.py +++ b/openep/analysis/analyse.py @@ -52,7 +52,7 @@ class ConductionVelocity: Attributes: _case (Case): Internal reference to the case data. values (np.ndarray): The calculated conduction velocity values. - points (np.ndarray): The points corresponding to the conduction velocity values. + centers (np.ndarray): The centers corresponding to the conduction velocity values. Parameters: case (Case): The case data from which conduction velocity is to be calculated. @@ -60,7 +60,7 @@ class ConductionVelocity: def __init__(self, case): self._case = case self._values = None - self._points = None + self._centers = None @property def values(self): @@ -70,11 +70,11 @@ def values(self): return self._values @property - def points(self): - if self._points is None: - raise ValueError('Before accessing ``conduction_velocity.points`` ' + def centers(self): + if self._centers is None: + raise ValueError('Before accessing ``conduction_velocity.centers`` ' 'run ``divergence.calculate_divergence()``') - return self._points + return self._centers def calculate_cv( self, @@ -109,14 +109,14 @@ def calculate_cv( Returns: tuple: A tuple containing two elements: - values (np.ndarray): An array of calculated conduction velocity values. - - points (np.ndarray): An array of points corresponding to the calculated values. + - centers (np.ndarray): An array of centers corresponding to the calculated values. Raises: ValueError: If the specified method is not among the available options. Example: >> cv = ConductionVelocity(case) - >> values, points = cv.calculate(method='plane_fitting', apply_scalar_field=True) + >> values, centers = cv.calculate(method='plane_fitting', apply_scalar_field=True) """ supported_cv_methods = { @@ -132,16 +132,16 @@ def calculate_cv( lat, bipolar_egm_pts = preprocess_lat_egm(self._case, include) cv_method = supported_cv_methods[method] - self._values, self._points = cv_method(bipolar_egm_pts, lat, **kwargs) + self._values, self._centers = cv_method(bipolar_egm_pts, lat, **kwargs) if apply_scalar_field: self._case.fields.conduction_velocity = interpolate_general_cloud_points_onto_surface( case=self._case, cloud_values=self.values, - cloud_points=self.points, + cloud_points=self.centers, ) - return self.values, self.points + return self.values, self.centers class Divergence: From b1f4809bfd8d075cd82ce4c115692654c86fa340 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Fri, 2 Feb 2024 09:43:29 +0000 Subject: [PATCH 27/31] deal with empty force, cv and divergence --- openep/data_structures/surface.py | 58 +++++++++++++++---------------- 1 file changed, 28 insertions(+), 30 deletions(-) diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index 77a05053..f2f31e11 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -20,6 +20,7 @@ from attr import attrs import numpy as np +import warnings __all__ = [] @@ -141,31 +142,28 @@ def extract_surface_data(surface_data): unipolar_voltage = None impedance = None force = None + elif surface_data['uni_imp_frc'].size == 2: + unipolar_voltage, impedance = surface_data['uni_imp_frc'].T.astype(float) + force = None + warnings.warn("force data was not detected in surface_data, force=None", UserWarning) else: unipolar_voltage, impedance, force = surface_data['uni_imp_frc'].T.astype(float) - if all(np.isnan(unipolar_voltage)): - unipolar_voltage = None - if all(np.isnan(impedance)): - impedance = None - if all(np.isnan(force)): - force = None - try: - thickness = surface_data['thickness'].astype(float) - except KeyError as e: - thickness = None + if isinstance(force, np.ndarray) and all(np.isnan(unipolar_voltage)): + unipolar_voltage = None + if isinstance(force, np.ndarray) and all(np.isnan(impedance)): + impedance = None + if isinstance(force, np.ndarray) and all(np.isnan(force)): + force = None - if isinstance(thickness, np.ndarray) and thickness.size == 0: - thickness = None + thickness = surface_data.get('thickness', None) + if isinstance(thickness, np.ndarray): + thickness = None if thickness.size == 0 else thickness.astype(float) # This is defined on a per-cell bases - try: - cell_region = surface_data['cell_region'].astype(int) - except KeyError as e: - cell_region = None - - if isinstance(cell_region, np.ndarray) and cell_region.size == 0: - cell_region = None + cell_region = surface_data.get('cell_region', None) + if isinstance(cell_region, np.ndarray): + cell_region = None if cell_region.size == 0 else cell_region.astype(int) # Fibre orientation are vectors defined on a per-cell basis try: @@ -185,22 +183,22 @@ def extract_surface_data(surface_data): transverse_fibres = None # Pacing site point ids (-1 is not pacing site) - try: - pacing_site = surface_data['pacing_site'].astype(int) - except KeyError as e: - pacing_site = None + pacing_site = surface_data.get('pacing_site', None) + if isinstance(pacing_site, np.ndarray): + pacing_site = None if pacing_site.size == 0 else pacing_site.astype(int) - if isinstance(pacing_site, np.ndarray) and pacing_site.size == 0: - pacing_site = None - try: - conduction_velocity = surface_data['signalMaps']['conduction_velocity_field']['value'].astype(float) - except KeyError as e: + 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']['value'].astype(float) - except KeyError as e: + 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 fields = Fields( From 9df2bb2aac1c48f1d34355e5d39f6532e992747b Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Fri, 2 Feb 2024 09:44:57 +0000 Subject: [PATCH 28/31] de-reference ecgNames from hdf5 during mat load --- openep/io/matlab.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openep/io/matlab.py b/openep/io/matlab.py index afe9a7a6..22be3156 100644 --- a/openep/io/matlab.py +++ b/openep/io/matlab.py @@ -51,6 +51,7 @@ def _visit_mat_v73_(file_pointer): 'userdata/electric/electrodeNames_uni', 'userdata/electric/names', 'userdata/electric/tags', + 'userdata/electric/ecgNames' } strings_to_decode = { From 503b7e250f758027f79caaa26da210deac7cbca9 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Fri, 2 Feb 2024 09:45:49 +0000 Subject: [PATCH 29/31] fix warn message typo --- openep/data_structures/surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index f2f31e11..db739d30 100644 --- a/openep/data_structures/surface.py +++ b/openep/data_structures/surface.py @@ -145,7 +145,7 @@ def extract_surface_data(surface_data): elif surface_data['uni_imp_frc'].size == 2: unipolar_voltage, impedance = surface_data['uni_imp_frc'].T.astype(float) force = None - warnings.warn("force data was not detected in surface_data, force=None", UserWarning) + warnings.warn("Force data was not detected in surface_data, force=None", UserWarning) else: unipolar_voltage, impedance, force = surface_data['uni_imp_frc'].T.astype(float) From 67c4b4414326949f85d5fd613e4dce5cf725cd13 Mon Sep 17 00:00:00 2001 From: vinush-vignes Date: Fri, 9 Feb 2024 12:48:51 +0000 Subject: [PATCH 30/31] added interpolation method kws --- openep/analysis/_conduction_velocity.py | 44 +++++++++++++++++++++++-- openep/analysis/analyse.py | 26 +++++++++++---- openep/case/case_routines.py | 1 + openep/converters/pyvista_converters.py | 2 +- 4 files changed, 63 insertions(+), 10 deletions(-) diff --git a/openep/analysis/_conduction_velocity.py b/openep/analysis/_conduction_velocity.py index 117ff16d..e251b0b5 100644 --- a/openep/analysis/_conduction_velocity.py +++ b/openep/analysis/_conduction_velocity.py @@ -160,10 +160,10 @@ def triangulation( Perform triangulation on electrogram data to calculate conduction velocities. Args: - bipolar_egm_pts (array-like): + bipolar_egm_pts (array): Array of bipolar electrogram points of size Nx3 - local_activation_time (array-like): + local_activation_time (array): Local activation times corresponding to bipolar_egm_pts. min_theta (float): @@ -291,16 +291,54 @@ def divergence( local_activation_time, collision_threshold=-1, focal_threshold=1, - output_binary_field=False + output_binary_field=False, + interpolation_kws=None, ): + """ + Calculate divergence of the conduction velocity (CV) field to identify regions of collision and + focal activation within cardiac tissue based on local activation times (LATs) and electrogram points. + + Args: + case (Case): The case object containing cardiac geometry and electrogram data. + + bipolar_egm_pts (np.ndarray): An array of coordinates for bipolar electrogram points. + + local_activation_time (np.ndarray): An array of local activation times corresponding to + bipolar electrogram points. + + collision_threshold (float, optional): Divergence value below which regions are considered + as collision fronts. Defaults to -1. + + focal_threshold (float, optional): Divergence value above which regions are considered as + focal activation fronts. Defaults to 1. + + output_binary_field (bool, optional): If True, the output divergence field is binary, with 1 + indicating regions meeting the collision or focal threshold, + and 0 otherwise. If False, the continuous divergence field + is returned. Defaults to False. + + interpolation_kws (dict, optional): Keyword arguments for interpolation function. + > interpolate_general_cloud_points_onto_surface(**interpolation_kws) + Defaults to None. + + Returns: + tuple: A tuple containing two elements: + - norm_cv_direction (np.ndarray): The normalized direction of conduction velocity + across the cardiac tissue surface. + - divergence (np.ndarray): The divergence of the activation direction field. This + can either be a continuous field or a binary field indicating focal and collision + regions based on the 'output_binary_field' parameter. + """ temp_mesh = case.create_mesh() basic_mesh = pv.PolyData(temp_mesh.points, temp_mesh.faces) + interpolation_kws = dict() if interpolation_kws is None else interpolation_kws interpolated_scalar = interpolate_general_cloud_points_onto_surface( case=case, cloud_values=local_activation_time, cloud_points=bipolar_egm_pts, + **interpolation_kws ) basic_mesh['LAT_scalar'] = interpolated_scalar diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py index 4fb50af2..ced65679 100644 --- a/openep/analysis/analyse.py +++ b/openep/analysis/analyse.py @@ -81,8 +81,9 @@ def calculate_cv( method='triangulation', include=None, apply_scalar_field=True, - **kwargs - ): + interpolation_kws=None, + **method_kwargs + ): """ Calculate the conduction velocity and interpolate points onto a mesh. @@ -104,7 +105,12 @@ def calculate_cv( apply_scalar_field (bool, optional): If True, applies the calculated conduction velocity values as a scalar field on the case object's mesh. Defaults to True. - **kwargs: Additional keyword arguments specific to the chosen calculation method. + + interpolation_kws (dict, optional): Keyword arguments for interpolation function. + > interpolate_general_cloud_points_onto_surface(**interpolation_kws) + Defaults to None. + + **method_kwargs: Additional keyword arguments specific to the chosen calculation method. Returns: tuple: A tuple containing two elements: @@ -128,17 +134,19 @@ def calculate_cv( if method.lower() not in supported_cv_methods: raise ValueError(f"`method` must be one of {supported_cv_methods.keys()}.") + 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) cv_method = supported_cv_methods[method] - self._values, self._centers = cv_method(bipolar_egm_pts, lat, **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( case=self._case, cloud_values=self.values, cloud_points=self.centers, + **interpolation_kws ) return self.values, self.centers @@ -180,7 +188,8 @@ def calculate_divergence( self, include=None, output_binary_field=False, - apply_scalar_field=True + apply_scalar_field=True, + interpolation_kws=None, ): """ Calculate the divergence of conduction velocity and optionally apply the result as a scalar field. @@ -199,6 +208,10 @@ def calculate_divergence( the case object. Defaults to True. + interpolation_kws (dict, optional): Keyword arguments for interpolation function. + > interpolate_general_cloud_points_onto_surface(**interpolation_kws) + Defaults to None. + Returns: tuple of (np.ndarray, np.ndarray): A tuple containing the divergence direction and divergence values. The direction is a 2D array of shape (N, 3), where N is the number @@ -218,7 +231,8 @@ def calculate_divergence( case=self._case, bipolar_egm_pts=bipolar_egm_pts, local_activation_time=lat, - output_binary_field=output_binary_field + output_binary_field=output_binary_field, + interpolation_kws=interpolation_kws, ) if apply_scalar_field: diff --git a/openep/case/case_routines.py b/openep/case/case_routines.py index 12fe6e68..69249c0d 100644 --- a/openep/case/case_routines.py +++ b/openep/case/case_routines.py @@ -76,6 +76,7 @@ class and associated keyword arguments to the `method` and `method_kws` 'interpolate_activation_time_onto_surface', 'interpolate_voltage_onto_surface', 'bipolar_from_unipolar_surface_points', + 'interpolate_general_cloud_points_onto_surface' ] diff --git a/openep/converters/pyvista_converters.py b/openep/converters/pyvista_converters.py index e670b8aa..92a6c893 100644 --- a/openep/converters/pyvista_converters.py +++ b/openep/converters/pyvista_converters.py @@ -93,7 +93,7 @@ def to_pyvista( mesh.point_data["LAT"] = case.fields.local_activation_time mesh.point_data["Impedance"] = case.fields.impedance mesh.point_data["Force"] = case.fields.force - mesh.point_data["Conduction Velocity"] = case.fields.conduction_velocity + mesh.point_data["Conduction velocity"] = case.fields.conduction_velocity mesh.point_data["Divergence"] = case.fields.cv_divergence return mesh From d225828ff30e569da42f2e17ead249ea0a23e30a Mon Sep 17 00:00:00 2001 From: Vinush Vigneswaran <150042692+vinush-vignes@users.noreply.github.com> Date: Sat, 14 Sep 2024 13:41:32 +0100 Subject: [PATCH 31/31] 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 | 2 +- openep/analysis/analyse.py | 60 ++++----- openep/data_structures/ablation.py | 8 +- openep/data_structures/arrows.py | 62 +++++++++ openep/data_structures/case.py | 3 + openep/data_structures/surface.py | 29 ++-- openep/draw/draw_routines.py | 40 ++++-- openep/io/readers.py | 89 +++++++++--- openep/io/writers.py | 208 +++++++++++++++++++++++------ 9 files changed, 382 insertions(+), 119 deletions(-) create mode 100644 openep/data_structures/arrows.py diff --git a/openep/__init__.py b/openep/__init__.py index e775d84c..e2d4fd82 100644 --- a/openep/__init__.py +++ b/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 diff --git a/openep/analysis/analyse.py b/openep/analysis/analyse.py index ced65679..84c2326d 100644 --- a/openep/analysis/analyse.py +++ b/openep/analysis/analyse.py @@ -17,7 +17,6 @@ # with this program (LICENSE.txt). If not, see """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, diff --git a/openep/data_structures/ablation.py b/openep/data_structures/ablation.py index 8f285968..d610c18a 100644 --- a/openep/data_structures/ablation.py +++ b/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) diff --git a/openep/data_structures/arrows.py b/openep/data_structures/arrows.py new file mode 100644 index 00000000..10ca39f3 --- /dev/null +++ b/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 diff --git a/openep/data_structures/case.py b/openep/data_structures/case.py index 57c64b01..50394065 100644 --- a/openep/data_structures/case.py +++ b/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): diff --git a/openep/data_structures/surface.py b/openep/data_structures/surface.py index db739d30..16711a31 100644 --- a/openep/data_structures/surface.py +++ b/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( diff --git a/openep/draw/draw_routines.py b/openep/draw/draw_routines.py index cb8d1065..8a56c9c0 100644 --- a/openep/draw/draw_routines.py +++ b/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,14 +77,15 @@ 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)] @@ -91,13 +93,31 @@ def draw_free_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 diff --git a/openep/io/readers.py b/openep/io/readers.py index eba6199c..4b9077e7 100644 --- a/openep/io/readers.py +++ b/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 diff --git a/openep/io/writers.py b/openep/io/writers.py index 854f2ff6..1f3c14c5 100644 --- a/openep/io/writers.py +++ b/openep/io/writers.py @@ -51,9 +51,10 @@ """ -import pathlib +import pandas as pd import numpy as np import scipy.io +import pathlib from openep.data_structures.ablation import Ablation from openep.data_structures.case import Case @@ -64,6 +65,8 @@ "export_openCARP", "export_openep_mat", "export_vtk", + "export_vtx", + "export_csv" ] @@ -72,6 +75,7 @@ def export_openCARP( prefix: str, scale_points: float = 1, export_transverse_fibres: bool = True, + export_pacing_site: bool = True, ): """Export mesh data from an OpenEP data to openCARP format. @@ -104,59 +108,58 @@ def export_openCARP( comments='', ) - # Save elements info + # 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) - elements = np.concatenate( - [ - cell_type[:, np.newaxis], - case.indices, - cell_region[:, np.newaxis], - ], - axis=1, - dtype=object, - ) + 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"), - elements, + combined_cell_data, fmt="%s %d %d %d %d", - header=str(n_triangles), + header=str(n_lines), comments='', ) - # Save fibres - n_fibre_vectors = 2 if export_transverse_fibres else 1 - fibres = np.zeros((n_triangles, n_fibre_vectors * 3), dtype=float) + # Save linear connections info + if case.arrows.linear_connections is not None: + 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) - if case.fields.longitudinal_fibres is not None: - fibres[:, :3] = case.fields.longitudinal_fibres - else: - fibres[:, 0] = 1 + 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) - if export_transverse_fibres: - if case.fields.transverse_fibres is not None: - fibres[:, 3:] = case.fields.transverse_fibres - else: - fibres[:, 3] = 1 + # 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", + ) - np.savetxt( - output_path.with_suffix('.lon'), - fibres, - fmt="%.6f", - header=str(n_fibre_vectors), - comments='', - ) + # Save fibres + if case.arrows.fibres is not None: + np.savetxt( + output_path.with_suffix('.lon'), + case.arrows.fibres, + fmt="%.6f", + comments='', + ) # Saving pacing sites if they exist - if case.fields.pacing_site is None: + if case.fields.pacing_site is None or not export_pacing_site: return # 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 @@ -169,6 +172,58 @@ def export_openCARP( ) +def export_vtx( + case: Case, + prefix: str, + pacing_site_internal_name: str +): + """ + Export vertex indices for a specified pacing site to a .vtx file. + + Args: + - case (Case): The case object containing the pacing site and landmark information. + - prefix (str): The file path prefix for saving the output .vtx file. + - pacing_site_internal_name (str): The tag of the pacing site to export. + + Returns: + - pathlib.Path: The path to the created .vtx file, or None if the pacing site is not found. + + Raises: + - IndexError: If the specified pacing site name is not found within the case landmarks. + """ + # TODO: Need to improve indexing method to find pacing site from landmarks, + # currently finds pacing sites using tag string. + + output_path = pathlib.Path(prefix).resolve() + + if case.fields.pacing_site is None: + return + + # This extracts the pacing sites from all landmarks by checking if the tag is of the form "Pacing site 1" + landmarks = case.electric.landmark_points + try: + landmark_index = np.where(landmarks.internal_names == pacing_site_internal_name)[0][0] + except IndexError: + raise IndexError(f"Pacing site: Expecting {landmarks.internal_names}, received \"{pacing_site_internal_name}\".") + + landmark_name = landmarks.names[landmark_index] + site_index = int(landmark_name.replace('Pacing site ', '')) + + pacing_site_points = np.nonzero(case.fields.pacing_site == site_index)[0] + n_points = pacing_site_points.size + + output_path_with_suffix = output_path.with_name(f'{output_path.name}_{pacing_site_internal_name}.vtx') + + np.savetxt( + output_path_with_suffix, + pacing_site_points, + header=f'{n_points}\nintra', + comments='', + fmt='%d', + ) + return output_path_with_suffix + + def export_openep_mat( case: Case, filename: str, @@ -195,13 +250,17 @@ def export_openep_mat( ) userdata['electric'] = _extract_electric_data(electric=case.electric) - userdata['electric'] = _add_electric_signal_props( - electric_data=userdata['electric'], - conduction_velocity=case.analyse.conduction_velocity, - divergence=case.analyse.divergence - ) - userdata['rf'] = _export_ablation_data(ablation=case.ablation) + if case.analyse.conduction_velocity.values is not None: + userdata['electric'] = _add_electric_signal_props( + electric_data=userdata['electric'], + conduction_velocity=case.analyse.conduction_velocity, + divergence=case.analyse.divergence + ) + + if case.ablation is not None: + userdata['rf'] = _export_ablation_data(ablation=case.ablation) + scipy.io.savemat( file_name=filename, mdict={'userdata': userdata}, @@ -236,6 +295,56 @@ def _add_surface_maps(surface_data, **kwargs): return surface_data +def export_csv( + system, + filename: str, + selections: dict, +): + + """Export data in CSV format. + + Saves selected field data. + + Args: + 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, + 'Unipolar voltage': case.fields.unipolar_voltage, + 'LAT': case.fields.local_activation_time, + 'Force': case.fields.force, + 'Impedance': case.fields.impedance, + 'Thickness': case.fields.thickness, + '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() + + for field_name, checked in selections.items(): + if checked: + field_data = available_exports.get(field_name) + if field_data is not None: + df[field_name] = pd.Series(field_data) + else: + df[field_name] = pd.NA + + df.to_csv(filename, index=False, encoding='utf-8') + + def _add_electric_signal_props(electric_data, **kwargs): conduction_velocity = kwargs.get('conduction_velocity') divergence = kwargs.get('divergence') @@ -255,7 +364,7 @@ def _add_electric_signal_props(electric_data, **kwargs): } signal_props['cvX'] = { 'name' : 'Conduction Velocity Coordinates', - 'value': conduction_velocity.points, + 'value': conduction_velocity.centers, } if divergence: @@ -474,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 \ No newline at end of file + 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