Python SimpleITK.ReadImage() Examples

def N4BiasFieldCorrection(src_path, dst_path):
        '''
        This function carry out BiasFieldCorrection for the files in a specific directory
        :param src_path: path of the source file
        :param dst_path: path of the target file
        :return:
        '''
        print("N4 bias correction runs.")
        inputImage = sitk.ReadImage(src_path)

        maskImage = sitk.OtsuThreshold(inputImage, 0, 1, 200)
        sitk.WriteImage(maskImage, dst_path)

        inputImage = sitk.Cast(inputImage, sitk.sitkFloat32)

        corrector = sitk.N4BiasFieldCorrectionImageFilter()

        # corrector.SetMaximumNumberOfIterations(10)

        output = corrector.Execute(inputImage, maskImage)
        sitk.WriteImage(output, dst_path)
        print("Finished N4 Bias Field Correction.....")

    # normalize the data(zero mean and unit variance) 

def __GetFeatureValuesEachModality(self, data_path, roi_path, key_name):
        if self.__is_ignore_tolerence:
            image = sitk.ReadImage(data_path)
            roi_image = sitk.ReadImage(roi_path)
            roi_image.CopyInformation(image)
            result = self.extractor.execute(image, roi_image)
        else:
            result = self.extractor.execute(data_path, roi_path)

        feature_names = []
        feature_values = []
        for feature_name, feature_value in zip(list(result.keys()), list(result.values())):
            if self.__is_ignore_diagnostic and 'diagnostics' in feature_name:
                continue
            feature_names.append(key_name + '_' + feature_name)
            feature_values.append(feature_value)
        return feature_names, feature_values 

def saveImagesFromDM(dataMatrix, outputPrefix, referenceImName):
    """Save 3D images from data matrix."""
    im_ref = sitk.ReadImage(referenceImName)
    im_ref_array = sitk.GetArrayFromImage(im_ref)  # get numpy array
    z_dim, x_dim, y_dim = im_ref_array.shape  # get 3D volume shape
    num_of_data = dataMatrix.shape[1]
    list_files = []
    for i in range(num_of_data):
        im = np.array(dataMatrix[:, i]).reshape(z_dim, x_dim, y_dim)
        img = sitk.GetImageFromArray(im)
        img.SetOrigin(im_ref.GetOrigin())
        img.SetSpacing(im_ref.GetSpacing())
        img.SetDirection(im_ref.GetDirection())
        fn = outputPrefix + str(i) + '.nrrd'
        list_files.append(fn)
        sitk.WriteImage(img, fn, True)
    del im_ref, im_ref_array
    return list_files 

def eyesize_noisy_test():

  downloader = imagedownloader.ImageDownloader()
  estimator = eyesize.EyeSize()

  image_name = "TralitusSaltrator.jpg"

  downloader.set_figshare_id("1066744")
  downloader.set_image_name(image_name)
  downloader.download()

  input_image = sitk.ReadImage(image_name)

  estimator.set_image(input_image)
  estimator.set_seed_point([204,400])

  eyes_segmented,radius_estimate = estimator.estimate()

  sitk.WriteImage(eyes_segmented,'SegmentedEye.png')

  expected_value = 85
  assert radius_estimate == expected_value, \
    "Problem with estimating radius: current: %s, expected: %s" % (radius_estimate, expected_value) 

def eyesize_basic_test():

  downloader = imagedownloader.ImageDownloader()
  estimator = eyesize.EyeSize()

  image_name = "TralitusSaltrator.jpg"

  downloader.set_figshare_id("1066744")
  downloader.set_image_name(image_name)
  downloader.download()

  input_image = sitk.ReadImage(image_name)

  estimator.set_image(input_image)
  estimator.set_seed_point([204,400])

  eyes_segmented,radius_estimate = estimator.estimate()

  sitk.WriteImage(eyes_segmented,'SegmentedEye.png')

  assert radius_estimate == 85 

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 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 getImages(testFilename, resultFilename):
    """Return the test and result images, thresholded and non-WMH masked."""
    testImage   = sitk.ReadImage(testFilename)
    resultImage = sitk.ReadImage(resultFilename)
    assert testImage.GetSize() == resultImage.GetSize()
    
    # Get meta data from the test-image, needed for some sitk methods that check this
    resultImage.CopyInformation(testImage)
    
    # Remove non-WMH from the test and result images, since we don't evaluate on that
    maskedTestImage = sitk.BinaryThreshold(testImage, 0.5,  1.5, 1, 0) # WMH == 1    
    nonWMHImage     = sitk.BinaryThreshold(testImage, 1.5,  2.5, 0, 1) # non-WMH == 2
    maskedResultImage = sitk.Mask(resultImage, nonWMHImage)
    
    # Convert to binary mask
    if 'integer' in maskedResultImage.GetPixelIDTypeAsString():
        bResultImage = sitk.BinaryThreshold(maskedResultImage, 1, 1000, 1, 0)
    else:
        bResultImage = sitk.BinaryThreshold(maskedResultImage, 0.5, 1000, 1, 0)
        
    return maskedTestImage, bResultImage 

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 __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 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 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 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_and_preprocess(mri_file):
    images = {}
    # t1
    images["T1"] = sitk.ReadImage(mri_file)

    properties_dict = {
        "spacing": images["T1"].GetSpacing(),
        "direction": images["T1"].GetDirection(),
        "size": images["T1"].GetSize(),
        "origin": images["T1"].GetOrigin()
    }

    for k in images.keys():
        images[k] = preprocess_image(images[k], is_seg=False, spacing_target=(1.5, 1.5, 1.5))

    properties_dict['size_before_cropping'] = images["T1"].shape

    imgs = []
    for seq in ['T1']:
        imgs.append(images[seq][None])
    all_data = np.vstack(imgs)
    print("image shape after preprocessing: ", str(all_data[0].shape))
    return all_data, properties_dict 

def showSlice(dataMatrix, title, color, subplotRow, referenceImName, slice_nr=-1):
    """Display a 2D slice within a given figure."""
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print "ShowSlice not supported - matplotlib not available"
        return
    im_ref = sitk.ReadImage(referenceImName)
    im_ref_array = sitk.GetArrayFromImage(im_ref)  # get numpy array
    z_dim, x_dim, y_dim = im_ref_array.shape  # get 3D volume shape
    if slice_nr == -1:
        slice_nr = z_dim / 2
    num_of_data = dataMatrix.shape[1]
    for i in range(num_of_data):
        plt.subplot2grid((3, num_of_data), (subplotRow, i))
        im = np.array(dataMatrix[:, i]).reshape(z_dim, x_dim, y_dim)
        implot = plt.imshow(np.fliplr(np.flipud(im[slice_nr, :, :])), color)
        plt.axis('off')
        plt.title(title + ' ' + str(i))
        # plt.colorbar()
    del im_ref, im_ref_array
    return 

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 showImageMidSlice(reference_im_fn, size_x=15, size_y=5):
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print "showImageMidSlice not supported - matplotlib not available"
        return
    im_ref = sitk.ReadImage(reference_im_fn)  # image in SITK format
    im_ref_array = sitk.GetArrayFromImage(im_ref)  # get numpy array
    z_dim, x_dim, y_dim = im_ref_array.shape  # get 3D volume shape
    # display reference image
    fig = plt.figure(figsize=(size_x,size_y))
    plt.subplot(131)
    implot = plt.imshow(np.flipud(im_ref_array[z_dim/2,:,:]),plt.cm.gray)
    plt.subplot(132)
    implot = plt.imshow(np.flipud(im_ref_array[:,x_dim/2,:]),plt.cm.gray)
    plt.subplot(133)
    implot = plt.imshow(np.flipud(im_ref_array[:,:,y_dim/2]),plt.cm.gray)
    plt.axis('off')
    plt.title('healthy atlas')
    fig.clf()
    del im_ref, im_ref_array
    return 

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 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_mri_series(mri_path):
    img_itk = sitk.ReadImage(mri_path)
    img_np = sitk.GetArrayFromImage(img_itk)
    # SimpleITK read image as (z, y, x), need to be transposed to (x, y, z)
    img_np = img_np.transpose((2, 1, 0)).astype('float')

    return torch.tensor(img_np).float() 

def eyesize_cleanplate_test():

  downloader = imagedownloader.ImageDownloader()
  estimator = eyesize.EyeSize()

  image_name = "TralitusSaltrator.jpg"

  inputs_dir = os.path.join(os.getcwd(), "inputs")
  outputs_dir = os.path.join(os.getcwd(), "outputs")
  
  image_path = os.path.join(inputs_dir, image_name)
  
  # Clear inputs folder
  if os.path.isdir(inputs_dir):
    shutil.rmtree(inputs_dir)
  os.mkdir(inputs_dir)
  
  # Clear outputs folder
  if os.path.isdir(outputs_dir):
    shutil.rmtree(outputs_dir)
  os.mkdir(outputs_dir)

  downloader.set_figshare_id("1066744")
  downloader.set_image_name(image_path)
  downloader.download()

  input_image = sitk.ReadImage(image_path)

  estimator.set_image(input_image)
  estimator.set_seed_point([204,400])

  eyes_segmented,radius_estimate = estimator.estimate()

  sitk.WriteImage(eyes_segmented,'SegmentedEye.png')

  expected_value = 85
  assert radius_estimate == expected_value, \
    "Problem with estimating radius: current: %s, expected: %s" % (radius_estimate, expected_value) 

def load_itk(filename):
    """
    load mhd files and normalization 0-255
    :param filename:
    :return:
    """
    rescalFilt = sitk.RescaleIntensityImageFilter()
    rescalFilt.SetOutputMaximum(255)
    rescalFilt.SetOutputMinimum(0)
    # Reads the image using SimpleITK
    itkimage = rescalFilt.Execute(sitk.Cast(sitk.ReadImage(filename), sitk.sitkFloat32))
    return itkimage 

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 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
        
        fileID='%d'%ind
        rate=1
        cubicCnt=cropCubic(mrimg,labelimg,fileID,dFA,step,rate)
        print '# of patches is ', cubicCnt 

def main():

    input_parser = InputArgparser(
        description="Multiply images. "
        "Pixel type is determined by first given image.",
    )

    input_parser.add_filenames(required=True)
    input_parser.add_output(required=True)
    input_parser.add_verbose(default=0)

    args = input_parser.parse_args()
    input_parser.print_arguments(args)

    if len(args.filenames) < 2:
        raise IOError("At least two images must be provided")

    out_sitk = sitk.ReadImage(args.filenames[0])
    for f in args.filenames[1:]:
        im_sitk = sitk.Cast(sitk.ReadImage(f), out_sitk.GetPixelIDValue())
        out_sitk = out_sitk * im_sitk

    dw.DataWriter.write_image(out_sitk, args.output)

    if args.verbose:
        args.filenames.insert(0, args.output)
        ph.show_niftis(args.filenames) 

def test_simulate_stacks_from_slices(self):

        cmd_args = []
        cmd_args.append("--filenames %s" % self.path_to_file)
        cmd_args.append("--dir-input-mc %s" %
                        os.path.join(
                            self.dir_data,
                            "recon_projections",
                            "motion_correction"))
        cmd_args.append("--reconstruction %s" % self.path_to_recon)
        cmd_args.append("--reconstruction-mask %s" % self.path_to_recon_mask)
        cmd_args.append("--copy-data 1")
        cmd_args.append("--suffix-mask %s" % self.suffix_mask)
        cmd_args.append("--dir-output %s" % self.dir_output)

        exe = os.path.abspath(simulate_stacks_from_reconstruction.__file__)
        cmd = "python %s %s" % (exe, (" ").join(cmd_args))
        self.assertEqual(ph.execute_command(cmd), 0)

        path_orig = os.path.join(self.dir_output, "%s.nii.gz" % self.filename)
        path_sim = os.path.join(
            self.dir_output, "Simulated_%s.nii.gz" % self.filename)

        path_orig_ref = os.path.join(self.dir_data,
                                     "recon_projections",
                                     "slices",
                                     "%s.nii.gz" % self.filename)
        path_sim_ref = os.path.join(self.dir_data,
                                    "recon_projections",
                                    "slices",
                                    "Simulated_%s.nii.gz" % self.filename)

        for res, ref in zip(
                [path_orig, path_sim], [path_orig_ref, path_sim_ref]):
            res_sitk = sitk.ReadImage(res)
            ref_sitk = sitk.ReadImage(ref)

            nda_diff = np.nan_to_num(
                sitk.GetArrayFromImage(res_sitk - ref_sitk))
            self.assertAlmostEqual(np.linalg.norm(
                nda_diff), 0, places=self.precision) 

def main(argv=None):
    if argv is None:
        argv = sys.argv

    if len(argv) != 5:
        print "Usage: python %s <CheckerboardImage> <OutlierFraction> <CorruptedImage> <LowRankImage>" % sys.argv[0]
        sys.exit(1)

    # outlier fraction
    p = float(sys.argv[2])

    # read image
    I = sitk.ReadImage(argv[1])
    # data for processing
    X = sitk.GetArrayFromImage(I)
    # number of pixel
    N = np.prod(X.shape)

    eps = np.round(np.random.uniform(-10, 10, 100))
    idx = np.random.random_integers(0, N-1, np.round(N*p))
    X.ravel()[idx] = np.array(200+eps, dtype=np.uint8)

    # write outlier image
    J = sitk.GetImageFromArray(X)
    sitk.WriteImage(J, sys.argv[3], True)

    # decompose X into L+S
    L, S, _, _, _, _ = ialm.recover(X)

    C = sitk.GetImageFromArray(np.asarray(L, dtype=np.uint8))
    sitk.WriteImage(C, sys.argv[4], True)

    # compute mean-square error and Frobenius norm
    print "MSE: %.4g" % np.sqrt(np.asmatrix((L-sitk.GetArrayFromImage(I))**2).sum())
    print "Frobenius-Norm: %.4g" % np.linalg.norm(L-sitk.GetArrayFromImage(I),ord='fro') 

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)