Python SimpleITK.GetImageFromArray() 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 CalDice(self, pred, gt, class_lbl=1):
        """
        Calculate dice score
        """
        # intersection = np.sum((pred == class_lbl) & (gt == class_lbl))
        # dice_4d = 2.0 *(intersection)/(np.sum(pred == class_lbl) + np.sum(gt == class_lbl))
        dices = []
        for i in range(pred.shape[3]):
            labelPred=sitk.GetImageFromArray(pred[:,:,:,i], isVector=False)
            labelTrue=sitk.GetImageFromArray(gt[:,:,:,i], isVector=False)
            dicecomputer=sitk.LabelOverlapMeasuresImageFilter()
            dicecomputer.Execute(labelTrue==class_lbl,labelPred==class_lbl)
            dice=dicecomputer.GetDiceCoefficient()
            dices.append(dice)
        # print (np.mean(dices), dice_4d)
        return np.mean(dices) 

def load_remove_save(input_file: str, output_file: str, for_which_classes: list,
                     minimum_valid_object_size: dict = None):
    # Only objects larger than minimum_valid_object_size will be removed. Keys in minimum_valid_object_size must
    # match entries in for_which_classes
    img_in = sitk.ReadImage(input_file)
    img_npy = sitk.GetArrayFromImage(img_in)
    volume_per_voxel = float(np.prod(img_in.GetSpacing(), dtype=np.float64))

    image, largest_removed, kept_size = remove_all_but_the_largest_connected_component(img_npy, for_which_classes,
                                                                                       volume_per_voxel,
                                                                                       minimum_valid_object_size)
    # print(input_file, "kept:", kept_size)
    img_out_itk = sitk.GetImageFromArray(image)
    img_out_itk = copy_geometry(img_out_itk, img_in)
    sitk.WriteImage(img_out_itk, output_file)
    return largest_removed, kept_size 

def convert_labels_back_to_BraTS_2018_2019_convention(input_folder: str, output_folder: str):
    """
    reads all prediction files (nifti) in the input folder, converts the labels back to BraTS convention and saves the
    result in output_folder
    :param input_folder:
    :param output_folder:
    :return:
    """
    maybe_mkdir_p(output_folder)
    nii = subfiles(input_folder, suffix='.nii.gz', join=False)
    for n in nii:
        a = sitk.ReadImage(join(input_folder, n))
        b = sitk.GetArrayFromImage(a)
        c = convert_labels_back_to_BraTS(b)
        d = sitk.GetImageFromArray(c)
        d.CopyInformation(a)
        sitk.WriteImage(d, join(output_folder, n)) 

def remove_all_but_the_two_largest_conn_comp(img_itk_file: str, file_out: str):
    """
    This was not used. I was just curious because others used this. Turns out this is not necessary for my networks
    """
    img_itk = sitk.ReadImage(img_itk_file)
    img_npy = sitk.GetArrayFromImage(img_itk)

    labelmap, num_labels = label((img_npy > 0).astype(int))

    if num_labels > 2:
        label_sizes = []
        for i in range(1, num_labels + 1):
            label_sizes.append(np.sum(labelmap == i))
        argsrt = np.argsort(label_sizes)[::-1] # two largest are now argsrt[0] and argsrt[1]
        keep_mask = (labelmap == argsrt[0] + 1) | (labelmap == argsrt[1] + 1)
        img_npy[~keep_mask] = 0
        new = sitk.GetImageFromArray(img_npy)
        new.CopyInformation(img_itk)
        sitk.WriteImage(new, file_out)
        print(os.path.basename(img_itk_file), num_labels, label_sizes)
    else:
        shutil.copy(img_itk_file, file_out) 

def save_image(img: np.ndarray, header: dict, output_file: str):
    """
    CAREFUL you need to restore_original_slice_orientation before saving!
    :param img:
    :param header:
    :return:
    """
    # reverse back
    img = reverse_axes(img)  # switch from zyx to xyz
    img_itk = sitk.GetImageFromArray(img)
    img_itk.SetSpacing(header['spacing'])
    img_itk.SetOrigin(header['origin'])
    if not isinstance(header['direction'], tuple):
        img_itk.SetDirection(header['direction'].flatten())
    else:
        img_itk.SetDirection(header['direction'])

    sitk.WriteImage(img_itk, output_file) 

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 export_segmentations_postprocess(indir, outdir):
    maybe_mkdir_p(outdir)
    niftis = subfiles(indir, suffix='nii.gz', join=False)
    for n in niftis:
        print("\n", n)
        identifier = str(n.split("_")[-1][:-7])
        outfname = join(outdir, "test-segmentation-%s.nii" % identifier)
        img = sitk.ReadImage(join(indir, n))
        img_npy = sitk.GetArrayFromImage(img)
        lmap, num_objects = label((img_npy > 0).astype(int))
        sizes = []
        for o in range(1, num_objects + 1):
            sizes.append((lmap == o).sum())
        mx = np.argmax(sizes) + 1
        print(sizes)
        img_npy[lmap != mx] = 0
        img_new = sitk.GetImageFromArray(img_npy)
        img_new.CopyInformation(img)
        sitk.WriteImage(img_new, outfname) 

def split_4d_nifti(filename, output_folder):
    img_itk = sitk.ReadImage(filename)
    dim = img_itk.GetDimension()
    file_base = filename.split("/")[-1]
    if dim == 3:
        shutil.copy(filename, join(output_folder, file_base[:-7] + "_0000.nii.gz"))
        return
    elif dim != 4:
        raise RuntimeError("Unexpected dimensionality: %d of file %s, cannot split" % (dim, filename))
    else:
        img_npy = sitk.GetArrayFromImage(img_itk)
        spacing = img_itk.GetSpacing()
        origin = img_itk.GetOrigin()
        direction = np.array(img_itk.GetDirection()).reshape(4,4)
        # now modify these to remove the fourth dimension
        spacing = tuple(list(spacing[:-1]))
        origin = tuple(list(origin[:-1]))
        direction = tuple(direction[:-1, :-1].reshape(-1))
        for i, t in enumerate(range(img_npy.shape[0])):
            img = img_npy[t]
            img_itk_new = sitk.GetImageFromArray(img)
            img_itk_new.SetSpacing(spacing)
            img_itk_new.SetOrigin(origin)
            img_itk_new.SetDirection(direction)
            sitk.WriteImage(img_itk_new, join(output_folder, file_base[:-7] + "_%04.0d.nii.gz" % i)) 

def loop(self):
    if self.trigger == "internal":
      if self.fps_label:
        while self.inputs[0].poll():
          self.camera.max_fps = self.inputs[0].recv()[self.fps_label]
      t,img = self.camera.read_image()
    elif self.trigger == "external":
      data = self.inputs[0].recv()  # wait for a signal
      if data is None:
        return
      t,img = self.camera.get_image()
    self.timer = time.time()
    if self.save_folder:
      if not sitk:
        raise IOError("[Camera] Cannot save image, sitk is not installed !")
      image = sitk.GetImageFromArray(img)
      sitk.WriteImage(image,
               self.save_folder + "img_%.6d_%.5f.tiff" % (
               self.loops, t-self.t0))
    self.loops += 1
    self.send([t-self.t0,img]) 

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 get_borders_mean(image, filter_otsu=True):
    # pylint: disable=bad-whitespace
    array = sitk.GetArrayViewFromImage(image)
    borders_tuple = (
        array[ 0,  :,  :],
        array[-1,  :,  :],
        array[ :,  0,  :],
        array[ :, -1,  :],
        array[ :,  :,  0],
        array[ :,  :, -1],
    )
    borders_flat = np.hstack([border.ravel() for border in borders_tuple])
    if not filter_otsu:
        return borders_flat.mean()
    borders_reshaped = borders_flat.reshape(1, 1, -1)
    borders_image = sitk.GetImageFromArray(borders_reshaped)
    otsu = sitk.OtsuThresholdImageFilter()
    otsu.Execute(borders_image)
    threshold = otsu.GetThreshold()
    values = borders_flat[borders_flat < threshold]
    if values.any():
        default_value = values.mean()
    else:
        default_value = borders_flat.mean()
    return default_value 

def __init__(self, input_image, 
                    input_mask=None, 
                    save_path=None, 
                    seq='Flair',
                    class_ = 'ET',
                    all_=True):
        
        self.input_image = input_image
        if not input_mask:
            self.input_mask = np.ones(tuple(list(self.input_image.shape)[:-1]))
        else: self.input_mask = input_mask
        
        self.img = sitk.GetImageFromArray(self.input_image)
        self.GT  = sitk.GetImageFromArray(self.input_mask)
        self.save_path = save_path
        self.seq = seq
        self.all_ = all_
        self.class_ = class_
        self.feat_dict = {} 

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 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 itk(self):

        # flip axis to
        df = Displacement(self.image.clone(), self.size, self.spacing, self.origin)
        df._reverse_axis()
        df.image = df.image.squeeze()
        df.image = df.image.cpu()

        if len(self.size) == 2:
            itk_displacement = sitk.GetImageFromArray(df.image.numpy(), isVector=True)
        elif len(self.size) == 3:
            itk_displacement = sitk.GetImageFromArray(df.image.numpy())

        itk_displacement.SetSpacing(spacing=self.spacing)
        itk_displacement.SetOrigin(origin=self.origin)
        return itk_displacement 

def ReadMhd_WriteNifty(impath,impath_new):

    '''
    input:param impath & impsth_new
    return: imsave which aim to save the format changed image to a new designated path;
    all functions used in this function are from SimpleITK
    sitk.ReadImage():aim to read the pending image in
    sitk.GetArrayFromImage():aim to get array from the image
    sitk.GetImageFromArray()&sitk.WriteImage():both achieve to write and save the new image as what we want

    '''

# main part of the function
    image=sitk.ReadImage(impath)
    image_arr=sitk.GetArrayFromImage(image)
    imnew=sitk.GetImageFromArray(image_arr)
    imsave=sitk.WriteImage(imnew,impath_new)
    return imsave

# an example of using the function
# read the new image in the format of Nifity 

def zero_one(image, mask=None, fill_value=0):
    """Normalizes an image by mapping the min to zero, and max to one."""

    img_array = sitk.GetArrayFromImage(image)
    img_array = img_array.astype(np.float32)

    msk_array = np.ones(img_array.shape)

    if mask is not None:
        msk_array = sitk.GetArrayFromImage(mask)

    min_value = np.min(img_array[msk_array>0])
    max_value = np.max(img_array[msk_array>0])

    img_array = (img_array - min_value) / (max_value - min_value)
    img_array[msk_array == 0] = fill_value

    image_normalised = sitk.GetImageFromArray(img_array)
    image_normalised.CopyInformation(image)

    return image_normalised 

def range_matching(image, mask=None, low_percentile=4, high_percentile=96, fill_value=0):
    """Normalizes an image by mapping the low_percentile to zero, and the high_percentile to one."""

    img_array = sitk.GetArrayFromImage(image)
    img_array = img_array.astype(np.float32)

    msk_array = np.ones(img_array.shape)

    if mask is not None:
        msk_array = sitk.GetArrayFromImage(mask)

    lo_p = np.percentile(img_array[msk_array>0], low_percentile)
    hi_p = np.percentile(img_array[msk_array>0], high_percentile)

    img_array = (img_array - lo_p) / (hi_p - lo_p)
    img_array[msk_array == 0] = fill_value

    image_normalised = sitk.GetImageFromArray(img_array)
    image_normalised.CopyInformation(image)

    return image_normalised 

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 zero_mean_unit_var(image, mask=None, fill_value=0):
    """Normalizes an image to zero mean and unit variance."""

    img_array = sitk.GetArrayFromImage(image)
    img_array = img_array.astype(np.float32)

    msk_array = np.ones(img_array.shape)

    if mask is not None:
        msk_array = sitk.GetArrayFromImage(mask)

    mean = np.mean(img_array[msk_array>0])
    std = np.std(img_array[msk_array>0])

    if std > 0:
        img_array = (img_array - mean) / std
        img_array[msk_array==0] = fill_value

    image_normalised = sitk.GetImageFromArray(img_array)
    image_normalised.CopyInformation(image)

    return image_normalised 

def average_surface_distance(predictions, labels, one_hot=False, unindexed_classes=0, spacing=[1, 1, 1]):
    def one_class_average_surface_distance(pred, lab):
        hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()
        batch = pred.shape[0]
        result = []
        for i in range(batch):
            pred_img = sitk.GetImageFromArray(pred[i].cpu().numpy())
            pred_img.SetSpacing(spacing)
            lab_img = sitk.GetImageFromArray(lab[i].cpu().numpy())
            lab_img.SetSpacing(spacing)
            hausdorff_distance_filter.Execute(pred_img, lab_img)
            result.append(hausdorff_distance_filter.GetAverageHausdorffDistance())
        return torch.tensor(np.asarray(result))

    return multi_class_score(one_class_average_surface_distance, predictions, labels, one_hot=one_hot,
                             unindexed_classes=unindexed_classes) 

def hausdorff_distance(predictions, labels, one_hot=False, unindexed_classes=0, spacing=[1, 1, 1]):
    def one_class_hausdorff_distance(pred, lab):
        hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()
        batch = pred.shape[0]
        result = []
        for i in range(batch):
            pred_img = sitk.GetImageFromArray(pred[i].cpu().numpy())
            pred_img.SetSpacing(spacing)
            lab_img = sitk.GetImageFromArray(lab[i].cpu().numpy())
            lab_img.SetSpacing(spacing)
            hausdorff_distance_filter.Execute(pred_img, lab_img)
            result.append(hausdorff_distance_filter.GetHausdorffDistance())
        return torch.tensor(np.asarray(result))

    return multi_class_score(one_class_hausdorff_distance, predictions, labels, one_hot=one_hot,
                             unindexed_classes=unindexed_classes) 

def run(self):
        mask_sitk = self._srr_mask_sitk

        # Smooth mask
        mask_sitk = sitk.SmoothingRecursiveGaussian(mask_sitk, self._sigma)

        # Binarize images given thresholds
        mask_sitk = sitk.BinaryThreshold(
            mask_sitk, lowerThreshold=self._lower, upperThreshold=self._upper)

        # Keep largest connected region only
        nda = sitk.GetArrayFromImage(mask_sitk)
        nda = tse.TemplateStackEstimator.get_largest_connected_region_mask(nda)

        self._mask_sitk = sitk.GetImageFromArray(nda)
        self._mask_sitk.CopyInformation(mask_sitk) 

def run(config_file_2d, config_file_3d, output_folder):
    cf_2d = imp.load_source("cf_2d", config_file_2d)
    cf_3d = imp.load_source("cf_3d", config_file_3d)

    dataset_base_dir = cf_2d.dataset_root_test
    results_folder_3D = os.path.join(cf_3d.results_dir, "test_predictions/")
    results_folder_2D = os.path.join(cf_2d.results_dir, "test_predictions/")
    patient_info = generate_patient_info(dataset_base_dir)

    if not os.path.isdir(output_folder):
        os.mkdir(output_folder)


    def resize_softmax_pred(softmax_output, new_shape, order=3):
        reshaped = np.zeros([len(softmax_output)] + list(new_shape), dtype=float)
        for i in range(len(softmax_output)):
            reshaped[i] = resize(softmax_output[i].astype(float), new_shape, order, mode="constant", cval=0, clip=True)
        return reshaped

    for patient in range(101, 151):
        for tpe in ['ed', 'es']:
            all_softmax = []
            raw_itk = sitk.ReadImage(os.path.join(dataset_base_dir, "patient%03.0d"%patient,
                                                  "patient%03.0d_frame%02.0d.nii.gz" %
                                                  (patient, patient_info[patient][tpe])))
            raw = sitk.GetArrayFromImage(raw_itk)
            for f in range(5):
                res_3d = np.load(os.path.join(results_folder_3D, "fold%d" % f, "patient%03.0d_3D_net.npz" % patient))
                res_2d = np.load(os.path.join(results_folder_2D, "fold%d" % f, "patient%03.0d_2D_net.npz" % patient))
                # resize softmax to original image size
                softmax_3d = resize_softmax_pred(res_3d[tpe], raw.shape, 3)
                softmax_2d = resize_softmax_pred(res_2d[tpe], raw.shape, 3)

                all_softmax += [softmax_3d[None], softmax_2d[None]]
            predicted_seg = postprocess_prediction(np.vstack(all_softmax).mean(0).argmax(0))

            itk_seg = sitk.GetImageFromArray(predicted_seg.astype(np.uint8))
            itk_seg.CopyInformation(raw_itk)
            sitk.WriteImage(itk_seg, os.path.join(output_folder, "patient%03.0d_%s.nii.gz" % (patient, tpe.upper()))) 

def _generate_identity_mask(self):
        shape = sitk.GetArrayFromImage(self.sitk).shape
        nda = np.ones(shape, dtype=np.uint8)

        binary_mask = sitk.GetImageFromArray(nda)
        binary_mask.CopyInformation(self.sitk)

        return binary_mask 

def run_dilation(self):

        time_start = ph.start_timing()

        dilater = sitk.BinaryDilateImageFilter()
        dilater.SetKernelType(eval("sitk.sitk" + self._dilation_kernel))
        dilater.SetKernelRadius(self._dilation_radius)

        if self._use_dilation_in_plane_only:

            shape = self._mask_sitk.GetSize()
            N_slices = shape[2]
            nda_mask = np.zeros(shape[::-1], dtype=np.uint8)

            for i in range(0, N_slices):
                slice_mask_sitk = self._mask_sitk[:, :, i:i + 1]
                mask_sitk = dilater.Execute(slice_mask_sitk)
                nda_mask[i, :, :] = sitk.GetArrayFromImage(mask_sitk)

            mask_sitk = sitk.GetImageFromArray(nda_mask)
            mask_sitk.CopyInformation(self._mask_sitk)
            self._mask_sitk = mask_sitk

        else:
            self._mask_sitk = dilater.Execute(self._mask_sitk)

        if self._stack is not None:
            self._stack = st.Stack.from_sitk_image(
                self._stack.sitk,
                image_sitk_mask=self._mask_sitk,
                filename=self._stack.get_filename(),
                slice_thickness=self._stack.get_slice_thickness(),
            )

        self._computational_time = ph.stop_timing(time_start) 

def _generate_identity_mask(self):
        shape = sitk.GetArrayFromImage(self.sitk).shape
        nda = np.ones(shape, dtype=np.uint8)

        binary_mask = sitk.GetImageFromArray(nda)
        binary_mask.CopyInformation(self.sitk)

        return binary_mask 

def sitk_new_blank_image(size, spacing, direction, origin, default_value=0.):
    image = sitk.GetImageFromArray(np.ones(size, dtype=np.float).T * default_value)
    image.SetSpacing(spacing)
    image.SetDirection(direction)
    image.SetOrigin(origin)
    return image 

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')