Python msct_image.Image类代码示例

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