【肺结节分割】luna数据mhd图片理解
just keep swimming
mhd格式图片在医疗领域很常见，本文重点介绍luna数据集中：
- mhd图片读取
- mhd图片可视化
- mhd图片坐标转换
1. Imports
import SimpleITK as sitk # 医疗图片处理包
import numpy as np
import csv
import os
from PIL import Image
import matplotlib.pyplot as plt
%matplotlib inline2.读取mhd图片
LUNA16数据集中，一个病例对应一个raw文件和一个mhd文件，raw文件包含图片数值信息，大小在50M~250M左右; mhd文件很小，包含图片其他信息，如：CT坐标原点，像素间距等。
img_path = '/home/dataset/medical/luna16/subset0/1.3.6.1.4.1.14519.5.2.1.6279.6001.108197895896446896160048741492.mhd'
cand_path = '/home/dataset/medical/luna16/CSVFILES/candidates.csv'
def load_itk_image(filename):
itkimage = sitk.ReadImage(filename)
numpyImage = sitk.GetArrayFromImage(itkimage)
numpyOrigin = np.array(list(itkimage.GetOrigin())) # CT原点坐标
numpySpacing = np.array(list(itkimage.GetSpacing())) # CT像素间隔
return numpyImage, numpyOrigin, numpySpacing
numpyImage, numpyOrigin, numpySpacing = load_itk_image(img_path)
print numpyImage.shape # 维度为(slice,w,h)
print numpyOrigin
print numpySpacing
'''
(119, 512, 512)
[-182.5 -190. -313.75]
[0.74218798 0.74218798 2.5 ]
'''3. 可视化CT切片
每个病例的CT切片数不一定，通常一个病例有上百个切片。
slice = 60
image = np.squeeze(numpyImage[slice, ...]) # if the image is 3d, the slice is integer
plt.imshow(image,cmap='gray')
plt.show()
4. 坐标转换
CT数据的坐标都是基于CT扫描仪的，坐标系单位为mm，每个病例的原点和体素间距都不一样；真实坐标系下，原点坐标统一为（0，0，0），体素间距为1。
结节的位置是CT扫坐标系（又称世界坐标系）下相对原点的mm值，需要将其转换到真实坐标轴位置。通过计算CT坐标系下结节位置在原点的相对位置，可以计算出真实坐标系下结节的位置，具体公式为：CT坐标系下的（结节位置-原点位置）/切片体素间距 = 真实坐标系下的（结节位置-（0，0，0））/1
4.1 坐标转换示例
例如，对病例1.3.6.1.4.1.14519.5.2.1.6279.6001.108197895896446896160048741492.mhd，其真实坐标系结节位置计算如下：
CT坐标系下
原点坐标：[-182.5 -190. -313.75]
像素间隔：[0.74 0.74 2.5 ]
结节位置：[-100.57, 67.26,-231.82]
原点与结节位置相对差：[ 81.93 , 257.26, 81.93 ]
真实坐标系下
结节位置 = CT坐标系下的（结节位置-原点位置）/切片体素间距=[ 81.93 , 257.26, 81.93 ]/[0.74218798 0.74218798 2.5] = [ 110.39 , 346.6 ,32.77]
4.2 坐标转换示例代码
'''read csv as a list of lists'''
def readCSV(filename):
lines = []
with open(filename, "rb") as f:
csvreader = csv.reader(f)
for line in csvreader:
lines.append(line)
return lines
'''convert world coordinate to real coordinate'''
def worldToVoxelCoord(worldCoord, origin, spacing):
stretchedVoxelCoord = np.absolute(worldCoord - origin)
voxelCoord = stretchedVoxelCoord / spacing
return voxelCoord
# 加载结节标注
anno_path = '/home/dataset/medical/luna16/CSVFILES/annotations.csv'
annos = readCSV(anno_path) # 共1186个结节标注
print(len(annos))
print(annos[0:3])
'''
1187
[['seriesuid', 'coordX', 'coordY', 'coordZ', 'diameter_mm'], ['1.3.6.1.4.1.14519.5.2.1.6279.6001.100225287222365663678666836860', '-128.6994211', '-175.3192718', '-298.3875064', '5.651470635'], ['1.3.6.1.4.1.14519.5.2.1.6279.6001.100225287222365663678666836860', '103.7836509', '-211.9251487', '-227.12125', '4.224708481']]
'''
# 获取一个结节标注
cand = annos[24]
print(cand)
'''
['1.3.6.1.4.1.14519.5.2.1.6279.6001.108197895896446896160048741492',
'-100.5679445',
'67.26051683',
'-231.816619',
'6.440878725']
'''
# 将世界坐标下肺结节标注转换为真实坐标系下的坐标标注
worldCoord = np.asarray([float(cand[1]),float(cand[2]),float(cand[3])])
voxelCoord = worldToVoxelCoord(worldCoord, numpyOrigin, numpySpacing)
print(voxelCoord)
'''
array([110.39259333, 346.62447366, 32.7733524 ])
'''5. 疑问
CT世界坐标系是怎样建立的，为什么原点不是（0，0，0）？
6. Reference
编辑于 2018-12-17 11:45