赞
踩
在来看看dicom.ImagePositionPatient中三个值代表的是哪个的那些轴注意看第一段代码中的x对应Origin[0],y对应Origin[1],z对应Origin[2].而Origin的获取是以下代码:Ps:这里面的一些道道给了我很大的困扰,这里在我脑袋清晰的时候将理解记录下来防止以后再次弄混.
先来看张图看看是怎么照CT:
所以先不看Z轴,只看X轴和Y轴这个面,看出X是列,而Y是行.这跟我们平时在图上画点很类似.图参考:https://www.cnblogs.com/h2zZhou/p/9072967.html
而dicom.ImagePositionPatient是直接从设备生成的dicom里读取的故而dicom.ImagePositionPatient里的三个值分别代表X轴,Y轴,Z轴的初始点.一定注意X是列,而Y是行.
来看看luna16标签数据里的xyz.在之前我的博文https://blog.csdn.net/qq_36401512/article/details/85163467用来将标签在二维上可视化.细节可以自己去看我们只看两段代码:
- x, y, z = int((nodules[idx, 0]-Origin[0])/SP[0]), int((nodules[idx, 1]-Origin[1])/SP[1]),
- int((nodules[idx, 2]-Origin[2])/SP[2])
- # 注意 y代表纵轴,x代表横轴
- data[max(0, y - radius):min(data.shape[0], y + radius),
- max(0, x - radius - pad):max(0, x - radius)] = 3000 # 竖线
- data[max(0, y - radius):min(data.shape[0], y + radius),
- min(data.shape[1], x + radius):min(data.shape[1], x + radius + pad)] = 3000 # 竖线
- data[max(0, y - radius - pad):max(0, y - radius),
- max(0, x - radius):min(data.shape[1], x + radius)] = 3000 # 横线
- data[min(data.shape[0], y + radius):min(data.shape[0], y + radius + pad),
- max(0, x - radius):min(data.shape[1], x + radius)] = 3000 # 横线
你可注意到x对应nodules[idx, 0],y对应nodules[idx,1],z对应nodules[idx,2].而后后面这段标出标签的方式就是赋一个很大值.但是有注释y代表纵轴,x代表横轴.在看上面那段代码第一行y相关的参数作为data赋值时的第一维度更加印证了上面的注释.
比如我们知道二维矩阵A,赋值取到点(x1,y1)代表A[x1,y1].而画图时想要标记出点(x1,y1)则需要如下代码:
- plt.ion()
- plt.scatter(y1, x1, marker='o', color='red', s=10, label='First')
- plt.imshow(A, cmap=plt.cm.bone)
如果还模棱两可,再来看看第一段代码中的x对应Origin[0],y对应Origin[1],z对应Origin[2].而Origin的获取是以下代码:
- filename='E:\\JLS\\dcm_data\\luna\\subset1\\1.3.6.1.4.1.14519.5.2.1.6279.6001.173106154739244262091404659845.mhd'
- itkimage = sitk.ReadImage(filename)#读取.mhd文件
- Origin=itkimage.GetOrigin()
而itkimage.GetOrigin()与这个CT序列的第一张dicom的ImagePositionPatient是一样的.这样就可以看出luna16的标签数据与dicom.ImagePositionPatient代表的轴是一样的,也是coordY代表列,coordX代表行.
PS:2019.11.6 博客的思路结论写了三遍,大体过程:是这样的,不反了反了是那样的,还是这样的.最后这次为了确保正确我还再次写了按照标签分别画出结果用luna16数据和目前在做项目的数据.具体代码如下:
- # coding=UTF-8
- import cv2
- import numpy as np
- from PIL import Image
- from scipy import misc
- import os
- import sys
- import cv2
- from skimage import measure, morphology, segmentation
- import matplotlib.pyplot as plt
-
- import SimpleITK as sitk
- import pandas
-
- 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 image
-
- def load_itk_image(filename):
- with open(filename) as f:
- contents = f.readlines()
- line = [k for k in contents if k.startswith('TransformMatrix')][0]
- transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
- transformM = np.round(transformM)
- if np.any(transformM != np.array([1, 0, 0, 0, 1, 0, 0, 0, 1])):
- isflip = True
- else:
- isflip = False
-
- itkimage = sitk.ReadImage(filename)
- numpyImage = sitk.GetArrayFromImage(itkimage)
-
- numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
- numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))
-
- return numpyImage, numpyOrigin, numpySpacing, isflip
-
-
- F_path='Tannotations.csv'
- L_path='annotations.csv'
- F='/FORNEWTEST/change/NewForTest'
- L='dcm_data/luna/subsets'
-
- annosF = np.array(pandas.read_csv(F_path))
- annosL = np.array(pandas.read_csv(L_path))
-
- nameF = annosF[0][0]
- nameL = annosL[55][0]
- coordXF,coordYF,coordZF=annosF[0][1],annosF[0][2],annosF[0][3]
- coordXL,coordYL,coordZL=annosL[55][1],annosL[55][2],annosL[55][3]
- print(nameF,nameL)
-
- for root, dirs, files in os.walk(F):
- for file in files:
- if (file[:-4] == nameF and file[-4:] == '.mhd'):
- sliceim, origin, spacing, isflip = load_itk_image(
- os.path.join(root, file))
- ZF=np.absolute(coordZF-origin[0])/spacing[0]
- YF=np.absolute(coordYF-origin[1])/spacing[1]
- XF = np.absolute(coordXF - origin[2]) / spacing[2]
- imgF=sliceim[ZF]
- plt.ion()
- plt.scatter(XF, YF, marker='o', color='red', s=10, label='First')
- plt.imshow(imgF, cmap=plt.cm.bone)
- plt.show()
-
- for root, dirs, files in os.walk(L):
- for file in files:
- if (file[:-4] == nameL and file[-4:] == '.mhd'):
- sliceim, origin, spacing, isflip = load_itk_image(os.path.join(root, file))
- ZL=int(np.absolute(coordZL-origin[0])/spacing[0])
- YL=np.absolute(coordYL-origin[1])/spacing[1]
- XL = np.absolute(coordXL - origin[2]) / spacing[2]
- imgL=sliceim[ZL]
- imgC=imgL[YL-20:YL+20,XL-20:XL+20]
- plt.ion()
- plt.scatter(XL, YL, marker='o', color='red', s=2, label='First')
- plt.imshow(imgL, cmap=plt.cm.bone)
- plt.show()
- plt.imshow(imgC, cmap=plt.cm.bone)
- plt.show()
现在思路终于十分清晰了.真是要眼见为实!
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。