海塞矩阵
#include "itkImage.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileWriter.h"
#include "itkNumericSeriesFileNames.h"
#include "itkGDCMImageIO.h"#include "itkShiftScaleImageFilter.h"
#include "itkHessianRecursiveGaussianImageFilter.h"
#include "itkHessian3DToVesselnessMeasureImageFilter.h"#include "itkRescaleIntensityImageFilter.h"
#include "itkShrinkImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"//二值化
#include "itkThresholdImageFilter.h"int main( int argc, char *argv[] )
{const char * inputImage = "D:\\se6\\ct3\\image%03d.dcm";const char * outputImage = "d:\\3.vtk";const unsigned int Dimension = 3;typedef short InputPixelType;typedef unsigned short OutputPixelType;typedef itk::Image< InputPixelType, Dimension > InputImageType;typedef itk::Image< OutputPixelType, Dimension > OutputImageType;typedef itk::NumericSeriesFileNames NameGeneratorType ;NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();nameGenerator->SetSeriesFormat(inputImage);nameGenerator->SetStartIndex(0);nameGenerator->SetEndIndex(213);nameGenerator->SetIncrementIndex(1);typedef itk::ImageSeriesReader<InputImageType> ReaderType ;ReaderType::Pointer reader = ReaderType::New();reader->SetImageIO(itk::GDCMImageIO::New());reader->SetFileNames(nameGenerator->GetFileNames());reader->Update();std::cout<<"reader updata"<<std::endl;typedef itk::ShrinkImageFilter <InputImageType, InputImageType> ShrinkImageFilterType;std::cout << "Original size: " << reader->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;ShrinkImageFilterType::Pointer shrinkFilter= ShrinkImageFilterType::New();shrinkFilter->SetInput(reader->GetOutput());shrinkFilter->SetShrinkFactor(0, 2); //shrinkFilter->SetShrinkFactor(1, 2); //shrinkFilter->SetShrinkFactor(2, 2); //shrinkFilter->Update();typedef itk::ThresholdImageFilter <InputImageType> ThresholdImageFilterType;ThresholdImageFilterType::Pointer thresholdFilter= ThresholdImageFilterType::New();thresholdFilter->SetInput(shrinkFilter->GetOutput());thresholdFilter->ThresholdOutside(-224, 676);thresholdFilter->SetOutsideValue(0);typedef itk::HessianRecursiveGaussianImageFilter< InputImageType >HessianFilterType;HessianFilterType::Pointer hessianFilter = HessianFilterType::New();hessianFilter->SetInput( thresholdFilter->GetOutput() );hessianFilter->SetSigma(1.0 ); typedef itk::Hessian3DToVesselnessMeasureImageFilter< OutputPixelType >VesselnessMeasureFilterType;VesselnessMeasureFilterType::Pointer vesselnessFilter = VesselnessMeasureFilterType::New();vesselnessFilter->SetInput( hessianFilter->GetOutput() );vesselnessFilter->SetAlpha1( 0.1 );vesselnessFilter->SetAlpha2( 2.0 );typedef itk::ImageFileWriter< OutputImageType > WriterType;WriterType::Pointer writer = WriterType::New();writer->SetInput( vesselnessFilter->GetOutput() );writer->SetFileName( outputImage );try{writer->Update();}catch( itk::ExceptionObject & error ){std::cerr << "Error: " << error << std::endl;return EXIT_FAILURE;}return EXIT_SUCCESS;
}