Compute_volume_ef

from read_data import * import glob from scipy.ndimage.morphology import binary_erosion, binary_dilation, generate_binary_structure import numpy as np from skimage.transform import resize from scipy.stats import pearsonr import cv2 import math from scipy.stats import pearsonr import matplotlib.pyplot as plt from datetime import date,datetime def get_mae(records_real, records_predict): """ MAE """ if len(records_real) == len(records_predict): return sum([abs(x - y) for x, y in zip(records_real, records_predict)]) / len(records_real) else: return None def diameter_lenght(ch, spacing_info): ch_endo = (ch == 1).astype(np.int16) #the area of left ventricle ch_la = (ch == 3).astype(np.int16) #the area of left atrium footprint = generate_binary_structure(2, 1) ch_endo_bottom = ch_endo & binary_dilation(ch_la, structure=footprint, iterations=1) index = np.where(ch_endo == 1) listofcoordinates = list(zip(index[0], index[1])) bottom_y_sum = np.sum(ch_endo_bottom, axis=0) # print(np.min(np.nonzero(bottom_y_sum)), np.max(np.nonzero(bottom_y_sum))) bottom_middle_y = (np.min(np.nonzero(bottom_y_sum)) + np.max(np.nonzero(bottom_y_sum))) / 2 # bottom_middle_x = ch_endo_bottom[:, np.array(bottom_y).astype(np.int16)].argmax() #k_bottom is the slope of botton line k_bottom=(np.max(np.nonzero(bottom_y_sum)) - np.min(np.nonzero(bottom_y_sum)))/\ (ch_endo_bottom[:, np.max(np.nonzero(bottom_y_sum)).astype(np.int16)].argmax()- ch_endo_bottom[:, np.min(np.nonzero(bottom_y_sum)).astype(np.int16)].argmax()+1e-10) b_bottom = np.max(np.nonzero(bottom_y_sum)) - \ k_bottom * ch_endo_bottom[:, np.max(np.nonzero(bottom_y_sum)).astype(np.int16)].argmax() bottom_middle_x=(bottom_middle_y-b_bottom)/k_bottom maxl = 0 for (a, b) in listofcoordinates: l = math.sqrt(math.pow(a - bottom_middle_x, 2) + math.pow(b - bottom_middle_y, 2)) if l > maxl: maxl = l (apex_x, apex_y) = (a, b) # apex cordis L = np.sqrt((np.power((bottom_middle_y.astype(np.float64) - apex_y), 2) * spacing_info[1]) + np.power(((bottom_middle_x.astype(np.float64) - apex_x) * spacing_info[0]), 2)) # line formed by the apex of the heart and the midpoint of the mitral valve k_long_axis = (bottom_middle_y - apex_y) / (bottom_middle_x - apex_x) b_long_axis = bottom_middle_y - k_long_axis * bottom_middle_x # k_radius = -1 / (k_long_axis+0.0000000001) #Straight line(cross_line) perpendicular to the longest diameter interval = (bottom_middle_x - apex_x) / 20 # 20 planes x_point_20 = np.arange(1, 21) * interval + apex_x y_point_20 = k_long_axis * x_point_20 + b_long_axis b_20 = y_point_20 - k_bottom * x_point_20 # endocardium_border ch_endo_border = ch_endo ^ binary_erosion(ch_endo, structure=footprint, iterations=1) ch_endo_border = ch_endo_border.astype(np.int16) - ch_endo_bottom.astype(np.int16) ch_endo_border_x_y = np.nonzero(ch_endo_border) y = k_long_axis * ch_endo_border_x_y[0] + b_long_axis right_border = y < ch_endo_border_x_y[1] right_border_x = ch_endo_border_x_y[0][right_border] right_border_y = ch_endo_border_x_y[1][right_border] left_border = y > ch_endo_border_x_y[1] left_border_x = ch_endo_border_x_y[0][left_border] left_border_y = ch_endo_border_x_y[1][left_border] # fig, axes=plt.subplots() # axes.scatter(left_border_x, left_border_y, linewidths=0.05) # axes.scatter(right_border_x, right_border_y, linewidths=0.05) # plt.savefig('./left_border.png') right_point_index = np.zeros((20, ), dtype=np.int64) left_point_index = np.zeros((20, ), dtype=np.int64) for i in range(20): right_point_index[i] = np.array((np.abs(k_bottom*right_border_x-right_border_y+b_20[i])).argmin()) #coordinate closest the cross_line # l=sorted(np.array((np.abs(k_bottom*right_border_x-right_border_y+b_20[i])))) # print(l[]) left_point_index[i] = np.array((np.abs(k_bottom*left_border_x-left_border_y+b_20[i])).argmin()) # right_index, left_index = right_point_index[i], left_point_index[i] # print(i, len(right_border_y), right_index, left_index) # print(right_border_x[right_index], left_border_x[left_index]) right_border_point_x = right_border_x[right_point_index] right_border_point_y = right_border_y[right_point_index] left_border_point_x = left_border_x[left_point_index] left_border_point_y = left_border_y[left_point_index] right_radius_lenght = np.sqrt((np.power(((right_border_point_y.astype(np.float64) - y_point_20)* spacing_info[1]), 2)) + np.power(((right_border_point_x.astype(np.float64) - x_point_20) * spacing_info[0]), 2)) left_radius_lenght = np.sqrt((np.power(((left_border_point_y.astype(np.float64) - y_point_20) * spacing_info[1]), 2)) + np.power(((left_border_point_x.astype(np.float64) - x_point_20) * spacing_info[0]), 2)) diameter_lenght = right_radius_lenght+left_radius_lenght return diameter_lenght, L def compute_volume(ch2, ch4, spacing): # row is x, col is y, left-upper is (0,0) diameter_ch2, L_ch2 = diameter_lenght(ch2, spacing) diameter_ch4, L_ch4 = diameter_lenght(ch4, spacing) # print(length_ch2, len(length_ch2), length_ch2.shape) volume = 3.1415926 / 4 * min(L_ch2, L_ch4) /20 * np.sum(diameter_ch2 * diameter_ch4) return round(volume/1000, 1) EDV_gt = np.zeros(450, ) ESV_gt = np.zeros(450, ) EF_gt = np.zeros(450, ) EDV_segflow=np.zeros(450,) ESV_segflow=np.zeros(450,) EF_segflow=np.zeros(450,) # EF_gt = read_train_sequences_volume_EF_Info() # val_pat_id1 = np.load('/data/ch/CAMUS/skip_clstm/dataloader/'+ '/cv_val_{}fold_id.npy'.format(1)) # for k in range(1, 11): # val_pat_id1 = np.load('/data/ch/CAMUS/skip_clstm/dataloader/CV10_validation_set/CV_val_set_' + str(k) + '.npy') # # val_pat_id1 = np.load('/data/ch/CAMUS/skip_clstm/dataloader/' + '/cv_val_{}fold_id.npy'.format(k)) # print(k, val_pat_id1, val_pat_id1.shape) # val_pat_id = val_pat_id1.tolist() # if k==1: # val_pat_id.remove(43) # if k==2: # val_pat_id.remove(191) # # if k==3: # # val_pat_id.remove(37) # # val_pat_id.remove(57) #5fold # if k==6: # val_pat_id.remove(434) # val_pat_id.remove(285) # if k==5: # val_pat_id.remove(57) # val_pat_id.remove(52) # if k==8: # val_pat_id.remove(37) # val_pat_id.remove(41) # if k==7: # val_pat_id.remove(86) # val_pat_id.remove(271) # if k==9: # val_pat_id.remove(28) # if k==10: # val_pat_id.remove(72) # val_pat_id.remove(98) # # print(val_pat_id1, val_pat_id) # files = glob.glob( '/data/ch/CAMUS/skip_clstm/output/segFlow_seq_attention_4down_200sample/val_mask/1/30' + '{}_1.0_pred.npy'.format(432)) # print(val_pat_id) # for i in range(len(val_pat_id)): # pat_id=val_pat_id[i] # print("Start Computing Pat{}".format(pat_id)) # CH4edfile = glob.glob('/data/ch/CAMUS/skip_clstm/output/segFlow_sharednet_seq_cos_4down_1CV10_catflow/val/seg/{}/'.format(k) + '{}_CH4_pred.npy'.format(pat_id)) # CH4esfile = glob.glob('/data/ch/CAMUS/skip_clstm/output/segFlow_sharednet_seq_cos_4down_1CV10_catflow/val/seg/{}/'.format(k) + '{}_CH4_pred.npy'.format(pat_id)) # CH2edfile = glob.glob('/data/ch/CAMUS/skip_clstm/output/segFlow_sharednet_seq_cos_4down_1CV10_catflow/val/seg/{}/'.format(k) + '{}_CH2_pred.npy'.format(pat_id)) # CH2esfile = glob.glob('/data/ch/CAMUS/skip_clstm/output/segFlow_sharednet_seq_cos_4down_1CV10_catflow/val/seg/{}/'.format(k)+ '{}_CH2_pred.npy'.fo