DICOM的理解与学习2

nii格式数据的python库nibabel：https://ricardo_liu.gitbooks.io/nibabel/content/kuai-su-kai-shi.html

一、医学CT图像移除检查床

1、问题描述

    

 ​​​

 想要去除掉图像中的扫描床，只保留人体的扫描部分。

 2、解决思路

可能在这个例子中腐蚀和找最大连通区域的结果相似，但是两者的作用并不相同，比如，腐蚀是将人体与扫描床分离，不要黏在一起；而找最大联通区域的作用是在不连通的几个区域中，保留最大的区域。

3、代码实现

源代码来自：CT医学图像去床_qq_44976409的博客-CSDN博客

import matplotlib.pyplot as plt
import numpy as np
import pydicom
import cv2
import SimpleITK as sitk
from skimage.morphology import disk, rectangle, binary_dilation, binary_erosion, binary_closing, binary_opening, \
    rectangle, remove_small_objects

def read_dicom_data(file_name):
    file = sitk.ReadImage(file_name)
    data = sitk.GetArrayFromImage(file)
    data = np.squeeze(data, axis=0)
    data = np.int32(data)

    dicom_dataset = pydicom.dcmread(file_name)
    slice_location = dicom_dataset.SliceLocation         #获取层间距
    return data, data.shape[0], data.shape[1], slice_location

def window(img_data):
    img_data[img_data < -140] = -140
    img_data[img_data > 260] = 260

    return img_data

def find_max_region(mask_sel):
    contours, hierarchy = cv2.findContours(mask_sel, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
    # 找到最大区域并填充
    area = []
    for j in range(len(contours)):
        area.append(cv2.contourArea(contours[j]))
    max_idx = np.argmax(area)
    max_area = cv2.contourArea(contours[max_idx])
    for k in range(len(contours)):
        if k != max_idx:
            cv2.fillPoly(mask_sel, [contours[k]], 0)
    return mask_sel

#读入dicom图像
dcm_path = './1.dcm'
pixel_array, rows, columns, slice_location = read_dicom_data(dcm_path)
plt.subplot(241)
plt.title('Yuantu')
plt.imshow(pixel_array, cmap='gray')

#改变窗宽窗位
pixel_array = window(pixel_array)
imageData = (pixel_array - pixel_array.min()) * 255.0 / (pixel_array.max() - pixel_array.min())
imageData = np.uint8(imageData)
plt.subplot(242)
plt.title('BianChuang')
plt.imshow(imageData, cmap='gray')

#二值化
ret, binary = cv2.threshold(imageData,3,255,cv2.THRESH_BINARY)
#print(binary.shape)
#print(binary.dtype)
#binary = np.uint8(binary)      #cv2.findCoutours()函数接收的图像数据格式要是uint8的格式
plt.subplot(243)
plt.title('Erzhitu')
plt.imshow(binary, cmap='gray')

#腐蚀
selem = disk(5)
fushi = binary_erosion(binary, selem)
plt.subplot(244)
plt.title('Fushitu')
plt.imshow(fushi, cmap='gray')

#找最大连通区域
fushi = np.uint8(fushi)
max_region = find_max_region(fushi)
plt.subplot(245)
plt.title('Zuidaliantongtu')
plt.imshow(max_region, cmap='gray')

#膨胀
selem = disk(10)
pengzhang = binary_dilation(max_region, selem)
plt.subplot(246)
plt.title('pengzhangtu')
plt.imshow(pengzhang, cmap='gray')

'''
#开运算
selem = rectangle(1, 8)
Open = binary_opening(pengzhang, selem)
plt.subplot(248)
plt.title('result')
plt.imshow(Open, cmap='gray')
plt.show()
'''

#填充
pengzhang = np.uint8(pengzhang)
h, w = pengzhang.shape[:2]
#print(pengzhang.shape[:2])
mask_tp = np.zeros((h + 2, w + 2), np.uint8)
temp = pengzhang.copy()
temp2 = pengzhang.copy()
cv2.floodFill(temp, mask_tp, (1, 1), 255)
cv2.floodFill(temp2, mask_tp, (w - 2, h - 2), 255)
rt = cv2.bitwise_not(temp)

Tianchong = rt
Tianchong[rt > 0] = 255
plt.subplot(247)
plt.title('Tianchong')
plt.imshow(Tianchong, cmap='gray')

#处理
image_process = np.uint8((Tianchong / 255) * imageData)
plt.subplot(248)
plt.title('result')
plt.imshow(image_process, cmap='gray')
plt.show()

4、结果

 二、nii格式移除检查床

 移除图像中的检查床往往是为了更好地进行图像的后续处理，那么一张一张DICOM处理就会比较费力，所以就需要批量处理nii格式地图像。

import matplotlib.pyplot as plt
import numpy as np
import pydicom
import cv2
import SimpleITK as sitk
from skimage.morphology import disk, rectangle, binary_dilation, binary_erosion, binary_closing, binary_opening, \
    rectangle, remove_small_objects
import scipy
import nibabel as nib

'''
def read_nii_data(niipath):
    nii = nib.load(niipath)
    niiimg = nii.dataobj
    width, height, queue = niiimg.shape
    # img = niiimg[:, :, 0]
    niiimg = np.int32(niiimg)

    return niiimg,width,height,queue
'''

def read_nii_data(niipath):
    nii = sitk.ReadImage(niipath)
    niiimg = sitk.GetArrayFromImage(nii)
    queue, width, height = niiimg.shape
    # img = sitk.GetImageFromArray(niiimg[1])
    niiimg = np.int32(niiimg)

    return niiimg, width, height, queue

def window(window, img_data):
    if window == 'Lung':
        img_data[img_data < -1150] = -1150
        img_data[img_data > 350] = 350
    elif window == 'Pelvic':
        img_data[img_data < -140] = -140
        img_data[img_data > 260] = 260
    else:
        img_data[img_data < 0] = 0
        img_data[img_data > 80] = 80

    return img_data

def find_max_region(img_data):
    contours, hierarchy = cv2.findContours(img_data, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
    # 找到最大区域并填充
    area = []
    for j in range(len(contours)):
        area.append(cv2.contourArea(contours[j]))
    max_idx = np.argmax(area)
    max_area = cv2.contourArea(contours[max_idx])
    for k in range(len(contours)):
        if k != max_idx:
            cv2.fillPoly(img_data, [contours[k]], 0)

    return img_data

def Remove_Bed(img_data):

    # 二值化
    ret, binary = cv2.threshold(img_data, 3, 255, cv2.THRESH_BINARY)
    # print(binary.shape)
    # print(binary.dtype)
    # binary = np.uint8(binary)      #cv2.findCoutours()函数接收的图像数据格式要是uint8的格式

    # 腐蚀
    selem = disk(5)
    fushi = binary_erosion(binary, selem)

    # 找最大连通区域
    fushi = np.uint8(fushi)
    max_region = find_max_region(fushi)

    # 膨胀
    selem = disk(8)
    pengzhang = binary_dilation(max_region, selem)

    pengzhang = np.uint8(pengzhang)
    h, w = pengzhang.shape[:2]
    # print(pengzhang.shape[:2])
    mask_tp = np.zeros((h + 2, w + 2), np.uint8)
    temp = pengzhang.copy()
    temp2 = pengzhang.copy()
    cv2.floodFill(temp, mask_tp, (1, 1), 255)
    cv2.floodFill(temp2, mask_tp, (w - 2, h - 2), 255)
    rt = cv2.bitwise_not(temp)

    Tianchong = rt
    Tianchong[rt > 0] = 255

    image_process = np.uint8((Tianchong / 255) * img_data)
    return image_process

niipath = './0.nii'
pixel_array, rows, columns, queue = read_nii_data(niipath)
pixel_array = window("Lung", pixel_array)
imageData = (pixel_array - pixel_array.min()) * 255.0 / (pixel_array.max() - pixel_array.min())
imageData = np.uint8(imageData)

for i in range(queue):
    img = imageData[i, :, :]
  
    img = Remove_Bed(img)
    img = np.fliplr(img)        #img左右翻转
  
    imageData[i, :, :] = img

#写入对应的nii中
out = sitk.GetImageFromArray(imageData)
sitk.WriteImage(out, '00.nii')

 读取结果保存的nii格式图像（下图是用上面代码被注释掉的那段nibabel包读取的，且没有倒数第五行的翻转代码），应该是最后两行SimpleITK包读取时改变了图像矩阵，nibabel包是正常读取：

 下面是用SimpleITK包读取的，使用这个包会让图像矩阵发生变化，可以通过翻转或者旋转解决。

#一些零碎知识点

#添加维度和删除维度
np.expend_dim(img,axis=0)    #添加一个维度
#img的shape为（m,n,c）
#np.expand_dims(a, axis=1)将得到shape为(m, 1, n, c)的新数组，新数组中的元素与原数组a完全相同。
#np.expand_dims(a, axis=2)将得到shape为(m, n, 1, c)的新数组，新数组中的元素与原数组a完全相同。
#np.expand_dims(a, axis=3)将得到shape为(m, n, c, 1)的新数组，新数组中的元素与原数组a完全相同。

np.squeeze(img, axis=0)     #删除一个维度
#img的shape为（1,m,1）
#np.squeeze(data, axis=0)将得到shape为(m, 1)的新数组，新数组中的元素与原数组a完全相同。


#三维数组的切片
#用nibabel包读取的nii图像是(512,512,n)
niipath = './0.nii'
nii = nib.load(niipath)
niiimg = nii.dataobj
width, height, queue = niiimg.shape
img = niiimg[:, :, 0]

#用SimpleITK包读取的nii图像是(n,512,512)
nii = sitk.ReadImage(niipath)
niiimg = sitk.GetArrayFromImage(nii)
queue, width, height = niiimg.shape
img = niiimg[0, :, :]