Python Image.copy方法代码示例

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
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

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]

#.........这里部分代码省略.........
        image_output = Image(self.image_input)
        image_output.data *= 0
        coordinates_input = self.image_input.getNonZeroCoordinates()

        # for all points with non-zeros neighbors, force the neighbors to 0
        for coord in coordinates_input:
            image_output.data[:,:,coord.z-width:coord.z+width] = offset + gap * coord.value

        return image_output

    def plan_ref(self):
        """
        This function generate a plan in the reference space for each label present in the input image
        """
        image_output = Image(self.image_ref)
        image_output.data *= 0
        coordinates_input = self.image_input.getNonZeroCoordinates()

        # for all points with non-zeros neighbors, force the neighbors to 0
        for coord in coordinates_input:
            image_output.data[:, :, coord.z] = coord.value

        return image_output
    def cubic_to_point(self):
        """
        This function calculates the center of mass of each group of labels and returns a file of same size with only a label by group at the center of mass.
        It is to be used after applying homothetic warping field to a label file as the labels will be dilated.
        :return:
        """
        from scipy import ndimage
        from numpy import array
        data = self.image_input.data

        image_output = self.image_input.copy()
        data_output = image_output.data
        data_output *= 0
        nx = image_output.data.shape[0]
        ny = image_output.data.shape[1]
        nz = image_output.data.shape[2]
        print '.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz)

        z_centerline = [iz for iz in range(0, nz, 1) if data[:,:,iz].any() ]
        nz_nonz = len(z_centerline)
        if nz_nonz==0 :
            print '\nERROR: Label file is empty'
            sys.exit()
        x_centerline = [0 for iz in range(0, nz_nonz, 1)]
        y_centerline = [0 for iz in range(0, nz_nonz, 1)]
        print '\nGet center of mass for each slice of the label file ...'
        for iz in xrange(len(z_centerline)):
            x_centerline[iz], y_centerline[iz] = ndimage.measurements.center_of_mass(array(data[:,:,z_centerline[iz]]))

        ## Calculate mean coordinate according to z for each cube of labels:
        list_cube_labels_x = [[]]
        list_cube_labels_y = [[]]
        list_cube_labels_z = [[]]
        count = 0
        for i in range(nz_nonz-1):
            # Make a list of group of slices that contains a non zero value
            if z_centerline[i] - z_centerline[i+1] == -1:
                # Verify if the value has already been recovered and add if not
                #If the group is empty add first value do not if it is not empty as it will copy it for a second time
                if len(list_cube_labels_z[count]) == 0 :#or list_cube_labels[count][-1] != z_centerline[i]:
                    list_cube_labels_z[count].append(z_centerline[i])
                    list_cube_labels_x[count].append(x_centerline[i])
                    list_cube_labels_y[count].append(y_centerline[i])

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
def vertebral_detection(fname, fname_seg, contrast):

    shift_AP = 14  # shift the centerline on the spine in mm default : 17 mm
    size_AP = 3  # mean around the centerline in the anterior-posterior direction in mm
    size_RL = 3  # mean around the centerline in the right-left direction in mm
    verbose = param.verbose

    if verbose:
        import matplotlib.pyplot as plt

    # open anatomical volume
    img = Image(fname)
    # orient to RPI
    img.change_orientation()
    # get dimension
    nx, ny, nz, nt, px, py, pz, pt = img.dim


    #==================================================
    # Compute intensity profile across vertebrae
    #==================================================

    shift_AP = shift_AP * py
    size_AP = size_AP * py
    size_RL = size_RL * px

    # orient segmentation to RPI
    run('sct_orientation -i ' + fname_seg + ' -s RPI')
    # smooth segmentation/centerline
    path_centerline, file_centerline, ext_centerline = extract_fname(fname_seg)
    x, y, z, Tx, Ty, Tz = smooth_centerline(path_centerline + file_centerline + '_RPI' + ext_centerline)

    # build intensity profile along the centerline
    I = np.zeros((len(y), 1))

    #  mask where intensity profile will be taken
    if verbose == 2:
        mat = img.copy()
        mat.data = np.zeros(mat.dim)

    for iz in range(len(z)):
        # define vector orthogonal to the cord in RL direction
        P1 = np.array([1, 0, -Tx[iz]/Tz[iz]])
        P1 = P1/np.linalg.norm(P1)
        # define vector orthogonal to the cord in AP direction
        P2 = np.array([0, 1, -Ty[iz]/Tz[iz]])
        P2 = P2/np.linalg.norm(P2)
        # define X and Y coordinates of the voxels to extract intensity profile from
        indexRL = range(-np.int(round(size_RL)), np.int(round(size_RL)))
        indexAP = range(0, np.int(round(size_AP)))+np.array(shift_AP)
        # loop over coordinates of perpendicular plane
        for i_RL in indexRL:
            for i_AP in indexAP:
                i_vect = np.round(np.array([x[iz], y[iz], z[iz]])+P1*i_RL+P2*i_AP)
                i_vect = np.minimum(np.maximum(i_vect, 0), np.array([nx, ny, nz])-1)  # check if index stays in image dimension
                I[iz] = I[iz] + img.data[i_vect[0], i_vect[1], i_vect[2]]

                # create a mask with this perpendicular plane
                if verbose == 2:
                    mat.data[i_vect[0], i_vect[1], i_vect[2]] = 1

    if verbose == 2:
        mat.file_name = 'mask'
        mat.save()

    # Detrending Intensity
    start_centerline_y = y[0]
    X = np.where(I == 0)
    mask2 = np.ones((len(y), 1), dtype=bool)
    mask2[X, 0] = False

    # low pass filtering
    import scipy.signal
    frequency = 2/pz
    Wn = 0.1/frequency
    N = 2              #Order of the filter
    #    b, a = scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba')
    b, a = scipy.signal.iirfilter(N, Wn, rp=None, rs=None, btype='high', analog=False, ftype='bessel', output='ba')
    I_detrend = scipy.signal.filtfilt(b, a, I[:, 0], axis=-1, padtype='constant', padlen=None)
    I_detrend = I_detrend/(np.amax(I_detrend))


    #==================================================
    # step 1 : Find the First Peak
    #==================================================
    if contrast == 't1':
        I_detrend2 = np.diff(I_detrend)
    elif contrast == 't2':
        space = np.linspace(-10/pz, 10/pz, round(21/pz), endpoint=True)
        pattern = (np.sinc((space*pz)/20)) ** 20
        I_corr = scipy.signal.correlate(-I_detrend.squeeze().squeeze()+1,pattern,'same')
        b, a = scipy.signal.iirfilter(N, Wn, rp=None, rs=None, btype='high', analog=False, ftype='bessel', output='ba')
        I_detrend2 = scipy.signal.filtfilt(b, a, I_corr, axis=-1, padtype='constant', padlen=None)

    I_detrend2[I_detrend2 < 0.2] = 0
    ind_locs = np.squeeze(scipy.signal.argrelextrema(I_detrend2, np.greater))

    # remove peaks that are too closed
    locsdiff = np.diff(z[ind_locs])
    ind = locsdiff > 10
#.........这里部分代码省略.........

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]

#.........这里部分代码省略.........
    yct = int(round(nyt/2))  # direction AP

    # define mean distance (in voxel) between adjacent discs: [C1/C2 -> C2/C3], [C2/C3 -> C4/C5], ..., [L1/L2 -> L2/L3]
    centerline_level = data_disc_template[xct, yct, :]
    # attribute value to each disc. Starts from max level, then decrease.
    min_level = centerline_level[centerline_level.nonzero()].min()
    max_level = centerline_level[centerline_level.nonzero()].max()
    list_disc_value_template = range(min_level, max_level)
    # add disc above top one
    list_disc_value_template.insert(int(0), min_level - 1)
    printv('\nDisc values from template: ' + str(list_disc_value_template), verbose)
    # get diff to find transitions (i.e., discs)
    diff_centerline_level = np.diff(centerline_level)
    # get disc z-values
    list_disc_z_template = diff_centerline_level.nonzero()[0].tolist()
    list_disc_z_template.reverse()
    printv('Z-values for each disc: ' + str(list_disc_z_template), verbose)
    list_distance_template = (
        np.diff(list_disc_z_template) * (-1)).tolist()  # multiplies by -1 to get positive distances
    printv('Distances between discs (in voxel): ' + str(list_distance_template), verbose)

    # if automatic mode, find C2/C3 disc
    if init_disc == [] and initc2 == 'auto':
        printv('\nDetect C2/C3 disk...', verbose)
        zrange = range(0, nz)
        ind_c2 = list_disc_value_template.index(2)
        z_peak = compute_corr_3d(src=data, target=data_template, x=xc, xshift=0, xsize=param.size_RL_initc2, y=yc, yshift=param.shift_AP_initc2, ysize=param.size_AP_initc2, z=0, zshift=param.shift_IS_initc2, zsize=param.size_IS_initc2, xtarget=xct, ytarget=yct, ztarget=list_disc_z_template[ind_c2], zrange=zrange, verbose=verbose, save_suffix='_initC2', gaussian_weighting=True, path_output=path_output)
        init_disc = [z_peak, 2]

    # if manual mode, open viewer for user to click on C2/C3 disc
    if init_disc == [] and initc2 == 'manual':
        from sct_viewer import ClickViewer
        # reorient image to SAL to be compatible with viewer
        im_input_SAL = im_input.copy()
        im_input_SAL.change_orientation('SAL')
        viewer = ClickViewer(im_input_SAL, orientation_subplot=['sag', 'ax'])
        viewer.number_of_slices = 1
        pz = 1
        viewer.gap_inter_slice = int(10 / pz)
        viewer.calculate_list_slices()
        viewer.help_url = 'https://sourceforge.net/p/spinalcordtoolbox/wiki/sct_label_vertebrae/attachment/label_vertebrae_viewer.png'
        # start the viewer that ask the user to enter a few points along the spinal cord
        mask_points = viewer.start()
        if mask_points:
            # create the mask containing either the three-points or centerline mask for initialization
            mask_filename = sct.add_suffix(fname, "_mask_viewer")
            sct.run("sct_label_utils -i " + fname + " -create " + mask_points + " -o " + mask_filename, verbose=False)
        else:
            sct.printv('\nERROR: the viewer has been closed before entering all manual points. Please try again.', verbose, type='error')
        # assign new init_disc_z value, which corresponds to the first vector of mask_points. Note, we need to substract from nz due to SAL orientation: in the viewer, orientation is S-I while in this code, it is I-S.
        init_disc = [nz-int(mask_points.split(',')[0]), 2]

    # display init disc
    if verbose == 2:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        plt.matshow(np.mean(data[xc-param.size_RL:xc+param.size_RL, :, :], axis=0).transpose(), fignum=50, cmap=plt.cm.gray, clim=[0, 800], origin='lower')
        plt.title('Anatomical image')
        plt.autoscale(enable=False)  # to prevent autoscale of axis when displaying plot
        plt.figure(50), plt.scatter(yc + param.shift_AP_visu, init_disc[0], c='yellow', s=50)
        plt.text(yc + param.shift_AP_visu + 4, init_disc[0], str(init_disc[1]) + '/' + str(init_disc[1] + 1),
                 verticalalignment='center', horizontalalignment='left', color='pink', fontsize=15), plt.draw()
        # plt.ion()  # enables interactive mode

    # FIND DISCS

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
def interpolate_im_to_ref(im_input, im_input_sc, new_res=0.3, sq_size_size_mm=22.5, interpolation_mode=3):
    nx, ny, nz, nt, px, py, pz, pt = im_input.dim

    im_input_sc = im_input_sc.copy()
    im_input= im_input.copy()

    # keep only spacing and origin in qform to avoid rotation issues
    input_qform = im_input.hdr.get_qform()
    for i in range(4):
        for j in range(4):
            if i!=j and j!=3:
                input_qform[i, j] = 0

    im_input.hdr.set_qform(input_qform)
    im_input.hdr.set_sform(input_qform)
    im_input_sc.hdr = im_input.hdr

    sq_size = int(sq_size_size_mm/new_res)
    # create a reference image : square of ones
    im_ref = Image(np.ones((sq_size, sq_size, 1), dtype=np.int), dim=(sq_size, sq_size, 1, 0, new_res, new_res, pz, 0), orientation='RPI')

    # copy input qform matrix to reference image
    im_ref.hdr.set_qform(im_input.hdr.get_qform())
    im_ref.hdr.set_sform(im_input.hdr.get_sform())

    # set correct header to reference image
    im_ref.hdr.set_data_shape((sq_size, sq_size, 1))
    im_ref.hdr.set_zooms((new_res, new_res, pz))

    # save image to set orientation to RPI (not properly done at the creation of the image)
    fname_ref = 'im_ref.nii.gz'
    im_ref.setFileName(fname_ref)
    im_ref.save()
    im_ref = set_orientation(im_ref, 'RPI', fname_out=fname_ref)

    # set header origin to zero to get physical coordinates of the center of the square
    im_ref.hdr.as_analyze_map()['qoffset_x'] = 0
    im_ref.hdr.as_analyze_map()['qoffset_y'] = 0
    im_ref.hdr.as_analyze_map()['qoffset_z'] = 0
    im_ref.hdr.set_sform(im_ref.hdr.get_qform())
    im_ref.hdr.set_qform(im_ref.hdr.get_qform())
    [[x_square_center_phys, y_square_center_phys, z_square_center_phys]] = im_ref.transfo_pix2phys(coordi=[[int(sq_size / 2), int(sq_size / 2), 0]])

    list_interpolate_images = []
    # iterate on z dimension of input image
    for iz in range(nz):
        # copy reference image: one reference image per slice
        im_ref_slice_iz = im_ref.copy()

        # get center of mass of SC for slice iz
        x_seg, y_seg = (im_input_sc.data[:, :, iz] > 0).nonzero()
        x_center, y_center = np.mean(x_seg), np.mean(y_seg)
        [[x_center_phys, y_center_phys, z_center_phys]] = im_input_sc.transfo_pix2phys(coordi=[[x_center, y_center, iz]])

        # center reference image on SC for slice iz
        im_ref_slice_iz.hdr.as_analyze_map()['qoffset_x'] = x_center_phys - x_square_center_phys
        im_ref_slice_iz.hdr.as_analyze_map()['qoffset_y'] = y_center_phys - y_square_center_phys
        im_ref_slice_iz.hdr.as_analyze_map()['qoffset_z'] = z_center_phys
        im_ref_slice_iz.hdr.set_sform(im_ref_slice_iz.hdr.get_qform())
        im_ref_slice_iz.hdr.set_qform(im_ref_slice_iz.hdr.get_qform())

        # interpolate input image to reference image
        im_input_interpolate_iz = im_input.interpolate_from_image(im_ref_slice_iz, interpolation_mode=interpolation_mode, border='nearest')
        # reshape data to 2D if needed
        if len(im_input_interpolate_iz.data.shape) == 3:
            im_input_interpolate_iz.data = im_input_interpolate_iz.data.reshape(im_input_interpolate_iz.data.shape[:-1])
        # add slice to list
        list_interpolate_images.append(im_input_interpolate_iz)

    return list_interpolate_images

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]

#.........这里部分代码省略.........
            if type_process == 'vert-continuous':
                self.output_image.save('float32')
            elif type_process != 'plan_ref':
                self.output_image.save('minimize_int')
            else:
                self.output_image.save()

    def add(self, value):
        """
        This function add a specified value to all non-zero voxels.
        """
        image_output = Image(self.image_input, self.verbose)
        # image_output.data *= 0
        coordinates_input = self.image_input.getNonZeroCoordinates()

        # for all points with non-zeros neighbors, force the neighbors to 0
        for i, coord in enumerate(coordinates_input):
            image_output.data[int(coord.x), int(coord.y), int(coord.z)] = image_output.data[int(coord.x), int(coord.y), int(coord.z)] + float(value)
        return image_output


    def create_label(self, add=False):
        """
        Create an image with labels listed by the user.
        This method works only if the user inserted correct coordinates.

        self.coordinates is a list of coordinates (class in msct_types).
        a Coordinate contains x, y, z and value.
        If only one label is to be added, coordinates must be completed with '[]'
        examples:
        For one label:  object_define=ProcessLabels( fname_label, coordinates=[coordi]) where coordi is a 'Coordinate' object from msct_types
        For two labels: object_define=ProcessLabels( fname_label, coordinates=[coordi1, coordi2]) where coordi1 and coordi2 are 'Coordinate' objects from msct_types
        """
        image_output = self.image_input.copy()
        if not add:
            image_output.data *= 0

        # loop across labels
        for i, coord in enumerate(self.coordinates):
            # display info
            sct.printv('Label #' + str(i) + ': ' + str(coord.x) + ',' + str(coord.y) + ',' + str(coord.z) + ' --> ' +
                       str(coord.value), 1)
            image_output.data[int(coord.x), int(coord.y), int(coord.z)] = coord.value

        return image_output


    def cross(self):
        """
        create a cross.
        :return:
        """
        output_image = Image(self.image_input, self.verbose)
        nx, ny, nz, nt, px, py, pz, pt = Image(self.image_input.absolutepath).dim

        coordinates_input = self.image_input.getNonZeroCoordinates()
        d = self.cross_radius  # cross radius in pixel
        dx = d / px  # cross radius in mm
        dy = d / py

        # clean output_image
        output_image.data *= 0

        cross_coordinates = self.get_crosses_coordinates(coordinates_input, dx, self.image_ref, self.dilate)

        for coord in cross_coordinates:

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
class Vesselness:
    def __init__(self, fname):
        self.image = Image(fname)

        self.pretreated_im = self.pretreat()

        self.result_im = self.image.copy()
        self.result_im.file_name = self.image.file_name + '_vessel_mask'
        self.result_im.path = './'

        self.hess = hessian(self.pretreated_im.data)
        self.vessel_mask = self.compute_vessel_mask()

        self.result_im.data = self.vessel_mask
        self.result_im.save()

    def pretreat(self):
        import scipy.ndimage.filters as scp_filters
        '''
        status, orientation = sct.run('sct_orientation.py -i ' + self.image.path + self.image.file_name + self.image.ext)
        orientation = orientation[4:7]
        if orientation != 'RPI':
            sct.run('sct_orientation.py -i ' + self.image.path + self.image.file_name + self.image.ext + ' -s RPI')
        fname = self.image.file_name[:-7] + '_RPI.nii.gz'
        '''
        fname = self.image.file_name

        pretreated = self.image.copy()
        pretreated.file_name = fname

        nx, ny, nz, nt, px, py, pz, pt = sct.get_dimension(pretreated.file_name)
        sc_size = 3  # in mm
        sc_npix = ((sc_size/px) + (sc_size/py))/2.0
        pretreated.data = scp_filters.gaussian_filter(pretreated.data, sigma=sc_npix)
        pretreated.file_name = pretreated.file_name + '_gaussian'
        pretreated.save()
        return pretreated


    def compute_vessel_mask(self):
        print 'COMPUTE VESSEL MASK'
        vessel_mask = np.zeros(self.pretreated_im.data.shape)
        print 'Image shape : ', str(self.pretreated_im.data.shape)
        for x in range(self.pretreated_im.data.shape[0]):
            print '----------> x =', x
            now_x = time.time()
            for y in range(self.pretreated_im.data.shape[1]):
                for z in range(self.pretreated_im.data.shape[2]):
                    # H = self.hess[x][y][z]
                    H = self.get_hessian(x, y, z)
                    vals, vects = np.linalg.eig(H)
                    l1, l2, l3 = sort_eigen_vals(vals)
                    vessel_mask[x][y][z] = line_filter(l1, l2, l3)
            print 'x ' + str(x) + ' --> in ' + str(time.time() - now_x) + ' sec'
        return vessel_mask

    def get_hessian(self, x, y, z):
        Hxx = self.hess[0][0][x][y][z]
        Hxy = self.hess[0][1][x][y][z]
        Hxz = self.hess[0][2][x][y][z]
        Hyx = self.hess[1][0][x][y][z]
        Hyy = self.hess[1][1][x][y][z]
        Hyz = self.hess[1][2][x][y][z]
        Hzx = self.hess[2][0][x][y][z]
        Hzy = self.hess[2][1][x][y][z]
        Hzz = self.hess[2][2][x][y][z]
        H = [[Hxx, Hxy, Hxz],
             [Hyx, Hyy, Hyz],
             [Hzx, Hzy, Hzz]]
        return H

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
def amu_processing(data_path):
    """
    get a segmentation image of the spinal cord an of the graymatter from a three level mask

    :param data_path: path to the data

    :return:
    """
    im_ext = '.nii'
    if data_path[-1] == '/':
        data_path = data_path[:-1]
    original_path = os.path.abspath('.')
    os.chdir(data_path)
    for subject_dir in os.listdir('.'):
        subject_path = data_path + '/' + subject_dir
        if os.path.isdir(subject_dir):
            os.chdir(subject_dir)
            sc_seg_list = []
            gm_seg_list = []
            im_list = []
            for file_name in os.listdir('.'):
                ext = sct.extract_fname(file_name)[2]
                if 'mask' in file_name and ext != '.hdr':
                    mask_im = Image(file_name)

                    sc_seg_im = mask_im.copy()
                    sc_seg_im.file_name = sct.extract_fname(file_name)[1][:-5] + '_manual_sc_seg'
                    sc_seg_im.ext = '.nii.gz'
                    sc_seg_im.data = (sc_seg_im.data > 1).astype(int)
                    # sc_seg_im = Image(param=sc_seg, absolutepath=subject_path + '/' + sct.extract_fname(file_name)[1][:-5] + '_manual_sc_seg.nii.gz')
                    # sc_seg_im.orientation = 'RPI'
                    sc_seg_im.save()
                    sc_seg_list.append(sc_seg_im.file_name + sc_seg_im.ext)

                    gm_seg_im = mask_im.copy()
                    gm_seg_im.file_name = sct.extract_fname(file_name)[1][:-5] + '_manual_gm_seg'
                    gm_seg_im.ext = '.nii.gz'
                    gm_seg_im.data = (gm_seg_im.data > 2).astype(int)
                    # gm_seg_im = Image(param=gm_seg, absolutepath=subject_path + '/' + sct.extract_fname(file_name)[1][:-5] + '_manual_gm_seg.nii.gz')
                    # gm_seg_im.orientation = 'RPI'
                    gm_seg_im.save()
                    gm_seg_list.append(gm_seg_im.file_name + gm_seg_im.ext)

                    im_list.append(file_name[:2] + im_ext)

            # merging the slice images into a 3D image
            im_list.reverse()
            gm_seg_list.reverse()
            sc_seg_list.reverse()
            cmd_merge = 'fslmerge -z '
            im_name = subject_dir + '_im.nii.gz '
            cmd_merge_im = cmd_merge + im_name
            gmseg_name = subject_dir + '_manual_gmseg.nii.gz '
            cmd_merge_gm_seg = cmd_merge + gmseg_name
            scseg_name = subject_dir + '_manual_scseg.nii.gz '
            cmd_merge_sc_seg = cmd_merge + scseg_name

            for im_i, gm_i, sc_i in zip(im_list, gm_seg_list, sc_seg_list):
                cmd_merge_im += im_i + ' '
                cmd_merge_gm_seg += gm_i + ' '
                cmd_merge_sc_seg += sc_i + ' '

            sct.run(cmd_merge_im)
            sct.run(cmd_merge_gm_seg)
            sct.run(cmd_merge_sc_seg)

            # creating a level image
            level_label = {0: '', 1: 'C1', 2: 'C2', 3: 'C3', 4: 'C4', 5: 'C5', 6: 'C6', 7: 'C7', 8: 'T1', 9: 'T2', 10: 'T3', 11: 'T4', 12: 'T5', 13: 'T6'}
            level_dat = np.zeros((mask_im.data.shape[0], mask_im.data.shape[1], len(im_list)))
            for i, im_i_name in enumerate(im_list):
                level_dat.T[:][:][i] = get_key_from_val(level_label, im_i_name[:2].upper())
            Image(param=level_dat, absolutepath=subject_dir + '_levels.nii.gz').save()

            # resampling
            resample_image(im_name)
            resample_image(gmseg_name, binary=True, thr=0.45)
            resample_image(scseg_name, binary=True, thr=0.55)

            # organizing data
            sct.run('mkdir original_data/')
            sct.run('mkdir extracted_data/')
            sct.run('mkdir 3d_data/')
            sct.run('mkdir 3d_resampled_data/')

            for file_name in os.listdir('.'):
                if 'mask' in file_name:
                    sct.run('mv ' + file_name + ' original_data/')
                elif 'manual' in file_name and 'G1' not in file_name:
                    sct.run('mv ' + file_name + ' extracted_data/')
                elif 'resampled.nii' in file_name:
                    sct.run('mv ' + file_name + ' 3d_resampled_data/')
                elif 'G1' in file_name:
                    sct.run('mv ' + file_name + ' 3d_data/')
                elif not os.path.isdir(os.path.abspath('.') + '/' + file_name):
                    sct.run('mv ' + file_name + ' original_data/')

            os.chdir('..')
    os.chdir(original_path)

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
def vanderbilt_processing(data_path):
    """
    get a segmentation image of the spinal cord an of the graymatter from a three level mask

    :param data_path: path to the data

    :return:
    """
    im_ext = '.nii.gz'
    if data_path[-1] == '/':
        data_path = data_path[:-1]
    original_path = os.path.abspath('.')
    os.chdir(data_path)
    for subject_dir in os.listdir('.'):
        if os.path.isdir(subject_dir):
            os.chdir(subject_dir)
            sc_seg_list = []
            gm_seg_list = []
            im_list = []
            for file_name in os.listdir('.'):
                if 'seg' in file_name:
                    mask_im = Image(file_name)

                    sc_seg_im = mask_im.copy()
                    sc_seg_im.file_name = sct.extract_fname(file_name)[1][:-4] + '_manual_sc_seg'
                    sc_seg_im.ext = '.nii.gz'
                    sc_seg_im.data = (sc_seg_im.data > 0).astype(int)
                    sc_seg_im.save()
                    sc_seg_list.append(sc_seg_im.file_name + sc_seg_im.ext)

                    gm_seg_im = mask_im.copy()
                    gm_seg_im.file_name = sct.extract_fname(file_name)[1][:-4] + '_manual_gm_seg'
                    gm_seg_im.ext = '.nii.gz'
                    gm_seg_im.data = (gm_seg_im.data > 1).astype(int)
                    gm_seg_im.save()
                    gm_seg_list.append(gm_seg_im.file_name + gm_seg_im.ext)

                    im_list.append(file_name[:17] + im_ext)

            # merging the slice images into a 3D image
            cmd_merge = 'fslmerge -z '
            im_name = subject_dir + '_im.nii.gz'
            cmd_merge_im = cmd_merge + im_name
            gmseg_name = subject_dir + '_manual_gmseg.nii.gz'
            cmd_merge_gm_seg = cmd_merge + gmseg_name
            scseg_name = subject_dir + '_manual_scseg.nii.gz'
            cmd_merge_sc_seg = cmd_merge + scseg_name

            for im_i, gm_i, sc_i in zip(im_list, gm_seg_list, sc_seg_list):
                cmd_merge_im += ' ' + im_i
                cmd_merge_gm_seg += ' ' + gm_i
                cmd_merge_sc_seg += ' ' + sc_i

            sct.run(cmd_merge_im)
            sct.run(cmd_merge_gm_seg)
            sct.run(cmd_merge_sc_seg)

            label_slices = [im_slice.split('_')[-1][2:4] for im_slice in im_list]
            i_slice_to_level = {0: 6, 1: 6, 2: 6, 3: 6, 4: 6, 5: 5, 6: 5, 7: 5, 8: 5, 9: 5, 10: 4, 11: 4, 12: 4, 13: 4, 14: 4, 15: 3, 16: 3, 17: 3, 18: 3, 19: 3, 20: 3, 21: 2, 22: 2, 23: 2, 24: 2, 25: 1, 26: 1, 27: 1, 28: 1, 29: 1}

            level_dat = np.zeros((mask_im.data.shape[0], mask_im.data.shape[1], len(im_list)))
            for i, l_slice in enumerate(label_slices):
                i_slice = int(l_slice) - 1
                level_dat.T[:][:][i] = i_slice_to_level[i_slice]
            Image(param=level_dat, absolutepath=subject_dir + '_levels.nii.gz').save()

            # resampling
            resample_image(im_name)
            resample_image(gmseg_name, binary=True)
            resample_image(scseg_name, binary=True)

            # organizing data
            sct.run('mkdir original_data/')
            sct.run('mkdir extracted_data/')
            sct.run('mkdir 3d_data/')
            sct.run('mkdir 3d_resampled_data/')
            sct.run('mkdir dic_data/')

            for file_name in os.listdir('.'):
                if '_manual_gm_seg' in file_name and 'sl' in file_name:
                    sct.run('cp ' + file_name + ' dic_data/')
                    sct.run('cp ' + file_name[:-21] + '.nii.gz dic_data/')
            for file_name in os.listdir('.'):
                if 'sl' in file_name and 'manual' not in file_name:
                    sct.run('mv ' + file_name + ' original_data/')
                elif 'manual' in file_name and 'sl' in file_name:
                    sct.run('mv ' + file_name + ' extracted_data/')
                elif 'resampled.nii' in file_name:
                    sct.run('mv ' + file_name + ' 3d_resampled_data/')
                elif '_sl' not in file_name and not os.path.isdir(os.path.abspath('.') + '/' + file_name) or 'level' in file_name:
                    sct.run('mv ' + file_name + ' 3d_data/')
                elif not os.path.isdir(os.path.abspath('.') + '/' + file_name):
                    sct.run('mv ' + file_name + ' original_data/')

            os.chdir('..')
    os.chdir(original_path)

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]

#.........这里部分代码省略.........
        # clean output_image
        output_image.data *= 0

        cross_coordinates = self.get_crosses_coordinates(coordinates_input, dx, self.image_ref, self.dilate)

        for coord in cross_coordinates:
            output_image.data[round(coord.x), round(coord.y), round(coord.z)] = coord.value

        return output_image
    # >>>

    def plan(self, width, offset=0, gap=1):
        """
        This function creates a plan of thickness="width" and changes its value with an offset and a gap between labels.
        """
        image_output = Image(self.image_input, self.verbose)
        image_output.data *= 0
        coordinates_input = self.image_input.getNonZeroCoordinates()

        # for all points with non-zeros neighbors, force the neighbors to 0
        for coord in coordinates_input:
            image_output.data[:,:,coord.z-width:coord.z+width] = offset + gap * coord.value

        return image_output

    def plan_ref(self):
        """
        This function generate a plan 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[:, :, coord.z] = -coord.value #PB: takes the int value of coord.value
        for coord in coordinates_input_pos:
            image_output.data[:, :, coord.z] = coord.value

        return image_output

    def cubic_to_point(self):
        """
        This function calculates the center of mass of each group of labels and returns a file of same size with only a
         label by group at the center of mass of this group.
        It is to be used after applying homothetic warping field to a label file as the labels will be dilated.
        Be careful: this algorithm computes the center of mass of voxels with same value, if two groups of voxels with
         the same value are present but separated in space, this algorithm will compute the center of mass of the two
         groups together.
        :return: image_output
        """
        from scipy import ndimage

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import copy [as 别名]
class MultiLabelRegistration:
    def __init__(
        self,
        fname_gm,
        fname_wm,
        path_template,
        fname_warp_template2target,
        param=None,
        fname_warp_target2template=None,
        apply_warp_template=0,
    ):
        if param is None:
            self.param = Param()
        else:
            self.param = param
        self.im_gm = Image(fname_gm)
        self.im_wm = Image(fname_wm)
        self.path_template = sct.slash_at_the_end(path_template, 1)
        if "MNI-Poly-AMU_GM.nii.gz" in os.listdir(self.path_template + "template/"):
            self.im_template_gm = Image(self.path_template + "template/MNI-Poly-AMU_GM.nii.gz")
            self.im_template_wm = Image(self.path_template + "template/MNI-Poly-AMU_WM.nii.gz")
            self.template = "MNI-Poly-AMU"
        else:
            self.im_template_gm = Image(self.path_template + "template/PAM50_gm.nii.gz")
            self.im_template_wm = Image(self.path_template + "template/PAM50_wm.nii.gz")
            self.template = "PAM50"

        # Previous warping fields:
        self.fname_warp_template2target = fname_warp_template2target
        self.fname_warp_target2template = fname_warp_target2template

        # new warping fields:
        self.fname_warp_template2gm = ""
        self.fname_wwarp_gm2template = ""

        # temporary fix - related to issue #871
        self.apply_warp_template = apply_warp_template

    def register(self):
        # accentuate separation WM/GM
        self.im_gm = thr_im(self.im_gm, 0.01, self.param.thr)
        self.im_wm = thr_im(self.im_wm, 0.01, self.param.thr)
        self.im_template_gm = thr_im(self.im_template_gm, 0.01, self.param.thr)
        self.im_template_wm = thr_im(self.im_template_wm, 0.01, self.param.thr)

        ## create multilabel images:
        # copy GM images to keep header information
        im_automatic_ml = self.im_gm.copy()
        im_template_ml = self.im_template_gm.copy()

        # create multi-label segmentation with GM*200 + WM*100 (100 and 200 encoded in self.param.gap)
        im_automatic_ml.data = self.param.gap[1] * self.im_gm.data + self.param.gap[0] * self.im_wm.data
        im_template_ml.data = (
            self.param.gap[1] * self.im_template_gm.data + self.param.gap[0] * self.im_template_wm.data
        )

        # set new names
        fname_automatic_ml = "multilabel_automatic_seg.nii.gz"
        fname_template_ml = "multilabel_template_seg.nii.gz"
        im_automatic_ml.setFileName(fname_automatic_ml)
        im_template_ml.setFileName(fname_template_ml)

        # Create temporary folder and put files in it
        tmp_dir = sct.tmp_create()

        path_gm, file_gm, ext_gm = sct.extract_fname(fname_gm)
        path_warp_template2target, file_warp_template2target, ext_warp_template2target = sct.extract_fname(
            self.fname_warp_template2target
        )

        convert(fname_gm, tmp_dir + file_gm + ext_gm)
        convert(fname_warp_template, tmp_dir + file_warp_template2target + ext_warp_template2target, squeeze_data=0)
        if self.fname_warp_target2template is not None:
            path_warp_target2template, file_warp_target2template, ext_warp_target2template = sct.extract_fname(
                self.fname_warp_target2template
            )
            convert(
                self.fname_warp_target2template,
                tmp_dir + file_warp_target2template + ext_warp_target2template,
                squeeze_data=0,
            )

        os.chdir(tmp_dir)
        # save images
        im_automatic_ml.save()
        im_template_ml.save()

        # apply template2image warping field
        if self.apply_warp_template == 1:
            fname_template_ml_new = sct.add_suffix(fname_template_ml, "_r")
            sct.run(
                "sct_apply_transfo -i "
                + fname_template_ml
                + " -d "
                + fname_automatic_ml
                + " -w "
                + file_warp_template2target
                + ext_warp_template2target
                + " -o "
                + fname_template_ml_new
#.........这里部分代码省略.........