Python SimpleITK.GetArrayFromImage() Examples

def preprocess_img(inputfile, output_preprocessed, zooms):
    img = nib.load(inputfile)
    data = img.get_data()
    affine = img.affine
    zoom = img.header.get_zooms()[:3]
    data, affine = reslice(data, affine, zoom, zooms, 1)
    data = np.squeeze(data)
    data = np.pad(data, [(0, 256 - len_) for len_ in data.shape], "constant")

    data_sub = data - gaussian_filter(data, sigma=1)
    img = sitk.GetImageFromArray(np.copy(data_sub))
    img = sitk.AdaptiveHistogramEqualization(img)
    data_clahe = sitk.GetArrayFromImage(img)[:, :, :, None]
    data = np.concatenate((data_clahe, data[:, :, :, None]), 3)
    data = (data - np.mean(data, (0, 1, 2))) / np.std(data, (0, 1, 2))
    assert data.ndim == 4, data.ndim
    assert np.allclose(np.mean(data, (0, 1, 2)), 0.), np.mean(data, (0, 1, 2))
    assert np.allclose(np.std(data, (0, 1, 2)), 1.), np.std(data, (0, 1, 2))
    data = np.float32(data)

    img = nib.Nifti1Image(data, affine)
    nib.save(img, output_preprocessed) 

def pp_patient(self, path):

        pid = path.split('/')[-1]
        img = sitk.ReadImage(os.path.join(path, '{}_ct_scan.nrrd'.format(pid)))
        img_arr = sitk.GetArrayFromImage(img)
        print('processing {} with GT(s) {}, spacing {} and img shape {}.'.format(
            pid, " and ".join(self.cf.gts_to_produce), img.GetSpacing(), img_arr.shape))
        img_arr = resample_array(img_arr, img.GetSpacing(), self.cf.target_spacing)
        img_arr = np.clip(img_arr, -1200, 600)
        #img_arr = (1200 + img_arr) / (600 + 1200) * 255  # a+x / (b-a) * (c-d) (c, d = new)
        img_arr = img_arr.astype(np.float32)
        img_arr = (img_arr - np.mean(img_arr)) / np.std(img_arr).astype('float16')

        df = pd.read_csv(os.path.join(self.cf.root_dir, 'characteristics.csv'), sep=';')
        df = df[df.PatientID == pid]

        np.save(os.path.join(self.cf.pp_dir, '{}_img.npy'.format(pid)), img_arr)
        if 'single_annotator' in self.cf.gts_to_produce or 'sa' in self.cf.gts_to_produce:
            self.produce_sa_gt(path, pid, df, img.GetSpacing(), img_arr.shape)
        if 'merged' in self.cf.gts_to_produce:
            self.produce_merged_gt(path, pid, df, img.GetSpacing(), img_arr.shape) 

def _add_slice_contribution(slice, coverage_sitk):

        #
        slice_sitk = sitk.Image(slice.sitk)
        spacing = np.array(slice_sitk.GetSpacing())
        spacing[-1] = slice.get_slice_thickness()
        slice_sitk.SetSpacing(spacing)

        contrib_nda = sitk.GetArrayFromImage(slice_sitk)
        contrib_nda[:] = 1
        contrib_sitk = sitk.GetImageFromArray(contrib_nda)
        contrib_sitk.CopyInformation(slice_sitk)

        coverage_sitk += sitk.Resample(
            contrib_sitk,
            coverage_sitk,
            sitk.Euler3DTransform(),
            sitk.sitkNearestNeighbor,
            0,
            coverage_sitk.GetPixelIDValue(),
        )

        return coverage_sitk 

def set_reconstruction(self, reconstruction):
        self._reconstruction = reconstruction

        # Extract information ready to use for itk image conversion operations
        self._reconstruction_shape = sitk.GetArrayFromImage(
            self._reconstruction.sitk).shape

        # Compute total amount of voxels of x:
        self._N_voxels_recon = np.array(
            self._reconstruction.sitk.GetSize()).prod()

    #
    # Set regularization parameter for Tikhonov regularization
    # \date       2017-07-25 15:15:54+0100
    #
    # \param      self   The object
    # \param      alpha  regularization parameter, scalar
    #
    # \return     { description_of_the_return_value }
    # 

def estimate(self):

### "segmented"
      color_region_growing = sitk.VectorConfidenceConnectedImageFilter()
      color_region_growing.SetNumberOfIterations(4)
      color_region_growing.SetMultiplier(5.3)
      color_region_growing.SetInitialNeighborhoodRadius(2)
      color_region_growing.SetReplaceValue(255)
      color_region_growing.AddSeed(self.seed_point)
      eyes_segmented = color_region_growing.Execute(self.input_image)

### "radius"
      distance_filter = sitk.SignedMaurerDistanceMapImageFilter()
      distance_filter.SetInsideIsPositive(True)
      distance_map = distance_filter.Execute(eyes_segmented)
      radius_estimate = np.amax(sitk.GetArrayFromImage(distance_map))

      return eyes_segmented,radius_estimate


### "overlay" 

def get_image(self):
    if self.frame == len(self.time_table):
      raise CrappyStop
    img_t = self.time_table[self.frame]
    img = sitk.GetArrayFromImage(sitk.ReadImage(self.img_dict[img_t]))
    t = time()
    delay = self.time_table[self.frame] - t + self.t0
    if delay > 0:
      if t > self.t0:
        sleep(delay)
      else:
        return img_t,img
    #else:
    #  print("Streamer is",-1000*delay,"ms late")
    self.frame += 1
    return img_t,img 

def preprocess_image(itk_image, is_seg=False, spacing_target=(1, 0.5, 0.5)):
    spacing = np.array(itk_image.GetSpacing())[[2, 1, 0]]
    image = sitk.GetArrayFromImage(itk_image).astype(float)

    assert len(image.shape) == 3, "The image has unsupported number of dimensions. Only 3D images are allowed"

    if not is_seg:
        if np.any([[i != j] for i, j in zip(spacing, spacing_target)]):
            image = resize_image(image, spacing, spacing_target).astype(np.float32)

        image -= image.mean()
        image /= image.std()
    else:
        new_shape = (int(np.round(spacing[0] / spacing_target[0] * float(image.shape[0]))),
                     int(np.round(spacing[1] / spacing_target[1] * float(image.shape[1]))),
                     int(np.round(spacing[2] / spacing_target[2] * float(image.shape[2]))))
        image = resize_segmentation(image, new_shape, 1)
    return image 

def load_nifty_volume_as_4d_array(filename):
    """Read a nifty image and return a dictionay storing data array, spacing and direction
    output['data_array'] 4d array with shape [C, D, H, W]
    output['spacing']    a list of spacing in z, y, x axis 
    output['direction']  a 3x3 matrix for direction
    """
    img_obj    = sitk.ReadImage(filename)
    data_array = sitk.GetArrayFromImage(img_obj)
    origin     = img_obj.GetOrigin()
    spacing    = img_obj.GetSpacing()
    direction  = img_obj.GetDirection()
    shape = data_array.shape
    if(len(shape) == 4):
        assert(shape[3] == 1) 
    elif(len(shape) == 3):
        data_array = np.expand_dims(data_array, axis = 0)
    else:
        raise ValueError("unsupported image dim: {0:}".format(len(shape)))
    output = {}
    output['data_array'] = data_array
    output['origin']     = origin
    output['spacing']    = (spacing[2], spacing[1], spacing[0])
    output['direction']  = direction
    return output 

def load_itk_image(filename):
    with open(filename) as f:
        contents = f.readlines()
        line = [k for k in contents if k.startswith('TransformMatrix')][0]
        transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
        transformM = np.round(transformM)
        if np.any( transformM!=np.array([1,0,0, 0, 1, 0, 0, 0, 1])):
            isflip = True
        else:
            isflip = False

    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)
     
    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))
     
    return numpyImage, numpyOrigin, numpySpacing, isflip 

def load_itk_image(filename):
    with open(filename) as f:
        contents = f.readlines()
        line = [k for k in contents if k.startswith('TransformMatrix')][0]
        transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
        transformM = np.round(transformM)
        if np.any( transformM!=np.array([1,0,0, 0, 1, 0, 0, 0, 1])):
            isflip = True
        else:
            isflip = False

    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)
    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))

    return numpyImage, numpyOrigin, numpySpacing, isflip 

def copy_BraTS_segmentation_and_convert_labels(in_file, out_file):
    # use this for segmentation only!!!
    # nnUNet wants the labels to be continuous. BraTS is 0, 1, 2, 4 -> we make that into 0, 1, 2, 3
    img = sitk.ReadImage(in_file)
    img_npy = sitk.GetArrayFromImage(img)

    uniques = np.unique(img_npy)
    for u in uniques:
        if u not in [0, 1, 2, 4]:
            raise RuntimeError('unexpected label')

    seg_new = np.zeros_like(img_npy)
    seg_new[img_npy == 4] = 3
    seg_new[img_npy == 2] = 1
    seg_new[img_npy == 1] = 2
    img_corr = sitk.GetImageFromArray(seg_new)
    img_corr.CopyInformation(img)
    sitk.WriteImage(img_corr, out_file) 

def create_tensor_image_from_itk_image(itk_image, dtype=th.float32, device='cpu'):

    # transform image in a unit direction
    image_dim = itk_image.GetDimension()
    if image_dim == 2:
        itk_image.SetDirection(sitk.VectorDouble([1, 0, 0, 1]))
    else:
        itk_image.SetDirection(sitk.VectorDouble([1, 0, 0, 0, 1, 0, 0, 0, 1]))

    image_spacing = itk_image.GetSpacing()
    image_origin = itk_image.GetOrigin()

    np_image = np.squeeze(sitk.GetArrayFromImage(itk_image))
    image_size = np_image.shape

    # adjust image spacing vector size if image contains empty dimension
    if len(image_size) != image_dim:
        image_spacing = image_spacing[0:len(image_size)]

    tensor_image = th.tensor(np_image, dtype=dtype, device=device).unsqueeze(0).unsqueeze(0)


    return Image(tensor_image, image_size, image_spacing, image_origin) 

def __init__ (self, path, uid=None):
        VolumeBase.__init__(self)
        self.uid = uid
        self.path = path
        if not os.path.exists(self.path):
            raise Exception('data not found for uid %s at %s' % (uid, self.path))
        pass
        #self.thumb_path = os.path.join(DATA_DIR, 'thumb', uid)
        # load path
        itkimage = itk.ReadImage(self.path)
        self.HU = (1.0, 0.0)
        self.images = itk.GetArrayFromImage(itkimage).astype(np.float32)
        #print type(self.images), self.images.dtype
        self.origin = list(reversed(itkimage.GetOrigin()))
        self.spacing = list(reversed(itkimage.GetSpacing()))
        self.view = AXIAL
        _, a, b = self.spacing
        #self.anno = LUNA_ANNO.get(uid, None)
        assert a == b
        # sanity check
        pass 

def load_itk_image(filename):
    with open(filename) as f:
        contents = f.readlines()
        line = [k for k in contents if k.startswith('TransformMatrix')][0]
        transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
        transformM = np.round(transformM)
        if np.any(transformM != np.array([1, 0, 0, 0, 1, 0, 0, 0, 1])):
            isflip = True
        else:
            isflip = False

    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)

    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))

    return numpyImage, numpyOrigin, numpySpacing, isflip 

def produceRandomlyTranslatedImage(image, label):
    sitkImage = sitk.GetImageFromArray(image, isVector=False)
    sitklabel = sitk.GetImageFromArray(label, isVector=False)

    itemindex = np.where(label > 0)
    randTrans = (0,np.random.randint(-np.min(itemindex[1])/2,(image.shape[1]-np.max(itemindex[1]))/2),np.random.randint(-np.min(itemindex[0])/2,(image.shape[0]-np.max(itemindex[0]))/2))
    translation = sitk.TranslationTransform(3, randTrans)

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(sitkImage)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(0)
    resampler.SetTransform(translation)

    outimgsitk = resampler.Execute(sitkImage)
    outlabsitk = resampler.Execute(sitklabel)

    outimg = sitk.GetArrayFromImage(outimgsitk)
    outimg = outimg.astype(dtype=float)

    outlbl = sitk.GetArrayFromImage(outlabsitk) > 0
    outlbl = outlbl.astype(dtype=float)

    return outimg, outlbl 

def load_itk_image(filename):
    with open(filename) as f:
        contents = f.readlines()
        line = [k for k in contents if k.startswith('TransformMatrix')][0]
        transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
        transformM = np.round(transformM)
        if np.any( transformM!=np.array([1,0,0, 0, 1, 0, 0, 0, 1])):
            isflip = True
        else:
            isflip = False

    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)
     
    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))
     
    return numpyImage, numpyOrigin, numpySpacing,isflip 

def read_images(self, path, affix):
        t1 = sitk.GetArrayFromImage(sitk.ReadImage(glob.glob(os.path.join(path, '*_t1' + affix))[0]))
        t1ce = sitk.GetArrayFromImage(sitk.ReadImage(glob.glob(os.path.join(path, '*_t1ce' + affix))[0]))
        t2 = sitk.GetArrayFromImage(sitk.ReadImage(glob.glob(os.path.join(path, '*_t2' + affix))[0]))
        flair = sitk.GetArrayFromImage(sitk.ReadImage(glob.glob(os.path.join(path, '*_flair' + affix))[0]))
        # scale to 0 to 1
        t1 = (t1 - np.amin(t1)) / (np.amax(t1) - np.amin(t1))
        t1ce = (t1ce - np.amin(t1ce)) / (np.amax(t1ce) - np.amin(t1ce))
        t2 = (t2 - np.amin(t2)) / (np.amax(t2) - np.amin(t2))
        flair = (flair - np.amin(flair)) / (np.amax(flair) - np.amin(flair))
        # scale to 0 mean, 1 std
        #t1 = (t1 - np.mean(t1)) / np.std(t1)
        #t1ce = (t1ce - np.mean(t1ce)) / np.std(t1ce)
        #t2 = (t2 - np.mean(t2)) / np.std(t2)
        #flair = (flair - np.mean(flair)) / np.std(flair)
        images = np.stack((t1, t1ce, t2, flair), axis=-1).astype(np.float32)
        return images 

def __init__ (self, uid):
        CaseBase.__init__(self)
        self.uid = uid
        self.path = os.path.join(LUNA_DIR_LOOKUP[uid], uid + '.mhd')
        if not os.path.exists(self.path):
            raise Exception('data not found for uid %s at %s' % (uid, self.path))
        pass
        #self.thumb_path = os.path.join(DATA_DIR, 'thumb', uid)
        # load path
        itkimage = itk.ReadImage(self.path)
        self.HU = (1.0, 0.0)
        self.images = itk.GetArrayFromImage(itkimage).astype(np.float32)
        #print type(self.images), self.images.dtype
        self.origin = list(reversed(itkimage.GetOrigin()))
        self.spacing = list(reversed(itkimage.GetSpacing()))
        self.view = AXIAL
        _, a, b = self.spacing
        self.anno = LUNA_ANNO.get(uid, None)
        assert a == b
        # sanity check
        pass 

def sitk_to_nib(image: sitk.Image) -> Tuple[np.ndarray, np.ndarray]:
    data = sitk.GetArrayFromImage(image).transpose()
    spacing = np.array(image.GetSpacing())
    direction = np.array(image.GetDirection())
    origin = image.GetOrigin()
    if len(direction) == 9:
        rotation = direction.reshape(3, 3)
    elif len(direction) == 4:  # ignore first dimension if 2D (1, 1, H, W)
        rotation_2d = direction.reshape(2, 2)
        rotation = np.eye(3)
        rotation[1:3, 1:3] = rotation_2d
        spacing = 1, *spacing
        origin = 0, *origin
    rotation = np.dot(FLIP_XY, rotation)
    rotation_zoom = rotation * spacing
    translation = np.dot(FLIP_XY, origin)
    affine = np.eye(4)
    affine[:3, :3] = rotation_zoom
    affine[:3, 3] = translation
    return data, affine 

def getLargestConnectedComponent_2D(itk_img):
    data = np.uint8(sitk.GetArrayFromImage(itk_img))
    for i in range(data.shape[0]):
        c,n = snd.label(data[i])
        sizes = snd.sum(data[i], c, range(n+1))
        mask_size = sizes < (max(sizes))
        remove_voxels = mask_size[c]
        c[remove_voxels] = 0
        c[np.where(c!=0)]=1
        data[i][np.where(c==0)] = 0
    return sitk.GetImageFromArray(data) 

def main():
    path='/shenlab/lab_stor3/dongnie/3T7T-Data/'
    saveto='/shenlab/lab_stor3/dongnie/3T7T-Data/'
   
    ids=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
    for ind in ids:
        datafilename='S%d/3t.hdr'%ind #provide a sample name of your filename of data here
        datafn=os.path.join(path,datafilename)
        labelfilename='S%d/7t.hdr'%ind  # provide a sample name of your filename of ground truth here
        labelfn=os.path.join(path,labelfilename)
        imgOrg=sitk.ReadImage(datafn)
        mrimg=sitk.GetArrayFromImage(imgOrg)
       	mu=np.mean(mrimg)
       	maxV, minV=np.percentile(mrimg, [99 ,25])
       	#mrimg=mrimg
       	mrimg=(mrimg-mu)/(maxV-minV)

        labelOrg=sitk.ReadImage(labelfn)
        labelimg=sitk.GetArrayFromImage(labelOrg) 
       	mu=np.mean(labelimg)
       	maxV, minV=np.percentile(labelimg, [99 ,25])
      	#labelimg=labelimg
       	labelimg=(labelimg-mu)/(maxV-minV)
        #you can do what you want here for for your label img
        #imgOrg=sitk.ReadImage(gtfn)
        #gtMat=sitk.GetArrayFromImage(imgOrg)
        prefn='s%d_3t.nii.gz'%ind
        preVol=sitk.GetImageFromArray(mrimg)
        sitk.WriteImage(preVol,prefn)
        outfn='s%d_7t.nii.gz'%ind
        preVol=sitk.GetImageFromArray(labelimg)
        sitk.WriteImage(preVol,outfn)
 
        fileID='%d'%ind
        rate=1
        #cubicCnt=cropCubic(mrimg,labelimg,fileID,dFA,step,rate)
        #print '# of patches is ', cubicCnt 

def getLargestConnectedComponent(itk_img):
    data = np.uint8(sitk.GetArrayFromImage(itk_img))
    c,n = snd.label(data)
    sizes = snd.sum(data, c, range(n+1))
    mask_size = sizes < (max(sizes))
    remove_voxels = mask_size[c]
    c[remove_voxels] = 0
    c[np.where(c!=0)]=1
    data[np.where(c==0)] = 0
    return sitk.GetImageFromArray(data) 

def apply_warp(x, warp_field, fill_mode='reflect',
               interpolator=sitk.sitkLinear,
               fill_constant=0):
    # Expand deformation field (and later the image), padding for the largest
    # deformation
    warp_field_arr = sitk.GetArrayFromImage(warp_field)
    max_deformation = np.max(np.abs(warp_field_arr))
    pad = np.ceil(max_deformation).astype(np.int32)
    warp_field_padded_arr = pad_image(warp_field_arr, pad_amount=pad,
                                      mode='nearest')
    warp_field_padded = sitk.GetImageFromArray(warp_field_padded_arr,
                                               isVector=True)

    # Warp x, one filter slice at a time
    x_warped = np.zeros(x.shape, dtype=np.float32)
    warp_filter = sitk.WarpImageFilter()
    warp_filter.SetInterpolator(interpolator)
    warp_filter.SetEdgePaddingValue(np.min(x).astype(np.double))
    for i, image in enumerate(x):
        x_tmp = np.zeros(image.shape, dtype=image.dtype)
        for j, channel in enumerate(image):
            image_padded = pad_image(channel, pad_amount=pad, mode=fill_mode,
                                     constant=fill_constant).T
            image_f = sitk.GetImageFromArray(image_padded)
            image_f_warped = warp_filter.Execute(image_f, warp_field_padded)
            image_warped = sitk.GetArrayFromImage(image_f_warped)
            x_tmp[j] = image_warped[pad:-pad, pad:-pad].T
        x_warped[i] = x_tmp
    return x_warped 

def apply_warp(x, warp_field, fill_mode='reflect',
               interpolator=None,
               fill_constant=0, rows_idx=1, cols_idx=2):
    '''Apply an spling warp field on an image'''
    import SimpleITK as sitk
    if interpolator is None:
        interpolator = sitk.sitkLinear
    # Expand deformation field (and later the image), padding for the largest
    # deformation
    warp_field_arr = sitk.GetArrayFromImage(warp_field)
    max_deformation = np.max(np.abs(warp_field_arr))
    pad = np.ceil(max_deformation).astype(np.int32)
    warp_field_padded_arr = pad_image(warp_field_arr, pad_amount=pad,
                                      mode='nearest')
    warp_field_padded = sitk.GetImageFromArray(warp_field_padded_arr,
                                               isVector=True)

    # Warp x, one filter slice at a time
    pattern = [el for el in range(0, x.ndim) if el not in [rows_idx, cols_idx]]
    pattern += [rows_idx, cols_idx]
    inv_pattern = [pattern.index(el) for el in range(x.ndim)]
    x = x.transpose(pattern)  # batch, channel, ...
    x_shape = list(x.shape)
    x = x.reshape([-1] + x_shape[2:])  # *, r, c
    warp_filter = sitk.WarpImageFilter()
    warp_filter.SetInterpolator(interpolator)
    warp_filter.SetEdgePaddingValue(np.min(x).astype(np.double))
    for i in range(x.shape[0]):
        bc_pad = pad_image(x[i], pad_amount=pad, mode=fill_mode,
                           constant=fill_constant).T
        bc_f = sitk.GetImageFromArray(bc_pad)
        bc_f_warped = warp_filter.Execute(bc_f, warp_field_padded)
        bc_warped = sitk.GetArrayFromImage(bc_f_warped)
        x[i] = bc_warped[pad:-pad, pad:-pad].T
    x = x.reshape(x_shape)  # unsquash
    x = x.transpose(inv_pattern)
    return x 

def analyze_lesion(self, pid, nodule_id):
        """print unique seg and counts of nodule nodule_id of patient pid.
        """
        nodule_id = nodule_id.lstrip("0")
        nodule_id_paths = [ii for ii in os.listdir(os.path.join(self.cf.raw_data_dir, pid)) if '.nii' in ii]
        nodule_id_paths = [ii for ii in nodule_id_paths if ii.split('_')[2].lstrip("0")==nodule_id]
        assert len(nodule_id_paths)==1
        nodule_path = nodule_id_paths[0]

        roi = sitk.ReadImage(os.path.join(self.cf.raw_data_dir, pid, nodule_path))
        roi_arr = sitk.GetArrayFromImage(roi).astype(np.uint8)

        print("pid {}, nodule {}, unique seg & counts: {}".format(pid, nodule_id, np.unique(roi_arr, return_counts=True)))
        return 

def test_leave_one_out(patient=0, flair=True, t1=True, full=True, first5=True, aug=True, verbose=False):
    if patient < 20: dir = 'raw/Utrecht/'
    elif patient < 40: dir = 'raw/Singapore/'
    else: dir = 'raw/GE3T/'
    dirs = os.listdir(dir)
    dirs.sort()
    dir += dirs[patient%20]
    FLAIR_image = sitk.ReadImage(dir + '/pre/FLAIR.nii.gz')
    T1_image = sitk.ReadImage(dir + '/pre/T1.nii.gz')
    FLAIR_array = sitk.GetArrayFromImage(FLAIR_image)
    T1_array = sitk.GetArrayFromImage(T1_image)
    if patient < 40: imgs_test = Utrecht_preprocessing(FLAIR_array, T1_array)
    else: imgs_test = GE3T_preprocessing(FLAIR_array, T1_array)
    if not flair: imgs_test = imgs_test[..., 1:2].copy()
    if not t1: imgs_test = imgs_test[..., 0:1].copy()
    img_shape = (rows_standard, cols_standard, flair+t1)
    model = get_unet(img_shape, first5)
    model_path = 'models/'
#if you want to test three models ensemble, just do like this: pred = (pred_1+pred_2+pred_3)/3
    model.load_weights(model_path + str(patient) + '.h5')
    pred = model.predict(imgs_test, batch_size=1, verbose=verbose)
    pred[pred > 0.5] = 1.
    pred[pred <= 0.5] = 0.
    if patient < 40: original_pred = Utrecht_postprocessing(FLAIR_array, pred)
    else: original_pred = GE3T_postprocessing(FLAIR_array, pred)
    filename_resultImage = model_path + str(patient) + '.nii.gz'
    sitk.WriteImage(sitk.GetImageFromArray(original_pred), filename_resultImage )
    filename_testImage = os.path.join(dir + '/wmh.nii.gz')
    testImage, resultImage = getImages(filename_testImage, filename_resultImage)
    dsc = getDSC(testImage, resultImage)
    avd = getAVD(testImage, resultImage) 
    h95 = getHausdorff(testImage, resultImage)
    recall, f1 = getLesionDetection(testImage, resultImage)
    return dsc, h95, avd, recall, f1 

def getLesionDetection(testImage, resultImage):    
    """Lesion detection metrics, both recall and F1."""
    
    # Connected components will give the background label 0, so subtract 1 from all results
    ccFilter = sitk.ConnectedComponentImageFilter()    
    ccFilter.SetFullyConnected(True)
    
    # Connected components on the test image, to determine the number of true WMH.
    # And to get the overlap between detected voxels and true WMH
    ccTest = ccFilter.Execute(testImage)    
    lResult = sitk.Multiply(ccTest, sitk.Cast(resultImage, sitk.sitkUInt32))
    
    ccTestArray = sitk.GetArrayFromImage(ccTest)
    lResultArray = sitk.GetArrayFromImage(lResult)
    
    # recall = (number of detected WMH) / (number of true WMH)
    recall = float(len(np.unique(lResultArray)) - 1) / (len(np.unique(ccTestArray)) - 1)
    
    # Connected components of results, to determine number of detected lesions
    ccResult = ccFilter.Execute(resultImage)
    ccResultArray = sitk.GetArrayFromImage(ccResult)
    
    # precision = (number of detected WMH) / (number of all detections)
    precision = float(len(np.unique(lResultArray)) - 1) / float(len(np.unique(ccResultArray)) - 1)
    
    f1 = 2.0 * (precision * recall) / (precision + recall)
    
    return recall, f1 

def getHausdorff(testImage, resultImage):
    """Compute the Hausdorff distance."""
        
    # Edge detection is done by ORIGINAL - ERODED, keeping the outer boundaries of lesions. Erosion is performed in 2D
    eTestImage   = sitk.BinaryErode(testImage, (1,1,0) )
    eResultImage = sitk.BinaryErode(resultImage, (1,1,0) )
    
    hTestImage   = sitk.Subtract(testImage, eTestImage)
    hResultImage = sitk.Subtract(resultImage, eResultImage)    
    
    hTestArray   = sitk.GetArrayFromImage(hTestImage)
    hResultArray = sitk.GetArrayFromImage(hResultImage)   
        
    # Convert voxel location to world coordinates. Use the coordinate system of the test image
    # np.nonzero   = elements of the boundary in numpy order (zyx)
    # np.flipud    = elements in xyz order
    # np.transpose = create tuples (x,y,z)
    # testImage.TransformIndexToPhysicalPoint converts (xyz) to world coordinates (in mm)
    testCoordinates   = np.apply_along_axis(testImage.TransformIndexToPhysicalPoint, 1, np.transpose( np.flipud( np.nonzero(hTestArray) )).astype(int) )
    resultCoordinates = np.apply_along_axis(testImage.TransformIndexToPhysicalPoint, 1, np.transpose( np.flipud( np.nonzero(hResultArray) )).astype(int) )
            
    # Use a kd-tree for fast spatial search
    def getDistancesFromAtoB(a, b):    
        kdTree = scipy.spatial.KDTree(a, leafsize=100)
        return kdTree.query(b, k=1, eps=0, p=2)[0]
    
    # Compute distances from test to result; and result to test
    dTestToResult = getDistancesFromAtoB(testCoordinates, resultCoordinates)
    dResultToTest = getDistancesFromAtoB(resultCoordinates, testCoordinates)    
    
    return max(np.percentile(dTestToResult, 95), np.percentile(dResultToTest, 95)) 

def apply_bet(img, bet, out_fname):
    img_itk = sitk.ReadImage(img)
    img_npy = sitk.GetArrayFromImage(img_itk)
    img_bet = sitk.GetArrayFromImage(sitk.ReadImage(bet))
    img_npy[img_bet == 0] = 0
    out = sitk.GetImageFromArray(img_npy)
    out.CopyInformation(img_itk)
    sitk.WriteImage(out, out_fname) 

def getDSC(testImage, resultImage):    
    """Compute the Dice Similarity Coefficient."""
    testArray   = sitk.GetArrayFromImage(testImage).flatten()
    resultArray = sitk.GetArrayFromImage(resultImage).flatten()
    
    # similarity = 1.0 - dissimilarity
    return 1.0 - scipy.spatial.distance.dice(testArray, resultArray)