【ITK】最小路径提取

// 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"
#include<itkNiftiImageIOFactory.h>
//#include<itkThresholdImageFilter.h>

#define NUM_THREADS 4

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< ImageType >;
	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 >;
	itk::NiftiImageIOFactory::RegisterOneFactory();
	//using FilterType = itk::ThresholdImageFilter<ImageType>;

	try
	{
		// Get filename arguments
		ReaderType::Pointer reader = ReaderType::New();
		reader->SetFileName("test.nii");
		reader->Update();
		typedef itk::RescaleIntensityImageFilter<ImageType, 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();

		ImageType::Pointer speed = rescaleFilter->GetOutput();
		std::cout << "GetLargestPossibleRegion:" << speed->GetLargestPossibleRegion() << std::endl;
		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);
		
		// 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;//世界坐标的三维点，不是索引index坐标
		//this is 3d image seed points

		//start[0] = -7.86; start[1] = -11.46; start[2] = -160.77;
		//end[0] = 3.65; end[1] = -22.97; end[2] = -207.65;

		start[0] = 25.14; start[1] = -12.87; start[2] = -210.15;
		end[0] = 40.17; end[1] = -10.06; end[2] = -209.32;

		// 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(): "<< path->GetVertexList()->Size() << std::endl;
			//for (int pathI = 0; pathI < path->GetVertexList()->Size(); pathI++)
			//{
			//	std::cout <<pathI<<" "<< path->GetVertexList()->at(pathI) << std::endl;
			//}
			
			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");
		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);
}