医学图像处理案例（十二）——最小路径提取算法-CFANZ编程社区

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

1、最小路径提取算法

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

医学图像处理案例（十二）——最小路径提取算法_#include

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[]){// Typedefsconstexpr unsigned int Dimension = 3;using PixelType = float;using OutputPixelType = unsigned char;usingusingusingusingusingusingusingusingusing
try  {// Get filename arguments    ReaderType::Pointer reader = ReaderType::New();"10.nii.gz");    reader->Update();typedef      rescaleIntensityType;    rescaleIntensityType::Pointer rescaleFilter = rescaleIntensityType::New();    rescaleFilter->SetInput(reader->GetOutput());0.0);1.0);    rescaleFilter->SetNumberOfThreads(NUM_THREADS);    rescaleFilter->Update();
// Read speed function
    ImageType::Pointer speed = rescaleFilter->GetOutput();    speed->DisconnectPipeline();
// Create interpolatorusing    InterpolatorType::Pointer interp = InterpolatorType::New();
// Create cost function    PathFilterType::CostFunctionType::Pointer cost =      PathFilterType::CostFunctionType::New();    cost->SetInterpolator(interp);
// Create optimizerusing    OptimizerType::Pointer optimizer = OptimizerType::New();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);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 points0] = 210; end[1] = 234; end[2] = 1;0] = 191; start[1] = 178; start[2] = 169;0] = 218; way0[1] = 210; way0[2] = 192;// Add path informationusing    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 pathfor (unsigned int i = 0; i < pathFilter->GetNumberOfOutputs(); i++)    {// Get the path      PathType::Pointer path = pathFilter->GetOutput(i);
// Check path is validif (path->GetVertexList()->Size() == 0)      {std::cout << "WARNING: Path " << (i + 1) << " contains no points!" << std::endl;continue;      }
// Iterate path and convert to imagestd::cout << path->GetVertexList()->Size() << std::endl;for (int pathI = 0; pathI < path->GetVertexList()->Size(); pathI++)      {std::cout      }PathIteratorType it(output, path);for      {        it.Set(itk::NumericTraits<OutputPixelType>::max());      }    }
// Write output    WriterType::Pointer writer = WriterType::New();"MinimalPath.nii.gz");    writer->SetInput(output);    writer->Update();  }catch  {std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;    getchar();return  }}
int main(int argc, char* argv[]){return}