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

    // General includes#include #include // 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      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;    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::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::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);}

    pip install itk-minimalpathextraction