大赛网址：天池医疗AI大赛[第一季]
AI的比赛是第一次参加，就当作一次实战练习了。由于比赛规则，这里也只能记录一些想法，不能将一些实际的实现方法以代码或者其他形式发布出来，请谅解。由于组委会规定，所以神经网络只能使用Caffe套件。
OpenCV @2017/04/19
官方教程第一步是进行图像分割，图像分割一般使用OpenCV来做，下面进行OpenCV的研究。
对于C#比较熟悉，所以暂时考虑在Win平台上进行研究。
在Windows上的安装
http://opencv.org/releases.html
选择Win Pack，安装其实就是一个解压缩的过程。
在这个路径中保存着所需要的可执行文件
你的解压路径\opencv\build\x64\vc14\bin
这个可执行文件需要添加到系统的Path里面去，以后其他程序可以直接调用。
EnguCV
EmguCV是一个C#用的OpenCV包装
(不知道是不是一定要安装OpenCV，可能这个东西可以脱离OpenCV直接用)
Install-Package EmguCV即可将这个包导入到项目中。
using System;
using Emgu.CV;
using Emgu.CV.CvEnum;
using Emgu.CV.Structure;
namespace OpenCV
{
class Program
{
static void Main(string[] args)
{
String win1 = "Test Window"; //The name of the window
CvInvoke.NamedWindow(win1); //Create the window using the specific name
Mat img = new Mat(200, 400, DepthType.Cv8U, 3); //Create a 3 channel image of 400x200
img.SetTo(new Bgr(255, 0, 0).MCvScalar); // set it to Blue color
//Draw "Hello, world." on the image using the specific font
CvInvoke.PutText(
img,
"Hello, world",
new System.Drawing.Point(10, 80),
FontFace.HersheyComplex,
1.0,
new Bgr(0, 255, 0).MCvScalar);
CvInvoke.Imshow(win1, img); //Show the image
CvInvoke.WaitKey(0); //Wait for the key pressing event
CvInvoke.DestroyWindow(win1); //Destroy the window if key is pressed
}
}
}
看一下是否能够运行，出现一个窗体如下所示。代码暂时不去理解了。
Emgu的HelloWorld
SimpleITK@2017/04/27
本次比赛的CT数据格式为MHD文件，由于采样用设备不同，每个文件的大小也不同。最小的在50M，最大的在250M左右。每个MHD（1K的小文件）和ZRAW文件为一套数据。整个比赛的文件大小约为50G，计算资源的多少是胜出的关键因素。
MHD的内容如下所示，保存着CT图片的一些摘要信息
ObjectType = Image
NDims = 3
BinaryData = True
BinaryDataByteOrderMSB = False
CompressedData = True
CompressedDataSize = 88265161
TransformMatrix = 1 0 0 0 1 0 0 0 1
Offset = -211.621 -358.121 -99.4
CenterOfRotation = 0 0 0
AnatomicalOrientation = RAI
ElementSpacing = 0.757813 0.757813 1
DimSize = 512 512 264
ElementType = MET_SHORT
ElementDataFile = LKDS-00012.zraw
注释版：
ObjectType = Image 目标类型是图像
NDims = 3 三维数据
BinaryData = True 二进制数据
BinaryDataByteOrderMSB = False 不是big endian（Unix格式），是windows linux格式二进制
CompressedData = True 压缩数据
CompressedDataSize = 91327856 压缩数据大小
TransformMatrix = 1 0 0 0 1 0 0 0 1 传输矩阵100,010,001 x,y,z
Offset = -229 -131 -506.4 原点坐标
CenterOfRotation = 0 0 0 中心旋转：无
AnatomicalOrientation = RAI 患者体位
ElementSpacing = 0.894531 0.894531 1 像素间隔x,y,z
DimSize = 512 512 365 数据大小，x，y, z
ElementType = MET_SHORT 数据类型，16bit整数
ElementDataFile = LKDS-00010.zraw 数据存储的文件名
读取mhd文件可以参考以下这篇文章：
[天池医疗AI大赛[第一季]：肺部结节智能诊断] simpleITK读取mhd的demo(4行,zraw和mhd放一起)
import SimpleITK as sitk
file_name = 'E:\MHD\LKDS-00012.mhd'
img = sitk.ReadImage(file_name)
img_array = sitk.GetArrayFromImage(img)
目录最好是纯英文的，不然可能出现无法读取的问题
下面根据这篇文章尝试将结节图像提取出来（这里只是进行一个假象的固定结节的保存），每个结节的采样文件大概是3M大小。
[天池医疗AI大赛[第一季]：肺部结节智能诊断] 参赛新人向导（一）数据读取和可视化--数据长什么样（含代码）
def process_image(file_name,output_path,nodule):
itk_img = sitk.ReadImage(file_name)
# load the data once
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)
# go through all nodes (why just the biggest?)
node_x = nodule.node_x
node_y = nodule.node_y
node_z = nodule.node_z
diam = nodule.diam
# just keep 3 slices
imgs = np.ndarray([3,height,width],dtype=np.float32)
masks = np.ndarray([3,height,width],dtype=np.uint8)
center = np.array([node_x, node_y, node_z]) # nodule center
v_center = SITKlib.worldToVoxel(center,origin,spacing) # nodule center in voxel space (still x,y,z ordering)
for i, i_z in enumerate(np.arange(int(v_center[2])-1, int(v_center[2])+2).clip(0, num_z-1)): # clip prevents going out of bounds in Z
mask = SITKlib.make_mask(center, diam, i_z*spacing[2]+origin[2],width, height, spacing, origin)
masks[i] = mask
imgs[i] = img_array[i_z]
np.save(os.path.join(output_path,"images.npy"),imgs)
np.save(os.path.join(output_path,"masks.npy"),masks)
SITKlib.show_img(imgs,masks)
接下来分析一下每部分的代码
读图
itk_img = sitk.ReadImage(file_name)
# load the data once
img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
num_z, height, width = img_array.shape #height X width constitute the transverse plane
使用SimpleITK的ReadImage函数可以获得CT图像。这个图像应该不可以直接进行可视化操作。
然后将这个图像放入一个数组中保存起来，注意这里的顺序，Z，Y，X。
数据准备
origin：表示CT图像最边缘的坐标
sapcing：真实世界和像素的比例关系
然后我们假定有一个结节，其位置和半径由 node_x，node_y，node_z，diam决定
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)
# go through all nodes (why just the biggest?)
node_x = nodule.node_x
node_y = nodule.node_y
node_z = nodule.node_z
diam = nodule.diam
真实世界坐标和CT像素坐标的转换
这里因为只保留和结节相关的3层，所以，z方向上固定为3.
然后比较核心的计算是这样的：
(center-origin)：结节的位置和图像边缘的位置做比较，就知道相差多少距离。
例如这里的Origin位置是：-211.621， -358.121， -99.4 （这个可以在mhd定义的地方看到）
结节的位置是：-100.56 ，67.26， -31.81
则相差的位置大约为 [ 111.061, 425.381, 67.59 ]
然后根据实际世界坐标和像素变换关系（ElementSpacing = 0.757813 ， 0.757813 ，1）
可以知道结节中心的位置所在的像素为：[ 147， 561， 68]。
这里请关注Z方向编号是68
# just keep 3 slices
imgs = np.ndarray([3,height,width],dtype=np.float32)
masks = np.ndarray([3,height,width],dtype=np.uint8)
center = np.array([node_x, node_y, node_z]) # nodule center
v_center = SITKlib.worldToVoxel(center,origin,spacing) # nodule center in voxel space (still x,y,z ordering)
def worldToVoxel(worldCoord, origin, spacing):
'''
结节世界坐标转成图像坐标
worldCoord:结节世界坐标
origin:边缘坐标
spacing：比例尺
'''
voxelCoord = np.absolute(worldCoord - origin)
voxelCoord = voxelCoord / spacing
voxelCoord = voxelCoord.astype(np.int32)
return voxelCoord
借用网络上的一张图说明问题：
原文：http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/03_Image_Details.html
坐标转换
保存图片
下面则将67，68，69，三个层面（每个层面是1个像素）的图像进行保存，整个层面的保存在"images.npy"里面。
结节部分的图片保存在"masks.npy"里面。
for i, i_z in enumerate(np.arange(int(v_center[2])-1, int(v_center[2])+2).clip(0, num_z-1)): # clip prevents going out of bounds in Z
mask = make_mask(center, diam, i_z*spacing[2]+origin[2],width, height, spacing, origin)
masks[i] = mask
imgs[i] = img_array[i_z]
np.save(os.path.join(output_path,"images.npy"),imgs)
np.save(os.path.join(output_path,"masks.npy"),masks)
可视化
彩色切片：原始图像
黑白切片：灰度图像
节点： 将节点区域使用白色展示在黑色背景中，作为节点切片的Mask层
节点切片：仅将结节表示在图像中，并且进行灰度处理（为什么非节点地方不是黑色的？？）
** imgs[i]*masks[i] ，mask里面只有 1 （白色，结节） 、0（黑色） 两种值 **
for i in range(len(imgs)):
print ("图片的第 %d 层" % i)
fig,ax = plt.subplots(2,2,figsize=[8,8])
ax[0,0].imshow(imgs[i])
ax[0,0].set_title(u'彩色切片')
ax[0,1].imshow(imgs[i],cmap='gray')
ax[0,1].set_title(u'黑白切片')
ax[1,0].imshow(masks[i],cmap='gray')
ax[1,0].set_title(u'节点')
ax[1,1].imshow(imgs[i]*masks[i],cmap='gray')
ax[1,1].set_title(u'节点切片')
plt.show()
print ('\n\n')
可视化
源代码（截至2017/04/28）
Solution.py
# coding: UTF-8
import os
import numpy as np
import SimpleITK as sitk
import matplotlib.pyplot as plt
import SITKlib
def process_image(file_name,output_path,nodule):
itk_img = sitk.ReadImage(file_name)
# load the data once
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)
# go through all nodes (why just the biggest?)
node_x = nodule.node_x
node_y = nodule.node_y
node_z = nodule.node_z
diam = nodule.diam
# just keep 3 slices
imgs = np.ndarray([3,height,width],dtype=np.float32)
masks = np.ndarray([3,height,width],dtype=np.uint8)
center = np.array([node_x, node_y, node_z]) # nodule center
v_center = SITKlib.worldToVoxel(center,origin,spacing) # nodule center in voxel space (still x,y,z ordering)
for i, i_z in enumerate(np.arange(int(v_center[2])-1, int(v_center[2])+2).clip(0, num_z-1)): # clip prevents going out of bounds in Z
mask = SITKlib.make_mask(center, diam, i_z*spacing[2]+origin[2],width, height, spacing, origin)
masks[i] = mask
imgs[i] = img_array[i_z]
np.save(os.path.join(output_path,"images.npy"),imgs)
np.save(os.path.join(output_path,"masks.npy"),masks)
SITKlib.show_img(imgs,masks)
def main():
#Global Setting
mhd_file_name = 'E:\\Huge CT Data\\train_subset00\\LKDS-00001.mhd'
output_path = "E:\\Huge CT Data"
nodule = SITKlib.Nodule();
nodule.node_x = -76.4498793983
nodule.node_y = -49.5405710363
nodule.node_z = 229.5
nodule.diam = 14.1804045239
process_image(mhd_file_name,output_path,nodule)
if __name__ == "__main__":
main()
SITKlib.py
# coding: UTF-8
import os
import numpy as np
import SimpleITK as sitk
import matplotlib.pyplot as plt
class Nodule:
def __init__(self):
self.node_x = 0
self.node_y = 0
self.node_z = 0
self.diam = 0
def worldToVoxel(worldCoord, origin, spacing):
'''
结节世界坐标转成图像坐标
worldCoord:结节世界坐标
origin:边缘坐标
spacing：比例尺
'''
voxelCoord = np.absolute(worldCoord - origin)
voxelCoord = voxelCoord / spacing
voxelCoord = voxelCoord.astype(np.int32)
return voxelCoord
def show_img(imgs,masks):
for i in range(len(imgs)):
print ("图片的第 %d 层" % i)
fig,ax = plt.subplots(2,2,figsize=[8,8])
ax[0,0].imshow(imgs[i])
ax[0,0].set_title(u'彩色切片')
ax[0,1].imshow(imgs[i],cmap='gray')
ax[0,1].set_title(u'黑白切片')
ax[1,0].imshow(masks[i],cmap='gray')
ax[1,0].set_title(u'节点')
ax[1,1].imshow(imgs[i]*masks[i],cmap='gray')
ax[1,1].set_title(u'节点切片')
plt.show()
print ('\n\n')
#raw_input("hit enter to cont : ")
#只显示结节
def make_mask(center,diam,z,width,height,spacing,origin):
'''
Center : 圆的中心 px -- list of coordinates x,y,z
diam : 圆的直径 px -- diameter
widthXheight : pixel dim of image
spacing = mm/px conversion rate np array x,y,z
origin = x,y,z mm np.array
z = z position of slice in world coordinates mm
'''
mask = np.zeros([height,width])
# 0's everywhere except nodule swapping x,y to match img
#convert to nodule space from world coordinates
# Defining the voxel range in which the nodule falls
v_center = (center-origin)/spacing
v_diam = int(diam/spacing[0]+5)
v_xmin = np.max([0,int(v_center[0]-v_diam)-5])
v_xmax = np.min([width-1,int(v_center[0]+v_diam)+5])
v_ymin = np.max([0,int(v_center[1]-v_diam)-5])
v_ymax = np.min([height-1,int(v_center[1]+v_diam)+5])
v_xrange = range(v_xmin,v_xmax+1)
v_yrange = range(v_ymin,v_ymax+1)
# Convert back to world coordinates for distance calculation
x_data = [x*spacing[0]+origin[0] for x in range(width)]
y_data = [x*spacing[1]+origin[1] for x in range(height)]
# Fill in 1 within sphere around nodule
for v_x in v_xrange:
for v_y in v_yrange:
p_x = spacing[0]*v_x + origin[0]
p_y = spacing[1]*v_y + origin[1]
if np.linalg.norm(center-np.array([p_x,p_y,z]))<=diam:
mask[int((p_y-origin[1])/spacing[1]),int((p_x-origin[0])/spacing[0])] = 1.0
return(mask)
def matrix2int16(matrix):
'''
matrix must be a numpy array NXN
Returns uint16 version
'''
m_min= np.min(matrix)
m_max= np.max(matrix)
matrix = matrix-m_min
return(np.array(np.rint( (matrix-m_min)/float(m_max-m_min) * 65535.0),dtype=np.uint16))
#
# Helper function to get rows in data frame associated
# with each file
def get_filename(file_list, case):
for f in file_list:
if case in f:
return(f)
#
# The locations of the nodes
def normalize(image, MIN_BOUND=-1000.0, MAX_BOUND=400.0):
"""数据标准化"""
image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
image[image > 1] = 1.
image[image < 0] = 0.
return image
#---数据标准化
def set_window_width(image, MIN_BOUND=-1000.0, MAX_BOUND=400.0):
"""设置窗宽"""
image[image > MAX_BOUND] = MAX_BOUND
image[image < MIN_BOUND] = MIN_BOUND
return image
#---设置窗宽
网友评论