一、數(shù)據(jù)格式

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

后綴為 .dcm。

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

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

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格式

一個(gè)raw通常有幾百兆，對(duì)應(yīng)的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的數(shù)組不是直接的像素值，而是相對(duì)于CT掃描中原點(diǎn)位置的差值，需要做進(jìn)一步轉(zhuǎn)換。

1.3 查看CT掃描文件軟件一個(gè)開源免費(fèi)的查看軟件 mango

二 dicom格式數(shù)據(jù)處理過程

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

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

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

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

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

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

# Load the scans in given folder path

def load_scan(path):

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

#對(duì)一個(gè)病人的所有slices進(jìn)行排序，x指的是一個(gè)slice。slice里面有好多屬性，

#有一個(gè)是ImagePositionPatient.按照他的這個(gè)屬性進(jìn)行對(duì)這些slices排序，方便我們組三維rendering。

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

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

for s in slices:

s.SliceThickness = slice_thickness #切片厚度

return slices

2.4 灰度值轉(zhuǎn)換為HU單元首先去除灰度值為-2000的pixl_array(pixl_array包含了真實(shí)數(shù)據(jù)),CT掃描邊界之外的灰度值固定為-2000(dicom和mhd都是這個(gè)值)。第一步設(shè)定這些值為0，當(dāng)前對(duì)應(yīng)為空氣(值為0).

回到HU單元，乘以rescale比率并加上intercept(存儲(chǔ)在掃描面的元數(shù)據(jù)中)。(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 重采樣不同掃描面的像素尺寸，粗細(xì)粒度是不同的，這不利于我們進(jìn)行CNN任務(wù)，我們可以使用同構(gòu)采樣。

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

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

# 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) #返回浮點(diǎn)數(shù)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') #使用所請(qǐng)求順序的樣條插值來縮放數(shù)組。

return image, new_spacing

# 現(xiàn)在重新取樣病人的像素，將其映射到一個(gè)同構(gòu)分辨率 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推進(jìn)立方體算法來查找3D體積數(shù)據(jù)中的曲面。

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) #創(chuàng)建3Dpoly

face_color = [0.5, 0.5, 1]

mesh.set_facecolor(face_color) #設(shè)置顏色

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

# 調(diào)用函數(shù)

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

2.6 輸出一個(gè)病人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 數(shù)據(jù)標(biāo)準(zhǔn)化處理歸一化處理：當(dāng)前的值范圍是[-1024,2000]。而任意大于400的值并不是處理肺結(jié)節(jié)需要考慮，因?yàn)樗鼈兌际遣煌瓷涿芏认碌墓穷^。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.

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

PIXEL_MEAN = 0.25

def zero_center(image):

image = image - PIXEL_MEAN

return image

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

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

3.1 載入必要的包

import SimpleITK as sitk

import numpy as np

import csv

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

import pandas as pd

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

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

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

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

#

# Helper function to get rows in data frame associated

# with each file

def get_filename(case):

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

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) #調(diào)用get_filename函數(shù)，并函數(shù)參數(shù)為df_node["seriesuid"]

df_node = df_node.dropna() #將所有含有nan項(xiàng)的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格式數(shù)據(jù)的值一直在尋找MHD格式數(shù)據(jù)的處理方法，對(duì)于dicom格式的CT有很多論文根據(jù)其HU值域可以輕易地分割肺、骨頭、血液等，但是對(duì)于MHD沒有這樣的參考。從LUNA16論壇得到的解釋是，LUNA16的MHD數(shù)據(jù)已經(jīng)轉(zhuǎn)換為HU值了，不需要再使用slope和intercept來做rescale變換了。此論壇主題下，有人提出MHD格式?jīng)]有提供pixel spacing(mm) 和 slice thickness(mm) ，而標(biāo)準(zhǔn)文件annotation.csv文件中結(jié)節(jié)的半徑和坐標(biāo)都是mm單位，最后確認(rèn)的是MHD格式文件中只保留了體素尺寸以及坐標(biāo)原點(diǎn)位置，沒有保存slice thickness。即，dicom才是原始數(shù)據(jù)格式。

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

坐標(biāo)體系變換如下：

相應(yīng)的代碼處理如下：

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的標(biāo)注CSV文件中標(biāo)注了結(jié)節(jié)中心的X,Y,Z軸坐標(biāo)，但是實(shí)際取值的時(shí)候取的是Z軸最后三層的數(shù)組(img_array)。

下述代碼只提取了包含結(jié)節(jié)的最后三個(gè)slice的數(shù)據(jù)，代碼參考自 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 查看節(jié)點(diǎn)以下代碼用于查看原始CT和結(jié)節(jié)mask。其實(shí)就是用matplotlib打印上一步存儲(chǔ)的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格式數(shù)據(jù)差不多，腐蝕膨脹、連通區(qū)域標(biāo)記等。

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

LUT全稱look up table，實(shí)際上就是一張像素灰度值的映射表，它將實(shí)際采樣到的像素灰度值經(jīng)過一定的變換如閾值、反轉(zhuǎn)、二值化、對(duì)比度調(diào)整、線性變換等，變成了另外一 個(gè)與之對(duì)應(yīng)的灰度值，這樣可以起到突出圖像的有用信息，增強(qiáng)圖像的光對(duì)比度的作用。

---------------------------------另外，我還整理了一份知乎萬贊的程序員學(xué)習(xí)大禮包，包括視頻教程、項(xiàng)目源碼、必看書籍、開發(fā)工具等【點(diǎn)擊此處一鍵取走】

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。