[Insight-users] combining output from segmentation with label map attributes using Gaetan Lehmann's paper: "Label object representation and manipulation with ITK"

Tue Dec 8 14:38:16 EST 2009

Hello,

I have segmented the ventricles from an image using a connected threshold
image filter.
I then tried to convert my outputted segmentation to a label map filter and
then to get the label's attributes.

I am using parts of the example code attribute_values.cxx that is described
in section 9.4 of
Gaetan Lehmann's paper: "Label object representation and manipulation with
ITK"

My problem is that my code is telling me that the number of labels is 0.
How do I apply a label to my segmented ventricles?

Below is my code:

/*
to run type:
./ConnectedThresholdImageFilter correctedsubject5.dcm outsubject5.dcm 103
142 95 17100 17300
*/

#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif

#include "itkConnectedThresholdImageFilter.h"

#include "itkImage.h"
#include "itkCastImageFilter.h"

#include "itkCurvatureFlowImageFilter.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkGDCMImageIO.h"

#include "itkVersion.h"

#include "itkOrientedImage.h"
#include "itkMinimumMaximumImageFilter.h"

#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkNumericSeriesFileNames.h"

#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"

#include "itkResampleImageFilter.h"
#include "itkShiftScaleImageFilter.h"

#include "itkIdentityTransform.h"
#include "itkLinearInterpolateImageFunction.h"

#include <itksys/SystemTools.hxx>

#include "gdcm/src/gdcmFile.h"
#include "gdcm/src/gdcmUtil.h"

#include <string>

//added per attribute_values.cxx
#include "itkImageFileReader.h"

#include "itkShapeLabelObject.h"
//#include "itkStatisticsLabelObject.h"

#include "itkLabelMap.h"

#include "itkBinaryImageToShapeLabelMapFilter.h"
//#include "itkBinaryImageToStatisticsLabelMapFilter.h"

//end of added code

int main( int argc, char *argv[])
{
  if( argc < 7 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " inputImage  outputImage seedX seedY seedZ lowerThreshold
upperThreshold" << std::endl;

    return 1;
    }

  typedef   float           InternalPixelType;

  const     unsigned int    Dimension = 3;

  typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;

  typedef signed short OutputPixelType;

  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
  typedef itk::Image< float, Dimension > OutputImageType2;
  typedef itk::CastImageFilter< InternalImageType, OutputImageType >
    CastingFilterType;
  CastingFilterType::Pointer caster = CastingFilterType::New();

  const    unsigned int    ImageDimension = 3;
  typedef  signed short    PixelType;

  typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
  typedef itk::Image< float, ImageDimension >  FloatImageType;

  typedef  itk::ImageFileReader< FixedImageType > ReaderType;
  typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
  typedef  itk::ImageFileWriter<  FloatImageType  > WriterType2;

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();
  WriterType2::Pointer writer2 = WriterType2::New();

  typedef itk::GDCMImageIO           ImageIOTypefixed;
  ImageIOTypefixed::Pointer gdcmImageIOfixed = ImageIOTypefixed::New();
  reader->SetImageIO( gdcmImageIOfixed );

  typedef itk::GDCMImageIO           ImageIOTypefixed2;
  ImageIOTypefixed2::Pointer gdcmImageIOfixed2 = ImageIOTypefixed2::New();

  reader->SetFileName( argv[1] );

  reader->Update();

  typedef itk::CurvatureFlowImageFilter< InternalImageType,
InternalImageType >
    CurvatureFlowImageFilterType;

  CurvatureFlowImageFilterType::Pointer smoothing =
                         CurvatureFlowImageFilterType::New();

  typedef itk::ConnectedThresholdImageFilter< InternalImageType,
                                    InternalImageType > ConnectedFilterType;

  ConnectedFilterType::Pointer connectedThreshold =
ConnectedFilterType::New();

  typedef signed short InputAPixelType;
  typedef float OutputBPixelType;

  typedef itk::Image< InputAPixelType, 3 > InputAImageType;
  typedef itk::Image< OutputBPixelType, 3 > OutputBImageType;

  typedef itk::CastImageFilter< InputAImageType, OutputBImageType >
CastFilterType;

  CastFilterType::Pointer castFilter = CastFilterType::New();

  castFilter->SetInput( reader->GetOutput() );

  connectedThreshold->SetInput( castFilter->GetOutput() );

  caster->SetInput( connectedThreshold->GetOutput() );

  smoothing->SetNumberOfIterations( 20 ); //was 5
  smoothing->SetTimeStep( 0.125 );

  const InternalPixelType lowerThreshold = atof( argv[6] );
  const InternalPixelType upperThreshold = atof( argv[7] );

  connectedThreshold->SetLower(  lowerThreshold  );
  connectedThreshold->SetUpper(  upperThreshold  );

  connectedThreshold->SetReplaceValue( 255 );

  InternalImageType::IndexType  index;

  index[0] = atoi( argv[3] );
  index[1] = atoi( argv[4] );

  //added
  index[2] = atoi( argv[5] );

  std::cout << index << std::endl;

  // Software Guide : BeginCodeSnippet
  connectedThreshold->SetSeed( index );

  //obtain a 5 x 5 bounding region of seeds
  int ii, jj, kk;

  ii = index[0];
  jj = index[1];
  kk = index[2];

  for (int i = ii; i < ii + 5; i++)
    for (int j = jj; j < jj + 5; j++)
      for (int k = kk; k < kk + 5; k++)
    {

      index[0] = i;
      index[1] = j;
      index[2] = k;
      connectedThreshold->AddSeed( index );
    }

  for (int i = ii; i > ii - 5; i--)
    for (int j = jj; j > jj - 5; j--)
      for (int k = kk; k > kk - 5; k--)
    {

      index[0] = i;
      index[1] = j;
      index[2] = k;
      connectedThreshold->AddSeed( index );
    }

  connectedThreshold->Print(std::cout,17100);

  typedef itk::MetaDataDictionary DictionaryType;

  DictionaryType inputdict = reader->GetMetaDataDictionary();

  writer->SetMetaDataDictionary( inputdict );

  writer->SetFileName( argv[2] );

  //added per attribute_values.cxx

  // define the object type. Here the ShapeLabelObject type
  // is chosen in order to read some attribute related to the shape
  // of the objects (by opposition to the content of the object, with
  // the StatisticsLabelObejct).
  typedef unsigned long LabelType;

  typedef itk::ShapeLabelObject< LabelType, 3 > LabelObjectType;

  typedef itk::LabelMap< LabelObjectType > LabelMapType;

  // convert the image in a collection of objects
  typedef itk::BinaryImageToLabelMapFilter< InternalImageType, LabelMapType
> ConverterType;

  ConverterType::Pointer converter = ConverterType::New();

  converter->SetInput( connectedThreshold->GetOutput() );

  //converter->SetForegroundValue( atoi(argv[2]) );
  //converter->SetForegroundValue( 255 );

  // valuate the attributes with the dedicated filter: ShapeLabelMapFilter

  typedef itk::ShapeLabelMapFilter< LabelMapType > ShapeFilterType;
  ShapeFilterType::Pointer shape = ShapeFilterType::New();
  shape->SetInput( converter->GetOutput() );

  shape->Update();

  // then we can read the attribute values we're interested in. The
BinaryImageToShapeLabelMapFilter
  // produce consecutive labels, so we can use a for loop and
GetLabelObject() method to retrieve
  // the label objects. If the labels are not consecutive, the
GetNthLabelObject() method must be
  // use instead of GetLabelObject(), or an iterator on the label object
container of the label map.
  LabelMapType::Pointer labelMap = shape->GetOutput();

  std::cout << "Number of label objects = " <<
labelMap->GetNumberOfLabelObjects() << std::endl;

  for( unsigned int label=1; label<=labelMap->GetNumberOfLabelObjects();
label++ )
    {
    // we don't need a SmartPointer of the label object here, because the
reference is kept in
    // in the label map.
    const LabelObjectType * labelObject = labelMap->GetLabelObject( label );
    std::cout << label << "\t" << labelObject->GetPhysicalSize() << "\t" <<
labelObject->GetCentroid() << std::endl;
    }

  writer->SetInput( caster->GetOutput() );

  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }

  std::cout << "output from reader->GetOutput()->GetDirection()" <<
std::endl;
  std::cout << reader->GetOutput()->GetDirection() << std::endl;

  std::cout << "output from castFilter->GetOutput()->GetDirection()" <<
std::endl;
  std::cout << castFilter->GetOutput()->GetDirection() << std::endl;

  std::cout << "output from connectedThreshold->GetOutput()->GetDirection()"
<< std::endl;
  std::cout << connectedThreshold->GetOutput()->GetDirection() << std::endl;

  std::cout << "output from caster->GetOutput()->GetDirection()" <<
std::endl;
  std::cout << caster->GetOutput()->GetDirection() << std::endl;

  return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091208/e54fb541/attachment-0001.htm>