大賽網址：天池醫(yī)療AI大賽[第一季]
AI的比賽是第一次參加，就當作一次實戰(zhàn)練習了。由于比賽規(guī)則，這里也只能記錄一些想法，不能將一些實際的實現(xiàn)方法以代碼或者其他形式發(fā)布出來，請諒解。由于組委會規(guī)定，所以神經網絡只能使用Caffe套件。
OpenCV @2017/04/19
官方教程第一步是進行圖像分割，圖像分割一般使用OpenCV來做，下面進行OpenCV的研究。
對于C#比較熟悉，所以暫時考慮在Win平臺上進行研究。
在Windows上的安裝
http://opencv.org/releases.html
選擇Win Pack，安裝其實就是一個解壓縮的過程。
在這個路徑中保存著所需要的可執(zhí)行文件
你的解壓路徑\opencv\build\x64\vc14\bin
這個可執(zhí)行文件需要添加到系統(tǒng)的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
}
}
}
看一下是否能夠運行，出現(xiàn)一個窗體如下所示。代碼暫時不去理解了。
Emgu的HelloWorld
SimpleITK@2017/04/27
本次比賽的CT數(shù)據格式為MHD文件，由于采樣用設備不同，每個文件的大小也不同。最小的在50M，最大的在250M左右。每個MHD（1K的小文件）和ZRAW文件為一套數(shù)據。整個比賽的文件大小約為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 三維數(shù)據
BinaryData = True 二進制數(shù)據
BinaryDataByteOrderMSB = False 不是big endian（Unix格式），是windows linux格式二進制
CompressedData = True 壓縮數(shù)據
CompressedDataSize = 91327856 壓縮數(shù)據大小
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 數(shù)據大小，x，y, z
ElementType = MET_SHORT 數(shù)據類型，16bit整數(shù)
ElementDataFile = LKDS-00010.zraw 數(shù)據存儲的文件名
讀取mhd文件可以參考以下這篇文章：
[天池醫(yī)療AI大賽[第一季]：肺部結節(jié)智能診斷] 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)
目錄最好是純英文的，不然可能出現(xiàn)無法讀取的問題
下面根據這篇文章嘗試將結節(jié)圖像提取出來（這里只是進行一個假象的固定結節(jié)的保存），每個結節(jié)的采樣文件大概是3M大小。
[天池醫(yī)療AI大賽[第一季]：肺部結節(jié)智能診斷] 參賽新人向導（一）數(shù)據讀取和可視化--數(shù)據長什么樣（含代碼）
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函數(shù)可以獲得CT圖像。這個圖像應該不可以直接進行可視化操作。
然后將這個圖像放入一個數(shù)組中保存起來，注意這里的順序，Z，Y，X。
數(shù)據準備
origin：表示CT圖像最邊緣的坐標
sapcing：真實世界和像素的比例關系
然后我們假定有一個結節(jié)，其位置和半徑由 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像素坐標的轉換
這里因為只保留和結節(jié)相關的3層，所以，z方向上固定為3.
然后比較核心的計算是這樣的：
(center-origin)：結節(jié)的位置和圖像邊緣的位置做比較，就知道相差多少距離。
例如這里的Origin位置是：-211.621， -358.121， -99.4 （這個可以在mhd定義的地方看到）
結節(jié)的位置是：-100.56 ，67.26， -31.81
則相差的位置大約為 [ 111.061, 425.381, 67.59 ]
然后根據實際世界坐標和像素變換關系（ElementSpacing = 0.757813 ， 0.757813 ，1）
可以知道結節(jié)中心的位置所在的像素為：[ 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):
'''
結節(jié)世界坐標轉成圖像坐標
worldCoord:結節(jié)世界坐標
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"里面。
結節(jié)部分的圖片保存在"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)
可視化
彩色切片：原始圖像
黑白切片：灰度圖像
節(jié)點： 將節(jié)點區(qū)域使用白色展示在黑色背景中，作為節(jié)點切片的Mask層
節(jié)點切片：僅將結節(jié)表示在圖像中，并且進行灰度處理（為什么非節(jié)點地方不是黑色的？？）
** imgs[i]*masks[i] ，mask里面只有 1 （白色，結節(jié)） 、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'節(jié)點')
ax[1,1].imshow(imgs[i]*masks[i],cmap='gray')
ax[1,1].set_title(u'節(jié)點切片')
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):
'''
結節(jié)世界坐標轉成圖像坐標
worldCoord:結節(jié)世界坐標
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'節(jié)點')
ax[1,1].imshow(imgs[i]*masks[i],cmap='gray')
ax[1,1].set_title(u'節(jié)點切片')
plt.show()
print ('\n\n')
#raw_input("hit enter to cont : ")
#只顯示結節(jié)
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):
"""數(shù)據標準化"""
image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
image[image > 1] = 1.
image[image < 0] = 0.
return image
#---數(shù)據標準化
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
#---設置窗寬
