import numpy as np
import os, sys, shutil, time, re
import h5py
import skimage.morphology as morph
from tqdm import tqdm
import matplotlib.pyplot as plt
from matplotlib import animation
import time
import pickle
# For ROI extraction
import skimage.transform
from scipy.fftpack import fftn, ifftn
from skimage.feature import peak_local_max, canny
from skimage.transform import hough_circle
# Nifti processing
import nibabel as nib
from collections import OrderedDict
# print sys.path
# sys.path.append("..")
import errno
np.random.seed(42)
# Helper functions
## Heart Metrics
def heart_metrics(seg_3Dmap, voxel_size, classes=[3, 1, 2]):
"""
Compute the volumes of each classes
"""
# Loop on each classes of the input images
volumes = []
for c in classes:
# Copy the gt image to not alterate the input
seg_3Dmap_copy = np.copy(seg_3Dmap)
seg_3Dmap_copy[seg_3Dmap_copy != c] = 0
# Clip the value to compute the volumes
seg_3Dmap_copy = np.clip(seg_3Dmap_copy, 0, 1)
# Compute volume
volume = seg_3Dmap_copy.sum() * np.prod(voxel_size) / 1000.
volumes += [volume]
return volumes
def ejection_fraction(ed_vol, es_vol):
"""
Calculate ejection fraction
"""
stroke_vol = ed_vol - es_vol
return (np.float(stroke_vol)/np.float(ed_vol))*100
def myocardialmass(myocardvol):
"""
Specific gravity of heart muscle (1.05 g/ml)
"""
return myocardvol*1.05
def imshow(*args,**kwargs):
""" Handy function to show multiple plots in on row, possibly with different cmaps and titles
Usage:
imshow(img1, title="myPlot")
imshow(img1,img2, title=['title1','title2'])
imshow(img1,img2, cmap='hot')
imshow(img1,img2,cmap=['gray','Blues']) """
cmap = kwargs.get('cmap', 'gray')
title= kwargs.get('title','')
if len(args)==0:
raise ValueError("No images given to imshow")
elif len(args)==1:
plt.title(title)
plt.imshow(args[0], interpolation='none')
else:
n=len(args)
if type(cmap)==str:
cmap = [cmap]*n
if type(title)==str:
title= [title]*n
plt.figure(figsize=(n*5,10))
for i in range(n):
plt.subplot(1,n,i+1)
plt.title(title[i])
plt.imshow(args[i], cmap[i])
plt.show()
def plot_roi(data4D, roi_center, roi_radii):
"""
Do the animation of full heart volume
"""
x_roi_center, y_roi_center = roi_center[0], roi_center[1]
x_roi_radius, y_roi_radius = roi_radii[0], roi_radii[1]
print ('nslices', data4D.shape[2])
zslices = data4D.shape[2]
tframes = data4D.shape[3]
slice_cnt = 0
for slice in [data4D[:,:,z,:] for z in range(zslices)]:
outdata = np.swapaxes(np.swapaxes(slice[:,:,:], 0,2), 1,2)
roi_mask = np.zeros_like(outdata[0])
roi_mask[x_roi_center - x_roi_radius:x_roi_center + x_roi_radius,
y_roi_center - y_roi_radius:y_roi_center + y_roi_radius] = 1
outdata[:, roi_mask > 0.5] = 0.8 * outdata[:, roi_mask > 0.5]
outdata[:, roi_mask > 0.5] = 0.8 * outdata[:, roi_mask > 0.5]
fig = plt.figure(1)
fig.canvas.set_window_title('slice_No' + str(slice_cnt))
slice_cnt+=1
def init_out():
im.set_data(outdata[0])
def animate_out(i):
im.set_data(outdata[i])
return im
im = fig.gca().imshow(outdata[0], cmap='gray')
anim = animation.FuncAnimation(fig, animate_out, init_func=init_out, frames=tframes, interval=50)
anim.save('Cine_MRI_SAX_%d.mp4'%slice_cnt, fps=50, extra_args=['-vcodec', 'libx264'])
plt.show()
def plot_4D(data4D):
"""
Do the animation of full heart volume
"""
print ('nslices', data4D.shape[2])
zslices = data4D.shape[2]
tframes = data4D.shape[3]
slice_cnt = 0
for slice in [data4D[:,:,z,:] for z in range(zslices)]:
outdata = np.swapaxes(np.swapaxes(slice[:,:,:], 0,2), 1,2)
fig = plt.figure(1)
fig.canvas.set_window_title('slice_No' + str(slice_cnt))
slice_cnt+=1
def init_out():
im.set_data(outdata[0])
def animate_out(i):
im.set_data(outdata[i])
return im
im = fig.gca().imshow(outdata[0], cmap='gray')
anim = animation.FuncAnimation(fig, animate_out, init_func=init_out, frames=tframes, interval=50)
plt.show()
def multilabel_split(image_tensor):
"""
image_tensor : Batch * H * W
Split multilabel images and return stack of images
Returns: Tensor of shape: Batch * H * W * n_class (4D tensor)
# TODO: Be careful: when using this code: labels need to be
defined, explictly before hand as this code does not handle
missing labels
So far, this function is okay as it considers full volume for
finding out unique labels
"""
labels = np.unique(image_tensor)
batch_size = image_tensor.shape[0]
out_shape = image_tensor.shape + (len(labels),)
image_tensor_4D = np.zeros(out_shape, dtype='uint8')
for i in xrange(batch_size):
cnt = 0
shape =image_tensor.shape[1:3] + (len(labels),)
temp = np.ones(shape, dtype='uint8')
for label in labels:
temp[...,cnt] = np.where(image_tensor[i] == label, temp[...,cnt], 0)
cnt += 1
image_tensor_4D[i] = temp
return image_tensor_4D
def save_data(data, filename, out_path):
out_filename = os.path.join(out_path, filename)
with open(out_filename, 'wb') as f:
pickle.dump(data, f, protocol=pickle.HIGHEST_PROTOCOL)
print ('saved to %s' % out_filename)
def load_pkl(path):
with open(path, 'rb') as f:
obj = pickle.load(f)
return obj
### Stratified Sampling of data
# Refer:
# http://www.echopedia.org/wiki/Left_Ventricular_Dimensions
# https://www.creatis.insa-lyon.fr/Challenge/acdc/databases.html
# https://en.wikipedia.org/wiki/Body_surface_area
# 30 normal subjects - NOR
NORMAL = 'NOR'
# 30 patients with previous myocardial infarction
# (ejection fraction of the left ventricle lower than 40% and several myocardial segments with abnormal contraction) - MINF
MINF = 'MINF'
# 30 patients with dilated cardiomyopathy
# (diastolic left ventricular volume >100 mL/m2 and an ejection fraction of the left ventricle lower than 40%) - DCM
DCM = 'DCM'
# 30 patients with hypertrophic cardiomyopathy
# (left ventricular cardiac mass high than 110 g/m2,
# several myocardial segments with a thickness higher than 15 mm in diastole and a normal ejecetion fraction) - HCM
HCM = 'HCM'
# 30 patients with abnormal right ventricle (volume of the right ventricular
# cavity higher than 110 mL/m2 or ejection fraction of the rigth ventricle lower than 40%) - RV
RV = 'RV'
def copy(src, dest):
"""
Copy function
"""
try:
shutil.copytree(src, dest, ignore=shutil.ignore_patterns())
except OSError as e:
# If the error was caused because the source wasn't a directory
if e.errno == errno.ENOTDIR:
shutil.copy(src, dest)
else:
print('Directory not copied. Error: %s' % e)
def read_patient_cfg(path):
"""
Reads patient data in the cfg file and returns a dictionary
"""
patient_info = {}
with open(os.path.join(path, 'Info.cfg')) as f_in:
for line in f_in:
l = line.rstrip().split(": ")
patient_info[l[0]] = l[1]
return patient_info
def group_patient_cases(src_path, out_path, force=False):
""" Group the patient data according to cardiac pathology"""
cases = sorted(next(os.walk(src_path))[1])
dest_path = os.path.join(out_path, 'Patient_Groups')
if force:
shutil.rmtree(dest_path)
if os.path.exists(dest_path):
return dest_path
os.makedirs(dest_path)
os.mkdir(os.path.join(dest_path, NORMAL))
os.mkdir(os.path.join(dest_path, MINF))
os.mkdir(os.path.join(dest_path, DCM))
os.mkdir(os.path.join(dest_path, HCM))
os.mkdir(os.path.join(dest_path, RV))
for case in cases:
full_path = os.path.join(src_path, case)
copy(full_path, os.path.join(dest_path,\
read_patient_cfg(full_path)['Group'], case))
def generate_train_validate_test_set(src_path, dest_path):
"""
Split the data into 70:15:15 for train-validate-test set
arg: path: input data path
"""
SPLIT_TRAIN = 0.7
SPLIT_VALID = 0.15
dest_path = os.path.join(dest_path,'dataset')
if os.path.exists(dest_path):
shutil.rmtree(dest_path)
os.makedirs(os.path.join(dest_path, 'train_set'))
os.makedirs(os.path.join(dest_path, 'validation_set'))
os.makedirs(os.path.join(dest_path, 'test_set'))
# print (src_path)
groups = next(os.walk(src_path))[1]
for group in groups:
group_path = next(os.walk(os.path.join(src_path, group)))[0]
patient_folders = next(os.walk(group_path))[1]
np.random.shuffle(patient_folders)
train_ = patient_folders[0:int(SPLIT_TRAIN*len(patient_folders))]
valid_ = patient_folders[int(SPLIT_TRAIN*len(patient_folders)):
int((SPLIT_TRAIN+SPLIT_VALID)*len(patient_folders))]
test_ = patient_folders[int((SPLIT_TRAIN+SPLIT_VALID)*len(patient_folders)):]
for patient in train_:
folder_path = os.path.join(group_path, patient)
copy(folder_path, os.path.join(dest_path, 'train_set', patient))
for patient in valid_:
folder_path = os.path.join(group_path, patient)
copy(folder_path, os.path.join(dest_path, 'validation_set', patient))
for patient in test_:
folder_path = os.path.join(group_path, patient)
copy(folder_path, os.path.join(dest_path, 'test_set', patient))
# Fourier-Hough Transform Based ROI Extraction
def extract_roi_fft(data4D, pixel_spacing, minradius_mm=15, maxradius_mm=45, kernel_width=5,
center_margin=8, num_peaks=10, num_circles=20, radstep=2):
"""
Returns center and radii of ROI region in (i,j) format
"""
# Data shape:
# radius of the smallest and largest circles in mm estimated from the train set
# convert to pixel counts
pixel_spacing_X, pixel_spacing_Y, _,_ = pixel_spacing
minradius = int(minradius_mm / pixel_spacing_X)
maxradius = int(maxradius_mm / pixel_spacing_Y)
ximagesize = data4D.shape[0]
yimagesize = data4D.shape[1]
zslices = data4D.shape[2]
tframes = data4D.shape[3]
xsurface = np.tile(range(ximagesize), (yimagesize, 1)).T
ysurface = np.tile(range(yimagesize), (ximagesize, 1))
lsurface = np.zeros((ximagesize, yimagesize))
allcenters = []
allaccums = []
allradii = []
for slice in range(zslices):
ff1 = fftn([data4D[:,:,slice, t] for t in range(tframes)])
fh = np.absolute(ifftn(ff1[1, :, :]))
fh[fh < 0.1 * np.max(fh)] = 0.0
image = 1. * fh / np.max(fh)
# find hough circles and detect two radii
edges = canny(image, sigma=3)
hough_radii = np.arange(minradius, maxradius, radstep)
# print hough_radii
hough_res = hough_circle(edges, hough_radii)
if hough_res.any():
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract num_peaks circles
peaks = peak_local_max(h, num_peaks=num_peaks)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
# Keep the most prominent num_circles circles
sorted_circles_idxs = np.argsort(accums)[::-1][:num_circles]
for idx in sorted_circles_idxs:
center_x, center_y = centers[idx]
allcenters.append(centers[idx])
allradii.append(radii[idx])
allaccums.append(accums[idx])
brightness = accums[idx]
lsurface = lsurface + brightness * np.exp(
-((xsurface - center_x) ** 2 + (ysurface - center_y) ** 2) / kernel_width ** 2)
lsurface = lsurface / lsurface.max()
# select most likely ROI center
roi_center = np.unravel_index(lsurface.argmax(), lsurface.shape)
# determine ROI radius
roi_x_radius = 0
roi_y_radius = 0
for idx in range(len(allcenters)):
xshift = np.abs(allcenters[idx][0] - roi_center[0])
yshift = np.abs(allcenters[idx][1] - roi_center[1])
if (xshift <= center_margin) & (yshift <= center_margin):
roi_x_radius = np.max((roi_x_radius, allradii[idx] + xshift))
roi_y_radius = np.max((roi_y_radius, allradii[idx] + yshift))
if roi_x_radius > 0 and roi_y_radius > 0:
roi_radii = roi_x_radius, roi_y_radius
else:
roi_radii = None
return roi_center, roi_radii
# Stddev-Hough Transform Based ROI Extraction
def extract_roi_stddev(data4D, pixel_spacing, minradius_mm=15, maxradius_mm=45, kernel_width=5,
center_margin=8, num_peaks=10, num_circles=20, radstep=2):
"""
Returns center and radii of ROI region in (i,j) format
"""
# Data shape:
# radius of the smallest and largest circles in mm estimated from the train set
# convert to pixel counts
pixel_spacing_X, pixel_spacing_Y, _,_ = pixel_spacing
minradius = int(minradius_mm / pixel_spacing_X)
maxradius = int(maxradius_mm / pixel_spacing_Y)
ximagesize = data4D.shape[0]
yimagesize = data4D.shape[1]
zslices = data4D.shape[2]
tframes = data4D.shape[3]
xsurface = np.tile(range(ximagesize), (yimagesize, 1)).T
ysurface = np.tile(range(yimagesize), (ximagesize, 1))
lsurface = np.zeros((ximagesize, yimagesize))
allcenters = []
allaccums = []
allradii = []
for slice in range(zslices):
ff1 = np.array([data4D[:,:,slice, t] for t in range(tframes)])
fh = np.std(ff1, axis=0)
fh[fh < 0.1 * np.max(fh)] = 0.0
image = 1. * fh / np.max(fh)
# find hough circles and detect two radii
edges = canny(image, sigma=3)
hough_radii = np.arange(minradius, maxradius, radstep)
# print hough_radii
hough_res = hough_circle(edges, hough_radii)
if hough_res.any():
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract num_peaks circles
peaks = peak_local_max(h, num_peaks=num_peaks)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
# Keep the most prominent num_circles circles
sorted_circles_idxs = np.argsort(accums)[::-1][:num_circles]
for idx in sorted_circles_idxs:
center_x, center_y = centers[idx]
allcenters.append(centers[idx])
allradii.append(radii[idx])
allaccums.append(accums[idx])
brightness = accums[idx]
lsurface = lsurface + brightness * np.exp(
-((xsurface - center_x) ** 2 + (ysurface - center_y) ** 2) / kernel_width ** 2)
lsurface = lsurface / lsurface.max()
# select most likely ROI center
roi_center = np.unravel_index(lsurface.argmax(), lsurface.shape)
# determine ROI radius
roi_x_radius = 0
roi_y_radius = 0
for idx in range(len(allcenters)):
xshift = np.abs(allcenters[idx][0] - roi_center[0])
yshift = np.abs(allcenters[idx][1] - roi_center[1])
if (xshift <= center_margin) & (yshift <= center_margin):
roi_x_radius = np.max((roi_x_radius, allradii[idx] + xshift))
roi_y_radius = np.max((roi_y_radius, allradii[idx] + yshift))
if roi_x_radius > 0 and roi_y_radius > 0:
roi_radii = roi_x_radius, roi_y_radius
else:
roi_radii = None
return roi_center, roi_radii
class Dataset(object):
def __init__(self, directory, subdir):
# type: (object, object) -> object
self.patient_data = {}
self.directory = directory
self.name = subdir
def _filename(self, file):
return os.path.join(self.directory, self.name, file)
def load_nii(self, img_path):
"""
Function to load a 'nii' or 'nii.gz' file, The function returns
everyting needed to save another 'nii' or 'nii.gz'
in the same dimensional space, i.e. the affine matrix and the header
Parameters
----------
img_path: string
String with the path of the 'nii' or 'nii.gz' image file name.
Returns
-------
Three element, the first is a numpy array of the image values,
the second is the affine transformation of the image, and the
last one is the header of the image.
"""
nimg = nib.load(self._filename(img_path))
return nimg.get_data(), nimg.affine, nimg.header
def read_patient_info_data(self):
"""
Reads patient data in the cfg file from patient folder
using Info.cfg
"""
print (self._filename('Info.cfg'))
with open(self._filename('Info.cfg')) as f_in:
for line in f_in:
l = line.rstrip().split(": ")
self.patient_data[l[0]] = l[1]
def read_patient_data(self, mode='train', roi_detect=True):
"""
Reads patient data in the cfg file and returns a dictionary and
extract End diastole and End Systole image from patient folder
using Info.cfg
"""
self.read_patient_info_data()
# Read patient Number
m = re.match("patient(\d{3})", self.name)
patient_No = int(m.group(1))
# Read Diastole frame Number
ED_frame_No = int(self.patient_data['ED'])
ed_img = "patient%03d_frame%02d.nii.gz" %(patient_No, ED_frame_No)
ed, affine, hdr = self.load_nii(ed_img)
# Read Systole frame Number
ES_frame_No = int(self.patient_data['ES'])
es_img = "patient%03d_frame%02d.nii.gz" %(patient_No, ES_frame_No)
es, _, _ = self.load_nii(es_img)
# Save Images:
self.patient_data['ED_VOL'] = ed
self.patient_data['ES_VOL'] = es
# Header Info for saving
header_info ={'affine':affine, 'hdr': hdr}
self.patient_data['header'] = header_info
if mode == 'reader':
# Read a particular volume number in 4D image
img_4d_name = "patient%03d_4d.nii.gz"%patient_No
# Load data
img_4D, _, hdr = self.load_nii(img_4d_name)
self.patient_data['4D'] = img_4D
ed_gt, _, _ = self.load_nii("patient%03d_frame%02d_gt.nii.gz" %(patient_No, ED_frame_No))
es_gt, _, _ = self.load_nii("patient%03d_frame%02d_gt.nii.gz" %(patient_No, ES_frame_No))
ed_lv, ed_rv, ed_myo = heart_metrics(ed_gt, hdr.get_zooms())
es_lv, es_rv, es_myo = heart_metrics(es_gt, hdr.get_zooms())
ef_lv = ejection_fraction(ed_lv, es_lv)
ef_rv = ejection_fraction(ed_rv, es_rv)
heart_param = {'EDV_LV': ed_lv, 'EDV_RV': ed_rv, 'ESV_LV': es_lv, 'ESV_RV': es_rv,
'ED_MYO': ed_myo, 'ES_MYO': es_myo, 'EF_LV': ef_lv, 'EF_RV': ef_rv}
self.patient_data['HP'] = heart_param
self.patient_data['ED_GT'] = ed_gt
self.patient_data['ES_GT'] = es_gt
return
if mode == 'train':
ed_gt, _, _ = self.load_nii("patient%03d_frame%02d_gt.nii.gz" %(patient_No, ED_frame_No))
es_gt, _, _ = self.load_nii("patient%03d_frame%02d_gt.nii.gz" %(patient_No, ES_frame_No))
ed_lv, ed_rv, ed_myo = heart_metrics(ed_gt, hdr.get_zooms())
es_lv, es_rv, es_myo = heart_metrics(es_gt, hdr.get_zooms())
ef_lv = ejection_fraction(ed_lv, es_lv)
ef_rv = ejection_fraction(ed_rv, es_rv)
heart_param = {'EDV_LV': ed_lv, 'EDV_RV': ed_rv, 'ESV_LV': es_lv, 'ESV_RV': es_rv,
'ED_MYO': ed_myo, 'ES_MYO': es_myo, 'EF_LV': ef_lv, 'EF_RV': ef_rv}
self.patient_data['HP'] = heart_param
self.patient_data['ED_GT'] = ed_gt
self.patient_data['ES_GT'] = es_gt
if mode == 'tester':
# Read a particular volume number in 4D image
img_4d_name = "patient%03d_4d.nii.gz"%patient_No
# Load data
img_4D, _, hdr = self.load_nii(img_4d_name)
self.patient_data['4D'] = img_4D
if roi_detect:
# Read a particular volume number in 4D image
img_4d_name = "patient%03d_4d.nii.gz"%patient_No
# Load data
img_4D, _, hdr = self.load_nii(img_4d_name)
c, r = extract_roi_stddev(img_4D, hdr.get_zooms())
self.patient_data['roi_center'], self.patient_data['roi_radii']=c,r
self.patient_data['4D'] = img_4D
# print c, r
# plot_roi(img_4D, c,r)
def convert_nii_np(data_path, mode, roi_detect):
"""
Prepare a dictionary of dataset and save it as numpy file
"""
patient_fulldata = OrderedDict()
print (data_path)
patient_folders = next(os.walk(data_path))[1]
for patient in tqdm(sorted(patient_folders)):
# print (patient)
dset = Dataset(data_path, patient)
dset.read_patient_data(mode=mode, roi_detect=roi_detect)
patient_fulldata[dset.name] = dset.patient_data
return patient_fulldata
if __name__ == '__main__':
start_time = time.time()
# Path to ACDC training database
complete_data_path = '../../ACDC_DataSet/training'
dest_path = '../../processed_acdc_dataset'
group_path = '../../processed_acdc_dataset/Patient_Groups'
# Training dataset
train_dataset = '../../processed_acdc_dataset/dataset/train_set'
validation_dataset = '../../processed_acdc_dataset/dataset/validation_set'
test_dataset = '../../processed_acdc_dataset/dataset/test_set'
out_path_train = '../../processed_acdc_dataset/pickled/full_data'
hdf5_out_path = '../../processed_acdc_dataset/hdf5_files'
#Final Test dataset
final_testing_dataset = '../../ACDC_DataSet/testing'
out_path_test = '../../processed_acdc_dataset/pickled/final_test'
# First perform stratified sampling
group_patient_cases(complete_data_path, dest_path)
generate_train_validate_test_set(group_path, dest_path)
print("---Time taken to stratify the dataset %s seconds ---" % (time.time() - start_time))
print ('ROI->ED->ES train dataset')
if not os.path.exists(out_path_train):
os.makedirs(out_path_train)
os.makedirs(out_path_test)
train_dataset = convert_nii_np(train_dataset, mode='train', roi_detect=True)
save_data(train_dataset, 'train_set.pkl', out_path_train)
print("---Processing Training dataset %s seconds ---" % (time.time() - start_time))
validation_dataset = convert_nii_np(validation_dataset, mode='train', roi_detect=True)
save_data(validation_dataset, 'validation_set.pkl', out_path_train)
print("---Processing Training dataset %s seconds ---" % (time.time() - start_time))
test_dataset = convert_nii_np(test_dataset, mode='train', roi_detect=True)
save_data(test_dataset, 'test_set.pkl', out_path_train)
print("---Processing Training dataset %s seconds ---" % (time.time() - start_time))
print ('ROI->ED->ES test dataset')
final_test_dataset = convert_nii_np(final_testing_dataset, mode='test', roi_detect=True)
save_data(final_test_dataset, 'final_testing_data.pkl', out_path_test)
print("---Processing final testing dataset %s seconds ---" % (time.time() - start_time))
# Generate 2D HDF5 files
modes = ['train_set', 'validation_set', 'test_set']
for mode in modes:
if os.path.exists(os.path.join(hdf5_out_path, mode)):
shutil.rmtree(os.path.join(hdf5_out_path, mode))
os.makedirs(os.path.join(hdf5_out_path, mode))
patient_data = load_pkl(os.path.join(out_path_train, mode+'.pkl'))
for patient_id in tqdm(patient_data.keys()):
# print (patient_id)
_id = patient_id[-3:]
n_slices = patient_data[patient_id]['ED_VOL'].shape[2]
# print (n_slices)
for slice in range(n_slices):
# ED frames
group = patient_data[patient_id]['Group']
slice_str ='_%02d_'%slice
roi_center = (patient_data[patient_id]['roi_center'][1], patient_data[patient_id]['roi_center'][0])
hp = h5py.File(os.path.join(hdf5_out_path, mode, 'P_'+_id+'_ED'+slice_str+group+'.hdf5'),'w')
hp.create_dataset('image', data=patient_data[patient_id]['ED_VOL'][:,:,slice].T)
hp.create_dataset('label', data=patient_data[patient_id]['ED_GT'][:,:,slice].T)
hp.create_dataset('roi_center', data=roi_center)
hp.create_dataset('roi_radii', data=patient_data[patient_id]['roi_radii'])
hp.create_dataset('pixel_spacing', data=patient_data[patient_id]['header']['hdr'].get_zooms())
hp.close()
# ES frames
hp = h5py.File(os.path.join(hdf5_out_path, mode, 'P_'+_id+'_ES'+slice_str+group+'.hdf5'),'w')
hp.create_dataset('image', data=patient_data[patient_id]['ES_VOL'][:,:,slice].T)
hp.create_dataset('label', data=patient_data[patient_id]['ES_GT'][:,:,slice].T)
hp.create_dataset('roi_center', data=roi_center)
hp.create_dataset('roi_radii', data=patient_data[patient_id]['roi_radii'])
hp.create_dataset('pixel_spacing', data=patient_data[patient_id]['header']['hdr'].get_zooms())
hp.close()
print("---Time taken to generate hdf5 files %s seconds ---" % (time.time() - start_time))
