Python Image.change_orientation方法代码示例

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

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

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

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

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
def set_orientation(fname_in, orientation, fname_out, inversion=False):
    if not inversion:
        sct.run('isct_orientation3d -i '+fname_in+' -orientation '+orientation+' -o '+fname_out, 0)
    else:
        from msct_image import Image
        input_image = Image(fname_in)
        input_image.change_orientation(orientation, True)
        input_image.setFileName(fname_out)
        input_image.save()
    # return full path
    return os.path.abspath(fname_out)

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
 def output_debug_file(self, img, data, file_name):
     """
     This method writes a nifti file that corresponds to a step in the algorithm for easy debug.
     The new nifti file uses the header from the the image passed as parameter
     :param data: data to be written to file
     :param file_name: filename...
     :return: None
     """
     if self.verbose == 2:
         current_folder = os.getcwd()
         # os.chdir(self.path_tmp)
         try:
             img = Image(img)
             img.data = data
             img.change_orientation(self.raw_orientation)
             img.file_name = file_name
             img.save()
         except Exception, e:
             print e

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
def label_discs(fname_seg_labeled, verbose=1):
    """
    Label discs from labaled_segmentation
    :param fname_seg_labeld: fname of the labeled segmentation
    :param verbose:
    :return:
    """
    # open labeled segmentation
    im_seg_labeled = Image(fname_seg_labeled)
    orientation_native = im_seg_labeled.change_orientation('RPI')
    nx, ny, nz = im_seg_labeled.dim[0], im_seg_labeled.dim[1], im_seg_labeled.dim[2]
    data_disc = np.zeros([nx, ny, nz])
    vertebral_level_previous = np.max(im_seg_labeled.data)
    # loop across z
    for iz in range(nz):
        # get 2d slice
        slice = im_seg_labeled.data[:, :, iz]
        # check if at least one voxel is non-zero
        if np.any(slice):
            slice_one = np.copy(slice)
            # set all non-zero values to 1
            slice_one[slice.nonzero()] = 1
            # compute center of mass
            cx, cy = [int(x) for x in np.round(center_of_mass(slice_one)).tolist()]
            # retrieve vertebral level
            vertebral_level = slice[cx, cy]
            # if smaller than previous level, then labeled as a disc
            if vertebral_level < vertebral_level_previous:
                # label disc
                # print 'iz='+iz+', disc='+vertebral_level
                data_disc[cx, cy, iz] = vertebral_level
            # update variable
            vertebral_level_previous = vertebral_level
    # save disc labeled file
    im_seg_labeled.file_name += '_disc'
    im_seg_labeled.data = data_disc
    im_seg_labeled.change_orientation(orientation_native)
    im_seg_labeled.save()

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
def create_label_z(fname_seg, z, value):
    """
    Create a label at coordinates x_center, y_center, z
    :param fname_seg: segmentation
    :param z: int
    :return: fname_label
    """
    fname_label = 'labelz.nii.gz'
    nii = Image(fname_seg)
    orientation_origin = nii.change_orientation('RPI')  # change orientation to RPI
    nx, ny, nz, nt, px, py, pz, pt = nii.dim  # Get dimensions
    # find x and y coordinates of the centerline at z using center of mass
    from scipy.ndimage.measurements import center_of_mass
    x, y = center_of_mass(nii.data[:, :, z])
    x, y = int(round(x)), int(round(y))
    nii.data[:, :, :] = 0
    nii.data[x, y, z] = value
    # dilate label to prevent it from disappearing due to nearestneighbor interpolation
    from sct_maths import dilate
    nii.data = dilate(nii.data, [3])
    nii.setFileName(fname_label)
    nii.change_orientation(orientation_origin)  # put back in original orientation
    nii.save()
    return fname_label

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
def main(args=None):

    # initializations
    initz = ''
    initcenter = ''
    initc2 = 'auto'
    param = Param()

    # check user arguments
    if not args:
        args = sys.argv[1:]

    # Get parser info
    parser = get_parser()
    arguments = parser.parse(sys.argv[1:])
    fname_in = arguments["-i"]
    fname_seg = arguments['-s']
    contrast = arguments['-c']
    path_template = sct.slash_at_the_end(arguments['-t'], 1)
    # if '-o' in arguments:
    #     file_out = arguments["-o"]
    # else:
    #     file_out = ''
    if '-ofolder' in arguments:
        path_output = sct.slash_at_the_end(os.path.abspath(arguments['-ofolder']), slash=1)
    else:
        path_output = sct.slash_at_the_end(os.path.abspath(os.curdir), slash=1)
    if '-initz' in arguments:
        initz = arguments['-initz']
    if '-initcenter' in arguments:
        initcenter = arguments['-initcenter']
    # if user provided text file, parse and overwrite arguments
    if '-initfile' in arguments:
        # open file
        file = open(arguments['-initfile'], 'r')
        initfile = ' '+file.read().replace('\n', '')
        arg_initfile = initfile.split(' ')
        for i in xrange(len(arg_initfile)):
            if arg_initfile[i] == '-initz':
                initz = [int(x) for x in arg_initfile[i+1].split(',')]
            if arg_initfile[i] == '-initcenter':
                initcenter = int(arg_initfile[i+1])
    if '-initc2' in arguments:
        initc2 = 'manual'
    if '-param' in arguments:
        param.update(arguments['-param'][0])
    verbose = int(arguments['-v'])
    remove_tmp_files = int(arguments['-r'])
    denoise = int(arguments['-denoise'])
    laplacian = int(arguments['-laplacian'])

    # if verbose, import matplotlib
    # if verbose == 2:
        # import matplotlib.pyplot as plt

    # create temporary folder
    printv('\nCreate temporary folder...', verbose)
    path_tmp = tmp_create(verbose=verbose)
    # path_tmp = '/Users/julien/Dropbox/documents/processing/20160813_wang/t12/tmp.160814213032_725693/'

    # Copying input data to tmp folder
    printv('\nCopying input data to tmp folder...', verbose)
    run('sct_convert -i '+fname_in+' -o '+path_tmp+'data.nii')
    run('sct_convert -i '+fname_seg+' -o '+path_tmp+'segmentation.nii.gz')

    # Go go temp folder
    os.chdir(path_tmp)

    # create label to identify disc
    printv('\nCreate label to identify disc...', verbose)
    initauto = False
    if initz:
        create_label_z('segmentation.nii.gz', initz[0], initz[1])  # create label located at z_center
    elif initcenter:
        # find z centered in FOV
        nii = Image('segmentation.nii.gz')
        nii.change_orientation('RPI')  # reorient to RPI
        nx, ny, nz, nt, px, py, pz, pt = nii.dim  # Get dimensions
        z_center = int(round(nz/2))  # get z_center
        create_label_z('segmentation.nii.gz', z_center, initcenter)  # create label located at z_center
    else:
        initauto = True
        # printv('\nERROR: You need to initialize the disc detection algorithm using one of these two options: -initz, -initcenter\n', 1, 'error')

    # Straighten spinal cord
    printv('\nStraighten spinal cord...', verbose)
    # check if warp_curve2straight and warp_straight2curve already exist (i.e. no need to do it another time)
    if os.path.isfile('../warp_curve2straight.nii.gz') and os.path.isfile('../warp_straight2curve.nii.gz') and os.path.isfile('../straight_ref.nii.gz'):
        # if they exist, copy them into current folder
        sct.printv('WARNING: Straightening was already run previously. Copying warping fields...', verbose, 'warning')
        shutil.copy('../warp_curve2straight.nii.gz', 'warp_curve2straight.nii.gz')
        shutil.copy('../warp_straight2curve.nii.gz', 'warp_straight2curve.nii.gz')
        shutil.copy('../straight_ref.nii.gz', 'straight_ref.nii.gz')
        # apply straightening
        sct.run('sct_apply_transfo -i data.nii -w warp_curve2straight.nii.gz -d straight_ref.nii.gz -o data_straight.nii')
    else:
        run('sct_straighten_spinalcord -i data.nii -s segmentation.nii.gz -r 0 -qc 0')

    # resample to 0.5mm isotropic to match template resolution
    printv('\nResample to 0.5mm isotropic...', verbose)
#.........这里部分代码省略.........

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
    warping_fields = [WarpingField(filename) for filename in warping_fields_filename]

    filenames_output = []

    # path_tmp = tmp_create()
    # tmp_copy_nifti(input_file, path_tmp, 'raw.nii')
    # run('cp '+warping_fields_filename[0]+' '+path_tmp)
    # chdir(path_tmp)
    run('mkdir images')
    run('mkdir niftis')
    while True:
        try:
            warping_fields[0].num_of_frames = number_of_frames
            image_output_iter, iteration = warping_fields[0].next()
            image_output_iter.save()
            filename_warp = image_output_iter.path + image_output_iter.file_name + image_output_iter.ext
            filename_output = "niftis/tmp.warped_image_" + str(iteration - 1) + image_output_iter.ext
            run("sct_apply_transfo -i " + input_file + " -d " + reference_image + " -w " + filename_warp +
                " -o " + filename_output)
            result=Image(filename_output)
            result.change_orientation()

            toimage(result.data[int(result.data.shape[0]/2)].squeeze(), cmin=0.0).save('images/'+extract_fname(filename_output)[1]+'.jpg')
            filenames_output.append(filename_output)
        except ValueError:
            printv('\nError during warping field generation...', 1, 'error')
        except StopIteration:
            printv('\nFinished iterations.')
            break

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
class ProcessLabels(object):
    def __init__(self, fname_label, fname_output=None, fname_ref=None, cross_radius=5, dilate=False,
                 coordinates=None, verbose=1, vertebral_levels=None, value=None):
        self.image_input = Image(fname_label, verbose=verbose)
        self.image_ref = None
        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.vertebral_levels = vertebral_levels
        self.dilate = dilate
        self.coordinates = coordinates
        self.verbose = verbose
        self.value = value

    def process(self, type_process):
        # for some processes, change orientation of input image to RPI
        change_orientation = False
        if type_process in ['vert-body', 'vert-disc', 'vert-continuous']:
            # get orientation of input image
            input_orientation = self.image_input.orientation
            # change orientation
            self.image_input.change_orientation('RPI')
            change_orientation = True
        if type_process == 'add':
            self.output_image = self.add(self.value)
        if type_process == 'cross':
            self.output_image = self.cross()
        if type_process == 'plan':
            self.output_image = self.plan(self.cross_radius, 100, 5)
        if type_process == 'plan_ref':
            self.output_image = self.plan_ref()
        if type_process == 'increment':
            self.output_image = self.increment_z_inverse()
        if type_process == 'disks':
            self.output_image = self.labelize_from_disks()
        if type_process == 'MSE':
            self.MSE()
            self.fname_output = None
        if type_process == 'remove':
            self.output_image = self.remove_label()
        if type_process == 'remove-symm':
            self.output_image = self.remove_label(symmetry=True)
        if type_process == 'centerline':
            self.extract_centerline()
        if type_process == 'create':
            self.output_image = self.create_label()
        if type_process == 'create-add':
            self.output_image = self.create_label(add=True)
        if type_process == 'display-voxel':
            self.display_voxel()
            self.fname_output = None
        if type_process == 'diff':
            self.diff()
            self.fname_output = None
        if type_process == 'dist-inter':  # second argument is in pixel distance
            self.distance_interlabels(5)
            self.fname_output = None
        if type_process == 'cubic-to-point':
            self.output_image = self.cubic_to_point()
        if type_process == 'vert-body':
            self.output_image = self.label_vertebrae(self.vertebral_levels)
        # if type_process == 'vert-disc':
        #     self.output_image = self.label_disc(self.vertebral_levels)
        # if type_process == 'label-vertebrae-from-disks':
        #     self.output_image = self.label_vertebrae_from_disks(self.vertebral_levels)
        if type_process == 'vert-continuous':
            self.output_image = self.continuous_vertebral_levels()

        # save the output image as minimized integers
        if self.fname_output is not None:
            self.output_image.setFileName(self.fname_output)
            if change_orientation:
                self.output_image.change_orientation(input_orientation)
            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
#.........这里部分代码省略.........

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

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [as 别名]
    def execute(self):
        print 'Execution of the SCAD algorithm'

        vesselness_file_name = "imageVesselNessFilter.nii.gz"
        raw_file_name = "raw.nii"

        if self.debug:
            import matplotlib.pyplot as plt # import for debug purposes

        # create tmp and copy input
        path_tmp = sct.tmp_create()
        sct.tmp_copy_nifti(self.input_image.absolutepath, path_tmp, raw_file_name)

        if self.vesselness_provided:
            sct.run('cp '+vesselness_file_name+' '+path_tmp+vesselness_file_name)
        os.chdir(path_tmp)

        # get input image information
        img = Image(raw_file_name)

        # save original orientation and change image to RPI
        self.raw_orientation = img.change_orientation()

        # get body symmetry
        sym = SymmetryDetector(raw_file_name, self.contrast, crop_xy=1)
        self.raw_symmetry = sym.execute()

        # vesselness filter
        if not self.vesselness_provided:
            sct.run('sct_vesselness -i '+raw_file_name+' -t ' + self._contrast)

        # load vesselness filter data and perform minimum path on it
        img = Image(vesselness_file_name)
        raw_orientation = img.change_orientation()
        self.minimum_path_data, self.J1_min_path, self.J2_min_path = get_minimum_path(img.data, invert=1, debug=1, smooth_factor=1)

        # Apply an exponent to the minimum path
        self.minimum_path_powered = np.power(self.minimum_path_data, self.minimum_path_exponent)

        # Saving in Image since smooth_minimal_path needs pixel dimensions
        img.data = self.minimum_path_powered

        # smooth resulting minimal path
        self.smoothed_min_path = smooth_minimal_path(img)

        # normalise symmetry values between 0 and 1
        normalised_symmetry = equalize_array_histogram(self.raw_symmetry)

        # multiply normalised symmetry data with the minimum path result
        self.spine_detect_data = np.multiply(self.smoothed_min_path.data, normalised_symmetry)

        # extract the centerline from the minimal path image
        centerline_with_outliers = get_centerline(self.spine_detect_data, self.spine_detect_data.shape)
        img.data = centerline_with_outliers
        img.change_orientation()
        img.file_name = "centerline_with_outliers"
        img.save()

        # use a b-spline to smooth out the centerline
        x, y, z, dx, dy, dz = smooth_centerline("centerline_with_outliers.nii.gz")

        # save the centerline
        centerline_dim = img.dim
        img.data = np.zeros(centerline_dim)
        for i in range(0, np.size(x)-1):
            img.data[int(x[i]), int(y[i]), int(z[i])] = 1

        img.change_orientation(raw_orientation)
        img.file_name = "centerline"
        img.save()

        # copy back centerline
        os.chdir('../')
        sct.tmp_copy_nifti(path_tmp + 'centerline.nii.gz',self.input_image.path,self.input_image.file_name+'_centerline'+self.input_image.ext)

        if self.rm_tmp_file == 1:
            import shutil
            shutil.rmtree(path_tmp)

        if self.produce_output:
            self.produce_output_files()

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import change_orientation [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 change_orientation [as 别名]
def scadMRValidation(algorithm, isPython=False, verbose=True):
    if not isinstance(algorithm, str) or not algorithm:
        print 'ERROR: You must provide the name of your algorithm as a string.'
        usage()

    import time
    import sct_utils as sct

    # creating a new folder with the experiment
    path_experiment = 'scad-experiment.'+algorithm+'.'+time.strftime("%y%m%d%H%M%S")
    #status, output = sct.run('mkdir '+path_experiment, verbose)

    # copying images from "data" folder into experiment folder
    sct.copyDirectory('data', path_experiment)

    # Starting validation
    os.chdir(path_experiment)
    # t1
    os.chdir('t1/')
    for subject_dir in os.listdir('./'):
        if os.path.isdir(subject_dir):
            os.chdir(subject_dir)

            # creating list of images and corresponding manual segmentation
            list_images = dict()
            for file_name in os.listdir('./'):
                if not 'manual_segmentation' in file_name:
                    for file_name_corr in os.listdir('./'):
                        if 'manual_segmentation' in file_name_corr and sct.extract_fname(file_name)[1] in file_name_corr:
                            list_images[file_name] = file_name_corr

            # running the proposed algorithm on images in the folder and analyzing the results
            for image, image_manual_seg in list_images.items():
                print image
                path_in, file_in, ext_in = sct.extract_fname(image)
                image_output = file_in+'_centerline'+ext_in
                if ispython:
                    try:
                        eval(algorithm+'('+image+', t1, verbose='+str(verbose)+')')
                    except Exception as e:
                        print 'Error during spinal cord detection on line {}:'.format(sys.exc_info()[-1].tb_lineno)
                        print 'Subject: t1/'+subject_dir+'/'+image
                        print e
                        sys.exit(2)
                else:
                    cmd = algorithm+' -i '+image+' -t t1'
                    if verbose:
                        cmd += ' -v'
                    status, output = sct.run(cmd, verbose=verbose)
                    if status != 0:
                        print 'Error during spinal cord detection on Subject: t1/'+subject_dir+'/'+image
                        print output
                        sys.exit(2)

                # analyzing the resulting centerline
                from msct_image import Image
                manual_segmentation_image = Image(image_manual_seg)
                manual_segmentation_image.change_orientation()
                centerline_image = Image(image_output)
                centerline_image.change_orientation()

                from msct_types import Coordinate
                # coord_manseg = manual_segmentation_image.getNonZeroCoordinates()
                coord_centerline = centerline_image.getNonZeroCoordinates()

                # check if centerline is in manual segmentation
                result_centerline_in = True
                for coord in coord_centerline:
                    if manual_segmentation_image.data[coord.x, coord.y, coord.z] == 0:
                        result_centerline_in = False
                        print 'failed on slice #' + str(coord.z)
                        break
                if result_centerline_in:
                    print 'OK: Centerline is inside manual segmentation.'
                else:
                    print 'FAIL: Centerline is outside manual segmentation.'


                # check the length of centerline compared to manual segmentation
                # import sct_process_segmentation as sct_seg
                # length_manseg = sct_seg.compute_length(image_manual_seg)
                # length_centerline = sct_seg.compute_length(image_output)
                # if length_manseg*0.9 <= length_centerline <= length_manseg*1.1:
                #     print 'OK: Length of centerline correspond to length of manual segmentation.'
                # else:
                #     print 'FAIL: Length of centerline does not correspond to length of manual segmentation.'
            os.chdir('..')