1、使用SimpleITK对齐图像
在看voxelmorph的代码，看到图像对齐部分，记录一下。
下面是从voxelmorph项目中截取的一段保存图像的函数。
函数输入分别是：配准后的图像、固定图像、要将配准图像保存的名字。
将图像对齐的操作需要将对齐的图像的原点、方向、间距设置成与 被对齐的图像一致。

1 def save_image(img, ref_img, name):
2     img = sitk.GetImageFromArray(img[0, 0, ...].cpu().detach().numpy())
3     img.SetOrigin(ref_img.GetOrigin())
4     img.SetDirection(ref_img.GetDirection())
5     img.SetSpacing(ref_img.GetSpacing())
6     sitk.WriteImage(img, os.path.join(args.result_dir, name))

这里只是做了属性的修改，并没有真正的对齐。

2、重采样并对齐图像

 1 def align_seg_with_raw_nrrd(dcm, seg):
 2     # Just for labelmap .... because of nearestNeighour interpolator
 3     resampler = sitk.ResampleImageFilter()
 4     resampler.SetReferenceImage(dcm)
 5     resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
 6     resampler.SetInterpolator(sitk.sitkNearestNeighbor)
 7     seg_new = resampler.Execute(seg)
 8     return seg_new
 9 
10 if __name__ == '__main__':
11     raw_dcm_file = r"path to your nrrd raw data"
12     seg_dcm_file= r"path to your seg file"
13     dcm = sitk.ReadImage(raw_dcm_file)
14     seg = sitk.ReadImage(seg_dcm_file)
15     print(raw_dcm_file, seg_new_file)
16     seg_new = align_seg_with_raw_nrrd(dcm, seg)

我要实现的目的是：将原始stl切下与原始数据重合的部分，并且用一个与原始数据一样的尺寸存放，并且对应的空间坐标一致（原点、spacing、direction）。

如下输入原始数据：

 原始的stl数据，绿色部分：

 对齐后的效果如下：
stl对齐后的volume信息如下

 对齐后的渲染效果如下

3、中间走过的弯路

多设置参数，得到错误的结果。

4、尺寸不一致时改变方向并对齐

  1 # -*- coding : UTF-8 -*-
  2 # @file   : resample_seg.py
  3 # @Time   : 2021/10/14 0014 19:58
  4 # @Author : wmz
  5 
  6 import os
  7 import SimpleITK as sitk
  8 from glob import glob
  9 import numpy as np
 10 
 11 
 12 def chnm_stdraw2seglabel(stdrawfile):
 13     seglabelfile = stdrawfile
 14     seglabelfile = seglabelfile.replace("_knee", "")
 15     if right_flag:
 16         if femur_flag:
 17             seglabelfile = 'Seg_' + seglabelfile[:-5] + '_femur_right.nrrd'
 18         else:
 19             seglabelfile = 'Seg_' + seglabelfile[:-5] + '_hip_right.nrrd'
 20     else:
 21         if femur_flag:
 22             seglabelfile = 'Seg_' + seglabelfile[:-5] + '_femur_left.nrrd'
 23         else:
 24             seglabelfile = 'Seg_' + seglabelfile[:-5] + '_hip_left.nrrd'
 25     return seglabelfile
 26 
 27 
 28 def align_seg_with_raw_dcm(dcm, seg):
 29     # Just for labelmap .... because of nearestNeighour interpolator
 30     # 读取文件的size和spacing信息
 31     outsize = [0, 0, 0]
 32     inputsize = seg.GetSize()
 33     inputspacing = seg.GetSpacing()
 34     input_origin = seg.GetOrigin()
 35     dcm_size = dcm.GetSize()
 36     dcm_spacing = dcm.GetSpacing()
 37     dcm_origin = dcm.GetOrigin()
 38     direction = dcm.GetDirection()
 39     seg_direction = seg.GetDirection()
 40 
 41     transform = sitk.Transform()
 42     transform.SetIdentity()
 43     # 计算改变spacing后的size，用物理尺寸/体素的大小
 44     outsize[0] = round(inputsize[0] * inputspacing[0] / dcm_spacing[0])
 45     outsize[1] = round(inputsize[1] * inputspacing[1] / dcm_spacing[1])
 46     # outsize[2] = round(inputsize[2] * inputspacing[2] / dcm_spacing[2]) # 原大尺寸，弃用
 47     outsize[2] = round((dcm_size[2] + dcm_origin[2] - input_origin[2])/dcm_spacing[2])
 48 
 49     resampler = sitk.ResampleImageFilter()
 50     resampler.SetReferenceImage(dcm)
 51     resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
 52     resampler.SetInterpolator(sitk.sitkNearestNeighbor)
 53     resampler.SetOutputOrigin(seg.GetOrigin())
 54     resampler.SetOutputSpacing(dcm.GetSpacing())
 55     resampler.SetOutputDirection(seg.GetDirection())
 56     resampler.SetSize(outsize)
 57     seg_new = resampler.Execute(seg)
 58     return seg_new
 59 
 60 
 61 def align_label2stdraw(stdrawfile, labelfile):
 62     std_img = sitk.ReadImage(stdrawfile)
 63     std_img_size = std_img.GetSize()
 64     std_img_direction = std_img.GetDirection()
 65     label_img = sitk.ReadImage(labelfile)
 66     new_seg = align_seg_with_raw_dcm(std_img, label_img)
 67     new_seg_origin = new_seg.GetOrigin()
 68     new_seg_direct = new_seg.GetDirection()
 69     new_seg_spacing = new_seg.GetSpacing()
 70     new_seg_size = new_seg.GetSize()
 71     new_seg_data = sitk.GetArrayFromImage(new_seg)
 72     new_seg_data = np.flip(new_seg_data, 1)
 73     new_seg_data = np.flip(new_seg_data, 2)
 74     stdlabel_img = sitk.GetImageFromArray(new_seg_data)
 75     flip_seg_origin = [new_seg_origin[0] + new_seg_direct[0]*new_seg_spacing[0]*new_seg_size[0],
 76                        new_seg_origin[1] + new_seg_direct[4]*new_seg_spacing[1]*new_seg_size[1], new_seg_origin[2]]
 77     stdlabel_img.SetOrigin(flip_seg_origin)
 78     stdlabel_img.SetSpacing(new_seg_spacing)
 79     stdlabel_img.SetDirection(std_img_direction)
 80 
 81     name = labelfile.split("\\")[-1]
 82     sitk.WriteImage(stdlabel_img, os.path.join(dst_dir, name))
 83 
 84 
 85 if __name__ == '__main__':
 86     std_full_dir = r'F:\dataset\align_data\std_raw'
 87     label_dir = r"F:\dataset\align_data\Seg_volume"
 88     dst_dir = r"F:\dataset\align_data\output"
 89     right_flag = True
 90     femur_flag = False
 91 
 92     files = glob(std_full_dir + '/*.nrrd')
 93     for indx, file in enumerate(files):
 94         print("processing", indx+1, "of", len(files))
 95         print("processing file:", file)
 96         # file = r"D:\Work\Data\Patient_Data\train_data\knee_femur\left\NanFang_03021_knee.nrrd"
 97         std_raw_file = file.split('\\')[-1]
 98         seg_label_file = chnm_stdraw2seglabel(std_raw_file)
 99         org_label = os.path.join(label_dir, seg_label_file)
100         if not os.path.exists(org_label):
101             print(org_label, "file not exist !")
102             continue
103         align_label2stdraw(file, org_label)

关键部分是下面这几句：

1      new_seg_data = np.flip(new_seg_data, 1)
2     new_seg_data = np.flip(new_seg_data, 2)
3     stdlabel_img = sitk.GetImageFromArray(new_seg_data)
4     flip_seg_origin = [new_seg_origin[0] + new_seg_direct[0]*new_seg_spacing[0]*new_seg_size[0],
5                        new_seg_origin[1] + new_seg_direct[4]*new_seg_spacing[1]*new_seg_size[1], new_seg_origin[2]]
6     stdlabel_img.SetOrigin(flip_seg_origin)
7     stdlabel_img.SetSpacing(new_seg_spacing)
8     stdlabel_img.SetDirection(std_img_direction)