[ITK-users] Vessel Segmentation

[Nissim Hadar]

Hello,

I am trying to use the /HessianToObjectnessMeasureImageFilter/ to segment
brain tissue.  As this code is part of a larger system, I am providing the
image as an array of /shorts/.  Other segmentation methods, such as Otsu,
flag the segmentation, so the result is a binary image (e.g. 255 and 0).  I
do not understand the results I am getting with the Hessian method.

My code is as follows:

*void HessianSegmentation(
	// Input data
	short			*inputVoxels,	const int		dim,
	const double    xSpacing,		const double	ySpacing,	const double	zSpacing,

	// Output data
	short			*outputVoxels
) {
	try
	{
		const   unsigned int Dimension = 3;

		typedef double                              InternalPixelType;
		typedef Image<InternalPixelType, Dimension> InternalImageType;
		
		typedef short                               VoxelPixelType;
		typedef Image<VoxelPixelType, Dimension>    VoxelImageType;

		// Image Importer - converts linear array of pixels in short format to an
Image
		//
		typedef ImportImageFilter<VoxelPixelType, Dimension> ImportFilterType;
		ImportFilterType::Pointer importFilter = ImportFilterType::New();

		ImportFilterType::IndexType start;
		start.Fill(0);
	
		ImportFilterType::SizeType size;
		size[0] = size[1] = size[2] = dim;
	
		ImportFilterType::RegionType region;
		region.SetIndex(start);
		region.SetSize (size );
		importFilter->SetRegion(region);

		const SpacePrecisionType origin[Dimension] = {0, 0, 0};
		importFilter->SetOrigin(origin);

		const SpacePrecisionType spacing[Dimension] = {xSpacing, ySpacing,
zSpacing};
		importFilter->SetSpacing(spacing);

		importFilter->SetImportPointer(inputVoxels, dim * dim * dim, false);

		// Caster from short to float (input is shorts, detection needs floats)
		typedef CastImageFilter<VoxelImageType, InternalImageType>
CastFilterTypeVoxelToInternal;		
		CastFilterTypeVoxelToInternal::Pointer casterVoxelToInternal =
CastFilterTypeVoxelToInternal::New();

		casterVoxelToInternal->	SetInput(importFilter->GetOutput());

		// Set up Hessian filter
		//
		typedef SymmetricSecondRankTensor<double, Dimension> HessianPixelType;
		typedef Image<HessianPixelType, Dimension>           HessianImageType;
		typedef HessianToObjectnessMeasureImageFilter<HessianImageType,
InternalImageType> ObjectnessFilterType;
		ObjectnessFilterType::Pointer objectnessFilter =
ObjectnessFilterType::New();
		
		objectnessFilter->SetBrightObject(true);
		objectnessFilter->SetScaleObjectnessMeasure(false);
		objectnessFilter->SetAlpha(0.5);
		objectnessFilter->SetBeta (0.5);
		objectnessFilter->SetGamma(5.0);

		typedef MultiScaleHessianBasedMeasureImageFilter<InternalImageType,
HessianImageType, InternalImageType> MultiScaleEnhancementFilterType;
		MultiScaleEnhancementFilterType::Pointer multiScaleEnhancementFilter =
MultiScaleEnhancementFilterType::New();
		multiScaleEnhancementFilter->SetHessianToMeasureFilter(objectnessFilter);
		multiScaleEnhancementFilter->SetSigmaStepMethodToLogarithmic();

		multiScaleEnhancementFilter->SetSigmaMinimum(0.2);
		multiScaleEnhancementFilter->SetSigmaMaximum(2.0);
		multiScaleEnhancementFilter->SetNumberOfSigmaSteps(10);

		multiScaleEnhancementFilter->SetInput(casterVoxelToInternal->GetOutput());

		// Setup rescaler
		//
		typedef RescaleIntensityImageFilter<InternalImageType, InternalImageType>
RescaleFilterType;
		RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
		rescaleFilter->SetOutputMinimum(0);
		rescaleFilter->SetOutputMaximum(2000);

		rescaleFilter->	SetInput(multiScaleEnhancementFilter->GetOutput());

		// Cast back from  float to short
		typedef CastImageFilter<InternalImageType, VoxelImageType>
CastFilterTypeInternalToVoxel;		
		CastFilterTypeInternalToVoxel::Pointer casterInternalToVoxel =
CastFilterTypeInternalToVoxel::New();

		casterInternalToVoxel->SetInput(rescaleFilter->GetOutput());
		
		// Run pipeline
		//
		casterInternalToVoxel->Update();

		// Extract results
		//
		VoxelImageType *outputImage = casterInternalToVoxel->GetOutput();
		InternalImageType::RegionType outputRegion =
outputImage->GetBufferedRegion();

		typedef itk::ImageRegionConstIterator<VoxelImageType> IteratorType;
		IteratorType it(outputImage, outputRegion);

		VoxelPixelType *ptr = outputVoxels;

		it.GoToBegin();
		while(!it.IsAtEnd())
		{
			*ptr = it.Get();
			++it;
			++ptr;
		}
	}
	catch (const std::exception& ex)
	{
		std::string s = ex.what();
		bool error = true;			// put breakpoint here to debug
	}
}*

The input data is a block of voxels from CT/MRI images.

Thank you very much
Cheers



--
View this message in context: http://itk-users.7.n7.nabble.com/Vessel-Segmentation-tp35113.html
Sent from the ITK - Users mailing list archive at Nabble.com.