def load_level(list_slices_target, fname_level):
    verbose = 1
    path_level, file_level, ext_level = extract_fname(fname_level)

    #  ####### Check if the level file is an image or a text file
    # Level file is an image
    if ext_level in ['.nii', '.nii.gz']:
        im_level = Image(fname_level)
        im_level.change_orientation('IRP')

        list_level = []
        list_med_level = []
        for slice_level in im_level.data:
            try:
                # vertebral level of the slice
                l = np.mean(slice_level[slice_level > 0])
                # median of the vertebral level of the slice: if all voxels are int, med will be an int.
                med = np.median(slice_level[slice_level > 0])
                # change med in int if it is an int
                med = int(med) if int(med)==med else med
            except Exception, e:
                printv('WARNING: ' + str(e) + '\nNo level label found. Level will be set to 0 for this slice', verbose, 'warning')
                l = 0
                med = 0
            list_level.append(l)
            list_med_level.append(med)

        # if all median of level are int for all slices : consider level as int
        if all([isinstance(med, int) for med in list_med_level]):
            # level as int are placed in the middle of each vertebra (that's why there is a "+0.5")
            list_level = [int(round(l))+0.5 for l in list_level]

def set_orientation(im, orientation, data_inversion=False, filename=False, fname_out=''):
    """
    Set orientation on image
    :param im: either Image object or file name. Carefully set param filename.
    :param orientation:
    :param data_inversion:
    :param filename:
    :return:
    """

    if fname_out:
        pass
    elif filename:
        path, fname, ext = extract_fname(im)
        fname_out = fname+'_'+orientation+ext
    else:
        fname_out = im.file_name+'_'+orientation+im.ext

    if not data_inversion:
        from sct_utils import run
        if filename:
            run('isct_orientation3d -i '+im+' -orientation '+orientation+' -o '+fname_out, 0)
            im_out = fname_out
        else:
            run('isct_orientation3d -i '+im.absolutepath+' -orientation '+orientation+' -o '+fname_out, 0)
            im_out = Image(fname_out)
    else:
        im_out = im.copy()
        im_out.change_orientation(orientation, True)
        im_out.setFileName(fname_out)
    return im_out

def multicomponent_merge(fname_list):
    from numpy import zeros, reshape
    # WARNING: output multicomponent is not optimal yet, some issues may be related to the use of this function

    im_0 = Image(fname_list[0])
    new_shape = list(im_0.data.shape)
    if len(new_shape) == 3:
        new_shape.append(1)
    new_shape.append(len(fname_list))
    new_shape = tuple(new_shape)

    data_out = zeros(new_shape)
    for i, fname in enumerate(fname_list):
        im = Image(fname)
        dat = im.data
        if len(dat.shape) == 2:
            data_out[:, :, 0, 0, i] = dat.astype('float32')
        elif len(dat.shape) == 3:
            data_out[:, :, :, 0, i] = dat.astype('float32')
        elif len(dat.shape) == 4:
            data_out[:, :, :, :, i] = dat.astype('float32')
        del im
        del dat
    im_out = im_0.copy()
    im_out.data = data_out.astype('float32')
    im_out.hdr.set_intent('vector', (), '')
    im_out.setFileName(im_out.file_name+'_multicomponent'+im_out.ext)
    return im_out

    def remove_label(self, symmetry=False):
        """
        Compare two label images and remove any labels in input image that are not in reference image.
        The symmetry option enables to remove labels from reference image that are not in input image
        """
        # image_output = Image(self.image_input.dim, orientation=self.image_input.orientation, hdr=self.image_input.hdr, verbose=self.verbose)
        image_output = Image(self.image_input, verbose=self.verbose)
        image_output.data *= 0  # put all voxels to 0

        result_coord_input, result_coord_ref = self.remove_label_coord(self.image_input.getNonZeroCoordinates(coordValue=True),
                                                                       self.image_ref.getNonZeroCoordinates(coordValue=True), symmetry)

        for coord in result_coord_input:
            image_output.data[int(coord.x), int(coord.y), int(coord.z)] = int(round(coord.value))

        if symmetry:
            # image_output_ref = Image(self.image_ref.dim, orientation=self.image_ref.orientation, hdr=self.image_ref.hdr, verbose=self.verbose)
            image_output_ref = Image(self.image_ref, verbose=self.verbose)
            for coord in result_coord_ref:
                image_output_ref.data[int(coord.x), int(coord.y), int(coord.z)] = int(round(coord.value))
            image_output_ref.setFileName(self.fname_output[1])
            image_output_ref.save('minimize_int')

            self.fname_output = self.fname_output[0]

        return image_output

def savePredictions(predictions, path_output, list_images, segmentation_image_size):
    number_of_images = len(list_images)
    predictions = numpy.reshape(predictions, [number_of_images, segmentation_image_size, segmentation_image_size, NUM_LABELS])
    predictions = predictions[:, :, :, 1]
    for i, pref in enumerate(predictions):
        im_pred = Image(pref)
        im_pred.setFileName(path_output+sct.add_suffix(list_images[i], '_pred'))
        im_pred.save()

    def continuous_vertebral_levels(self):
        """
        This function transforms the vertebral levels file from the template into a continuous file.
        Instead of having integer representing the vertebral level on each slice, a continuous value that represents
        the position of the slice in the vertebral level coordinate system.
        The image must be RPI
        :return:
        """
        im_input = Image(self.image_input, self.verbose)
        im_output = Image(self.image_input, self.verbose)
        im_output.data *= 0

        # 1. extract vertebral levels from input image
        #   a. extract centerline
        #   b. for each slice, extract corresponding level
        nx, ny, nz, nt, px, py, pz, pt = im_input.dim
        from sct_straighten_spinalcord import smooth_centerline
        x_centerline_fit, y_centerline_fit, z_centerline_fit, x_centerline_deriv, y_centerline_deriv, z_centerline_deriv = smooth_centerline(self.image_input, algo_fitting='nurbs', verbose=0)
        value_centerline = np.array([im_input.data[int(x_centerline_fit[it]), int(y_centerline_fit[it]), int(z_centerline_fit[it])] for it in range(len(z_centerline_fit))])

        # 2. compute distance for each vertebral level --> Di for i being the vertebral levels
        vertebral_levels = {}
        for slice_image, level in enumerate(value_centerline):
            if level not in vertebral_levels:
                vertebral_levels[level] = slice_image

        length_levels = {}
        for level in vertebral_levels:
            indexes_slice = np.where(value_centerline == level)
            length_levels[level] = np.sum([math.sqrt(((x_centerline_fit[indexes_slice[0][index_slice + 1]] - x_centerline_fit[indexes_slice[0][index_slice]])*px)**2 +
                                                     ((y_centerline_fit[indexes_slice[0][index_slice + 1]] - y_centerline_fit[indexes_slice[0][index_slice]])*py)**2 +
                                                     ((z_centerline_fit[indexes_slice[0][index_slice + 1]] - z_centerline_fit[indexes_slice[0][index_slice]])*pz)**2)
                                           for index_slice in range(len(indexes_slice[0]) - 1)])

        # 2. for each slice:
        #   a. identify corresponding vertebral level --> i
        #   b. calculate distance of slice from upper vertebral level --> d
        #   c. compute relative distance in the vertebral level coordinate system --> d/Di
        continuous_values = {}
        for it, iz in enumerate(z_centerline_fit):
            level = value_centerline[it]
            indexes_slice = np.where(value_centerline == level)
            indexes_slice = indexes_slice[0][indexes_slice[0] >= it]
            distance_from_level = np.sum([math.sqrt(((x_centerline_fit[indexes_slice[index_slice + 1]] - x_centerline_fit[indexes_slice[index_slice]]) * px * px) ** 2 +
                                                    ((y_centerline_fit[indexes_slice[index_slice + 1]] - y_centerline_fit[indexes_slice[index_slice]]) * py * py) ** 2 +
                                                    ((z_centerline_fit[indexes_slice[index_slice + 1]] - z_centerline_fit[indexes_slice[index_slice]]) * pz * pz) ** 2)
                                          for index_slice in range(len(indexes_slice) - 1)])
            continuous_values[iz] = level + 2.0 * distance_from_level / float(length_levels[level])

        # 3. saving data
        # for each slice, get all non-zero pixels and replace with continuous values
        coordinates_input = self.image_input.getNonZeroCoordinates()
        im_output.changeType('float32')
        # for all points in input, find the value that has to be set up, depending on the vertebral level
        for i, coord in enumerate(coordinates_input):
            im_output.data[int(coord.x), int(coord.y), int(coord.z)] = continuous_values[coord.z]

        return im_output

 def next(self):
     if self.iteration <= self.num_of_frames:
         result = Image(self)
         print "Iteration #" + str(self.iteration)
         result.data *= float(self.iteration) / float(self.num_of_frames)
         result.file_name = "tmp."+result.file_name+"_" + str(self.iteration)
         self.iteration += 1
         return result, self.iteration
     else:
         raise StopIteration()

    def plan_ref(self):
        """
        Generate a plane in the reference space for each label present in the input image
        """

        image_output = Image(self.image_ref, self.verbose)
        image_output.data *= 0

        image_input_neg = Image(self.image_input, self.verbose).copy()
        image_input_pos = Image(self.image_input, self.verbose).copy()
        image_input_neg.data *=0
        image_input_pos.data *=0
        X, Y, Z = (self.image_input.data< 0).nonzero()
        for i in range(len(X)):
            image_input_neg.data[X[i], Y[i], Z[i]] = -self.image_input.data[X[i], Y[i], Z[i]] # in order to apply getNonZeroCoordinates
        X_pos, Y_pos, Z_pos = (self.image_input.data> 0).nonzero()
        for i in range(len(X_pos)):
            image_input_pos.data[X_pos[i], Y_pos[i], Z_pos[i]] = self.image_input.data[X_pos[i], Y_pos[i], Z_pos[i]]

        coordinates_input_neg = image_input_neg.getNonZeroCoordinates()
        coordinates_input_pos = image_input_pos.getNonZeroCoordinates()

        image_output.changeType('float32')
        for coord in coordinates_input_neg:
            image_output.data[:, :, int(coord.z)] = -coord.value #PB: takes the int value of coord.value
        for coord in coordinates_input_pos:
            image_output.data[:, :, int(coord.z)] = coord.value

        return image_output

def concat_data(fname_in_list, dim, pixdim=None):
    """
    Concatenate data
    :param im_in_list: list of images.
    :param dim: dimension: 0, 1, 2, 3.
    :param pixdim: pixel resolution to join to image header
    :return im_out: concatenated image
    """
    # WARNING: calling concat_data in python instead of in command line causes a non understood issue (results are different with both options)
    from numpy import concatenate, expand_dims, squeeze

    dat_list = []
    data_concat_list = []

    # check if shape of first image is smaller than asked dim to concatenate along
    data0 = Image(fname_in_list[0]).data
    if len(data0.shape) <= dim:
        expand_dim = True
    else:
        expand_dim = False

    for i, fname in enumerate(fname_in_list):
        # if there is more than 100 images to concatenate, then it does it iteratively to avoid memory issue.
        if i != 0 and i % 100 == 0:
            data_concat_list.append(concatenate(dat_list, axis=dim))
            im = Image(fname)
            dat = im.data
            if expand_dim:
                dat = expand_dims(dat, dim)
            dat_list = [dat]
            del im
            del dat
        else:
            im = Image(fname)
            dat = im.data
            if expand_dim:
                dat = expand_dims(dat, dim)
            dat_list.append(dat)
            del im
            del dat
    if data_concat_list:
        data_concat_list.append(concatenate(dat_list, axis=dim))
        data_concat = concatenate(data_concat_list, axis=dim)
    else:
        data_concat = concatenate(dat_list, axis=dim)
    # write file
    im_out = Image(fname_in_list[0]).copy()
    im_out.data = data_concat
    im_out.setFileName(im_out.file_name+'_concat'+im_out.ext)

    if pixdim is not None:
        im_out.hdr['pixdim'] = pixdim

    return im_out

    def get_im_from_list(self, data):
        im = Image(data)
        # set pix dimension
        im.hdr.structarr['pixdim'][1] = self.param_data.axial_res
        im.hdr.structarr['pixdim'][2] = self.param_data.axial_res
        # set the correct orientation
        im.setFileName('im_to_orient.nii.gz')
        im.save()
        im = set_orientation(im, 'IRP')
        im = set_orientation(im, 'PIL', data_inversion=True)

        return im

    def __init__(self, fname_label, fname_output=None, fname_ref=None, cross_radius=5, dilate=False,
                 coordinates=None, verbose='1'):
        self.image_input = Image(fname_label)

        if fname_ref is not None:
            self.image_ref = Image(fname_ref)

        self.fname_output = fname_output
        self.cross_radius = cross_radius
        self.dilate = dilate
        self.coordinates = coordinates
        self.verbose = verbose

def register_seg(seg_input, seg_dest):
    """Slice-by-slice registration by translation of two segmentations.

    For each slice, we estimate the translation vector by calculating the difference of position of the two centers of
    mass.
    The segmentations can be of different sizes but the output segmentation must be smaller than the input segmentation.

    input:
        seg_input: name of moving segmentation file (type: string)
        seg_dest: name of fixed segmentation file (type: string)

    output:
        x_displacement: list of translation along x axis for each slice (type: list)
        y_displacement: list of translation along y axis for each slice (type: list)

    """
    seg_input_img = Image(seg_input)
    seg_dest_img = Image(seg_dest)
    seg_input_data = seg_input_img.data
    seg_dest_data = seg_dest_img.data

    x_center_of_mass_input = [0 for i in range(seg_dest_data.shape[2])]
    y_center_of_mass_input = [0 for i in range(seg_dest_data.shape[2])]
    print "\nGet center of mass of the input segmentation for each slice (corresponding to a slice in the output segmentation)..."  # different if size of the two seg are different
    # TO DO: select only the slices corresponding to the output segmentation
    coord_origin_dest = seg_dest_img.transfo_pix2phys([[0, 0, 0]])
    [[x_o, y_o, z_o]] = seg_input_img.transfo_phys2pix(coord_origin_dest)
    for iz in xrange(seg_dest_data.shape[2]):
        x_center_of_mass_input[iz], y_center_of_mass_input[iz] = ndimage.measurements.center_of_mass(
            array(seg_input_data[:, :, z_o + iz])
        )

    x_center_of_mass_output = [0 for i in range(seg_dest_data.shape[2])]
    y_center_of_mass_output = [0 for i in range(seg_dest_data.shape[2])]
    print "\nGet center of mass of the output segmentation for each slice ..."
    for iz in xrange(seg_dest_data.shape[2]):
        x_center_of_mass_output[iz], y_center_of_mass_output[iz] = ndimage.measurements.center_of_mass(
            array(seg_dest_data[:, :, iz])
        )

    x_displacement = [0 for i in range(seg_input_data.shape[2])]
    y_displacement = [0 for i in range(seg_input_data.shape[2])]
    print "\nGet displacement by voxel..."
    for iz in xrange(seg_dest_data.shape[2]):
        x_displacement[iz] = -(
            x_center_of_mass_output[iz] - x_center_of_mass_input[iz]
        )  # WARNING: in ITK's coordinate system, this is actually Tx and not -Tx
        y_displacement[iz] = (
            y_center_of_mass_output[iz] - y_center_of_mass_input[iz]
        )  # This is Ty in ITK's and fslview' coordinate systems

    return x_displacement, y_displacement

    def execute(self):
        """
        This method executes the symmetry detection
        :return: returns the symmetry data
        """
        img = Image(self.input_image)
        raw_orientation = img.change_orientation()
        data = np.squeeze(img.data)
        dim = data.shape
        section_length = dim[1]/self.nb_sections

        result = np.zeros(dim)

        for i in range(0, self.nb_sections):
            if (i+1)*section_length > dim[1]:
                y_length = (i+1)*section_length - ((i+1)*section_length - dim[1])
                result[:, i*section_length:i*section_length + y_length, :] = symmetry_detector_right_left(data[:, i*section_length:i*section_length + y_length, :],  cropped_xy=self.crop_xy)
            sym = symmetry_detector_right_left(data[:, i*section_length:(i+1)*section_length, :], cropped_xy=self.crop_xy)
            result[:, i*section_length:(i+1)*section_length, :] = sym

        result_image = Image(img)
        if len(result_image.data) == 4:
            result_image.data = result[:,:,:,np.newaxis]
        else:
            result_image.data = result

        result_image.change_orientation(raw_orientation)

        return result_image.data

    def __init__(self, fname_label, fname_output=None, fname_ref=None, cross_radius=5, dilate=False,
                 coordinates=None, verbose=1):
        self.image_input = Image(fname_label, verbose=verbose)

        if fname_ref is not None:
            self.image_ref = Image(fname_ref, verbose=verbose)

        if isinstance(fname_output, list):
            if len(fname_output) == 1:
                self.fname_output = fname_output[0]
            else:
                self.fname_output = fname_output
        else:
            self.fname_output = fname_output
        self.cross_radius = cross_radius
        self.dilate = dilate
        self.coordinates = coordinates
        self.verbose = verbose

def get_minimum_path_nii(fname):
    from msct_image import Image
    data=Image(fname)
    vesselness_data = data.data
    raw_orient=data.change_orientation()
    data.data=get_minimum_path(data.data, invert=1)
    data.change_orientation(raw_orient)
    data.file_name += '_minimalpath'
    data.save()

    def crop_from_mask_with_background(self):
        from numpy import asarray, einsum
        image_in = Image(self.input_filename)
        data_array = asarray(image_in.data)
        data_mask = asarray(Image(self.mask).data)
        assert data_array.shape == data_mask.shape

        # Element-wise matrix multiplication:
        new_data = None
        dim = len(data_array.shape)
        if dim == 3:
            new_data = einsum('ijk,ijk->ijk', data_mask, data_array)
        elif dim == 2:
            new_data = einsum('ij,ij->ij', data_mask, data_array)

        if self.background != 0:
            from sct_maths import get_data_or_scalar
            data_background = get_data_or_scalar(str(self.background), data_array)
            data_mask_inv = data_mask.max() - data_mask
            if dim == 3:
                data_background = einsum('ijk,ijk->ijk', data_mask_inv, data_background)
            elif dim == 2:
                data_background = einsum('ij,ij->ij', data_mask_inv, data_background)
            new_data += data_background

        # set image out
        image_in.setFileName(self.output_filename)
        image_in.data = new_data
        image_in.save()

def concat_data(fname_in, fname_out, dim):
    """
    Concatenate data
    :param fname_in: list of file names.
    :param fname_out:
    :param dim: dimension: 0, 1, 2, 3.
    :return: none
    """
    # create empty list
    list_data = []

    # loop across files
    for i in range(len(fname_in)):
        # append data to list
        list_data.append(Image(fname_in[i]).data)

    # expand dimension of all elements in the list if necessary
    if dim > list_data[0].ndim-1:
        list_data = [expand_dims(i, dim) for i in list_data]
    # concatenate
    try:
        data_concat = concatenate(list_data, axis=dim)
    except Exception as e:
        sct.printv('\nERROR: Concatenation on line {}'.format(sys.exc_info()[-1].tb_lineno)+'\n'+str(e)+'\n', 1, 'error')

    # write file
    im = Image(fname_in[0])
    im.data = data_concat
    im.setFileName(fname_out)
    im.save()

def main(args = None):

    dim_list = ['x', 'y', 'z', 't']

    if not args:
        args = sys.argv[1:]

    # Get parser info
    parser = get_parser()
    arguments = parser.parse(sys.argv[1:])
    fname_in = arguments["-i"]
    fname_out = arguments["-o"]
    verbose = int(arguments['-v'])

    # Build fname_out
    if fname_out == '':
        path_in, file_in, ext_in = extract_fname(fname_in)
        fname_out = path_in+file_in+'_mean'+ext_in

    # Open file.
    nii = Image(fname_in)
    data = nii.data

    # run command
    if '-otsu' in arguments:
        param = arguments['-otsu']
        data_out = otsu(data, param)
    elif '-otsu_adap' in arguments:
        param = arguments['-otsu_adap']
        data_out = otsu_adap(data, param[0], param[1])
    elif '-otsu_median' in arguments:
        param = arguments['-otsu_median']
        data_out = otsu_median(data, param[0], param[1])
    elif '-thr' in arguments:
        param = arguments['-thr']
        data_out = threshold(data, param)
    elif '-percent' in arguments:
        param = arguments['-percent']
        data_out = perc(data, param)
    elif '-mean' in arguments:
        dim = dim_list.index(arguments['-mean'])
        data_out = compute_mean(data, dim)
    elif '-std' in arguments:
        dim = dim_list.index(arguments['-std'])
        data_out = compute_std(data, dim)
    elif '-dilate' in arguments:
        data_out = dilate(data, arguments['-dilate'])
    elif '-erode' in arguments:
        data_out = erode(data, arguments['-dilate'])
    else:
        printv('No process applied.', 1, 'warning')
        return

    # Write output
    nii.data = data_out
    nii.setFileName(fname_out)
    nii.save()

    # display message
    printv('Created file:\n--> '+fname_out+'\n', verbose, 'info')

    def __init__(self, target_fname, sc_seg_fname, t2_data=None, denoising=True):

        self.t2star = 't2star.nii.gz'
        self.sc_seg = 't2star_sc_seg.nii.gz'
        self.t2 = 't2.nii.gz'
        self.t2_seg = 't2_seg.nii.gz'
        self.t2_landmarks = 't2_landmarks.nii.gz'
        self.resample_to = 0.3

        sct.run('cp ../' + target_fname + ' ./' + self.t2star)
        sct.run('cp ../' + sc_seg_fname + ' ./' + self.sc_seg)

        nx, ny, nz, nt, self.original_px, self.original_py, pz, pt = sct.get_dimension(self.t2star)

        if round(self.original_px, 2) != self.resample_to or round(self.original_py, 2) != self.resample_to:
            self.t2star = resample_image(self.t2star, npx=self.resample_to, npy=self.resample_to)
            self.sc_seg = resample_image(self.sc_seg, binary=True, npx=self.resample_to, npy=self.resample_to)

        t2star_im = Image(self.t2star)
        if denoising:
            t2star_im.denoise_ornlm()
            t2star_im.save()
            self.t2star = t2star_im.file_name + t2star_im.ext
        '''
        status, t2_star_orientation = sct.run('sct_orientation -i ' + self.t2star)
        self.original_orientation = t2_star_orientation[4:7]
        '''
        self.original_orientation = t2star_im.orientation

        self.square_mask = crop_t2_star(self.t2star, self.sc_seg, box_size=75)

        self.treated_target = sct.extract_fname(self.t2star)[1] + '_seg_in_croped.nii.gz'

        self.level_fname = None
        if t2_data is not None:
            sct.run('cp ../' + t2_data[0] + ' ./' + self.t2)
            sct.run('cp ../' + t2_data[1] + ' ./' + self.t2_seg)
            sct.run('cp ../' + t2_data[2] + ' ./' + self.t2_landmarks)

            self.level_fname = compute_level_file(self.t2star, self.sc_seg, self.t2, self.t2_seg, self.t2_landmarks)

def label_segmentation(fname_seg, list_disc_z, list_disc_value, verbose=1):
    """
    Label segmentation image
    :param fname_seg: fname of the segmentation
    :param list_disc_z: list of z that correspond to a disc
    :param list_disc_value: list of associated disc values
    :param verbose:
    :return:
    """
    # open segmentation
    seg = Image(fname_seg)
    dim = seg.dim
    ny = dim[1]
    nz = dim[2]
    # open labeled discs
    im_discs = Image(fname_seg)
    # loop across z
    for iz in range(nz):
        # get index of the disc right above iz
        try:
            ind_above_iz = max([i for i in range(len(list_disc_z)) if list_disc_z[i] > iz])
        except ValueError:
            # if ind_above_iz is empty, attribute value 0
            vertebral_level = 0
        else:
            # assign vertebral level (add one because iz is BELOW the disk)
            vertebral_level = list_disc_value[ind_above_iz] + 1
            # print vertebral_level
        # get voxels in mask
        ind_nonzero = np.nonzero(seg.data[:, :, iz])
        seg.data[ind_nonzero[0], ind_nonzero[1], iz] = vertebral_level
        if verbose == 2:
            import matplotlib
            matplotlib.use('Agg')
            import matplotlib.pyplot as plt
            plt.figure(50)
            plt.scatter(int(round(ny/2)), iz, c=vertebral_level, vmin=min(list_disc_value), vmax=max(list_disc_value), cmap='prism', marker='_', s=200)
    # write file
    seg.file_name += '_labeled'
    seg.save()