一、數據格式

1.1 dicomDICOM是醫學圖像中的標準文件，這些文件包含了諸多元數據信息(比如像素尺寸)，此處以kaggle Data Science Bowl數據集為例：data-science-bowl-2017,數據列表如下：

后綴為 .dcm。

每個病人的一次掃描CT(scan)可能有幾十到一百多個dcm數據文件(slices)。可以使用 python的dicom包讀取，讀取示例代碼如下:

dicom.read_file('/data/lung_competition/stage1/7050f8141e92fa42fd9c471a8b2f50ce/498d16aa2222d76cae1da144ddc59a13.dcm')其pixl_array包含了真實數據。

slices = [dicom.read_file(os.path.join(folder_name,filename)) for filename in os.listdir(folder_name)]

slices = np.stack([s.pixel_array for s in slices])

1.2 mhd格式

一個raw通常有幾百兆，對應的mhd文件只有1kb。mhd文件需要借助python的SimpleITK包來處理。SimpleITK 示例代碼如下:

import SimpleITK as sitk

itk_img = sitk.ReadImage(img_file)

img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)

num_z, height, width = img_array.shape #heightXwidth constitute the transverse plane

origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm)

spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)

需要注意的是，SimpleITK的img_array的數組不是直接的像素值，而是相對于CT掃描中原點位置的差值，需要做進一步轉換。

1.3 查看CT掃描文件軟件一個開源免費的查看軟件 mango

二 dicom格式數據處理過程

2.1 處理思路首先，需要明白的是醫學掃描圖像其實是三維圖像，使用代碼讀取之后查看不同的切面的切片(slices),可以從不同軸切割。

如下圖展示了一個病人CT掃描中，其中部分切片slices：

其次，CT掃描圖是包含了所有組織的，如果直接去看，看不到任何有用的信息，需要做一些預處理，預處理中一個重要概念是仿射劑量，衡量單位為HU(Hounsfield Unit),下表是不同放射劑量對應的組織器官：

Hounsfield Unit = pixel_value * rescale_slope + rescale_intercept一般情況rescale slope = 1, intercept = -1024。

上表中肺部組織的HU數值為-500,但通常是大于這個值，比如-320、-400。挑選出這些區域，然后做其他變換抽取出肺部像素點。

2.2 先載入必要的包

# -*- coding:utf-8 -*-

'''

this script is used for basic process of lung 2017 in Data Science Bowl

'''

import glob

import os

import pandas as pd

import SimpleITK as sitk

import 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 dicom

import scipy.misc

import numpy as np

2.3 將厚度加入到元數據如下代碼是載入一個掃描面，包含了多個(slices)，我們僅簡化的將其存儲為python列表，數據集中每個目錄都是一個掃描集(一個病人)。有個元數據域丟失，即Z軸方向上的像素尺寸，也即切片的厚度，所幸，我們可以用其他值推測出來，并加入到元數據中。

# Load the scans in given folder path

def load_scan(path):

slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]

#對一個病人的所有slices進行排序，x指的是一個slice。slice里面有好多屬性，

#有一個是ImagePositionPatient.按照他的這個屬性進行對這些slices排序，方便我們組三維rendering。

#imageOrientationPatient表示的是當前圖像的第一行在空間中的三維方向向量與第一列的三維方向向量。

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) #SliceLocation：表示的圖像平面的相對位置。

for s in slices:

s.SliceThickness = slice_thickness #切片厚度

return slices

2.4 灰度值轉換為HU單元首先去除灰度值為-2000的pixl_array(pixl_array包含了真實數據),CT掃描邊界之外的灰度值固定為-2000(dicom和mhd都是這個值)。第一步設定這些值為0，當前對應為空氣(值為0).

回到HU單元，乘以rescale比率并加上intercept(存儲在掃描面的元數據中)。(Hounsfield Unit = pixel_value * rescale_slope + rescale_intercept).

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)

image = image.astype(np.int16)

# Set outside-of-scan pixels to 0

# The intercept is usually -1024, so air is approximately 0

image[image == -2000] = 0

# Convert to Hounsfield units (HU)

for slice_number in range(len(slices)):

intercept = slices[slice_number].RescaleIntercept #Intercept

slope = slices[slice_number].RescaleSlope #Rescale

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)可以查看病人的掃描HU分布值情況：

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()

2.5 重采樣不同掃描面的像素尺寸，粗細粒度是不同的，這不利于我們進行CNN任務，我們可以使用同構采樣。

一個掃描面的像素區間可能是[2.5,0.5,0.5],即切片之間的距離為2.5mm。可能另外一個掃描面的范圍是[1.5,0.725,0.725]。這可能不利于自動分析。常見的處理方法是從全數據集中以固定的同構分辨率重新采樣，將所有的東西采樣為(1,1,1).

def resample(image, scan, new_spacing=[1,1,1]): # scan是load_scan函數返回的結果

# Determine current pixel spacing

spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))

spacing = np.array(list(spacing))

resize_factor = spacing / new_spacing

new_real_shape = image.shape * resize_factor

new_shape = np.round(new_real_shape) #返回浮點數x的四舍五入值。

real_resize_factor = new_shape / image.shape

new_spacing = spacing / real_resize_factor

image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest') #使用所請求順序的樣條插值來縮放數組。

return image, new_spacing

# 現在重新取樣病人的像素，將其映射到一個同構分辨率 1mm x1mm x1mm。

pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])使用matplotlib輸出肺部掃描的3D圖像方法。可能需要一兩分鐘。

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 = image.transpose(2,1,0) #將掃描件豎直放置

verts, faces = measure.marching_cubes(p, threshold) #Liner推進立方體算法來查找3D體積數據中的曲面。

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.1) #創建3Dpoly

face_color = [0.5, 0.5, 1]

mesh.set_facecolor(face_color) #設置顏色

ax.add_collection3d(mesh)

ax.set_xlim(0, p.shape[0])

ax.set_ylim(0, p.shape[1])

ax.set_zlim(0, p.shape[2])

plt.show()

# 調用函數

plot_3d(pix_resampled, 400)打印函數有個閾值(threshold)參數，來打印特定的結構，比如tissue或者骨頭。400是一個僅僅打印骨頭的閾值(HU對照表)，如下圖：

2.6 輸出一個病人scans中所有的slices

def plot_ct_scan(scan):

'''

plot a few more images of the slices

:param scan:

:return:

'''

f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(50, 50))

for i in range(0, scan.shape[0], 5):

plots[int(i / 20), int((i % 20) / 5)].axis('off')

plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone)此方法的效果示例如下:

2.7 數據標準化處理歸一化處理：當前的值范圍是[-1024,2000]。而任意大于400的值并不是處理肺結節需要考慮，因為它們都是不同反射密度下的骨頭。LUNA16競賽中常用來做歸一化處理的閾值集是-1000和400.以下代碼：

MIN_BOUND = -1000.0

MAX_BOUND = 400.0

def normalize(image):

image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)

image[image>1] = 1.

image[image<0] = 0.

return image0值中心化：簡單來說就是所有像素值減去均值。LUNA16競賽中的均值大約是0.25.

不要對每一張圖像做零值中心化(此處像是在kernel中完成的)CT掃描器返回的是校準后的精確HU計量。不會出現普通圖像中會出現某些圖像低對比度和明亮度的情況

PIXEL_MEAN = 0.25

def zero_center(image):

image = image - PIXEL_MEAN

return image

三 mhd格式數據處理過程mhd的數據只是格式與dicom不一樣，其實質包含的都是病人的掃描，處理MHD需要借助SimpleITK這個包，處理思路詳情可以參考Data Science Bowl2017的toturail Data Science Bowl 2017.需要注意的是MHD格式的數據沒有HU值，它的值域范圍與dicom很不同。

我們以LUNA2016年的數據處理流程為例。參考代碼為: LUNA2016數據切割.

3.1 載入必要的包

import SimpleITK as sitk

import numpy as np

import csv

from glob import glob #用它可以查找符合自己目的的文件

import pandas as pd

# glob方法返回所有匹配的文件路徑列表(list)；該方法需要一個參數用來指定匹配的路徑字符串，

# 其返回的文件名只包括當前目錄里的文件名，不包括子文件夾里的文件。

file_list=glob(luna_subset_path+"*.mhd")

#####################

#

# Helper function to get rows in data frame associated

# with each file

def get_filename(case):

# 如果你想要為一個定義在函數外的變量，那么你就得告訴Python這個變量名不是局部的，而是 全局 的。

global file_list

for f in file_list:

if case in f:

return(f)

#

# The locations of the nodes

df_node = pd.read_csv(luna_path+"annotations.csv")

df_node["file"] = df_node["seriesuid"].apply(get_filename) #調用get_filename函數，并函數參數為df_node["seriesuid"]

df_node = df_node.dropna() #將所有含有nan項的row刪除

#####

#

# Looping over the image files

#

fcount = 0

for img_file in file_list:

print "Getting mask for image file %s" % img_file.replace(luna_subset_path,"")

mini_df = df_node[df_node["file"]==img_file] #get all nodules associate with file

if len(mini_df)>0: # some files may not have a nodule--skipping those

biggest_node = np.argsort(mini_df["diameter_mm"].values)[-1] # just using the biggest node

node_x = mini_df["coordX"].values[biggest_node]

node_y = mini_df["coordY"].values[biggest_node]

node_z = mini_df["coordZ"].values[biggest_node]

diam = mini_df["diameter_mm"].values[biggest_node]

3.2 LUNA16的MHD格式數據的值一直在尋找MHD格式數據的處理方法，對于dicom格式的CT有很多論文根據其HU值域可以輕易地分割肺、骨頭、血液等，但是對于MHD沒有這樣的參考。從LUNA16論壇得到的解釋是，LUNA16的MHD數據已經轉換為HU值了，不需要再使用slope和intercept來做rescale變換了。此論壇主題下，有人提出MHD格式沒有提供pixel spacing(mm) 和 slice thickness(mm) ，而標準文件annotation.csv文件中結節的半徑和坐標都是mm單位，最后確認的是MHD格式文件中只保留了體素尺寸以及坐標原點位置，沒有保存slice thickness。即，dicom才是原始數據格式。

3.4 坐標體系變換MHD值的坐標體系是體素，以mm為單位(dicom的值是GV灰度值)。結節的位置是CT scanner坐標軸里面相對原點的mm值，需要將其轉換到真實坐標軸位置，可以使用SimpleITK包中的 GetOrigin() GetSpacing()。圖像數據是以512x512數組的形式給出的。

坐標體系變換如下：

相應的代碼處理如下：

itk_img = sitk.ReadImage(img_file)

img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)

center = np.array([node_x,node_y,node_z]) # nodule center

origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm)

spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)

# np.rint(a) 各元素四舍五入

v_center = np.rint((center-origin)/spacing) # nodule center in voxel space (still x,y,z ordering)在LUNA16的標注CSV文件中標注了結節中心的X,Y,Z軸坐標，但是實際取值的時候取的是Z軸最后三層的數組(img_array)。

下述代碼只提取了包含結節的最后三個slice的數據，代碼參考自 LUNA_mask_extraction.py

i = 0

for i_z in range(int(v_center[2])-1,int(v_center[2])+2):

mask = make_mask(center,diam,i_z*spacing[2]+origin[2],width,height,spacing,origin)

masks[i] = mask

imgs[i] = matrix2int16(img_array[i_z])

i+=1

np.save(output_path+"images_%d.npy" % (fcount) ,imgs)

np.save(output_path+"masks_%d.npy" % (fcount) ,masks)

3.5 查看節點以下代碼用于查看原始CT和結節mask。其實就是用matplotlib打印上一步存儲的npy文件。

import matplotlib.pyplot as plt

imgs = np.load(output_path+'images_0.npy')

masks = np.load(output_path+'masks_0.npy')

for i in range(len(imgs)):

print "image %d" % i

fig,ax = plt.subplots(2,2,figsize=[8,8])

ax[0,0].imshow(imgs[i],cmap='gray')

ax[0,1].imshow(masks[i],cmap='gray')

ax[1,0].imshow(imgs[i]*masks[i],cmap='gray')

plt.show()

raw_input("hit enter to cont : ")接下來的處理和DICOM格式數據差不多，腐蝕膨脹、連通區域標記等。

參考信息灰度值是pixel value經過重重LUT轉換得到的用來進行顯示的值，而這個轉換過程是不可逆的，也就是說，灰度值無法轉換為ct值。只能根據窗寬窗位得到一個大概的范圍。 pixel value經過modality lut得到Hu，但是懷疑pixelvalue的讀取出了問題。dicom文件中存在(0028，0106)(0028，0107)兩個tag，分別是最大最小pixel value，可以用來檢驗你讀取的pixel value 矩陣是否正確。

LUT全稱look up table，實際上就是一張像素灰度值的映射表，它將實際采樣到的像素灰度值經過一定的變換如閾值、反轉、二值化、對比度調整、線性變換等，變成了另外一 個與之對應的灰度值，這樣可以起到突出圖像的有用信息，增強圖像的光對比度的作用。

---------------------------------另外，我還整理了一份知乎萬贊的程序員學習大禮包，包括視頻教程、項目源碼、必看書籍、開發工具等【點擊此處一鍵取走】

推薦閱讀：干貨 | 共享免費資源整理(上)：學習資源篇?mp.weixin.qq.com干貨 | 共享免費資源整理(下)：程序員篇?mp.weixin.qq.com

總結

以上是生活随笔為你收集整理的图像处理中ct图的通道是多少_常见医疗扫描图像处理步骤的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。