今天将分享使用快速行进算法(FastMarching)对医学图像分割案例。
1、FastMarching简介
快速行进方法(FastMarching)是水平集演化方法的一种简化形式,其仅使用正速度项来控制微分方程,生成的水平集轮廓随着时间增长。在实际中,FastMarching算法可以看作是由速度图像控制的高级区域增长分割方法。该算法具体推导请参考原文连接。
2、使用SimpleITK函数来实现FastMarching分割算法
用FastMarching算法分割有5个步骤:(1)、首先使用各向异性扩散方法对输入图像进行平滑处理;(2)、其次对平滑后的图像进行梯度计算,生成边缘图像,在梯度计算过程中可调节高斯sigma参数,来控制水平集减速到接近边缘;(3)、然后使用逻辑回归(Sigmoid)函数对边缘图像进行线性变换,保证边界接区域近零,平坦区域接近1,回归可调参数有alpha和beta;(4)、接着手动设置置FastMarching算法的初始种子点和起始值,该种子点是水平集的起始位置。FastMarching的输出是时间跨度图,表示传播的水平集面到达的时间;(5)、最后通过阈值方法将FastMarching结果限制在水平集面传播区域而形成分割的区域。
该例子既可以在C++中使用,也可以在Python中使用,下面将给出C++和Python的使用例子代码。
C++代码:
*=========================================================================
// This example is based on ITK's FastMarchingImageFilter.cxx example
#include <SimpleITK.h>
#include <iostream>
#include <string>
#include <cstdlib>
namespace sitk = itk::simple;
int main(int argc, char *argv[])
{
if ( argc < 10 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " inputImage outputImage seedX seedY";
std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingTime" << std::endl;
return EXIT_FAILURE;
}
const std::string inputFilename(argv[1]);
const std::string outputFilename(argv[2]);
unsigned int seedPosition[2];
seedPosition[0] = atoi( argv[3] );
seedPosition[1] = atoi( argv[4] );
const double sigma = atof( argv[5] );
const double alpha = atof( argv[6] );
const double beta = atof( argv[7] );
const double timeThreshold = atof( argv[8] );
const double stoppingTime = atof( argv[9] );
sitk::Image inputImage = sitk::ReadImage( inputFilename, sitk::sitkFloat32 );
sitk::CurvatureAnisotropicDiffusionImageFilter smoothing;
smoothing.SetTimeStep( 0.125 );
smoothing.SetNumberOfIterations( 5 );
smoothing.SetConductanceParameter( 9.0 );
sitk::Image smoothingOutput = smoothing.Execute( inputImage );
sitk::GradientMagnitudeRecursiveGaussianImageFilter gradientMagnitude;
gradientMagnitude.SetSigma( sigma );
sitk::Image gradientMagnitudeOutput = gradientMagnitude.Execute( smoothingOutput );
sitk::SigmoidImageFilter sigmoid;
sigmoid.SetOutputMinimum( 0.0 );
sigmoid.SetOutputMaximum( 1.0 );
sigmoid.SetAlpha( alpha );
sigmoid.SetBeta( beta );
sitk::Image sigmoidOutput = sigmoid.Execute( gradientMagnitudeOutput );
sitk::FastMarchingImageFilter fastMarching;
std::vector< unsigned int > trialPoint(3);
trialPoint[0] = seedPosition[0];
trialPoint[1] = seedPosition[1];
trialPoint[2] = 0u; // Seed Value
fastMarching.AddTrialPoint( trialPoint );
fastMarching.SetStoppingValue(stoppingTime);
sitk::Image fastmarchingOutput = fastMarching.Execute( sigmoidOutput );
sitk::BinaryThresholdImageFilter thresholder;
thresholder.SetLowerThreshold( 0.0 );
thresholder.SetUpperThreshold( timeThreshold );
thresholder.SetOutsideValue( 0 );
thresholder.SetInsideValue( 255 );
sitk::Image result = thresholder.Execute(fastmarchingOutput);
sitk::WriteImage(result, outputFilename);
return 0;
}
Python代码:
#!/usr/bin/env python
from __future__ import print_function
import SimpleITK as sitk
import sys
import os
if len(sys.argv) < 10:
print("Usage: {0} <inputImage> <outputImage> <seedX> <seedY> <Sigma> <SigmoidAlpha> <SigmoidBeta> <TimeThreshold>".format(sys.argv[0]))
sys.exit(1)
inputFilename = sys.argv[1]
outputFilename = sys.argv[2]
seedPosition = (int(sys.argv[3]), int(sys.argv[4]))
sigma = float(sys.argv[5])
alpha = float(sys.argv[6])
beta = float(sys.argv[7])
timeThreshold = float(sys.argv[8])
stoppingTime = float(sys.argv[9])
inputImage = sitk.ReadImage(inputFilename, sitk.sitkFloat32)
print(inputImage)
smoothing = sitk.CurvatureAnisotropicDiffusionImageFilter()
smoothing.SetTimeStep(0.125)
smoothing.SetNumberOfIterations(5)
smoothing.SetConductanceParameter(9.0)
smoothingOutput = smoothing.Execute(inputImage)
gradientMagnitude = sitk.GradientMagnitudeRecursiveGaussianImageFilter()
gradientMagnitude.SetSigma(sigma)
gradientMagnitudeOutput = gradientMagnitude.Execute(smoothingOutput)
sigmoid = sitk.SigmoidImageFilter()
sigmoid.SetOutputMinimum(0.0)
sigmoid.SetOutputMaximum(1.0)
sigmoid.SetAlpha(alpha)
sigmoid.SetBeta(beta)
sigmoid.DebugOn()
sigmoidOutput = sigmoid.Execute(gradientMagnitudeOutput)
fastMarching = sitk.FastMarchingImageFilter()
seedValue = 0
trialPoint = (seedPosition[0], seedPosition[1], seedValue)
fastMarching.AddTrialPoint(trialPoint)
fastMarching.SetStoppingValue(stoppingTime)
fastMarchingOutput = fastMarching.Execute(sigmoidOutput)
thresholder = sitk.BinaryThresholdImageFilter()
thresholder.SetLowerThreshold(0.0)
thresholder.SetUpperThreshold(timeThreshold)
thresholder.SetOutsideValue(0)
thresholder.SetInsideValue(255)
result = thresholder.Execute(fastMarchingOutput)
sitk.WriteImage(result, outputFilename)
3、FastMarching分割效果
在MRI脑部图像上进行脑室、灰质和白质的分割测试,如图所示依次是MRI原始图像,左脑室分割结果,右脑室分割结果,白质分割结果,灰质分割结果。分割算法的参数典型设置:sigma=0.5,alpha=-0.3,beta=2.0, timeThreshold=200, stoppingTime=210,seedPosition为期望分割区域内任意点坐标即可。
如果碰到任何问题,随时留言,我会尽量去回答的。