Python Image.data[vox_coord[0],vox_coord[1],vox_coord[2]]方法代码示例

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import data[vox_coord[0],vox_coord[1],vox_coord[2]] [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
#.........这里部分代码省略.........