首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >医学图像处理案例(十二)——最小路径提取算法

医学图像处理案例(十二)——最小路径提取算法

作者头像
医学处理分析专家
发布2020-06-29 15:43:12
1.5K0
发布2020-06-29 15:43:12
举报

今天将分享人体血管两点间最小路径提取案例。

1、最小路径提取算法

最小路径提取算法在很多领域都有广泛应用,医学图像分析,机器人导航等。2008年来自昆士兰科技大学的Dan Mueller开源了基于Fast Marching方式的最小路径提取算法,原理:利用Fast Marching到达函数T的梯度是与波前正交的事实来求解仅有一个的局部最小值,这也是全局最小值。通过从给定种子(路径终点)反向传播到起点来提取最小路径。起点和终点是隐式嵌入在T中的,反向传播可以通过梯度下降和正阶梯度下降来实现。

2、使用ITK函数来实现最小路径提取算法

Dan Mueller写了基于ITK的最小路径提取算法,C++源码下载请见原文链接。该函数使用时需要有三个输入,(1)、有意义的速度函数来生成到达函数,一般速度函数是归一化(0-1)的原始图像;(2)、起点(一个),终点(一个)和航点(路径必须经过其附近,多个)组成的路径信息;(3)、 优化器,其沿着到达函数垂直于Fast Marching面方向进行优化。详情可以阅读原文链接中的文章。

该函数既可以在C++中使用,也可以在Python中使用,下面将给出C++使用例子,并给出如何在Python上安装。

C++代码:


// General includes
#include <string>
#include <iostream>

// ITK includes
#include "itkNumericTraits.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkPolyLineParametricPath.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkArrivalFunctionToPathFilter.h"
#include "itkSpeedFunctionToPathFilter.h"
#include "itkPathIterator.h"
#include "itkGradientDescentOptimizer.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkIterateNeighborhoodOptimizer.h"
#include "itkRescaleIntensityImageFilter.h"

int example_gradientdescent(int argc, char* argv[])
{
  // Typedefs
  constexpr unsigned int Dimension = 3;
  using PixelType = float;
  using OutputPixelType = unsigned char;
  using ImageType = itk::Image< PixelType, Dimension >;
  using ucharImagetype = itk::Image< OutputPixelType, Dimension >;
  using OutputImageType = itk::Image< OutputPixelType, Dimension >;
  using ReaderType = itk::ImageFileReader< ucharImagetype >;
  using WriterType = itk::ImageFileWriter< OutputImageType >;
  using PathType = itk::PolyLineParametricPath< Dimension >;
  using PathFilterType = itk::SpeedFunctionToPathFilter< ImageType, PathType >;
  using CoordRepType = PathFilterType::CostFunctionType::CoordRepType;
  using PathIteratorType = itk::PathIterator< OutputImageType, PathType >;

  try
  {
    // Get filename arguments
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName("10.nii.gz");
    reader->Update();
    typedef itk::RescaleIntensityImageFilter<ucharImagetype, ImageType>
      rescaleIntensityType;
    rescaleIntensityType::Pointer rescaleFilter = rescaleIntensityType::New();
    rescaleFilter->SetInput(reader->GetOutput());
    rescaleFilter->SetOutputMinimum(0.0);
    rescaleFilter->SetOutputMaximum(1.0);
    rescaleFilter->SetNumberOfThreads(NUM_THREADS);
    rescaleFilter->Update();

    // Read speed function

    ImageType::Pointer speed = rescaleFilter->GetOutput();
    speed->DisconnectPipeline();

    // Create interpolator
    using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType, CoordRepType>;
    InterpolatorType::Pointer interp = InterpolatorType::New();

    // Create cost function
    PathFilterType::CostFunctionType::Pointer cost =
      PathFilterType::CostFunctionType::New();
    cost->SetInterpolator(interp);

    // Create optimizer
    using OptimizerType = itk::GradientDescentOptimizer;
    OptimizerType::Pointer optimizer = OptimizerType::New();
    optimizer->SetNumberOfIterations(1000);
    std::cout << "GetLargestPossibleRegion:" << speed->GetLargestPossibleRegion() << std::endl;
    // Create path filter
    PathFilterType::Pointer pathFilter = PathFilterType::New();
    pathFilter->SetInput(speed);
    pathFilter->SetCostFunction(cost);
    pathFilter->SetOptimizer(optimizer);
    pathFilter->SetTerminationValue(2);

    // Setup path points
    PathFilterType::PointType start, end, way0;
    //this is 2d image seed points
    /*end[0] = 289; end[1] = 252;
    start[0] = 232; start[1] = 490;
    way0[0] = 194; way0[1] = 309;*/
    //this is 3d image seed points
    end[0] = 210; end[1] = 234; end[2] = 1;
    start[0] = 191; start[1] = 178; start[2] = 169;
    way0[0] = 218; way0[1] = 210; way0[2] = 192;
    // Add path information
    using PathInformationType = PathFilterType::PathInformationType;
    PathInformationType::Pointer info = PathInformationType::New();
    info->SetStartPoint(start);
    info->SetEndPoint(end);
    info->AddWayPoint(way0);
    pathFilter->AddPathInformation(info);

    // Compute the path
    pathFilter->Update();

    // Allocate output image
    OutputImageType::Pointer output = OutputImageType::New();
    output->SetRegions(speed->GetLargestPossibleRegion());
    output->SetSpacing(speed->GetSpacing());
    output->SetOrigin(speed->GetOrigin());
    output->Allocate();
    output->FillBuffer(itk::NumericTraits<OutputPixelType>::Zero);

    // Rasterize path
    for (unsigned int i = 0; i < pathFilter->GetNumberOfOutputs(); i++)
    {
      // Get the path
      PathType::Pointer path = pathFilter->GetOutput(i);

      // Check path is valid
      if (path->GetVertexList()->Size() == 0)
      {
        std::cout << "WARNING: Path " << (i + 1) << " contains no points!" << std::endl;
        continue;
      }

      // Iterate path and convert to image
      std::cout << path->GetVertexList()->Size() << std::endl;
      for (int pathI = 0; pathI < path->GetVertexList()->Size(); pathI++)
      {
        std::cout << path->GetVertexList()->at(pathI);
      }
      PathIteratorType it(output, path);
      for (it.GoToBegin(); !it.IsAtEnd(); ++it)
      {
        it.Set(itk::NumericTraits<OutputPixelType>::max());
      }
    }

    // Write output
    WriterType::Pointer writer = WriterType::New();
    writer->SetFileName("MinimalPath.nii.gz");
    writer->SetInput(output);
    writer->Update();
  }
  catch (itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    getchar();
    return EXIT_FAILURE;
  }
}

int main(int argc, char* argv[])
{
  return example_gradientdescent(argc, argv);
}

Python安装指令:

pip install itk-minimalpathextraction

3、最小路径提取效果(2D和3D)

分别在DSA血管和CTA血管上进行了测试,如图所示是血管原始图像,图中红色点目标就是最小路径提取结果。

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

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-10-25,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 最新医学影像技术 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档