查看病人的扫描HU值分布情况（DICOM文件格式）

文件格式）"/>

查看病人的扫描HU值分布情况（DICOM文件格式）

以一个病人的腹部淋巴结数据为例，查看病人的扫描HU值分布情况

# __author: Y
# date: 2020/1/8
# -*- coding:utf-8 -*-# 载入相应的包
import glob
import os
import pandas as pd
import SimpleITK as sitkimport numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import skimage, os
from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing
from skimage.measure import label,regionprops, perimeter
from skimage.morphology import binary_dilation, binary_opening
from skimage.filters import roberts, sobel
from skimage import measure, feature
from skimage.segmentation import clear_border
from skimage import data
from scipy import ndimage as ndi
import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import pydicom
import scipy.misc
import numpy as np# 载入一个扫描面，包含了多个切片（slices）。数据集中每个目录都是一个扫描面集（一个病人）
# 有个元数据域丢失，即z轴方向上的像素尺寸，也即切片的厚度，可以用其他值推测出来，并加入到元数据中。
# Load the scans in given folder path
# path是数据集的存放路径，根据需要换为自己数据集的存放路径
# path='E:/Python/Practice/ImageProcessing/ABD_LYMPH_001/09-14-2014-ABDLYMPH001-abdominallymphnodes-30274/abdominallymphnodes-26828'
def load_scan(path):# 获取切片slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]# 按ImagePositionPatient[2]排序，否则得到的扫描面是混乱无序的slices.sort(key=lambda x: int(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_thicknessreturn slices# 将图像的灰度值转换为CT值
def get_pixels_hu(slices):# 按顺序将切片的像素数组压入栈中image = np.stack([s.pixel_array for s in slices])# Convert to int16 (from sometimes int16),# should be possible as values should always be low enough (<32k)# 数据类型转换为int16类型的image = image.astype(np.int16)# Set outside-of-scan pixels to 0   将CT扫描边界之外的灰度值固定为-2000，对应空气为0# The intercept is usually -1024, so air is approximately 0image[image == -2000] = 0# Convert to Hounsfield units (HU)for slice_number in range(len(slices)):intercept = slices[slice_number].RescaleInterceptslope = slices[slice_number].RescaleSlopeif 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)# 查看病人的扫描HU值分布情况
INPUT_FOLDER = 'E:/Python/Practice/ImageProcessing/ABD_LYMPH_001/09-14-2014-ABDLYMPH001-abdominallymphnodes-30274/abdominallymphnodes-26828'
first_patient = load_scan(INPUT_FOLDER)
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()