Python Image.data[ind_nonzero[0],ind_nonzero[1],iz]方法代码示例

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import data[ind_nonzero[0],ind_nonzero[1],iz] [as 别名]
def label_segmentation(fname_seg, list_disc_z, list_disc_value, verbose=1):
    """
    Label segmentation image
    :param fname_seg: fname of the segmentation
    :param list_disc_z: list of z that correspond to a disc
    :param list_disc_value: list of associated disc values
    :param verbose:
    :return:
    """
    # open segmentation
    seg = Image(fname_seg)
    dim = seg.dim
    ny = dim[1]
    nz = dim[2]
    # open labeled discs
    im_discs = Image(fname_seg)
    # loop across z
    for iz in range(nz):
        # get index of the disc right above iz
        try:
            ind_above_iz = max([i for i in range(len(list_disc_z)) if list_disc_z[i] > iz])
        except ValueError:
            # if ind_above_iz is empty, attribute value 0
            vertebral_level = 0
        else:
            # assign vertebral level (add one because iz is BELOW the disk)
            vertebral_level = list_disc_value[ind_above_iz] + 1
            # print vertebral_level
        # get voxels in mask
        ind_nonzero = np.nonzero(seg.data[:, :, iz])
        seg.data[ind_nonzero[0], ind_nonzero[1], iz] = vertebral_level
        if verbose == 2:
            import matplotlib
            matplotlib.use('Agg')
            import matplotlib.pyplot as plt
            plt.figure(50)
            plt.scatter(int(round(ny/2)), iz, c=vertebral_level, vmin=min(list_disc_value), vmax=max(list_disc_value), cmap='prism', marker='_', s=200)
    # write file
    seg.file_name += '_labeled'
    seg.save()

# 需要导入模块: from msct_image import Image [as 别名]
# 或者: from msct_image.Image import data[ind_nonzero[0],ind_nonzero[1],iz] [as 别名]
def vertebral_detection(fname, fname_seg, init_disc, verbose):

    shift_AP = 17  # shift the centerline towards the spine (in mm).
    size_AP = 4  # window size in AP direction (=y) in mm
    size_RL = 7  # window size in RL direction (=x) in mm
    size_IS = 7  # window size in RL direction (=z) in mm
    searching_window_for_maximum = 5  # size used for finding local maxima
    thr_corr = 0.2  # disc correlation threshold. Below this value, use template distance.
    gaussian_std_factor = 5  # the larger, the more weighting towards central value. This value is arbitrary-- should adjust based on large dataset
    fig_anat_straight = 1  # handle for figure
    fig_pattern = 2  # handle for figure
    fig_corr = 3  # handle for figure
    # define mean distance between adjacent discs: C1/C2 -> C2/C3, C2/C3 -> C4/C5, ..., L1/L2 -> L2/L3.
    mean_distance = np.array([18, 16, 17.0000, 16.0000, 15.1667, 15.3333, 15.8333,   18.1667,   18.6667,   18.6667,
    19.8333,   20.6667,   21.6667,   22.3333,   23.8333,   24.1667,   26.0000,   28.6667,   30.5000,   33.5000,
    33.0000,   31.3330])


    if verbose == 2:
        import matplotlib.pyplot as plt
        plt.ion()  # enables interactive mode

    # open anatomical volume
    img = Image(fname)
    data = img.data

    # smooth data
    from scipy.ndimage.filters import gaussian_filter
    data = gaussian_filter(data, [3, 1, 0], output=None, mode="reflect")

    # 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

    # define z: vector of indices along spine
    z = range(nz)
    # define xc and yc (centered in the field of view)
    xc = int(round(nx/2))  # direction RL
    yc = int(round(ny/2))  # direction AP

    # display stuff
    if verbose == 2:
        plt.matshow(np.mean(data[xc-size_RL:xc+size_RL, :, :], axis=0).transpose(), fignum=fig_anat_straight, cmap=plt.cm.gray, origin='lower')
        plt.title('Anatomical image')
        plt.autoscale(enable=False)  # to prevent autoscale of axis when displaying plot
        plt.figure(fig_anat_straight), plt.scatter(yc+shift_AP, init_disc[0], c='y', s=50)  # display init disc
        plt.text(yc+shift_AP+4, init_disc[0], 'init', verticalalignment='center', horizontalalignment='left', color='yellow', fontsize=15), plt.draw()


    # FIND DISCS
    # ===========================================================================
    printv('\nDetect intervertebral discs...', verbose)
    # assign initial z and disc
    current_z = init_disc[0]
    current_disc = init_disc[1]
    # adjust to pix size
    mean_distance = mean_distance * pz
    mean_distance_real = np.zeros(len(mean_distance))
    # create list for z and disc
    list_disc_z = []
    list_disc_value = []
    # do local adjustment to be at the center of the disc
    printv('.. local adjustment to center disc', verbose)
    pattern = data[xc-size_RL:xc+size_RL+1, yc+shift_AP-size_AP:yc+shift_AP+size_AP+1, current_z-size_IS:current_z+size_IS+1]
    current_z = local_adjustment(xc, yc, current_z, current_disc, data, size_RL, shift_AP, size_IS, searching_window_for_maximum, verbose)
    if verbose == 2:
        plt.figure(fig_anat_straight), plt.scatter(yc+shift_AP, current_z, c='g', s=50)
        plt.text(yc+shift_AP+4, current_z, str(current_disc)+'/'+str(current_disc+1), verticalalignment='center', horizontalalignment='left', color='green', fontsize=15)
        # plt.draw()
    # append value to main list
    list_disc_z = np.append(list_disc_z, current_z).astype(int)
    list_disc_value = np.append(list_disc_value, current_disc).astype(int)
    # update initial value (used when switching disc search to inferior direction)
    init_disc[0] = current_z
    # define mean distance to next disc
    approx_distance_to_next_disc = int(round(mean_distance[current_disc]))
    # find_disc(data, current_z, current_disc, approx_distance_to_next_disc, direction)
    # loop until potential new peak is inside of FOV
    direction = 'superior'
    search_next_disc = True
    while search_next_disc:
        printv('Current disc: '+str(current_disc)+' (z='+str(current_z)+'). Direction: '+direction, verbose)
        # Get pattern centered at z = current_z
        pattern = data[xc-size_RL:xc+size_RL+1, yc+shift_AP-size_AP:yc+shift_AP+size_AP+1, current_z-size_IS:current_z+size_IS+1]
        pattern1d = pattern.ravel()
        # display pattern
        if verbose == 2:
            plt.figure(fig_pattern)
            plt.matshow(np.flipud(np.mean(pattern[:, :, :], axis=0).transpose()), fignum=fig_pattern, cmap=plt.cm.gray)
            plt.title('Pattern in sagittal averaged across R-L')
        # compute correlation between pattern and data within
        printv('.. approximate distance to next disc: '+str(approx_distance_to_next_disc)+' mm', verbose)
#.........这里部分代码省略.........

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

    shift_AP = 15  # shift the centerline towards the spine (in mm).
    size_AP = 3  # mean around the centerline in the anterior-posterior direction in mm
    size_RL = 5  # mean around the centerline in the right-left direction in mm
    verbose = param.verbose

    if verbose == 2:
        import matplotlib.pyplot as plt
        plt.ion()  # enables interactive mode

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

    # matshow(img.data[:, :, 100]), show()

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

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

    # define z: vector of indices along spine
    z = range(nz)

    # define xc and yc (centered in the field of view)
    xc = round(nx/2)  # direction RL
    yc = round(ny/2)  # direction AP
    I = np.zeros((nz, 1))
    data_masked = img.data
    data = img.data
    for iz in range(nz):
        vox_in_spine = np.mgrid[xc-size_RL:xc+size_RL+1, yc+shift_AP-size_AP:yc+shift_AP+size_AP+1]
        # average intensity within box in the spine (shifted from spinal cord)
        I[iz] = np.mean(img.data[vox_in_spine[0, :, :].ravel().astype(int),
                                 vox_in_spine[1, :, :].ravel().astype(int),
                                 iz])
        # just for visualization:
        data_masked[vox_in_spine[0, :, :].ravel().astype(int),
                    vox_in_spine[1, :, :].ravel().astype(int),
                    iz] = 0

    # Display mask
    if verbose == 2:
        plt.matshow(np.flipud(data[xc, :, :].transpose()), cmap=plt.cm.gray)
        plt.title('Anatomical image')
        plt.draw()
        plt.matshow(np.flipud(data_masked[xc, :, :].transpose()), cmap=plt.cm.gray)
        plt.title('Anatomical image with mask')
        plt.draw()


    # display intensity along spine
    if verbose == 2:
        plt.figure()
        plt.plot(I)
        plt.title('Averaged intensity within spine. x=0: most caudal.')
        plt.draw()

    # find local extrema
    from scipy.signal import argrelextrema
    peaks = argrelextrema(I, np.greater, order=10)[0]
    nb_peaks = len(peaks)
    printv('.. Number of peaks found: '+str(nb_peaks), verbose)

    if verbose == 2:
        plt.figure()
        plt.plot(I)
        plt.plot(peaks, I[peaks], 'ro')
        plt.draw()

    # LABEL PEAKS
    labeled_peaks = np.array(range(nb_peaks, 0, -1)).astype(int)
    # find peak index closest to user input
    peak_ind_closest = np.argmin(abs(peaks-init_disc[0]))
    # build vector of peak labels
    # labeled_peaks = np.array(range(nb_peaks))
    # add the difference between "peak_ind_closest" and the init_disc value
    labeled_peaks = init_disc[1] - labeled_peaks[peak_ind_closest] + labeled_peaks

    # REMOVE WRONG LABELS (ASSUMING NO PEAK IS VISIBLE ABOVE C2/C3 DISK)
    ind_true_labels = np.where(labeled_peaks>1)[0]
    peaks = peaks[ind_true_labels]
    labeled_peaks = labeled_peaks[ind_true_labels]

    # ADD C1 label (ASSUMING DISTANCE FROM THE ADULT TEMPLATE)
    distance_c1_c2 = 20.8300/pz  # in mm
    # check if C2 disk is there
    if np.min(labeled_peaks) == 2:
        printv('.. C2 disk is present. Adding C2 vertebrae based on template...')
        peaks = np.append(peaks, (np.max(peaks) + distance_c1_c2).astype(int))
        labeled_peaks = np.append(labeled_peaks, 1)
    printv('.. Labeled peaks: '+str(labeled_peaks[:-1]), verbose)

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