当前位置:   article > 正文

luna16标签数据里的xyz,以及CT的dicom.ImagePositionPatient里的三个值分别代表哪些轴的初始点

imagepositionpatient

在来看看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用来将标签在二维上可视化.细节可以自己去看我们只看两段代码:

  1. x, y, z = int((nodules[idx, 0]-Origin[0])/SP[0]), int((nodules[idx, 1]-Origin[1])/SP[1]),
  2. int((nodules[idx, 2]-Origin[2])/SP[2])
  1. # 注意 y代表纵轴,x代表横轴
  2. data[max(0, y - radius):min(data.shape[0], y + radius),
  3. max(0, x - radius - pad):max(0, x - radius)] = 3000  # 竖线
  4. data[max(0, y - radius):min(data.shape[0], y + radius),
  5. min(data.shape[1], x + radius):min(data.shape[1], x + radius + pad)] = 3000  # 竖线
  6. data[max(0, y - radius - pad):max(0, y - radius),
  7. max(0, x - radius):min(data.shape[1], x + radius)] = 3000  # 横线
  8. data[min(data.shape[0], y + radius):min(data.shape[0], y + radius + pad),
  9. 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)则需要如下代码:

  1. plt.ion()
  2. plt.scatter(y1, x1, marker='o', color='red', s=10, label='First')
  3. plt.imshow(A, cmap=plt.cm.bone)

如果还模棱两可,再来看看第一段代码中的x对应Origin[0],y对应Origin[1],z对应Origin[2].而Origin的获取是以下代码:

  1. filename='E:\\JLS\\dcm_data\\luna\\subset1\\1.3.6.1.4.1.14519.5.2.1.6279.6001.173106154739244262091404659845.mhd'
  2. itkimage = sitk.ReadImage(filename)#读取.mhd文件
  3. Origin=itkimage.GetOrigin()

而itkimage.GetOrigin()与这个CT序列的第一张dicom的ImagePositionPatient是一样的.这样就可以看出luna16的标签数据与dicom.ImagePositionPatient代表的轴是一样的,也是coordY代表列,coordX代表行.

 

 

PS:2019.11.6 博客的思路结论写了三遍,大体过程:是这样的,不反了反了是那样的,还是这样的.最后这次为了确保正确我还再次写了按照标签分别画出结果用luna16数据和目前在做项目的数据.具体代码如下:

  1. # coding=UTF-8
  2. import cv2
  3. import numpy as np
  4. from PIL import Image
  5. from scipy import misc
  6. import os
  7. import sys
  8. import cv2
  9. from skimage import measure, morphology, segmentation
  10. import matplotlib.pyplot as plt
  11. import SimpleITK as sitk
  12. import pandas
  13. MIN_BOUND = -1000.0
  14. MAX_BOUND = 400.0
  15. def normalize(image):
  16. image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
  17. image[image > 1] = 1.
  18. image[image < 0] = 0.
  19. return image
  20. def load_itk_image(filename):
  21. with open(filename) as f:
  22. contents = f.readlines()
  23. line = [k for k in contents if k.startswith('TransformMatrix')][0]
  24. transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
  25. transformM = np.round(transformM)
  26. if np.any(transformM != np.array([1, 0, 0, 0, 1, 0, 0, 0, 1])):
  27. isflip = True
  28. else:
  29. isflip = False
  30. itkimage = sitk.ReadImage(filename)
  31. numpyImage = sitk.GetArrayFromImage(itkimage)
  32. numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
  33. numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))
  34. return numpyImage, numpyOrigin, numpySpacing, isflip
  35. F_path='Tannotations.csv'
  36. L_path='annotations.csv'
  37. F='/FORNEWTEST/change/NewForTest'
  38. L='dcm_data/luna/subsets'
  39. annosF = np.array(pandas.read_csv(F_path))
  40. annosL = np.array(pandas.read_csv(L_path))
  41. nameF = annosF[0][0]
  42. nameL = annosL[55][0]
  43. coordXF,coordYF,coordZF=annosF[0][1],annosF[0][2],annosF[0][3]
  44. coordXL,coordYL,coordZL=annosL[55][1],annosL[55][2],annosL[55][3]
  45. print(nameF,nameL)
  46. for root, dirs, files in os.walk(F):
  47. for file in files:
  48. if (file[:-4] == nameF and file[-4:] == '.mhd'):
  49. sliceim, origin, spacing, isflip = load_itk_image(
  50. os.path.join(root, file))
  51. ZF=np.absolute(coordZF-origin[0])/spacing[0]
  52. YF=np.absolute(coordYF-origin[1])/spacing[1]
  53. XF = np.absolute(coordXF - origin[2]) / spacing[2]
  54. imgF=sliceim[ZF]
  55. plt.ion()
  56. plt.scatter(XF, YF, marker='o', color='red', s=10, label='First')
  57. plt.imshow(imgF, cmap=plt.cm.bone)
  58. plt.show()
  59. for root, dirs, files in os.walk(L):
  60. for file in files:
  61. if (file[:-4] == nameL and file[-4:] == '.mhd'):
  62. sliceim, origin, spacing, isflip = load_itk_image(os.path.join(root, file))
  63. ZL=int(np.absolute(coordZL-origin[0])/spacing[0])
  64. YL=np.absolute(coordYL-origin[1])/spacing[1]
  65. XL = np.absolute(coordXL - origin[2]) / spacing[2]
  66. imgL=sliceim[ZL]
  67. imgC=imgL[YL-20:YL+20,XL-20:XL+20]
  68. plt.ion()
  69. plt.scatter(XL, YL, marker='o', color='red', s=2, label='First')
  70. plt.imshow(imgL, cmap=plt.cm.bone)
  71. plt.show()
  72. plt.imshow(imgC, cmap=plt.cm.bone)
  73. plt.show()

现在思路终于十分清晰了.真是要眼见为实!

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/花生_TL007/article/detail/165570
推荐阅读
  

闽ICP备14008679号