医疗影像处理入门篇读入*.mhd数据&解析ElementSpacing参数

关于Dicom数据

使用pydicom包即可，网上有很多写的比较详细的，我就不搬运了，不了解的还请自行百度

读取*.mhd数据

• 有一些数据集提供的是mhd格式的数据，还有一些是*.nii.gz的数据，这些格式的数据可以使用SimpleITK包来读取，首先安装SimpleITK。

• mhd文件夹内会有一个同mhd同名的.raw文件，这个就是我们的数据文件，加载的时候只需要提供 *.mhd的文件路径即可，可以使用记事本打开查看：

  ObjectType = Image
  NDims = 3
  BinaryData = True
  BinaryDataByteOrderMSB = False
  CompressedData = True
  CompressedDataSize = 349063905
  TransformMatrix = 1 0 0 0 1 0 0 0 1
  Offset = -90.8574 -213.357 590.75
  CenterOfRotation = 0 0 0
  ElementSpacing = 0.285156 0.285156 0.75
  DimSize = 512 512 1262
  AnatomicalOrientation = ???
  ElementType = MET_SHORT
  ElementDataFile = Leg.raw
  

  其中比较重要的是TransformMatrix 、Offset 和ElementSpacing ，
  其中：
  TransformMatrix ：代表x y z的方向，用于判定CT的方向
  Offset ：CT的原点坐标位置，一般在做器官分割的时候这个参数用不上，但是后期工程落地的时候是必须考虑的！！！
  ElementSpacing ：是层厚，单位：mm/pixel，划重点这个参数不论是前期的算法建模还是后期的工程落地都非常重要。

• 读取代码

  def read_mhd(mhd_dir):
      import SimpleITK as sitk
      mhd_file = os.path.join(mhd_dir, 'leg.mhd')
      itkimage = sitk.ReadImage(mhd_file)
      ct_value= sitk.GetArrayFromImage(itkimage)  # 这里一定要注意，得到的是[z,y,x]格式    
      direction = itkimage.GetDirection()  # mhd文件中的TransformMatrix 
      origin = np.array(itkimage.GetOrigin())
      spacing = np.array(itkimage.GetSpacing())  # 文件中的ElementSpacing 
  

  使用sitk.GetArrayFromImage()就可以得到我们需要的HU值，然后就可以通过设置窗宽/窗位来显示我们感兴趣的部位了。唯一需要注意的是得到的数据格式是[z,y,x]的顺序，切记！！

• 显示图像
  刚才使用sitk.GetArrayFromImage()得到的是HU值，我们首先设置一下窗宽/窗位

  def set_window(cls, img_data, win_width, win_center):
      img_temp = img_data
      min = (2 * win_center - win_width) / 2.0 + 0.5
      max = (2 * win_center + win_width) / 2.0 + 0.5
      dFactor = 255.0 / (max - min)
  
      img_temp = ((img_temp - min) * dFactor).astype(np.int)
      img_temp = np.clip(img_temp, 0, 255)
      return img_temp
  
  

  ct_value_new = set_window(ct_value, window_center[0], window_center[1])
  # k取DimSize 中的值
  tar_img = ct_value_new[k, :, :] # 水平面
  cor_img = ct_value_new[:, k, :] # 冠状面
  sag_img = ct_value_new[:, :, k] # 矢状面
  

  这样就可以得到三个面的二维视图了。比如想得到全部的冠状面的视图，直接写个循环，其余视图以此类推。

  for k in range(ct_value.shape[1]):
      cor_gray_img = ct_value_new[:, k, :].astype(np.uint8)
      cv2.imwrite('你的路径', cor_gray_img )
  

  目前一切看起来好像都没什么问题，但是有些数据保存后会看着很别扭，有些看着还挺正常。比如
  
  猛一看，压根不知道这是什么鬼！这其实是右下肢的冠状面，天呐不说还以为是个靴子呢（/捂脸）
  正常的图像应该长这样：
  
  这样再看是不是正常了。为什么会出现这种情况？就是因为我们上面的代码都没有考虑spacing这个参数！！

• 使用ElementSpacing

  1. 前面介绍的时候说过这个值的单位是mm/pixel，所以比如你的水平面的宽度是512，这个值是0.285，那么就可以得到实际的宽度512*0.285 = 145.92mm。上面的代码忽略了一个事实，在得到的ct_value矩阵中，并不是同像素等价的，假如ElementSpacing 的值是[0.285156, 0.285156, 0.75]，则水平面上没问题，一旦到了冠状面或矢状面[H, W]，以上面的mhd文件为例，冠状面的尺寸是[1262, 512]，H方向上每个像素是0.75mm，而W方向上却是每个像素0.285156mm，单位不一样，强行组成了一张图片，自然就变的扭曲了(强扭的瓜还是不甜啊)。这里有点绕，如果是初次接触这些东西，需要多理解一下。

  2. 其实仔细对比一下变形的图片和正常的图片也能看出来，变形的图片明显在短边的方向被拉伸了。

  3. 正确的打开方式
     3.1 为了后期的工程落地不再折腾，建议前期预处理的时候统一间隔，比如全部统一成[1.0, 1.0, 1.0]，后期工程落地的时候以实际的物理坐标作为同软件组的对接标准即可。
     3.2 实现方式

     def resample(image, spacing, new_spacing=[1,1,1]):
         import scipy
         
         raw_shape = image.shape
         resize_factor = spacing / new_spacing
         new_real_shape = image.shape * resize_factor
         new_shape = np.round(new_real_shape).astype(np.int)   # 返回浮点数x的四舍五入值。
         real_resize_factor = new_shape / image.shape
         new_spacing = spacing / real_resize_factor
         image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')  # 使用所请求顺序的样条插值来缩放数组。
     image = resample(ct_value, spacing[::-1])  # spacing[::-1]是因为ct_value的顺序是[z,y,x]，而读入的spacing是[x,y,z]
     

     这样差值以后再去显示任何二维视图就全部正常了，但是不论是检测还是分割，网络的输出都是像素坐标，工程落地的时候切记需要根据需要变化到物理坐标哦！

最后放一张itk-snap读取显示的冠状面结果：