--- a
+++ b/code/image_utils.py
@@ -0,0 +1,473 @@
+import os
+import numpy as np
+import nibabel as nib
+from numpy.linalg import inv
+
+def rescale_intensity(image, thres=(1.0, 99.0)):
+    """ Rescale the image intensity to the range of [0, 1] """
+    val_l, val_h = np.percentile(image, thres)
+    image2 = image
+    image2[image < val_l] = val_l
+    image2[image > val_h] = val_h
+    image2 = (image2.astype(np.float32) - val_l) / (val_h - val_l + 1e-6)
+    return image2
+
+
+def imagePreprocessing(originalNii, data_dir, atlas_dir):
+     
+    os.system('headertool '
+              '{0} '
+              '{1}/lvsa_.nii.gz '
+              '-reset'
+              .format(originalNii, data_dir))
+
+    os.system('temporalalign '
+              '{0}/temporal.gipl '
+              '{1}/lvsa_.nii.gz '
+              '{1}/lvsa_.nii.gz ' 
+              '-St1 0 '
+              '-St2 0 '
+              '-Et1 1 '
+              '-Et2 1'
+              .format(atlas_dir, data_dir))
+
+    os.system('autocontrast '
+              '{0}/lvsa_.nii.gz '
+              '{0}/lvsa_.nii.gz'
+              .format(data_dir))
+        
+    os.system('cardiacphasedetection '
+              '{0}/lvsa_.nii.gz '
+              '{0}/lvsa_ED.nii.gz '
+              '{0}/lvsa_ES.nii.gz'
+              .format(data_dir))
+    
+    print('  Image preprocessing is done ...')
+
+
+def clearBaseManbrance(data_dir, output_name):
+
+    nim = nib.load('{0}/{1}'.format(data_dir, output_name))
+    seg = nim.get_data()
+    if seg.ndim == 4:
+        seg = np.squeeze(seg, axis=-1).astype(np.int16)
+        
+    # find the last slice that has RV component 
+    rv = (seg == 4).astype(np.int16)   
+    for i in reversed(range(seg.shape[2])):
+        if np.sum(rv[:,:,i]) > 0: 
+            break
+    lastAppearSlice = i 
+    
+    # create a 2D mask that average the last 12 slice of RV which is robust to noise
+    slicesUsed = 12
+    tmp = np.zeros((seg.shape[0], seg.shape[1]), dtype=np.int16)
+    for j in range(slicesUsed):
+        tmp = tmp + rv[:,:,lastAppearSlice-j]
+    tmp = (tmp > 0).astype(np.int16)   
+    
+    # create a volume mask used to remove the base manbrance caused by RV myocardium
+    mask = np.ones((seg.shape[2], seg.shape[0], seg.shape[1]), dtype=np.int16)
+    mask[lastAppearSlice-slicesUsed:,:,:] = 1 - tmp
+    mask = np.transpose(mask, axes=(1, 2, 0))
+    seg = np.multiply(seg, mask)
+    
+    # fill the gap/hole of RV 
+    flat = False
+    if flat:
+        seg = np.transpose(seg, axes=(2, 0, 1))
+        seg[lastAppearSlice-slicesUsed:lastAppearSlice-3,:,:] = seg[lastAppearSlice-slicesUsed:lastAppearSlice-3,:,:] + 4*tmp
+        seg = np.transpose(seg, axes=(1, 2, 0))
+    else:
+        seg[:,:,lastAppearSlice-slicesUsed:] = seg[:,:,lastAppearSlice-slicesUsed:] + 4*rv[:,:,lastAppearSlice-slicesUsed:]
+    
+    ###############################################################################
+    # find the last slice that has LV component 
+    lv = (seg == 1).astype(np.int16)   
+    for i in reversed(range(seg.shape[2])):
+        if np.sum(lv[:,:,i]) > 0: 
+            break
+    lastAppearSlice = i 
+    
+    # create a 2D mask that average the last 5 slice of RV which is robust to noise
+    slicesUsed = 10
+    tmp = np.zeros((seg.shape[0], seg.shape[1]), dtype=np.int16)
+    for j in range(slicesUsed):
+        tmp = tmp + lv[:,:,lastAppearSlice-j]
+    tmp = (tmp > 0).astype(np.int16)   
+    
+    # create a volume mask used to remove the manbrance caused by LV myocardium
+    mask = np.ones((seg.shape[2], seg.shape[0], seg.shape[1]), dtype=np.int16)
+    mask[lastAppearSlice-slicesUsed:,:,:] = 1 - tmp
+    mask = np.transpose(mask, axes=(1, 2, 0))
+    seg = np.multiply(seg, mask)
+    
+    # fill the gap/hole of LV 
+    flat = True
+    if flat:
+        seg = np.transpose(seg, axes=(2, 0, 1))
+        seg[lastAppearSlice-slicesUsed:lastAppearSlice-1,:,:] = seg[lastAppearSlice-slicesUsed:lastAppearSlice-1,:,:] + tmp
+        seg = np.transpose(seg, axes=(1, 2, 0))
+    else:
+        seg[:,:,lastAppearSlice-slicesUsed:] = seg[:,:,lastAppearSlice-slicesUsed:] + lv[:,:,lastAppearSlice-slicesUsed:]
+
+    # save the result
+    nim2 = nib.Nifti1Image(seg, nim.affine)
+    nim2.header['pixdim'] = nim.header['pixdim']
+    nib.save(nim2, '{0}/{1}'.format(data_dir, output_name))
+
+
+def removeSegsAboveBase(data_dir, output_name):
+
+    # Read segmentations
+    nim = nib.load('{0}/{1}'.format(data_dir, output_name))
+    seg = nim.get_data()
+    if seg.ndim == 4:
+        seg = np.squeeze(seg, axis=-1).astype(np.int16)
+    os.system('vtk2txt {0}/landmarks.vtk {0}/landmarks.txt'.format(data_dir)) 
+    
+    # convert txt file into matrix
+    file = open('{0}/landmarks.txt'.format(data_dir), 'r') 
+    A = file.read()
+    tmp = np.zeros(18)  
+    i, c = 0, 0
+    for p in range(len(A)):
+        if A[p] == ' ' or A[p] == '\n':
+            tmp[i] = np.float32(A[c:p])
+            i = i + 1
+            c = p + 1;
+    tmp = np.reshape(tmp,(6,3))        
+    landmarks = np.ones((6,4)) 
+    landmarks[:,:-1] = tmp
+    landmarks = np.transpose(landmarks).astype(np.float32)           
+    os.system('rm {0}/landmarks.txt'.format(data_dir))    
+       
+    # map world coordinate system to image pixel position
+    pixelsPositions = np.ceil(np.dot(inv(nim.affine),landmarks)).astype(np.int16) 
+    pixelsPositions = np.delete(pixelsPositions, (-1), axis=0)
+    pixelsPositions = np.transpose(pixelsPositions).astype(np.int16) 
+    
+    if output_name[9:11]=='ED':
+        seg[:,:,pixelsPositions[5,2]+1:] = 0
+    else:
+        seg[:,:,pixelsPositions[5,2]:] = 0
+    
+    nim2 = nib.Nifti1Image(seg, nim.affine)
+    nim2.header['pixdim'] = nim.header['pixdim']
+    nib.save(nim2, '{0}/{1}'.format(data_dir, output_name))
+    
+
+def formHighResolutionImg(subject_dir, fr): 
+    
+    os.system('resample ' 
+              '{0}/lvsa_{1}.nii.gz '
+              '{0}/lvsa_SR_{1}.nii.gz '
+              '-size 1.25 1.25 2'
+              .format(subject_dir, fr))
+        
+#    os.system('enlarge_image '
+#              '{0}/lvsa_SR_{1}.nii.gz '
+#              '{0}/lvsa_SR_{1}.nii.gz '
+#              '-z 20 '
+#              '-value 0'
+#              .format(subject_dir, fr))
+
+    
+def convertImageSegment(data_dir, fr):
+   
+    os.system('convert '
+              '{0}/seg_lvsa_SR_{1}.nii.gz '
+              '{0}/PHsegmentation_{1}.gipl'
+              .format(data_dir, fr))
+    
+    os.system('cp '
+              '{0}/lvsa_SR_{1}.nii.gz '
+              '{0}/lvsa_{1}_enlarged_SR.nii.gz'
+              .format(data_dir, fr))
+    
+    
+def outputVolumes(subject_dir, data_dir, subject, fr):
+     
+    os.system('rm '
+              '{0}/{1}_{2}.txt'
+              .format(data_dir, subject, fr))
+    
+    os.system('cardiacvolumecount '
+              '{0}/PHsegmentation_{3}.gipl 1 '
+              '-output '
+              '{1}/{2}_{3}.txt'
+              .format(subject_dir, data_dir, subject, fr))
+    
+    os.system('cardiacvolumecount '
+              '{0}/PHsegmentation_{3}.gipl 2 '
+              '-output '
+              '{1}/{2}_{3}.txt '
+              '-scale 1.05'
+              .format(subject_dir, data_dir, subject, fr))
+    
+    os.system('cardiacvolumecount '
+              '{0}/PHsegmentation_{3}.gipl 4 '
+              '-output '
+              '{1}/{2}_{3}.txt'
+              .format(subject_dir, data_dir, subject, fr))
+    
+    os.system('cardiacvolumecount '
+              '{0}/PHsegmentation_{3}.gipl 3 '
+              '-output '
+              '{1}/{2}_{3}.txt '
+              '-scale 1.05'
+              .format(subject_dir, data_dir, subject, fr))
+    
+    
+def moveVolumes(subject_dir, sizes_dir, fr):
+     
+    os.system('cp '
+              '{0}/lvsa_{2}.nii.gz '
+              '{1}/lvsa_{2}.nii.gz'
+              .format(subject_dir, sizes_dir, fr))
+    
+    os.system('rm '
+              '{0}/lvsa_{1}.nii.gz '
+              .format(subject_dir, fr))
+    
+    os.system('cp '
+              '{0}/seg_lvsa_{2}.nii.gz '
+              '{1}/2D_segmentation_{2}.nii.gz'
+              .format(subject_dir, sizes_dir, fr))
+    
+    os.system('rm '
+              '{0}/seg_lvsa_{1}.nii.gz '
+              .format(subject_dir, fr))
+    
+    os.system('cp '
+              '{0}/lvsa_SR_{2}.nii.gz '
+              '{1}/lvsa_{2}_SR.nii.gz'
+              .format(subject_dir, sizes_dir, fr))
+    
+    os.system('rm '
+              '{0}/lvsa_SR_{1}.nii.gz '
+              .format(subject_dir, fr))
+    
+    os.system('cp '
+              '{0}/seg_lvsa_SR_{2}.nii.gz '
+              '{1}/3D_segmentation_{2}.nii.gz'
+              .format(subject_dir, sizes_dir, fr))
+    
+    os.system('rm '
+              '{0}/seg_lvsa_SR_{1}.nii.gz '
+              .format(subject_dir, fr))
+
+    
+def refineFusionResults(data_dir, output_name, alfa): 
+    
+    ##########################################
+    os.system('binarize '
+              '{0}/{1} '
+              '{0}/tmps/hrt.nii.gz '
+              '1 4 255 0'
+              .format(data_dir, output_name))
+    
+    os.system('blur '
+              '{0}/tmps/hrt.nii.gz '
+              '{0}/tmps/hrt.nii.gz '
+              '{1}'
+              .format(data_dir, alfa))
+    
+    os.system('threshold '
+              '{0}/tmps/hrt.nii.gz '
+              '{0}/tmps/hrt.nii.gz '
+              '130'
+              .format(data_dir))
+    
+    ##########################################
+    os.system('binarize '
+              '{0}/{1} '
+              '{0}/tmps/rvendo.nii.gz '
+              '4 4 255 0'
+              .format(data_dir, output_name))
+    
+    os.system('blur '
+              '{0}/tmps/rvendo.nii.gz '
+              '{0}/tmps/rvendo.nii.gz '
+              '{1}'
+              .format(data_dir, alfa))
+    
+    os.system('threshold '
+              '{0}/tmps/rvendo.nii.gz '
+              '{0}/tmps/rvendo.nii.gz '
+              '130'
+              .format(data_dir))
+    
+    ##########################################
+    os.system('binarize '
+              '{0}/{1} '
+              '{0}/tmps/lvepi.nii.gz '
+              '1 2 255 0'
+              .format(data_dir, output_name))
+       
+    os.system('blur '
+              '{0}/tmps/lvepi.nii.gz '
+              '{0}/tmps/lvepi.nii.gz '
+              '{1}'
+              .format(data_dir, alfa))
+    
+    os.system('threshold '
+              '{0}/tmps/lvepi.nii.gz '
+              '{0}/tmps/lvepi.nii.gz '
+              '115'
+              .format(data_dir))
+    
+    ##########################################
+    os.system('binarize '
+              '{0}/{1} '
+              '{0}/tmps/lvendo.nii.gz '
+              '1 1 255 0'
+              .format(data_dir, output_name))
+    
+    os.system('blur '
+              '{0}/tmps/lvendo.nii.gz '
+              '{0}/tmps/lvendo.nii.gz '
+              '{1}'
+              .format(data_dir, alfa))
+    
+    os.system('threshold '
+              '{0}/tmps/lvendo.nii.gz '
+              '{0}/tmps/lvendo.nii.gz '
+              '130'
+              .format(data_dir))
+       
+    ##########################################
+    os.system('padding '
+              '{0}/tmps/hrt.nii.gz '
+              '{0}/tmps/hrt.nii.gz '
+              '{0}/tmps/hrt.nii.gz '
+              '1 3'
+              .format(data_dir))
+    
+    os.system('padding '
+              '{0}/tmps/hrt.nii.gz '
+              '{0}/tmps/rvendo.nii.gz '
+              '{0}/tmps/rvendo.nii.gz '
+              '1 4'
+              .format(data_dir))
+    
+    os.system('padding '
+              '{0}/tmps/rvendo.nii.gz '
+              '{0}/tmps/lvepi.nii.gz '
+              '{0}/tmps/lvepi.nii.gz '
+              '1 2'
+              .format(data_dir))
+    
+    os.system('padding '
+              '{0}/tmps/lvepi.nii.gz '
+              '{0}/tmps/lvendo.nii.gz '
+              '{0}/{1} '
+              '1 1'
+              .format(data_dir, output_name))
+    
+    
+def allAtlasShapeSelection(dataset_dir):
+   
+    atlases_list, landmarks_list = {}, {}
+    
+    for fr in ['ED', 'ES']:
+            
+        atlases_list[fr], landmarks_list[fr] = [], [] 
+        i = 0
+        for atlas in sorted(os.listdir(dataset_dir)): 
+            
+            atlas_dir = os.path.join(dataset_dir, atlas)
+            
+            if not os.path.isdir(atlas_dir):
+                
+                print('  {0} is not a valid atlas directory, Discard'.format(atlas_dir))
+                
+                continue 
+            
+            atlas_3D_shape = '{0}/PHsegmentation_{1}.nii.gz'.format(atlas_dir, fr)
+           
+            landmarks = '{0}/landmarks.vtk'.format(atlas_dir)
+        
+            if i < 400:
+                if os.path.exists(atlas_3D_shape) or os.path.exists(landmarks):
+                    
+                    atlases_list[fr] += [atlas_3D_shape]
+                
+                    landmarks_list[fr] += [landmarks]
+            else:
+                print(atlases_list)
+                break
+            
+            i = i + 1
+                 
+    return atlases_list, landmarks_list
+
+
+def topSimilarAtlasShapeSelection(atlases, atlas_landmarks, subject_landmarks, tmps_dir, dofs_dir, DLSeg, param_dir, topSimilarAltasNo):            
+    
+    landmarks = False
+    
+    nmi = []
+    
+    topSimilarAtlases_list = []
+    
+    atlasNo = len(atlases)
+    
+    os.system('rm {0}/shapenmi*.txt'.format(tmps_dir))
+    
+    for i in range(atlasNo):
+        
+        if landmarks:
+            
+            os.system('pareg '
+                      '{0} '
+                      '{1} '
+                      '-dofout {2}/shapelandmarks_{3}.dof.gz'
+                      .format(subject_landmarks, atlas_landmarks[i], dofs_dir, i))
+        else: 
+            
+            os.system('areg '
+                      '{0} '
+                      '{1} '
+                      '-dofout {2}/shapelandmarks_{4}.dof.gz '
+                      '-parin {3}/segareg.txt'
+                      .format(DLSeg, atlases[i], dofs_dir, param_dir, i))
+                      
+        os.system('cardiacimageevaluation '
+                  '{0} '
+                  '{1} '
+                  '-nbins_x 64 '
+                  '-nbins_y 64 '
+                  '-dofin {2}/shapelandmarks_{4}.dof.gz '
+                  '-output {3}/shapenmi_{4}.txt'
+                  .format(DLSeg, atlases[i], dofs_dir, tmps_dir, i))
+                        
+        if os.path.exists('{0}/shapenmi_{1}.txt'.format(tmps_dir, i)):
+
+            similarities = np.genfromtxt('{0}/shapenmi_{1}.txt'.format(tmps_dir, i))
+            nmi += [similarities[3]]
+            
+        else:
+            nmi += [0]
+    
+    if topSimilarAltasNo < atlasNo:
+        
+        sortedIndexes = np.array(nmi).argsort()[::-1]  
+        
+        savedInd = np.zeros(topSimilarAltasNo, dtype=int)
+        
+        for i in range(topSimilarAltasNo):
+            
+            topSimilarAtlases_list += [atlases[sortedIndexes[i]]]
+            
+            savedInd[i] = sortedIndexes[i]
+    else:
+        
+        savedInd = np.zeros(atlasNo, dtype=int)
+        
+        topSimilarAtlases_list = atlases
+        
+        savedInd = np.arange(atlasNo)
+        
+    return topSimilarAtlases_list, savedInd    
+