医学图像处理及配准代码整理

医学图像处理开源库

  • SimpleITK
    安装方式 pip install SimpleITK
  • MedPy
    安装方式 pip install MedPy

其中SimpleITK中包含了医学图像数据读取,输入,保存,配准,分割的模块,MedPy目前只介绍数据读入的模块(其他模块暂未使用)

医学图像数据读入:

import SimpleITK as sitk
import medpy.io

#数据读入:
image = sitk.ReadImage('image_path')
#此时读入的图像的格式是SimpleITK image的格式,包含了医学图像的信息,常用的有Spacing, Origin, Direction,Dimension等

image_array = sitk.GetArrayFromImage(image)
#若只需要使用图像的灰度信息,即只需要numpy格式的医学图像灰度值组成的数组

#获取医学图像的其他附加信息
image_spacing = image.GetSpacing()  #空间分辨率
image_origin = image.GetOrigin()  #原点
image_direction = image.GetDirection() #方向
image_dimension = image.GetDimension() #维度
image_size = image.GetSize() #尺寸
#ps:一般深度学习进行医学图像处理时,只会用到图像的灰度信息,也就是只把图像当成数组来处理,暂未发现使用附加信息的。
#传统对医学图像的处理中,会考虑图像的分辨率,角度,方向等信息

#数据保存
#SimpleITK保存数据需要数据为SimpleITK image的格式,如果需要将数组保存为医学图像的格式,需做以下处理:
image_save = sitk.GetImageFromArray(image_array)
#保存
sitk.WriteImage(image_save, 'save_path+name.nii.gz') #医学图像格式,.nii, .nii.gz, .mhd, .mha


#数据读入使用medpy
image_medpy, header = medpy.io("image_path")
#其中image_medpy是该医学图像灰度值组成的数组,其轴的方向默认被转化为(x, y, z)
#headr中包含了四个部分,direction, offset, sitkimage, spacing,建议获取附加信息,仍使用sitkimage中信息
#for example
image_spacing = header.sitkimage.GetSpacing()
image_direction = header.sitkimage.GetDirection()
image_origin = header.sitkimage.GetOrigin()
image_size = header.sitkimage.GetSize()
#medpy这个库用的不是很多,SimpleITK的库大于medpy,所以建议使用SimpleITK

#如果想要给医学图像重新设置分辨率该如何实现呢?
#关于空间分辨率spacing和图像尺寸shape的关系不做介绍,直接给出修改分辨率的代码,以将图像分辨率调整到1*1*1为例

def Normalization(img):
    ary = sitk.GetArrayFromImage(img) #转变为数组
    MaxValue = ary.max() 
    MinValue = ary.min()
    objmax = np.ones_like(ary) * MaxValue   # Return an array of ones with the same shape and type as a given array.
    objmin = np.ones_like(ary) * MinValue
    NormalizedimageArray = 255.0 * (ary - objmin) / (objmax - objmin)
    NormalizedImage = sitk.GetImageFromArray(NormalizedimageArray)
    NormalizedImage.SetOrigin(img.GetOrigin()) #为图像添加原点信息
    NormalizedImage.SetSpacing(img.GetSpacing()) #为图像添加空间分辨率信息
    NormalizedImage.SetDirection(img.GetDirection()) #为图像添加方向信息
    return NormalizedImage

def Resampling(img,lable = False):
    original_size = img.GetSize() #获取图像原始尺寸
    original_spacing = img.GetSpacing() #获取图像原始分辨率
    new_spacing = [1, 1, 1] #设置图像新的分辨率为1*1*1
    new_size = [int(round(original_size[0] * (original_spacing[0] /1))),
                int(round(original_size[1] * (original_spacing[1] / 1))),
                int(round(original_size[2] * (original_spacing[2] / 1)))] #计算图像在新的分辨率下尺寸大小
    resampleSliceFilter = sitk.ResampleImageFilter() #初始化
    if lable == False:
        Resampleimage = resampleSliceFilter.Execute(img, new_size, sitk.Transform(), sitk.sitkBSpline,
                                                img.GetOrigin(), new_spacing, img.GetDirection(), 0,
                                                img.GetPixelIDValue())
        ResampleimageArray = sitk.GetArrayFromImage(Resampleimage)
        ResampleimageArray[ResampleimageArray < 0] = 0 #将图中小于0的元素置为0
    else:# for label, should use sitk.sitkLinear to make sure the original and resampled label are the same!!!
        Resampleimage = resampleSliceFilter.Execute(img, new_size, sitk.Transform(), sitk.sitkLinear,
                                                    img.GetOrigin(), new_spacing, img.GetDirection(), 0,
                                                    img.GetPixelIDValue())
        ResampleimageArray = sitk.GetArrayFromImage(Resampleimage)


    NormalizedImage = sitk.GetImageFromArray(ResampleimageArray) #对处理后的图像的灰度值进行归一化处理(可不做处理)
    return NormalizedImage

医学图像配准工具包(Python调用)

  • SimpleElastix (需要自己编译源码,官方文档有误,自己可以尝试踩坑)
  • pyelastix (直接拷贝Github代码,基础环境还是需要自己配置,linux下使用需要修改部分代码)
  • SimpleITK (简单粗暴,但方法过多,需要自己优化参数)

将配准的方式分为以下几种

  1. 使用SimpleITK的方式进行配准,同时可以将moving对应的mask配准到fixed上
def registration(fixed, moving, moving_mask):
    '''

    :param fixed: SimpleITK image, which can get spacing ,origin, direction
    :param moving: SimpleITK image, which can get spac  ing ,origin, direction
    :param moving_mask: SimpleITK image, which can get spacing ,origin, direction
    :return: registered moving image and registered mask image
    '''
    R = sitk.ImageRegistrationMethod()

    R.SetMetricAsCorrelation()
    R.SetOptimizerAsRegularStepGradientDescent(4.0, .01, 200)
    # R.SetOptimizerAsGradientDescent(learningRate=1.0,
    #                                 numberOfIterations=300)
    R.SetInitialTransform(sitk.TranslationTransform(fixed.GetDimension()))  # get InitialTransform
    R.SetInterpolator(sitk.sitkLinear)  # Interpolator params

    outTx = R.Execute(fixed, moving)  # get transform params

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(fixed)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(0)
    resampler.SetTransform(outTx)

    out_movingImage = resampler.Execute(moving)

    moving_mask.SetOrigin(moving.GetOrigin())
    moving_mask.SetOffset(moving.GetOffset())
    moving_mask.SetDirection(moving.GetDirection())
    moving_mask.SetSpacing(moving.GetSpacing())

    out_mask = resampler.Execute(moving_mask)

    return out_movingImage, out_mask

图示举例:
fixedImage

movingImage

从上面两幅中可以看到,fixedImage和movingImage的尺寸和分辨率都不相同,且fixedImage没有mask。我们将movingImage配准到fixedImage上,同时将配准的信息作用到movingMask上,使得配准后的mask也可以作用在fixedImage上。可以使用上述封装API
效果图如下:
registered movingImage

movingMask on fixedImage

从图中可以看出来,movingImage经过配准后,维度,分辨率都与fixed保持了一致,同时经过配准后,movingMask也可以直接作用在fixedImage上。

  1. 若直接对数组进行配准,不考虑图像的其他信息,使用Pyelastix
import pyelastix
import numpy as np
import os
import SimpleITK as sitk

os.environ['LD_LIBRARY_PATH'] = '/data/elastix/lib'
def imageRegistration(fixedImage, movingImage, isRigid=True):
    if isRigid:
        params = pyelastix.get_default_params(type='RIGID')
    else:
        params = pyelastix.get_default_params(type='BSPLINE')
    params.MaximumNumberOfIterations = 200
    params.FinalGridSpacingInVoxels = 10
    params.NumberOfResolutions = 3

    outputMovingImage, field = pyelastix.register(movingImage.astype('float32'), fixedImage.astype('float32'), params, verbose=0)
    outputMovingImage = np.array(outputMovingImage)
    return outputMovingImage

上述代码分为刚性和非刚性,参数默认,可以解决大部分配准问题,通时返回的为配准后的数组。

代补充:

对图像产生随机形变
根据产生形变场作用到图像上,2D和3D图像方式不同

发表新评论