[Insight-users] Visualization problem in vtk after extracting vascular by level set segmentation

[Zhenqiu Feng]

Hi all,
	I need some advice. I’m triying to segment the artery from CT Dicom
datasets and visualize it in VTK. Thanks to the previous mailing list
archive, I noticed that the example “itkReadITKImage3DSegmentShowVTK.cxx” in
the InsightApplications is quite helpful. I followed the example, used
ConnectedThresholdImageFilter to segment the artery and visualized it in
vtk, it worked fine.(figure1 is the screenshot). 
But when I replaced region growing method with level set
(ShapeDetectionLevelSetFilter, ThresholdSegmentationLevelSetImageFilter),
problem came up, the segmented artery can’t display within
vtkImagePlaneWidget with 3d CT image in VTK, that is, the segmented artery
is out of vtkImagePlaneWidget(figure2 shows). I have excluded the reason of
dataset and have no idea why this happened.
Dose anyone has the similar experience and resolved this issue? Anyone knows
what’s the reason and what should I do?

http://itk-insight-users.2283740.n2.nabble.com/file/n7561343/figure1-segmented_artery_within_imageplanewidget.png
figure1-segmented_artery_within_imageplanewidget.png 
http://itk-insight-users.2283740.n2.nabble.com/file/n7561343/figure2-segmented_artery_out_of_vtkImagePlaneWidget.png
figure2-segmented_artery_out_of_vtkImagePlaneWidget.png 

The ITK version I am using is 3.2, the VTK version I am using is 5.8 and
compiled with MS visual c++ 2008. Here’s my code:


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkThresholdSegmentationLevelSetImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkVTKImageExport.h"

#include "vtkDataSet.h"
#include "vtkImageData.h"
#include "vtkPolyData.h"
#include "vtkImageShrink3D.h"
#include "vtkContourFilter.h"
#include "vtkSmoothPolyDataFilter.h"
#include "vtkDecimatePro.h"
#include "vtkImageImport.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkInteractorStyleTrackballCamera.h"
#include "vtkProperty.h"
#include "vtkCamera.h"
#include "vtkCellPicker.h"
#include "vtkImagePlaneWidget.h"
#include "vtkTextProperty.h"
#include "vtkSTLWriter.h"
#include "vtkSTLReader.h"

template <typename ITK_Exporter,typename VTK_Importer>   
void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer)   
{   

importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());   

importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());   
	importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());   
	importer->SetSpacingCallback(exporter->GetSpacingCallback());   
	importer->SetOriginCallback(exporter->GetOriginCallback());   
	importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());   

importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());   

importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());   
	importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());   
	importer->SetDataExtentCallback(exporter->GetDataExtentCallback());   
	importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());   
	importer->SetCallbackUserData(exporter->GetCallbackUserData());   
} 

int main(int argc, char * argv[])
{
	argv[1] = "D:/MedialImag/Data/ThresholdSegLevelSetTest";				//Input file
name
	argv[2] = "";					//Output file name 
	argv[3] = "3.0";				//Initial distance
	argv[4] = "1.0";				//Propagation Scaling
	argv[5] = "1.0";				//Curvature Scaling
	argv[6] = "600";				//Upper threshold
	argv[7] = "180";				//Lower threshold
	argv[8] = "0.02";				//RMS
	argv[9] = "2000";				//Number of Iteration
	argv[10] = "243";				//Seed X
	argv[11] = "374";				//seed y
	argv[12] = "5";					//seed z

	typedef float		PixelType;
	const unsigned int	Dimension = 3;
	typedef itk::Image< PixelType, Dimension >	ImageType;
	typedef itk::ImageSeriesReader< ImageType >	ReaderType;
	typedef itk::GDCMImageIO					ImageIOType;
	typedef itk::GDCMSeriesFileNames			NameGeneratorType;

	ImageIOType::Pointer gdcmIO = ImageIOType::New();
	NameGeneratorType::Pointer nameGenerator = NameGeneratorType::New();

	nameGenerator->SetInputDirectory( argv[1] );
	const ReaderType::FileNamesContainer & filenames = 
		nameGenerator->GetInputFileNames();

	unsigned int numberOfFilenames = filenames.size();
	std::cout << "Number of CT files:" << numberOfFilenames << std::endl;
	for(unsigned int i = 0; i < numberOfFilenames; i++)
	{
		std::cout << "filename #" << i << "=";
		std::cout << filenames[i] << std::endl;
	}

	ReaderType::Pointer reader = ReaderType::New();
	reader->SetImageIO( gdcmIO );
	reader->SetFileNames( filenames );

	try
	{
		std::cout << "Reading CT files..." << std::endl;
		reader->Update();
		std::cout << "Reading CT files...Completed" << std::endl;
	}
	catch(itk::ExceptionObject &excp)
	{
		std::cout << "Exception thrown when reading CT files." << std::endl;
		std::cout << excp << std::endl;
		return EXIT_FAILURE;
	}

	//Use fast marching to generate the initial level set
	typedef itk::FastMarchingImageFilter< ImageType, ImageType >
FastMarchingFilterType;
	FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();
	typedef FastMarchingFilterType::NodeContainer		NodeContainerType;
	typedef FastMarchingFilterType::NodeType			NodeType;
	NodeContainerType::Pointer seeds = NodeContainerType::New();

	ImageType::IndexType seedPosition1;
	seedPosition1[0] = atoi( argv[10] );
	seedPosition1[1] = atoi( argv[11] );
	seedPosition1[2] = atoi( argv[12] );

	NodeType node1;
	const double initialDistance = atof( argv[3] );
	const double seedValue = -initialDistance;
	node1.SetValue( seedValue );
	node1.SetIndex( seedPosition1 );
	
	seeds->Initialize();
	seeds->InsertElement(0, node1);

	fastMarching->SetTrialPoints( seeds );
	fastMarching->SetOutputSize(
reader->GetOutput()->GetBufferedRegion().GetSize());

	try
	{
		std::cout << "Fastmarching..." << std::endl;
		fastMarching->Update();
		std::cout << "Fastmarching... Completed" << std::endl;
	}
	catch (itk::ExceptionObject &excp)
	{
		std::cerr << "Exception thrown while computing initial level set!" <<
std::endl;
		std::cerr << excp << std::endl;
		return EXIT_FAILURE;
	}


	//Threshold Segmentation Level set
	typedef itk::ThresholdSegmentationLevelSetImageFilter< ImageType,
		ImageType>	ThresholdLeveSetFilterType;
	ThresholdLeveSetFilterType::Pointer thresholdSegmentation = 
		ThresholdLeveSetFilterType::New();

	thresholdSegmentation->SetPropagationScaling( atof(argv[4]) );
	thresholdSegmentation->SetCurvatureScaling( atof(argv[5]) );
	thresholdSegmentation->SetUpperThreshold( atof(argv[6]) );
	thresholdSegmentation->SetLowerThreshold( atof(argv[7]) );
	thresholdSegmentation->SetIsoSurfaceValue( 0.0 );
	thresholdSegmentation->SetMaximumRMSError( atof(argv[8]) );
	thresholdSegmentation->SetNumberOfIterations( atof(argv[9]) );
	
	thresholdSegmentation->SetInput( fastMarching->GetOutput() );
	thresholdSegmentation->SetFeatureImage( reader->GetOutput() );
	thresholdSegmentation->Update();

	typedef itk::BinaryThresholdImageFilter<ImageType,ImageType> ThresholdType;
	ThresholdType::Pointer thresholder = ThresholdType::New();

	thresholder->SetLowerThreshold( -1000 );
	thresholder->SetUpperThreshold( 0 );

	thresholder->SetOutsideValue( 0 );
	thresholder->SetInsideValue( 255 );

	thresholder->SetInput( thresholdSegmentation->GetOutput() );

	// Print out some useful information 
	std::cout << std::endl;
	std::cout << "Max. no. iterations: " <<
thresholdSegmentation->GetNumberOfIterations() << std::endl;
	std::cout << "Max. RMS error: " <<
thresholdSegmentation->GetMaximumRMSError() << std::endl;
	std::cout << std::endl;
	std::cout << "No. elpased iterations: " <<
thresholdSegmentation->GetElapsedIterations() << std::endl;
	std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() <<
std::endl;

	typedef itk::VTKImageExport <ImageType>			ExportFilterType;
	ExportFilterType::Pointer itkExporter1 = ExportFilterType::New();
	ExportFilterType::Pointer itkExporter2 = ExportFilterType::New();
	itkExporter1->SetInput( reader->GetOutput() );
	itkExporter2->SetInput( thresholder->GetOutput() );


	//Connect itk to vtk
	vtkImageImport *vtkImporter1 = vtkImageImport::New();
	vtkImageImport *vtkImporter2 = vtkImageImport::New();
	ConnectPipelines( itkExporter1, vtkImporter1 );
	ConnectPipelines( itkExporter2, vtkImporter2 );

	vtkImporter1->Update();

	ImageType::Pointer inputImage = reader->GetOutput();   
	ImageType::SizeType  size  = inputImage->GetBufferedRegion().GetSize(); 

	vtkRenderer *renderer = vtkRenderer::New();
	vtkRenderWindow *renWin = vtkRenderWindow::New();
	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
	vtkInteractorStyleTrackballCamera *trackballCameraStyle = 
		vtkInteractorStyleTrackballCamera::New();

	renWin->SetSize( 1270,950 );
	renWin->AddRenderer( renderer );
	iren->SetRenderWindow( renWin );
	iren->SetInteractorStyle( trackballCameraStyle );

	vtkCellPicker *picker = vtkCellPicker::New();
	picker->SetTolerance( 0.005 );

	vtkProperty *ipwProp = vtkProperty::New();

	vtkImagePlaneWidget *xImagePlaneWidget = vtkImagePlaneWidget::New();
	vtkImagePlaneWidget *yImagePlaneWidget = vtkImagePlaneWidget::New();
	vtkImagePlaneWidget *zImagePlaneWidget = vtkImagePlaneWidget::New();

	xImagePlaneWidget->DisplayTextOn();
	xImagePlaneWidget->SetInput(vtkImporter1->GetOutput());
	xImagePlaneWidget->SetPlaneOrientationToXAxes();
	xImagePlaneWidget->SetSliceIndex( size[0]/2 );
	xImagePlaneWidget->SetPicker( picker );
	xImagePlaneWidget->RestrictPlaneToVolumeOn();
	xImagePlaneWidget->SetKeyPressActivationValue('x');
	xImagePlaneWidget->GetPlaneProperty()->SetColor(1,1,0);
	xImagePlaneWidget->SetTexturePlaneProperty(ipwProp);
	xImagePlaneWidget->SetResliceInterpolateToNearestNeighbour();
	xImagePlaneWidget->GetTextProperty()->SetColor(1,0,0);
	xImagePlaneWidget->GetSelectedPlaneProperty()->SetColor(1,0,0);
	xImagePlaneWidget->GetCursorProperty()->SetColor(1,0,0);

	yImagePlaneWidget->DisplayTextOn();   
	yImagePlaneWidget->SetInput(vtkImporter1->GetOutput());   
	yImagePlaneWidget->SetPlaneOrientationToYAxes();   
	yImagePlaneWidget->SetSliceIndex(size[1]/2);   
	yImagePlaneWidget->SetPicker(picker);   
	yImagePlaneWidget->RestrictPlaneToVolumeOn();   
	yImagePlaneWidget->SetKeyPressActivationValue('y');   
	yImagePlaneWidget->GetPlaneProperty()->SetColor(0, 1, 0);   
	yImagePlaneWidget->SetTexturePlaneProperty(ipwProp);   
	yImagePlaneWidget->SetLookupTable(xImagePlaneWidget->GetLookupTable());  
	yImagePlaneWidget->GetTextProperty()->SetColor(1, 0, 0);
	yImagePlaneWidget->GetSelectedPlaneProperty()->SetColor(1, 0, 0);
	yImagePlaneWidget->GetCursorProperty()->SetColor(1, 0, 0);

	zImagePlaneWidget->DisplayTextOn();   
	zImagePlaneWidget->SetInput(vtkImporter1->GetOutput());   
	zImagePlaneWidget->SetPlaneOrientationToZAxes();   
	zImagePlaneWidget->SetSliceIndex(size[2]/2);   
	zImagePlaneWidget->SetPicker(picker);   
	zImagePlaneWidget->SetKeyPressActivationValue('z');   
	zImagePlaneWidget->GetPlaneProperty()->SetColor(0, 0, 1);   
	zImagePlaneWidget->SetTexturePlaneProperty(ipwProp);   
	zImagePlaneWidget->SetLookupTable(xImagePlaneWidget->GetLookupTable()); 
	zImagePlaneWidget->GetTextProperty()->SetColor(1, 0, 0);
	zImagePlaneWidget->GetSelectedPlaneProperty()->SetColor(1, 0, 0);
	zImagePlaneWidget->GetCursorProperty()->SetColor(1, 0, 0);

	xImagePlaneWidget->SetInteractor( iren );   
	xImagePlaneWidget->On();   

	yImagePlaneWidget->SetInteractor( iren );   
	yImagePlaneWidget->On();   

	zImagePlaneWidget->SetInteractor( iren );   
	zImagePlaneWidget->On();  

	renderer->SetBackground(0.4392, 0.5020, 0.5647);  
	
	vtkImageShrink3D *shrink = vtkImageShrink3D::New();
	shrink->SetInput( vtkImporter2->GetOutput() );
	shrink->SetShrinkFactors( 3,3,1 );

	vtkContourFilter *contour = vtkContourFilter::New();
	contour->SetInput( shrink->GetOutput() );
	contour->SetValue( 0,128 );

	vtkSmoothPolyDataFilter *smooth = vtkSmoothPolyDataFilter::New();
	smooth->SetInput( contour->GetOutput() );
	smooth->SetNumberOfIterations( 20 );
	smooth->BoundarySmoothingOn();
	smooth->SetFeatureAngle( 100 );
	smooth->SetEdgeAngle( 100 );

	vtkDecimatePro *deci = vtkDecimatePro::New();
	deci->SetInputConnection(smooth->GetOutputPort() );
	deci->SetTargetReduction( 0.3 );

	vtkSTLWriter * writeSTL = vtkSTLWriter::New();
	writeSTL->SetFileName("poly.stl");
	writeSTL->SetInput( deci->GetOutput() );
	writeSTL->Update();

	vtkPolyDataMapper *polyMapper = vtkPolyDataMapper::New();
	vtkActor		*polyActor = vtkActor::New();

	polyActor->SetMapper( polyMapper );
	polyMapper->SetInput( deci->GetOutput() );
	polyMapper->ScalarVisibilityOff();

	vtkProperty *property = vtkProperty::New();
	property->SetAmbient( 0.1 );
	property->SetDiffuse( 0.1 );
	property->SetSpecular( 0.5 );
	property->SetColor( 1.0, 0.0, 0.0 );
	property->SetLineWidth( 2.0 );
	property->SetRepresentationToSurface();

	polyActor->SetProperty( property );

	renderer->AddActor( polyActor );

	vtkCamera *camera = vtkCamera::New();
	camera->SetPosition( 0,-1,0 );
	camera->SetFocalPoint( 0,0,0 );
	camera->SetViewUp(0,0,1);
	camera->Azimuth(45);
	camera->Pitch(-10);
	camera->Roll(30);
	renderer->SetActiveCamera( camera );
	renderer->ResetCamera();


	renWin->Render();
	iren->Start();

	polyActor->Delete();   
	picker->Delete();   
	ipwProp->Delete();   
	vtkImporter1->Delete();   
	vtkImporter2->Delete();   
	xImagePlaneWidget->Delete();   
	yImagePlaneWidget->Delete();   
	zImagePlaneWidget->Delete();   
	contour->Delete();   
	property->Delete();   
	polyMapper->Delete();   
	renWin->Delete();   
	renderer->Delete();   
	iren->Delete();  
	shrink->Delete();   
	smooth->Delete();   
	deci->Delete(); 
}


--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/Visualization-problem-in-vtk-after-extracting-vascular-by-level-set-segmentation-tp7561343.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.