SimpleITK图像对齐
1、使用SimpleITK对齐图像
在看voxelmorph的代码,看到图像对齐部分,记录一下。
下面是从voxelmorph项目中截取的一段保存图像的函数。
函数输入分别是:配准后的图像、固定图像、要将配准图像保存的名字。
将图像对齐的操作需要将对齐的图像的原点、方向、间距设置成与 被对齐的图像一致。
def save_image(img, ref_img, name):
img = sitk.GetImageFromArray(img[0, 0, ...].cpu().detach().numpy())
img.SetOrigin(ref_img.GetOrigin())
img.SetDirection(ref_img.GetDirection())
img.SetSpacing(ref_img.GetSpacing())
sitk.WriteImage(img, os.path.join(args.result_dir, name))
这里只是做了属性的修改,并没有真正的对齐。
2、重采样并对齐图像
def align_seg_with_raw_nrrd(dcm, seg):
# Just for labelmap .... because of nearestNeighour interpolator
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(dcm)
resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
seg_new = resampler.Execute(seg)
return seg_new
if __name__ == '__main__':
raw_dcm_file = r"path to your nrrd raw data"
seg_dcm_file= r"path to your seg file"
dcm = sitk.ReadImage(raw_dcm_file)
seg = sitk.ReadImage(seg_dcm_file)
print(raw_dcm_file, seg_new_file)
seg_new = align_seg_with_raw_nrrd(dcm, seg)
我要实现的目的是:将原始stl切下与原始数据重合的部分,并且用一个与原始数据一样的尺寸存放,并且对应的空间坐标一致(原点、spacing、direction)。
如下输入原始数据:
原始的stl数据,绿色部分:
对齐后的效果如下:
stl对齐后的volume信息如下
对齐后的渲染效果如下
3、中间走过的弯路
多设置参数,得到错误的结果。
4、尺寸不一致时改变方向并对齐
# -*- coding : UTF-8 -*-
# @file : resample_seg.py
# @Time : 2021/10/14 0014 19:58
# @Author : wmz
import os
import SimpleITK as sitk
from glob import glob
import numpy as np
def chnm_stdraw2seglabel(stdrawfile):
seglabelfile = stdrawfile
seglabelfile = seglabelfile.replace("_knee", "")
if right_flag:
if femur_flag:
seglabelfile = 'Seg_' + seglabelfile[:-5] + '_femur_right.nrrd'
else:
seglabelfile = 'Seg_' + seglabelfile[:-5] + '_hip_right.nrrd'
else:
if femur_flag:
seglabelfile = 'Seg_' + seglabelfile[:-5] + '_femur_left.nrrd'
else:
seglabelfile = 'Seg_' + seglabelfile[:-5] + '_hip_left.nrrd'
return seglabelfile
def align_seg_with_raw_dcm(dcm, seg):
# Just for labelmap .... because of nearestNeighour interpolator
# 读取文件的size和spacing信息
outsize = [0, 0, 0]
inputsize = seg.GetSize()
inputspacing = seg.GetSpacing()
input_origin = seg.GetOrigin()
dcm_size = dcm.GetSize()
dcm_spacing = dcm.GetSpacing()
dcm_origin = dcm.GetOrigin()
direction = dcm.GetDirection()
seg_direction = seg.GetDirection()
transform = sitk.Transform()
transform.SetIdentity()
# 计算改变spacing后的size,用物理尺寸/体素的大小
outsize[0] = round(inputsize[0] * inputspacing[0] / dcm_spacing[0])
outsize[1] = round(inputsize[1] * inputspacing[1] / dcm_spacing[1])
# outsize[2] = round(inputsize[2] * inputspacing[2] / dcm_spacing[2]) # 原大尺寸,弃用
outsize[2] = round((dcm_size[2] + dcm_origin[2] - input_origin[2])/dcm_spacing[2])
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(dcm)
resampler.SetTransform(sitk.Transform(3, sitk.sitkIdentity))
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
resampler.SetOutputOrigin(seg.GetOrigin())
resampler.SetOutputSpacing(dcm.GetSpacing())
resampler.SetOutputDirection(seg.GetDirection())
resampler.SetSize(outsize)
seg_new = resampler.Execute(seg)
return seg_new
def align_label2stdraw(stdrawfile, labelfile):
std_img = sitk.ReadImage(stdrawfile)
std_img_size = std_img.GetSize()
std_img_direction = std_img.GetDirection()
label_img = sitk.ReadImage(labelfile)
new_seg = align_seg_with_raw_dcm(std_img, label_img)
new_seg_origin = new_seg.GetOrigin()
new_seg_direct = new_seg.GetDirection()
new_seg_spacing = new_seg.GetSpacing()
new_seg_size = new_seg.GetSize()
new_seg_data = sitk.GetArrayFromImage(new_seg)
new_seg_data = np.flip(new_seg_data, 1)
new_seg_data = np.flip(new_seg_data, 2)
stdlabel_img = sitk.GetImageFromArray(new_seg_data)
flip_seg_origin = [new_seg_origin[0] + new_seg_direct[0]*new_seg_spacing[0]*new_seg_size[0],
new_seg_origin[1] + new_seg_direct[4]*new_seg_spacing[1]*new_seg_size[1], new_seg_origin[2]]
stdlabel_img.SetOrigin(flip_seg_origin)
stdlabel_img.SetSpacing(new_seg_spacing)
stdlabel_img.SetDirection(std_img_direction)
name = labelfile.split("\\")[-1]
sitk.WriteImage(stdlabel_img, os.path.join(dst_dir, name))
if __name__ == '__main__':
std_full_dir = r'F:\dataset\align_data\std_raw'
label_dir = r"F:\dataset\align_data\Seg_volume"
dst_dir = r"F:\dataset\align_data\output"
right_flag = True
femur_flag = False
files = glob(std_full_dir + '/*.nrrd')
for indx, file in enumerate(files):
print("processing", indx+1, "of", len(files))
print("processing file:", file)
# file = r"D:\Work\Data\Patient_Data\train_data\knee_femur\left\NanFang_03021_knee.nrrd"
std_raw_file = file.split('\\')[-1]
seg_label_file = chnm_stdraw2seglabel(std_raw_file)
org_label = os.path.join(label_dir, seg_label_file)
if not os.path.exists(org_label):
print(org_label, "file not exist !")
continue
align_label2stdraw(file, org_label)
关键部分是下面这几句:
new_seg_data = np.flip(new_seg_data, 1)
new_seg_data = np.flip(new_seg_data, 2)
stdlabel_img = sitk.GetImageFromArray(new_seg_data)
flip_seg_origin = [new_seg_origin[0] + new_seg_direct[0]*new_seg_spacing[0]*new_seg_size[0],
new_seg_origin[1] + new_seg_direct[4]*new_seg_spacing[1]*new_seg_size[1], new_seg_origin[2]]
stdlabel_img.SetOrigin(flip_seg_origin)
stdlabel_img.SetSpacing(new_seg_spacing)
stdlabel_img.SetDirection(std_img_direction)