专栏首页最新医学影像技术医学图像处理案例(十三)——快速行进算法分割医学图像

医学图像处理案例(十三)——快速行进算法分割医学图像

今天将分享使用快速行进算法(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为期望分割区域内任意点坐标即可。

如果碰到任何问题,随时留言,我会尽量去回答的。

本文分享自微信公众号 - 最新医学影像技术(MedicalHealthNews),作者:最新医学影像技术

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-12-27

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 医学图像处理案例(二十)——医学图像处理案例代码详解

    在前面分享的医学图像处理案例中,给出了很多具体案例,但有些读者还是渴望可以深入分享案例代码详解。那么今天我将从骨骼分割,气管分割,肺组织分割,血管分割这四个具体...

    用户7498388
  • 医学图像处理教程(四)——医学图像去噪算法

    均值滤波是典型的线性滤波算法,它是指在图像上对目标像素给一个模板,该模板包括了其周围的临近像素(以目标像素为中心的周围8个像素,构成一个滤波模板,即去掉目标像素...

    用户7498388
  • BraTS18——多模态MR图像脑肿瘤分割挑战赛续4

    今天将继续分享从网络结构上进行改进提出SEVNet模型来分割脑肿瘤。为了方便大家学习理解整个分割流程,我将整个流程步骤进行了整理,并给出每个步骤的结果,希望对大...

    用户7498388
  • 论文阅读——利用Binary Hash Codes的深度图像检索

    这篇文章是阅读《Deep Learning of Binary Hash Codes for Fast Image Retrieval》后的总结,该文章提出了...

    zhaozhiyong
  • 腾讯出品!这个免费网盘小程序,让你快速「找私货」| 亲儿子 #17

    在 2017 年 1 月 16 日前,腾讯微云为每个帐号分配了非常大的网盘容量,很多人都用它来当照片、文件等的「仓库」。

    知晓君
  • 当年的Windows 手机系统,为何不如安卓系统受欢迎?

    微软作为PC端的绝对霸主一直想在移动端开辟属于自己的战场,中间做了很多次的尝试都以失败而告终,最后孤注一掷拿下诺基亚结果windows手机还是处于不温不火的状态...

    程序员互动联盟
  • 03 Jme3和Nifty1.4.2中文显示

    用JMonkey最大的问题就是中文问题,IDE不是中文没关系,反正可以迁移到Idea里,但是打包发布的项目以及Nifty做的GUI里没有中文就心塞塞了。好在找到...

    刘开心_1266679
  • SVN版本管理工具的使用

    https://tortoisesvn.net/downloads.html 下载网站

    半指温柔乐
  • git 发送patch给别人

    用户3765803
  • Shell常见命令实践

    版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/sinat_35512245/articl...

    大黄大黄大黄

扫码关注云+社区

领取腾讯云代金券