Python Image.transfo_pix2phys方法代码示例

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

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register_seg(seg_input, seg_dest):
    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]):
        print iz
        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]
        )  # strangely, this is the inverse of x_displacement when the same equation defines y_displacement
        y_displacement[iz] = y_center_of_mass_output[iz] - y_center_of_mass_input[iz]

    return x_displacement, y_displacement

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

#.........这里部分代码省略.........
            sct.printv('WARNING: Slice #' + str(iz) + ' is empty. It will be ignored.', verbose, 'warning')

    # regularize rotation
    if not poly == 0 and rot == 1:
        from msct_smooth import polynomial_fit
        angle_src_dest_regularized = polynomial_fit(z_nonzero, angle_src_dest[z_nonzero], poly)[0]
        # display
        if verbose == 2:
            plt.plot(180 * angle_src_dest[z_nonzero] / np.pi)
            plt.plot(180 * angle_src_dest_regularized / np.pi, 'r', linewidth=2)
            plt.grid()
            plt.xlabel('z')
            plt.ylabel('Angle (deg)')
            plt.savefig(path_qc+'register2d_centermassrot_regularize_rotation.png')
            plt.close()
        # update variable
        angle_src_dest[z_nonzero] = angle_src_dest_regularized

    # initialize warping fields
    # N.B. forward transfo is defined in destination space and inverse transfo is defined in the source space
    warp_x = np.zeros(data_dest.shape)
    warp_y = np.zeros(data_dest.shape)
    warp_inv_x = np.zeros(data_src.shape)
    warp_inv_y = np.zeros(data_src.shape)

    # construct 3D warping matrix
    for iz in z_nonzero:
        print str(iz)+'/'+str(nz)+'..',
        # get indices of x and y coordinates
        row, col = np.indices((nx, ny))
        # build 2xn array of coordinates in pixel space
        coord_init_pix = np.array([row.ravel(), col.ravel(), np.array(np.ones(len(row.ravel())) * iz)]).T
        # convert coordinates to physical space
        coord_init_phy = np.array(im_src.transfo_pix2phys(coord_init_pix))
        # get centermass coordinates in physical space
        centermass_src_phy = im_src.transfo_pix2phys([[centermass_src[iz, :].T[0], centermass_src[iz, :].T[1], iz]])[0]
        centermass_dest_phy = im_src.transfo_pix2phys([[centermass_dest[iz, :].T[0], centermass_dest[iz, :].T[1], iz]])[0]
        # build rotation matrix
        R = np.matrix(((cos(angle_src_dest[iz]), sin(angle_src_dest[iz])), (-sin(angle_src_dest[iz]), cos(angle_src_dest[iz]))))
        # build 3D rotation matrix
        R3d = np.eye(3)
        R3d[0:2, 0:2] = R
        # apply forward transformation (in physical space)
        coord_forward_phy = np.array(np.dot((coord_init_phy - np.transpose(centermass_dest_phy)), R3d) + np.transpose(centermass_src_phy))
        # apply inverse transformation (in physical space)
        coord_inverse_phy = np.array(np.dot((coord_init_phy - np.transpose(centermass_src_phy)), R3d.T) + np.transpose(centermass_dest_phy))
        # display rotations
        if verbose == 2 and not angle_src_dest[iz] == 0:
            # compute new coordinates
            coord_src_rot = coord_src[iz] * R
            coord_dest_rot = coord_dest[iz] * R.T
            # generate figure
            plt.figure('iz=' + str(iz) + ', angle_src_dest=' + str(angle_src_dest[iz]), figsize=(9, 9))
            # plt.ion()  # enables interactive mode (allows keyboard interruption)
            # plt.title('iz='+str(iz))
            for isub in [221, 222, 223, 224]:
                # plt.figure
                plt.subplot(isub)
                # ax = matplotlib.pyplot.axis()
                if isub == 221:
                    plt.scatter(coord_src[iz][:, 0], coord_src[iz][:, 1], s=5, marker='o', zorder=10, color='steelblue',
                                alpha=0.5)
                    pcaaxis = pca_src[iz].components_.T
                    pca_eigenratio = pca_src[iz].explained_variance_ratio_
                    plt.title('src')
                elif isub == 222:

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register2d_columnwise(fname_src, fname_dest, fname_warp='warp_forward.nii.gz', fname_warp_inv='warp_inverse.nii.gz', verbose=0, path_qc='./', smoothWarpXY=1):
    """
    Column-wise non-linear registration of segmentations. Based on an idea from Allan Martin.
    - Assumes src/dest are segmentations (not necessarily binary), and already registered by center of mass
    - Assumes src/dest are in RPI orientation.
    - Split along Z, then for each slice:
    - scale in R-L direction to match src/dest
    - loop across R-L columns and register by (i) matching center of mass and (ii) scaling.
    :param fname_src:
    :param fname_dest:
    :param fname_warp:
    :param fname_warp_inv:
    :param verbose:
    :return:
    """

    # initialization
    th_nonzero = 0.5  # values below are considered zero

    # for display stuff
    if verbose == 2:
        import matplotlib
        matplotlib.use('Agg')  # prevent display figure
        import matplotlib.pyplot as plt

    # Get image dimensions and retrieve nz
    sct.printv('\nGet image dimensions of destination image...', verbose)
    nx, ny, nz, nt, px, py, pz, pt = Image(fname_dest).dim
    sct.printv('  matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz), verbose)
    sct.printv('  voxel size:  '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm', verbose)

    # Split source volume along z
    sct.printv('\nSplit input volume...', verbose)
    from sct_image import split_data
    im_src = Image('src.nii')
    split_source_list = split_data(im_src, 2)
    for im in split_source_list:
        im.save()

    # Split destination volume along z
    sct.printv('\nSplit destination volume...', verbose)
    im_dest = Image('dest.nii')
    split_dest_list = split_data(im_dest, 2)
    for im in split_dest_list:
        im.save()

    # open image
    data_src = im_src.data
    data_dest = im_dest.data

    if len(data_src.shape) == 2:
        # reshape 2D data into pseudo 3D (only one slice)
        new_shape = list(data_src.shape)
        new_shape.append(1)
        new_shape = tuple(new_shape)
        data_src = data_src.reshape(new_shape)
        data_dest = data_dest.reshape(new_shape)

    # initialize forward warping field (defined in destination space)
    warp_x = np.zeros(data_dest.shape)
    warp_y = np.zeros(data_dest.shape)

    # initialize inverse warping field (defined in source space)
    warp_inv_x = np.zeros(data_src.shape)
    warp_inv_y = np.zeros(data_src.shape)

    # Loop across slices
    sct.printv('\nEstimate columnwise transformation...', verbose)
    for iz in range(0, nz):
        print str(iz)+'/'+str(nz)+'..',

        # PREPARE COORDINATES
        # ============================================================
        # get indices of x and y coordinates
        row, col = np.indices((nx, ny))
        # build 2xn array of coordinates in pixel space
        # ordering of indices is as follows:
        # coord_init_pix[:, 0] = 0, 0, 0, ..., 1, 1, 1..., nx, nx, nx
        # coord_init_pix[:, 1] = 0, 1, 2, ..., 0, 1, 2..., 0, 1, 2
        coord_init_pix = np.array([row.ravel(), col.ravel(), np.array(np.ones(len(row.ravel())) * iz)]).T
        # convert coordinates to physical space
        coord_init_phy = np.array(im_src.transfo_pix2phys(coord_init_pix))
        # get 2d data from the selected slice
        src2d = data_src[:, :, iz]
        dest2d = data_dest[:, :, iz]
        # julien 20161105
        #<<<
        # threshold at 0.5
        src2d[src2d < th_nonzero] = 0
        dest2d[dest2d < th_nonzero] = 0
        # get non-zero coordinates, and transpose to obtain nx2 dimensions
        coord_src2d = np.array(np.where(src2d > 0)).T
        coord_dest2d = np.array(np.where(dest2d > 0)).T
        # here we use 0.5 as threshold for non-zero value
        # coord_src2d = np.array(np.where(src2d > th_nonzero)).T
        # coord_dest2d = np.array(np.where(dest2d > th_nonzero)).T
        #>>>

        # SCALING R-L (X dimension)
        # ============================================================
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
    # Apply straightening to labels
    sct.printv('\nApply straightening to labels...', verbose)
    sct.run('sct_apply_transfo -i landmarks_rpi_cross3x3.nii.gz -o landmarks_rpi_cross3x3_straight.nii.gz -d segmentation_rpi_crop_straight.nii.gz -w warp_curve2straight.nii.gz -x nn')

    # Convert landmarks from FLOAT32 to INT
    sct.printv('\nConvert landmarks from FLOAT32 to INT...', verbose)
    sct.run('isct_c3d landmarks_rpi_cross3x3_straight.nii.gz -type int -o landmarks_rpi_cross3x3_straight.nii.gz')

    # Remove labels that do not correspond with each others.
    sct.printv('\nRemove labels that do not correspond with each others.', verbose)
    sct.run('sct_label_utils -t remove-symm -i landmarks_rpi_cross3x3_straight.nii.gz -o landmarks_rpi_cross3x3_straight.nii.gz,template_label_cross.nii.gz -r template_label_cross.nii.gz')

    # Estimate affine transfo: straight --> template (landmark-based)'
    sct.printv('\nEstimate affine transfo: straight anat --> template (landmark-based)...', verbose)
    # converting landmarks straight and curved to physical coordinates
    image_straight = Image('landmarks_rpi_cross3x3_straight.nii.gz')
    landmark_straight = image_straight.getNonZeroCoordinates(sorting='value')
    image_template = Image('template_label_cross.nii.gz')
    landmark_template = image_template.getNonZeroCoordinates(sorting='value')
    # Reorganize landmarks
    points_fixed, points_moving = [], []
    landmark_straight_mean = []
    for coord in landmark_straight:
        if coord.value not in [c.value for c in landmark_straight_mean]:
            temp_landmark = coord
            temp_number = 1
            for other_coord in landmark_straight:
                if coord.hasEqualValue(other_coord) and coord != other_coord:
                    temp_landmark += other_coord
                    temp_number += 1
            landmark_straight_mean.append(temp_landmark / temp_number)

    for coord in landmark_straight_mean:
        point_straight = image_straight.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        points_moving.append([point_straight[0][0], point_straight[0][1], point_straight[0][2]])
    for coord in landmark_template:
        point_template = image_template.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        points_fixed.append([point_template[0][0], point_template[0][1], point_template[0][2]])

    # Register curved landmarks on straight landmarks based on python implementation
    sct.printv('\nComputing rigid transformation (algo=translation-scaling-z) ...', verbose)
    import msct_register_landmarks
    (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = \
        msct_register_landmarks.getRigidTransformFromLandmarks(
            points_fixed, points_moving, constraints='translation-scaling-z', show=False)

    # writing rigid transformation file
    text_file = open("straight2templateAffine.txt", "w")
    text_file.write("#Insight Transform File V1.0\n")
    text_file.write("#Transform 0\n")
    text_file.write("Transform: FixedCenterOfRotationAffineTransform_double_3_3\n")
    text_file.write("Parameters: %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f\n" % (
        1.0/rotation_matrix[0, 0], rotation_matrix[0, 1],     rotation_matrix[0, 2],
        rotation_matrix[1, 0],     1.0/rotation_matrix[1, 1], rotation_matrix[1, 2],
        rotation_matrix[2, 0],     rotation_matrix[2, 1],     1.0/rotation_matrix[2, 2],
        translation_array[0, 0],   translation_array[0, 1],   -translation_array[0, 2]))
    text_file.write("FixedParameters: %.9f %.9f %.9f\n" % (points_moving_barycenter[0],
                                                           points_moving_barycenter[1],
                                                           points_moving_barycenter[2]))
    text_file.close()

    # Apply affine transformation: straight --> template
    sct.printv('\nApply affine transformation: straight --> template...', verbose)
    sct.run('sct_concat_transfo -w warp_curve2straight.nii.gz,straight2templateAffine.txt -d template.nii -o warp_curve2straightAffine.nii.gz')
    sct.run('sct_apply_transfo -i data_rpi.nii -o data_rpi_straight2templateAffine.nii -d template.nii -w warp_curve2straightAffine.nii.gz')
    sct.run('sct_apply_transfo -i segmentation_rpi.nii.gz -o segmentation_rpi_straight2templateAffine.nii.gz -d template.nii -w warp_curve2straightAffine.nii.gz -x linear')

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

#.........这里部分代码省略.........
    sct.printv('\nApply straightening to labels...', verbose)
    sct.run('sct_apply_transfo -i '+ftmp_label+' -o '+add_suffix(ftmp_label, '_straight')+' -d '+add_suffix(ftmp_seg, '_straight')+' -w warp_curve2straight.nii.gz -x nn')
    ftmp_label = add_suffix(ftmp_label, '_straight')

    # Create crosses for the template labels and get coordinates
    sct.printv('\nCreate a 15 mm cross for the template labels...', verbose)
    template_image = Image(ftmp_template_label)
    coordinates_input = template_image.getNonZeroCoordinates(sorting='value')
    # jcohenadad, issue #628 <<<<<
    # landmark_template = ProcessLabels.get_crosses_coordinates(coordinates_input, gapxy=15)
    landmark_template = coordinates_input
    # >>>>>
    if verbose == 2:
        # TODO: assign cross to image before saving
        template_image.setFileName(add_suffix(ftmp_template_label, '_cross'))
        template_image.save(type='minimize_int')

    # Create crosses for the input labels into straight space and get coordinates
    sct.printv('\nCreate a 15 mm cross for the input labels...', verbose)
    label_straight_image = Image(ftmp_label)
    coordinates_input = label_straight_image.getCoordinatesAveragedByValue()  # landmarks are sorted by value
    # jcohenadad, issue #628 <<<<<
    # landmark_straight = ProcessLabels.get_crosses_coordinates(coordinates_input, gapxy=15)
    landmark_straight = coordinates_input
    # >>>>>
    if verbose == 2:
        # TODO: assign cross to image before saving
        label_straight_image.setFileName(add_suffix(ftmp_label, '_cross'))
        label_straight_image.save(type='minimize_int')

    # Reorganize landmarks
    points_fixed, points_moving = [], []
    for coord in landmark_straight:
        point_straight = label_straight_image.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        points_moving.append([point_straight[0][0], point_straight[0][1], point_straight[0][2]])

    for coord in landmark_template:
        point_template = template_image.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        points_fixed.append([point_template[0][0], point_template[0][1], point_template[0][2]])

    # Register curved landmarks on straight landmarks based on python implementation
    sct.printv('\nComputing rigid transformation (algo=translation-scaling-z) ...', verbose)

    import msct_register_landmarks
    # for some reason, the moving and fixed points are inverted between ITK transform and our python-based transform.
    # and for another unknown reason, x and y dimensions have a negative sign (at least for translation and center of rotation).
    if verbose == 2:
        show_transfo = True
    else:
        show_transfo = False
    (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = msct_register_landmarks.getRigidTransformFromLandmarks(points_moving, points_fixed, constraints='translation-scaling-z', show=show_transfo)
    # writing rigid transformation file
    text_file = open("straight2templateAffine.txt", "w")
    text_file.write("#Insight Transform File V1.0\n")
    text_file.write("#Transform 0\n")
    text_file.write("Transform: AffineTransform_double_3_3\n")
    text_file.write("Parameters: %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f\n" % (
        rotation_matrix[0, 0], rotation_matrix[0, 1], rotation_matrix[0, 2],
        rotation_matrix[1, 0], rotation_matrix[1, 1], rotation_matrix[1, 2],
        rotation_matrix[2, 0], rotation_matrix[2, 1], rotation_matrix[2, 2],
        -translation_array[0, 0], -translation_array[0, 1], translation_array[0, 2]))
    text_file.write("FixedParameters: %.9f %.9f %.9f\n" % (-points_moving_barycenter[0],
                                                           -points_moving_barycenter[1],
                                                           points_moving_barycenter[2]))
    text_file.close()

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register_images(fname_source, fname_dest, mask='', paramreg=Paramreg(step='0', type='im', algo='Translation', metric='MI', iter='5', shrink='1', smooth='0', gradStep='0.5'),
                    ants_registration_params={'rigid': '', 'affine': '', 'compositeaffine': '', 'similarity': '', 'translation': '','bspline': ',10', 'gaussiandisplacementfield': ',3,0',
                                              'bsplinedisplacementfield': ',5,10', 'syn': ',3,0', 'bsplinesyn': ',1,3'}, remove_tmp_folder = 1):
    """Slice-by-slice registration of two images.

    We first split the 3D images into 2D images (and the mask if inputted). Then we register slices of the two images
    that physically correspond to one another looking at the physical origin of each image. The images can be of
    different sizes but the destination image must be smaller thant the input image. We do that using antsRegistration
    in 2D. Once this has been done for each slices, we gather the results and return them.
    Algorithms implemented: translation, rigid, affine, syn and BsplineSyn.
    N.B.: If the mask is inputted, it must also be 3D and it must be in the same space as the destination image.

    input:
        fname_source: name of moving image (type: string)
        fname_dest: name of fixed image (type: string)
        mask[optional]: name of mask file (type: string) (parameter -x of antsRegistration)
        paramreg[optional]: parameters of antsRegistration (type: Paramreg class from sct_register_multimodal)
        ants_registration_params[optional]: specific algorithm's parameters for antsRegistration (type: dictionary)

    output:
        if algo==translation:
            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)
        if algo==rigid:
            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)
            theta_rotation: list of rotation angle in radian (and in ITK's coordinate system) for each slice (type: list)
        if algo==affine or algo==syn or algo==bsplinesyn:
            creation of two 3D warping fields (forward and inverse) that are the concatenations of the slice-by-slice
            warps.
    """
    # Extracting names
    path_i, root_i, ext_i = sct.extract_fname(fname_source)
    path_d, root_d, ext_d = sct.extract_fname(fname_dest)

    # set metricSize
    if paramreg.metric == 'MI':
        metricSize = '32'  # corresponds to number of bins
    else:
        metricSize = '4'  # corresponds to radius (for CC, MeanSquares...)


    # Get image dimensions and retrieve nz
    print '\nGet image dimensions of destination image...'
    nx, ny, nz, nt, px, py, pz, pt = Image(fname_dest).dim
    print '.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz)
    print '.. voxel size:  '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm'

    # Define x and y displacement as list
    x_displacement = [0 for i in range(nz)]
    y_displacement = [0 for i in range(nz)]
    theta_rotation = [0 for i in range(nz)]

    # create temporary folder
    print('\nCreate temporary folder...')
    path_tmp = 'tmp.'+time.strftime("%y%m%d%H%M%S")
    sct.create_folder(path_tmp)
    print '\nCopy input data...'
    sct.run('cp '+fname_source+ ' ' + path_tmp +'/'+ root_i+ext_i)
    sct.run('cp '+fname_dest+ ' ' + path_tmp +'/'+ root_d+ext_d)
    if mask:
        sct.run('cp '+mask+ ' '+path_tmp +'/mask.nii.gz')

    # go to temporary folder
    os.chdir(path_tmp)

    # Split input volume along z
    print '\nSplit input volume...'
    from sct_split_data import split_data
    split_data(fname_source, 2, '_z')

    # Split destination volume along z
    print '\nSplit destination volume...'
    split_data(fname_dest, 2, '_z')

    # Split mask volume along z
    if mask:
        print '\nSplit mask volume...'
        split_data('mask.nii.gz', 2, '_z')

    im_dest_img = Image(fname_dest)
    im_input_img = Image(fname_source)
    coord_origin_dest = im_dest_img.transfo_pix2phys([[0,0,0]])
    coord_origin_input = im_input_img.transfo_pix2phys([[0,0,0]])
    coord_diff_origin = (asarray(coord_origin_dest[0]) - asarray(coord_origin_input[0])).tolist()
    [x_o, y_o, z_o] = [coord_diff_origin[0] * 1.0/px, coord_diff_origin[1] * 1.0/py, coord_diff_origin[2] * 1.0/pz]

    if paramreg.algo == 'BSplineSyN' or paramreg.algo == 'SyN' or paramreg.algo == 'Affine':
        list_warp_x = []
        list_warp_x_inv = []
        list_warp_y = []
        list_warp_y_inv = []
        name_warp_final = 'Warp_total' #if modified, name should also be modified in msct_register (algo slicereg2d_bsplinesyn and slicereg2d_syn)

    # loop across slices
    for i in range(nz):
        # set masking
        num = numerotation(i)
        num_2 = numerotation(int(num) + int(z_o))
        if mask:
#.........这里部分代码省略.........

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
#
# generate_warping_field('data_T2_RPI.nii.gz', x_disp_2_smooth, y_disp_2_smooth, fname='warping_field_im_trans.nii.gz')
# sct.run('sct_apply_transfo -i data_RPI_registered_reg1.nii.gz -d data_T2_RPI.nii.gz -w warping_field_im_trans.nii.gz -o data_RPI_registered_reg2.nii.gz -x spline')


f_1 = "/Users/tamag/data/data_template/independant_templates/Results_magma/t2_avg_RPI.nii.gz"
f_2 = "/Users/tamag/data/data_template/independant_templates/Results_magma/t1_avg.independent_RPI_reg1_unpad.nii.gz"
f_3 = "/Users/tamag/data/data_template/independant_templates/Results_magma/t1_avg.independent_RPI.nii.gz"

os.chdir("/Users/tamag/data/data_template/independant_templates/Results_magma")

im_1 = Image(f_1)
im_2 = Image(f_2)

data_1 = im_1.data

coord_test1 = [[1,1,1]]
coord_test = [[1,1,1],[2,2,2],[3,3,3]]

coordi_phys = im_1.transfo_pix2phys(coordi=coord_test)
coordi_pix = im_1.transfo_phys2pix(coordi = coordi_phys)
bla

# im_3 = nibabel.load(f_3)
# data_3 = im_3.get_data()
# hdr_3 = im_3.get_header()
#
# data_f = data_3 - laplace(data_3)
#
# img_f = nibabel.Nifti1Image(data_f, None, hdr_3)
# nibabel.save(img_f, "rehauss.nii.gz")

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register_images(
    im_input,
    im_dest,
    mask="",
    paramreg=Paramreg(
        step="0", type="im", algo="Translation", metric="MI", iter="5", shrink="1", smooth="0", gradStep="0.5"
    ),
    remove_tmp_folder=1,
):

    path_i, root_i, ext_i = sct.extract_fname(im_input)
    path_d, root_d, ext_d = sct.extract_fname(im_dest)
    path_m, root_m, ext_m = sct.extract_fname(mask)

    # set metricSize
    if paramreg.metric == "MI":
        metricSize = "32"  # corresponds to number of bins
    else:
        metricSize = "4"  # corresponds to radius (for CC, MeanSquares...)

    # initiate default parameters of antsRegistration transformation
    ants_registration_params = {
        "rigid": "",
        "affine": "",
        "compositeaffine": "",
        "similarity": "",
        "translation": "",
        "bspline": ",10",
        "gaussiandisplacementfield": ",3,0",
        "bsplinedisplacementfield": ",5,10",
        "syn": ",3,0",
        "bsplinesyn": ",3,32",
    }

    # Get image dimensions and retrieve nz
    print "\nGet image dimensions of destination image..."
    nx, ny, nz, nt, px, py, pz, pt = sct.get_dimension(im_dest)
    print ".. matrix size: " + str(nx) + " x " + str(ny) + " x " + str(nz)
    print ".. voxel size:  " + str(px) + "mm x " + str(py) + "mm x " + str(pz) + "mm"

    # Define x and y displacement as list
    x_displacement = [0 for i in range(nz)]
    y_displacement = [0 for i in range(nz)]
    theta_rotation = [0 for i in range(nz)]
    matrix_def = [0 for i in range(nz)]

    # create temporary folder
    print ("\nCreate temporary folder...")
    path_tmp = "tmp." + time.strftime("%y%m%d%H%M%S")
    sct.create_folder(path_tmp)
    print "\nCopy input data..."
    sct.run("cp " + im_input + " " + path_tmp + "/" + root_i + ext_i)
    sct.run("cp " + im_dest + " " + path_tmp + "/" + root_d + ext_d)
    if mask:
        sct.run("cp " + mask + " " + path_tmp + "/mask.nii.gz")

    # go to temporary folder
    os.chdir(path_tmp)

    # Split input volume along z
    print "\nSplit input volume..."
    sct.run(sct.fsloutput + "fslsplit " + im_input + " " + root_i + "_z -z")
    # file_anat_split = ['tmp.anat_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]

    # Split destination volume along z
    print "\nSplit destination volume..."
    sct.run(sct.fsloutput + "fslsplit " + im_dest + " " + root_d + "_z -z")
    # file_anat_split = ['tmp.anat_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]

    # Split mask volume along z
    if mask:
        print "\nSplit mask volume..."
        sct.run(sct.fsloutput + "fslsplit mask.nii.gz mask_z -z")
        # file_anat_split = ['tmp.anat_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]

    im_dest_img = Image(im_dest)
    im_input_img = Image(im_input)
    coord_origin_dest = im_dest_img.transfo_pix2phys([[0, 0, 0]])
    coord_origin_input = im_input_img.transfo_pix2phys([[0, 0, 0]])
    coord_diff_origin_z = coord_origin_dest[0][2] - coord_origin_input[0][2]
    [[x_o, y_o, z_o]] = im_input_img.transfo_phys2pix([[0, 0, coord_diff_origin_z]])

    # loop across slices
    for i in range(nz):
        # set masking
        num = numerotation(i)
        num_2 = numerotation(int(num) + z_o)
        if mask:
            masking = "-x mask_z" + num + ".nii"
        else:
            masking = ""

        cmd = (
            "isct_antsRegistration "
            "--dimensionality 2 "
            "--transform "
            + paramreg.algo
            + "["
            + paramreg.gradStep
            + ants_registration_params[paramreg.algo.lower()]
#.........这里部分代码省略.........

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register_landmarks(fname_src, fname_dest, dof, fname_affine="affine.txt", verbose=1, path_qc="./"):
    """
    Register two NIFTI volumes containing landmarks
    :param fname_src: fname of source landmarks
    :param fname_dest: fname of destination landmarks
    :param dof: degree of freedom. Separate with "_". Example: Tx_Ty_Tz_Rx_Ry_Sz
    :param fname_affine: output affine transformation
    :param verbose: 0, 1, 2
    :return:
    """
    from msct_image import Image

    # open src label
    im_src = Image(fname_src)
    # coord_src = im_src.getNonZeroCoordinates(sorting='value')  # landmarks are sorted by value
    coord_src = im_src.getCoordinatesAveragedByValue()  # landmarks are sorted by value
    # open dest labels
    im_dest = Image(fname_dest)
    # coord_dest = im_dest.getNonZeroCoordinates(sorting='value')
    coord_dest = im_dest.getCoordinatesAveragedByValue()
    # Reorganize landmarks

    points_src, points_dest = [], []
    for coord in coord_src:
        point_src = im_src.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        # convert NIFTI to ITK world coordinate
        # points_src.append([point_src[0][0], point_src[0][1], point_src[0][2]])
        points_src.append([-point_src[0][0], -point_src[0][1], point_src[0][2]])
    for coord in coord_dest:
        point_dest = im_dest.transfo_pix2phys([[coord.x, coord.y, coord.z]])
        # convert NIFTI to ITK world coordinate
        # points_dest.append([point_dest[0][0], point_dest[0][1], point_dest[0][2]])
        points_dest.append([-point_dest[0][0], -point_dest[0][1], point_dest[0][2]])

    # display
    sct.printv("Labels src: " + str(points_src), verbose)
    sct.printv("Labels dest: " + str(points_dest), verbose)
    sct.printv("Degrees of freedom (dof): " + dof, verbose)

    if len(coord_src) != len(coord_dest):
        raise Exception(
            "Error: number of source and destination landmarks are not the same, so landmarks cannot be paired."
        )

    # estimate transformation
    # N.B. points_src and points_dest are inverted below, because ITK uses inverted transformation matrices, i.e., src->dest is defined in dest instead of src.
    # (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = getRigidTransformFromLandmarks(points_dest, points_src, constraints=dof, verbose=verbose, path_qc=path_qc)
    (rotation_matrix, translation_array, points_moving_reg, points_moving_barycenter) = getRigidTransformFromLandmarks(
        points_src, points_dest, constraints=dof, verbose=verbose, path_qc=path_qc
    )
    # writing rigid transformation file
    # N.B. x and y dimensions have a negative sign to ensure compatibility between Python and ITK transfo
    text_file = open(fname_affine, "w")
    text_file.write("#Insight Transform File V1.0\n")
    text_file.write("#Transform 0\n")
    text_file.write("Transform: AffineTransform_double_3_3\n")
    text_file.write(
        "Parameters: %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f %.9f\n"
        % (
            rotation_matrix[0, 0],
            rotation_matrix[0, 1],
            rotation_matrix[0, 2],
            rotation_matrix[1, 0],
            rotation_matrix[1, 1],
            rotation_matrix[1, 2],
            rotation_matrix[2, 0],
            rotation_matrix[2, 1],
            rotation_matrix[2, 2],
            translation_array[0, 0],
            translation_array[0, 1],
            translation_array[0, 2],
        )
    )
    text_file.write(
        "FixedParameters: %.9f %.9f %.9f\n"
        % (points_moving_barycenter[0], points_moving_barycenter[1], points_moving_barycenter[2])
    )
    text_file.close()

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register_images(
    im_input,
    im_dest,
    mask="",
    paramreg=Paramreg(
        step="0", type="im", algo="Translation", metric="MI", iter="5", shrink="1", smooth="0", gradStep="0.5"
    ),
    ants_registration_params={
        "rigid": "",
        "affine": "",
        "compositeaffine": "",
        "similarity": "",
        "translation": "",
        "bspline": ",10",
        "gaussiandisplacementfield": ",3,0",
        "bsplinedisplacementfield": ",5,10",
        "syn": ",3,0",
        "bsplinesyn": ",1,3",
    },
    remove_tmp_folder=1,
):
    """Slice-by-slice registration of two images.

    We first split the 3D images into 2D images (and the mask if inputted). Then we register slices of the two images
    that physically correspond to one another looking at the physical origin of each image. The images can be of
    different sizes but the destination image must be smaller thant the input image. We do that using antsRegistration
    in 2D. Once this has been done for each slices, we gather the results and return them.
    Algorithms implemented: translation, rigid, affine, syn and BsplineSyn.
    N.B.: If the mask is inputted, it must also be 3D and it must be in the same space as the destination image.

    input:
        im_input: name of moving image (type: string)
        im_dest: name of fixed image (type: string)
        mask[optional]: name of mask file (type: string) (parameter -x of antsRegistration)
        paramreg[optional]: parameters of antsRegistration (type: Paramreg class from sct_register_multimodal)
        ants_registration_params[optional]: specific algorithm's parameters for antsRegistration (type: dictionary)

    output:
        if algo==translation:
            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)
        if algo==rigid:
            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)
            theta_rotation: list of rotation angle in radian (and in ITK's coordinate system) for each slice (type: list)
        if algo==affine or algo==syn or algo==bsplinesyn:
            creation of two 3D warping fields (forward and inverse) that are the concatenations of the slice-by-slice
            warps.
    """
    # Extracting names
    path_i, root_i, ext_i = sct.extract_fname(im_input)
    path_d, root_d, ext_d = sct.extract_fname(im_dest)

    # set metricSize
    if paramreg.metric == "MI":
        metricSize = "32"  # corresponds to number of bins
    else:
        metricSize = "4"  # corresponds to radius (for CC, MeanSquares...)

    # Get image dimensions and retrieve nz
    print "\nGet image dimensions of destination image..."
    nx, ny, nz, nt, px, py, pz, pt = sct.get_dimension(im_dest)
    print ".. matrix size: " + str(nx) + " x " + str(ny) + " x " + str(nz)
    print ".. voxel size:  " + str(px) + "mm x " + str(py) + "mm x " + str(pz) + "mm"

    # Define x and y displacement as list
    x_displacement = [0 for i in range(nz)]
    y_displacement = [0 for i in range(nz)]
    theta_rotation = [0 for i in range(nz)]

    # create temporary folder
    print ("\nCreate temporary folder...")
    path_tmp = "tmp." + time.strftime("%y%m%d%H%M%S")
    sct.create_folder(path_tmp)
    print "\nCopy input data..."
    sct.run("cp " + im_input + " " + path_tmp + "/" + root_i + ext_i)
    sct.run("cp " + im_dest + " " + path_tmp + "/" + root_d + ext_d)
    if mask:
        sct.run("cp " + mask + " " + path_tmp + "/mask.nii.gz")

    # go to temporary folder
    os.chdir(path_tmp)

    # Split input volume along z
    print "\nSplit input volume..."
    sct.run(sct.fsloutput + "fslsplit " + im_input + " " + root_i + "_z -z")

    # Split destination volume along z
    print "\nSplit destination volume..."
    sct.run(sct.fsloutput + "fslsplit " + im_dest + " " + root_d + "_z -z")

    # Split mask volume along z
    if mask:
        print "\nSplit mask volume..."
        sct.run(sct.fsloutput + "fslsplit mask.nii.gz mask_z -z")

    im_dest_img = Image(im_dest)
    im_input_img = Image(im_input)
    coord_origin_dest = im_dest_img.transfo_pix2phys([[0, 0, 0]])
    coord_origin_input = im_input_img.transfo_pix2phys([[0, 0, 0]])
#.........这里部分代码省略.........

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def register_seg(seg_input, seg_dest, verbose=1):
    """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 in voxel unit.
    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] * seg_dest_data.shape[2]
    y_center_of_mass_input = [0] * seg_dest_data.shape[2]
    sct.printv('\nGet center of mass of the input segmentation for each slice '
               '(corresponding to a slice in the output segmentation)...', verbose)  # different if size of the two seg are different
    # TODO: select only the slices corresponding to the output segmentation

    # grab physical coordinates of destination origin
    coord_origin_dest = seg_dest_img.transfo_pix2phys([[0, 0, 0]])

    # grab the voxel coordinates of the destination origin from the source image
    [[x_o, y_o, z_o]] = seg_input_img.transfo_phys2pix(coord_origin_dest)

    # calculate center of mass for each slice of the input image
    for iz in xrange(seg_dest_data.shape[2]):
        # starts from z_o, which is the origin of the destination image in the source image
        x_center_of_mass_input[iz], y_center_of_mass_input[iz] = ndimage.measurements.center_of_mass(array(seg_input_data[:, :, z_o + iz]))

    # initialize data
    x_center_of_mass_output = [0] * seg_dest_data.shape[2]
    y_center_of_mass_output = [0] * seg_dest_data.shape[2]

    # calculate center of mass for each slice of the destination image
    sct.printv('\nGet center of mass of the destination segmentation for each slice ...', verbose)
    for iz in xrange(seg_dest_data.shape[2]):
        try:
            x_center_of_mass_output[iz], y_center_of_mass_output[iz] = ndimage.measurements.center_of_mass(array(seg_dest_data[:, :, iz]))
        except Exception as e:
            sct.printv('WARNING: Exception error in msct_register_regularized during register_seg:', 1, 'warning')
            print 'Error on line {}'.format(sys.exc_info()[-1].tb_lineno)
            print e

    # calculate displacement in voxel space
    x_displacement = [0] * seg_input_data.shape[2]
    y_displacement = [0] * seg_input_data.shape[2]
    sct.printv('\nGet displacement by voxel...', verbose)
    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, None

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def validate_scad(folder_input, contrast):
    """
    Expecting folder to have the following structure :
    errsm_01:
    - t2
    -- errsm_01.nii.gz or t2.nii.gz
    --
    :param folder_input:
    :return:
    """
    from sct_get_centerline import ind2sub
    import time
    import math
    import numpy

    t0 = time.time()

    current_folder = os.getcwd()
    os.chdir(folder_input)

    try:
        patients = next(os.walk('.'))[1]
        overall_distance = {}
        max_distance = {}
        standard_deviation = 0
        overall_std = {}
        rmse = {}
        for i in patients:
            directory = i + "/" + str(contrast)
            try:
                os.chdir(directory)
            except Exception, e:
                print str(i)+" : "+contrast+" directory not found"

            try:
                if os.path.isfile(i+"_"+contrast+".nii.gz"):
                    raw_image = Image(i+"_"+contrast+".nii.gz")
                elif os.path.isfile(contrast+".nii.gz"):
                    raw_image = Image(contrast+".nii.gz")
                else:
                    raise Exception("Patient scan not found")

                if os.path.isfile(i+"_"+contrast+"_manual_segmentation.nii.gz"):
                    raw_orientation = raw_image.change_orientation()
                    scad = SCAD(raw_image, contrast=contrast, rm_tmp_file=1, verbose=1)
                    scad.execute()

                    manual_seg = Image(i+"_"+contrast+"_manual_segmentation.nii.gz")
                    manual_orientation = manual_seg.change_orientation()

                    from scipy.ndimage.measurements import center_of_mass
                    # find COM
                    iterator = range(manual_seg.data.shape[2])
                    com_x = [0 for ix in iterator]
                    com_y = [0 for iy in iterator]

                    for iz in iterator:
                        com_x[iz], com_y[iz] = center_of_mass(manual_seg.data[:, :, iz])

                    centerline_scad = Image(i+"_"+contrast+"_centerline.nii.gz")
                    # os.remove(i+"_"+contrast+"_centerline.nii.gz")

                    centerline_scad.change_orientation()
                    distance = {}
                    for iz in range(1, centerline_scad.data.shape[2]-1):
                        ind1 = np.argmax(centerline_scad.data[:, :, iz])
                        X,Y = ind2sub(centerline_scad.data[:, :, iz].shape,ind1)
                        com_phys = np.array(manual_seg.transfo_pix2phys([[com_x[iz], com_y[iz], iz]]))
                        scad_phys = np.array(centerline_scad.transfo_pix2phys([[X, Y, iz]]))
                        distance_magnitude = np.linalg.norm([com_phys[0][0]-scad_phys[0][0], com_phys[0][1]-scad_phys[0][1], 0])
                        if math.isnan(distance_magnitude):
                            print "Value is nan"
                        else:
                            distance[iz] = distance_magnitude

                    f = open(i+"_"+contrast+"_results.txt", 'w+')
                    f.write("Patient,Slice,Distance")
                    for key, value in distance.items():
                        f.write(i+","+str(key)+","+str(value))

                    standard_deviation = np.std(np.array(distance.values()))
                    average = sum(distance.values())/len(distance)
                    root_mean_square = np.sqrt(np.mean(np.square(distance.values())))

                    f.write("\nAverage : "+str(average))
                    f.write("\nStandard Deviation : "+str(standard_deviation))

                    f.close()

                    overall_distance[i] = average
                    max_distance[i] = max(distance.values())
                    overall_std[i] = standard_deviation
                    rmse[i] = root_mean_square

                else:
                    printv("Cannot find the manual segmentation", type="warning")

            except Exception, e:
                print e.message

#.........这里部分代码省略.........

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import transfo_pix2phys [as 别名]
def validate_scad(folder_input):
    """
    Expecting folder to have the following structure :
    errsm_01:
    - t2
    -- errsm_01.nii.gz or t2.nii.gz
    :param folder_input:
    :return:
    """
    current_folder = os.getcwd()
    os.chdir(folder_input)
    try:
        patients = next(os.walk('.'))[1]
        for i in patients:
            if i != "errsm_01" and i !="errsm_02":
                directory = i + "/t2"
                os.chdir(directory)
                try:
                    if os.path.isfile(i+"_t2.nii.gz"):
                        raw_image = Image(i+"_t2.nii.gz")
                    elif os.path.isfile("t2.nii.gz"):
                        raw_image = Image("t2.nii.gz")
                    else:
                        raise Exception("t2.nii.gz or "+i+"_t2.nii.gz file is not found")

                    raw_orientation = raw_image.change_orientation()
                    SCAD(raw_image, contrast="t2", rm_tmp_file=1, verbose=1).test_debug()

                    manual_seg = Image(i+"_t2_manual_segmentation.nii.gz")
                    manual_orientation = manual_seg.change_orientation()

                    from scipy.ndimage.measurements import center_of_mass
                    # find COM
                    iterator = range(manual_seg.data.shape[2])
                    com_x = [0 for ix in iterator]
                    com_y = [0 for iy in iterator]

                    for iz in iterator:
                        com_x[iz], com_y[iz] = center_of_mass(manual_seg.data[:, :, iz])
                    #raw_image.change_orientation(raw_orientation)
                    #manual_seg.change_orientation(manual_orientation)

                    centerline_scad = Image(i+"_t2_centerline.nii.gz")
                    os.remove(i+"_t2_centerline.nii.gz")

                    centerline_scad.change_orientation()
                    distance = []
                    for iz in range(centerline_scad.data.shape[2]):
                        ind1 = np.argmax(centerline_scad.data[:, :, iz])
                        X,Y = scad.ind2sub(centerline_scad.data[:, :, i].shape,ind1)
                        com_phys = centerline_scad.transfo_pix2phys([[com_x[iz], com_y[iz], iz]])
                        scad_phys = centerline_scad.transfo_pix2phys([[X, Y, iz]])
                        distance_magnitude = np.linalg.norm(com_phys-scad_phys)
                        distance.append(distance_magnitude)



                    os.chdir(folder_input)

                except Exception, e:
                    print e.message
                pass
    except Exception, e:
        print e.message