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)