今天将分享人体血管两点间最小路径提取案例。
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血管上进行了测试,如图所示是血管原始图像,图中红色点目标就是最小路径提取结果。
如果碰到任何问题,随时留言,我会尽量去回答的。