医学图像重采样resample记录(220825补充matlab resample)-20201102

医学影像重采样
利用simpleITK:

import SimpleITK as sitk
"""
resample
"""

def resampleVolume(outspacing,vol):
    """
    将体数据重采样的指定的spacing大小\n
    paras：
    outpacing：指定的spacing，例如[1,1,1]
    vol：sitk读取的image信息，这里是体数据\n
    return：重采样后的数据
    """
    outsize = [0,0,0]
    inputspacing = 0
    inputsize = 0
    inputorigin = [0,0,0]
    inputdir = [0,0,0]

    #读取文件的size和spacing信息
    
    inputsize = vol.GetSize()
    inputspacing = vol.GetSpacing()

    transform = sitk.Transform()
    transform.SetIdentity()
    #计算改变spacing后的size，用物理尺寸/体素的大小
    outsize[0] = int(inputsize[0]*inputspacing[0]/outspacing[0] + 0.5)
    outsize[1] = int(inputsize[1]*inputspacing[1]/outspacing[1] + 0.5)
    outsize[2] = int(inputsize[2]*inputspacing[2]/outspacing[2] + 0.5)

    #设定重采样的一些参数
    resampler = sitk.ResampleImageFilter()
    resampler.SetTransform(transform)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetOutputOrigin(vol.GetOrigin())
    resampler.SetOutputSpacing(outspacing)
    resampler.SetOutputDirection(vol.GetDirection())
    resampler.SetSize(outsize)
    newvol = resampler.Execute(vol)
    return newvol

    

def main():
    #读文件
    vol = sitk.Image(sitk.ReadImage("input.mha"))

    #重采样
    newvol = resampleVolume([1,1,1],vol)

    #写文件
    wriiter = sitk.ImageFileWriter()
    wriiter.SetFileName("output.mha")
    wriiter.Execute(newvol)

利用scipy插值

def resample(image, scan, new_spacing=[1,1,1]): # scan是load_scan函数返回的结果
    # Determine current pixel spacing
    spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))
    spacing = np.array(list(spacing))
    resize_factor = spacing / new_spacing
    new_real_shape = image.shape * resize_factor
    new_shape = np.round(new_real_shape)   #返回浮点数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') #使用所请求顺序的样条插值来缩放数组。
    return image, new_spacing
# 现在重新取样病人的像素，将其映射到一个同构分辨率 1mm x1mm x1mm。
pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])

220825补充Matlab差值：

imgs_single = double(imgs);
d = size(imgs_single);
scaleCoeff_z=thickness;
Xin = 1:1:d(1);
Yin = 1:1:d(2);
Zin = 1:1/scaleCoeff_z:d(3);
[Xout, Yout, Zout] = ndgrid(Xin, Yin, Zin);
rescaled_imgs = interpn(imgs_single, Xout, Yout, Zout);

以及

sz = size(imgs);
xg = 1:sz(1);
yg = 1:sz(2);
zg = 1:sz(3);
F = griddedInterpolant({xg,yg,zg},double(imgs));
xq = (1:1:sz(1))';
yq = (1:1:sz(2))';
zq = (1:1/ori_sliceness:sz(3))';
resample_imgs = int16(F({xq,yq,zq}));