[ITK-users] Vessel Segmentation

[Matt McCormick]

Hi Nissim,

It is unclear from your message what you do not understand about the results.

More information on the filters can be found in Doxygen [1] and The
ITK Software Guide [2].

HTH,
Matt

[1]  http://www.itk.org/Doxygen/html/classitk_1_1HessianToObjectnessMeasureImageFilter.html

[2]  http://itk.org/ItkSoftwareGuide.pdf

On Wed, Jan 14, 2015 at 5:40 PM, Nissim Hadar
<nhadar at surgicaltheater.net> wrote:
> 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.
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users