pythondicom图像分割_处理医疗影像的Python利器：
PyDicom
Pydicom是⼀个⽤于处理DICOM格式⽂件的Python包，可以处理包括如医学图像（CT等）、报告等。
国家环保
一年级下册口算Pydicom⽀持DICOM格式的读取：可以将dicom⽂件读⼊python结构，同时⽀持修改后的数据集可以再次写⼊DICOM格式⽂件。但需要注意，它不是被设计为查看图像，主要是⽤来操作DICOM⽂件的各种数据元素。
PyDicom的安装
⽀持PIP和Conda安装，因此⾮常⽅便。以Anaconda安装为例：
conda install pydicom --channel conda-forge
读取Dicom⽂件
导⼊PyDicom包
import pydicom
读取⼀个dicom⽂件并显⽰
ds = pydicom.dcmread(file)
plt.figure(figsize=(10, 10))
plt.imshow(ds.pixel_array, bone)
plt.show()
⼀个完整的CT图像预处理的例⼦
对于CT图像，通常以患者的⼀次拍摄为⼀个⽂件夹，⽂件夹下有⼀序列的dicom⽂件，每个⽂件称为⼀个切⽚（slice）。但是每个患者的情况不同，所以slice间的间距不同，并且可能slice的排序也不同，因此需要在训练数据前做预处理。
这⾥参考Kaggle中的处理肺部照⽚的⼀个Kernel，对使⽤PyDicom处理CT图像做⼀个介绍参考评论⾥的反馈，这⾥将Kernel的Link以卡⽚形式放出：Full Preprocessing Tutorial | Kaggle w
1、导⼊libraries，将所有患者的⽬录列出来：
%matplotlib inline
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import dicom
import os
import scipy.ndimage
import matplotlib.pyplot as plt
from skimage import measure, morphology
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
# 包含所有患者⽬录的根⽬录
INPUT_FOLDER = '../input/sample_images/'
patients = os.listdir(INPUT_FOLDER)
patients.sort()
2、扫描⼀个患者的⽬录，加载所有的切⽚，按切换的z⽅向排序切⽚，并获取切⽚厚度
def load_scan(path):
slices = [ad_file(path + '/' + s) for s in os.listdir(path)]
slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
try:
slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
except:
slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
for s in slices:
s.SliceThickness = slice_thickness
return slices
CT扫描中的测量单位是Hounsfield单位（HU），它是辐射密度的量度。 CT扫描仪经过精⼼校准，可以精确测量。 来⾃维基百科：
3、默认情况下，从DICOM⽂件中获得的值是HU这个单位。 需要解决这个问题。
某些扫描仪具有圆柱扫描边界，但输出图像为⽅形。 落在这些边界之外的像素获得固定值-2000。 第⼀步是将这些值设置为0，当前对应于air。 接下来，回到HU单位，乘以重新缩放斜率并添加截距（⽅便地存储在扫描的元数据中！）。
def get_pixels_hu(slices):
image = np.stack([s.pixel_array for s in slices])
# 转换为int16，int16是ok的，因为所有的数值都应该 <32k
image = image.astype(np.int16)
# 设置边界外的元素为0
image[image == -2000] = 0
# 转换为HU单位
for slice_number in range(len(slices)):
intercept = slices[slice_number].RescaleIntercept
slope = slices[slice_number].RescaleSlope
if slope != 1:
image[slice_number] = slope * image[slice_number].astype(np.float64)
image[slice_number] = image[slice_number].astype(np.int16)
image[slice_number] += np.int16(intercept)
return np.array(image, dtype=np.int16)骨结核症状
4、查看⼀个患者的图像：
first_patient = load_scan(INPUT_FOLDER + patients[0])
first_patient_pixels = get_pixels_hu(first_patient)
plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
plt.xlabel("Hounsfield Units (HU)")
plt.ylabel("Frequency")什么中药补肾
plt.show()
# 显⽰⼀个中间位置的切⽚
plt.imshow(first_patient_pixels[80], ay)
plt.show()
5、重新采样
CT 扫描可能的像素间距为[2.5, 0.5, 0.5]，代表着切⽚间的距离是2.5毫⽶。对于不同的扫描，切⽚距离可能不同，对于⾃动分析是⼀个问题。
常⽤的⽅法是将整个数据集重新采样为相同分辨率的切⽚。例如将所有切⽚采样为[1 1 1]毫⽶的间距。这样就可以使⽤3D⽹格，⽽⽆需担⼼切⽚厚度的不确定性。
def resample(image, scan, new_spacing=[1,1,1]):
# Determine current pixel spacing
spacing = np.array([scan[0].SliceThickness] + scan[0].PixelSpacing, dtype=np.float32)
resize_factor = spacing / new_spacing
new_real_shape = image.shape * resize_factor
new_shape = np.round(new_real_shape)
real_resize_factor = new_shape / image.shape
new_spacing = spacing / real_resize_factor
image = scipy.(image, real_resize_factor, mode='nearest')
return image, new_spacing
pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])
print("Shape before resampling\t", first_patient_pixels.shape)
print("Shape after resampling\t", pix_resampled.shape)
6、画3D图像
显⽰扫描3D图像，对数据有个直观的感受对处理数据是有⽤的。
不幸的是，将使⽤⽴⽅体为我们的3D对象创建⼀个近似⽹格，并使⽤matplotlib绘制它。
def plot_3d(image, threshold=-300):
# Position the scan upright,
# so the head of the patient would be at the top facing the camera
p = anspo(2,1,0)
verts, faces = measure.marching_cubes(p, threshold)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
# Fancy indexing: `verts[faces]` to generate a collection of triangles
桂圆酒mesh = Poly3DCollection(verts[faces], alpha=0.70)
face_color = [0.45, 0.45, 0.75]
mesh.t_facecolor(face_color)
ax.add_collection3d(mesh)
市场营销案例
ax.t_xlim(0, p.shape[0])
ax.t_ylim(0, p.shape[1])
ax.t_zlim(0, p.shape[2])
plt.show()
绘图函数采⽤的阈值可以被⽤来绘制某些结构，例如所有组织或仅⾻骼。 400是仅显⽰⾻骼的阈值（可以参见上⾯的Hounsfield单位
表）。
星星的守护7、肺部切割
为了降低问题的复杂度，我们先对肺部进⾏切割。
它涉及很多巧妙的步骤。由区域⽣长和形态学的⼀系列操作组成。在这种情况下，我们将仅使⽤连通分量分析。阈值图像（-320 HU是⼀个很好的阈值，但这种⽅法并不重要）
做连接组件，确定⼈周围的空⽓标签，在⼆进制图像中⽤1s填充
可选：对于扫描中的每个轴向切⽚，确定最⼤的固体连接组件（⼈体周围的⾝体+空⽓），并将其他组件设置为0。这样可以填充⾯罩中肺部的结构。
只保留最⼤的⽓袋（⼈体在这⾥和那⾥都有其他的⽓袋）。
def largest_label_volume(im, bg=-1):
vals, counts = np.unique(im, return_counts=True)
counts = counts[vals != bg]
vals = vals[vals != bg]
if len(counts) > 0:
return vals[np.argmax(counts)]
el:
return None
def gment_lung_mask(image, fill_lung_structures=True):
# not actually binary, but 1 and 2.
# 0 is treated as background, which we do not want
binary_image = np.array(image > -320, dtype=np.int8)+1
labels = measure.label(binary_image)
# Pick the pixel in the very corner to determine which label is air.
# Improvement: Pick multiple background labels from around the patient
# More resistant to "trays" on which the patient lays cutting the air
# around the person in half
background_label = labels[0,0,0]
#Fill the air around the person
binary_image[background_label == labels] = 2
# Method of filling the lung structures (that is superior to something like
# morphological closing)
if fill_lung_structures:
# For every slice we determine the largest solid structure
for i, axial_slice in enumerate(binary_image):
怎么让头发蓬松axial_slice = axial_slice - 1
labeling = measure.label(axial_slice)
l_max = largest_label_volume(labeling, bg=0)
if l_max is not None: #This slice contains some lung
binary_image[i][labeling != l_max] = 1
binary_image -= 1 #Make the image actual binary
binary_image = 1-binary_image # Invert it, lungs are now 1
# Remove other air pockets insided body
labels = measure.label(binary_image, background=0)
l_max = largest_label_volume(labels, bg=0)
if l_max is not None: # There are air pockets
binary_image[labels != l_max] = 0
return binary_image
gmented_lungs = gment_lung_mask(pix_resampled, Fal)
gmented_lungs_fill = gment_lung_mask(pix_resampled, True)
plot_3d(gmented_lungs, 0)
看起来不错，但是我们可以进⼀步，在肺内包含结构可能是⼀个好主意（因为结节是固体），我们不希望在肺部只有空⽓。
plot_3d(gmented_lungs_fill, 0)
我们还可以看出两者之间的差异。
看起来不错！
当你想使⽤这个时，记得⾸先在它上⾯应⽤扩张形态学操作（即⽤圆形内核）。这会在所有⽅向上扩展蒙版。仅肺部的空⽓+结构将不包含所有结节，特别是它会遗漏那些粘在肺部侧⾯的结节，它们经常出现在那⾥！所以扩⼤⾯具⼀点:)
对于某些边缘情况，此分段可能会失败。它依赖于患者体外的空⽓不与肺部空⽓相连的事实。如果患者进⾏了⽓管造⼝术，情况也可能并⾮如此，不知道这是否存在于数据集中。此外，特别是噪声图像（例如由于下图中的起搏器），这种⽅法也可能失败。相反，⾝体中的第⼆⼤⽓袋将被分割。您可以通过检查蒙版对应的图像分数来识别这⼀点，对于这种情况，这将是⾮常⼩的。然后，你可以⾸先使⽤⼏毫⽶⼤⼩的内核进⾏形态学关闭操作以关闭这些孔，之后它应该可以⼯作（或者更简单地说，不要对此图像使⽤蒙版）。
