[SimpleITK] Python으로 medical image processing하기 - (2) Registration 1

[SimpleITK] Python으로 medical image processing하기 - (1) Load, Resample와 이어진다.

SimpleITK library의 간단한 사용법 역시 이 글에서 확인할 수 있다.

SimpleITK library를 이용해 Python으로 resample 하기

 

 

SimpleITK를 이용해 Python으로 processing 하기

[SimpleITK] Python으로 medical image processing하기 - (1) Load, Resample

[SimpleITK] Python으로 medical image processing하기 - (2) Registration 1

[SimpleITK] Python으로 medical image processing하기 - (3) Registration 2

[SimpleITK] Python으로 medical image processing하기 - (4) Visualization

 

 

SimpleITK를 이용해 Resample하기에 이어서, Registration 하는 방법을 소개해 보려고 한다.

 

우선, Registration을 하는 방법에는 두 가지가 있다.

 

첫 번째는 parameter map을 설정해 학습을 통해 registration parameter을 찾는 방법과, (ElastixImageFilter)

두 번째는 parameter file을 불러오는 방법이다. (TransformixImageFilter)

 

이 글에서는 첫 번째 방법을 설명하고, 두 번째 방법은 다음 글에서 설명하도록 한다.

 

 

Registration

Registration은 fixed image에 moving image를 맞추는 processing이다.

 

위와 같이 비슷한 구조의 두 이미지가 있는데, 위치를 맞추고 싶을 경우 사용할 수 있다.

(left) registration 전의 두 image의 overlap / (right) registration 후 두 image의 overlap

두 image 중 하나를 fixed image로 두고, 다른 하나를 moving image로 설정 후

fixed image에 moving image를 맞추는 방식이다.

 

위 예시의 경우, 간단한 rigid registration (linear transform) 만으로 두 image를 맞출 수 있지만,

구조/형태를 바꾸는 affine registration 이나 B-spline registration 을 이용해야 하는 경우도 있다.

 

각 rigid/affine/b-spline 등에 대해서는 이 글에서는 설명하지 않기로 한다.

간단히 말하자면, rigid는 이미지의 이동이나 회전, 대칭 등만을 포함하는 linear transform이고,

affine/b-spline 은 구조의 변형을 포함하는 nonlinear transform이다.

 

보통 복잡한 image의 registration을 할 때, 한 가지 종류의 registration만 사용한다기보다는

rigid -> affine -> B-splline 등과 같이 여러 종류의 registration을 순차적으로 적용하는 것 같다.

 

 

SimpleITK를 이용한 Registration

- registration 역시, resample과 마찬가지로 다음 세 가지 단계로 진행할 수 있다.

 

1) 전처리를 수행할 객체를 생성

• sitk.ElastixImageFilter()은 registration을 수행해주는 객체이다.
• 이 객체에 원하는 parameter을 설정해줌으로써 rigid/affine/b-spline 등 다양한 종류의 registration을 수행할 수 있다.

 

2) 객체에 parameter 설정

•  크게 세 가지 step으로 나뉜다.

       (1) SetFixedImage : Fixed image로 사용할 sitk image 객체를 설정해준다.

       (2) SetMovingImage : Moving image로 사용할 sitk image 객체를 설정해준다.

       (3) SetParameterMap : Registration에 필요한 parameter map을 설정해준다. (parameter map은 registration에 사용할 parameter의 종류를 담고 있는 map으로, 아래에서 자세히 설명)

 

3) Execute() method를 이용해 전처리 수행

• Execute() method를 하면 설정한 fixed image, moving image, parameter map을 이용해 registration이 진행된다.
• 이후 GetResultImage() method를 통해 결과 image (fixed image에 맞춰진 moving image)를 가져올 수 있다.

 

+) parameter value를 file 형태로 저장하기

•  Registration은 learning-base이기 때문에 registration이 끝나면 parameter value가 생성된다.
• 이 값들을 저장하여 다음 번에 다시 동일한 registration을 수행할 수 있다.
• Parameter map이 동일해도, parameter value가 다르면 다른 결과가 나올 수 있으므로 동일한 결과를 원하면 parameter value를 저장해야 한다.

 

import SimpleITK as sitk

"""1) Set ElastixImageFilter"""
elastixImageFilter = sitk.ElastixImageFilter()

"""2) Set Parameters"""
elastixImageFilter.SetFixedImage(sitk.ReadImage('fixedImage.nii'))
elastixImageFilter.SetMovingImage(sitk.ReadImage('movingImage.nii'))
elastixImageFilter.SetParameterMap(sitk.GetDefaultParameterMap('translation')) #method

"""3) Execute"""
elastixImageFilter.Execute()
resultimage = elastixImageFilter.GetResultImage()

#save parameter map
sitk.WriteParameterFile(elastixImageFilter.GetTransformParameterMap()[0], 'save_path')

 

 

Parameter map

Registration에 필요한 parameter들을 묶어 parameter map의 형태로 ElastixImageFilter에 전달해 주어야 한다.

Parameter map을 불러오는 방법에는 두 가지가 있는데, 하나는 Default parameter map을 불러오는 것이고, 다른 하나는 Parameter Map file을 불러오는 것이다. 

 

1. Default parameter map 불러오기

GetDefaultParamerMap() method를 이용해 Default parameter map을 불러올 수 있다.

parametermap = sitk.GetDefaultParameterMap('translation')
elastixImageFilter.SetParameterMap(parametermap)

 

PrintParameterMap() method를 이용하면 parameter map의 parameter들을 확인할 수 있다.

Default pixel value (background를 채울 pixel value의 값), learning에 필요한 metric과 optimizer, interpolation 방법 등을 지정할 수 있다.

sitk.PrintParameterMap(sitk.GetDefaultParameterMap("translation"))

>>>
ParameterMap 0:
  (AutomaticParameterEstimation "true")
  (CheckNumberOfSamples "true")
  (DefaultPixelValue 0)
  (FinalBSplineInterpolationOrder 2)
  (FixedImagePyramid "FixedSmoothingImagePyramid")
  (FixedImagePyramidSchedule 8 8 8 4 4 4 2 2 2 1 1 1)
  (ImageSampler "RandomCoordinate")
  (Interpolator "LinearInterpolator")
  (MaximumNumberOfIterations 32)
  (MaximumNumberOfSamplingAttempts 8)
  (Metric "AdvancedMattesMutualInformation")
  (MovingImagePyramid "MovingSmoothingImagePyramid")
  (MovingImagePyramidSchedule 8 8 8 4 4 4 2 2 2 1 1 1)
  (NewSamplesEveryIteration "true")
  (NumberOfResolutions 4)
  (NumberOfSamplesForExactGradient 4096)
  (NumberOfSpatialSamples 4096)
  (Optimizer "AdaptiveStochasticGradientDescent")
  (Registration "MultiResolutionRegistration")
  (ResampleInterpolator "FinalBSplineInterpolator")
  (Resampler "DefaultResampler")
  (Transform "TranslationTransform")
  (WriteResultImage "true")

 

이 Default parameter map을 그대로 SetParameterMap() method에 전달할 수도 있고,

다음과 같이 일부를 바꿔 사용할 수도 있다.

parameterMap['MaximumNumberOfIterations'] = ['512']

 

각 parameter에 대한 설명은 아래에서 확인할 수 있다. (열심히 찾았지만 더 친절한 docs는 찾을 수가 없었다..)

https://elastix.lumc.nl/doxygen/parameter.html

elastix: Parameters

 

2. Parameter Map file 불러오기

여러 Parameter을 바꿔가며 실험하기 귀찮다.. 그렇다면 다른 사람이 미리 설정해 놓은 Parameter map의 file을 불러올 수도 있다.

논문에서 저자가 자신들이 사용한 parameter map을 공개하는 경우가 있는데, 이렇게 잘 tuning된 parameter map을 사용하면 성능이 잘 나올 가능성이 높다.

 

https://elastix.lumc.nl/modelzoo/

Model Zoo Elastix

위 링크는 다양한 parameter map들이 모여 있는 database이다.

위 링크에서 parameter map file을 다운받아 사용할 수 있다.

 

parameter map file을 사용해 parameter map을 설정하려면 sitk.ReadParameterFile() method를 이용하면 된다.

Default parameter map과 마찬가지로, 이 parametermap을 바로 SetParameterMap() method에 전달하거나 수정하여 사용하면 된다.

parametermap = sitk.ReadParameterFile('file_path')
elastixImageFilter.SetParameterMap(parametermap)

 

---

 

sitk.WriteParameterFile(elastixImageFilter.GetTransformParameterMap()[0], 'save_path')

Registration을 마친 후 생성된 Parameter value들을 file로 저장하면 동일한 registration을 다른 image에 적용할 수 있다.

이는 sitk.TransformixImageFilter() 을 이용하면 되는데, 이 방법은 다음 게시글에서 설명하도록 하겠다.

 

 

 

참고: https://simpleelastix.readthedocs.io/HelloWorld.html#registration-with-translation-transform

Hello World — SimpleElastix 0.1 documentation