医学图像处理案例（十三）——快速行进算法分割医学图像-CFANZ编程社区

今天将分享使用快速行进算法（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>
namespaceint main(int argc, char{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    }
const std::string inputFilename(argv[1]);const std::string outputFilename(argv[2]);
unsigned int seedPosition[2];0] = atoi( argv[3] );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;0.12559.0  sitk::Image smoothingOutput = smoothing.Execute( inputImage );
  sitk::GradientMagnitudeRecursiveGaussianImageFilter gradientMagnitude;  gradientMagnitude.SetSigma(  sigma  );  sitk::Image gradientMagnitudeOutput = gradientMagnitude.Execute( smoothingOutput );
  sitk::SigmoidImageFilter sigmoid;0.01.0  sigmoid.SetAlpha( alpha );  sigmoid.SetBeta(  beta  );  sitk::Image sigmoidOutput = sigmoid.Execute( gradientMagnitudeOutput );

  sitk::FastMarchingImageFilter fastMarching;
std::vector< unsigned int > trialPoint(3);

0] = seedPosition[0];1] = seedPosition[1];
2] = 0u; // Seed Value
  fastMarching.AddTrialPoint( trialPoint );
  fastMarching.SetStoppingValue(stoppingTime);
  sitk::Image fastmarchingOutput = fastMarching.Execute( sigmoidOutput );

  sitk::BinaryThresholdImageFilter thresholder;0.0  thresholder.SetUpperThreshold( timeThreshold  );0255
  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 sitkimport sysimport 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()
seedValuetrialPoint = (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为期望分割区域内任意点坐标即可。

医学图像处理案例（十三）——快速行进算法分割医学图像_#include