[Insight-users] ITK and VTK

Pietro Nardelli pie.nardelli at gmail.com
Tue May 29 09:42:51 EDT 2012


Hi Luis,

thank you for your help! I solved that error setting the output of the
region of interest filter egual to the input of the connected
threshold class. But now I have another problem: my method crop my
volume of interest but it doesn't segment anything. I tried different
fashions, but nothing changes. This is my code (I point out that I'm
writing this method for Slicer):

/*=========================================================================

  Program:   Slicer
  Language:  C++
  Module:    $HeadURL$
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) Brigham and Women's Hospital (BWH) All Rights Reserved.

  See License.txt or http://www.slicer.org/copyright/copyright.txt for details.

==========================================================================*/

#include "itkRegionOfInterestImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"
#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkPluginFilterWatcher.h"
#include "itkPluginUtilities.h"

#include "LungSegmentationCLP.h"

// Use an anonymous namespace to keep class types and function names
// from colliding when module is used as shared object module.  Every
// thing should be in an anonymous namespace except for the module
// entry point, e.g. main()
//
namespace
{

} // end of anonymous namespace

int main( int argc, char *argv[] )
{
  PARSE_ARGS;

  typedef signed short    InputPixelType;
  typedef signed short    cropOutputPixelType;
  typedef unsigned short  labelPixelType;
  const   unsigned int    Dimension = 3;

  typedef itk::Image<InputPixelType, Dimension>        cropInputImageType;
  typedef itk::Image<cropOutputPixelType, Dimension>   cropOutputImageType;
  typedef itk::Image<labelPixelType, Dimension>        LabelImageType;

  typedef itk::CastImageFilter<cropOutputImageType, LabelImageType>
	  CastingFilterType;
  CastingFilterType::Pointer       caster = CastingFilterType::New();

  typedef  itk::ImageFileReader<cropInputImageType>    cropReaderType;
  typedef  itk::ImageFileWriter<cropOutputImageType>   cropOutputWriterType;

  typedef  itk::ImageFileReader<cropOutputImageType>   segmentInputReaderType;
  typedef  itk::ImageFileWriter<LabelImageType>        labelWriterType;

  cropReaderType::Pointer          cropReader =	        cropReaderType::New();
  cropOutputWriterType::Pointer    cropWriter =
cropOutputWriterType::New();
  segmentInputReaderType::Pointer  tracheaReader =
segmentInputReaderType::New();
  labelWriterType::Pointer         tracheaVolume =      labelWriterType::New();

  typedef itk::RegionOfInterestImageFilter< cropInputImageType,
cropOutputImageType >
	  FilterType;
  FilterType::Pointer               cropFilter =        FilterType::New();
  itk::PluginFilterWatcher filterWatcher(cropFilter, "Cropping the
Trachea Volume", CLPProcessInformation);

  cropReader->SetFileName( inputVolume.c_str() );
  cropReader->Update();

  //split into 3 parts along x
  cropInputImageType::SizeType       sizeOfImage =
cropReader->GetOutput()->GetLargestPossibleRegion().GetSize();
  int XVoxelsNumber = sizeOfImage[0];
  int XNewDimension;
  XNewDimension = XVoxelsNumber / 3;

  cropOutputImageType::SizeType tracheaSize;
  tracheaSize[0] = XNewDimension - 100;//length of the region
  tracheaSize[1] = sizeOfImage[1];//width of the region
  tracheaSize[2] = sizeOfImage[2];//thick of the region

  cropOutputImageType::IndexType tracheaIndexPixel;
  tracheaIndexPixel[0] = XNewDimension + 51;//x position of starting pixel
  tracheaIndexPixel[1] = 0;//y position of starting pixel
  tracheaIndexPixel[2] = 0;//z position of starting pixel

  cropOutputImageType::RegionType tracheaDesiredRegion;
  tracheaDesiredRegion.SetSize(  tracheaSize  );
  tracheaDesiredRegion.SetIndex( tracheaIndexPixel );

  cropFilter->SetRegionOfInterest( tracheaDesiredRegion );

  typedef itk::ConnectedThresholdImageFilter<cropOutputImageType,
cropOutputImageType>
	  ConnectedFilterType;
  ConnectedFilterType::Pointer thresholdConnected = ConnectedFilterType::New();
  ConnectedFilterType::ConnectivityEnumType neighborhood =
ConnectedFilterType::FullConnectivity;
  itk::PluginFilterWatcher     watcher(thresholdConnected, "Segmenting",
                                        CLPProcessInformation, 0.5, 0.5);

  tracheaVolume->SetFileName( tracheaLabel.c_str() );

  //cropWriter->SetFileName( newVolume.c_str() );
  cropFilter->SetInput( cropReader->GetOutput() );
  //cropWriter->SetInput( cropFilter->GetOutput() );
  //cropWriter->SetUseCompression(1);

  //tracheaReader->SetFileName( cropWriter->GetFileName() );
  thresholdConnected->SetInput( cropFilter->GetOutput() );
  caster->SetInput( thresholdConnected->GetOutput() );
  tracheaVolume->SetInput( caster->GetOutput() );
  tracheaVolume->SetUseCompression(1);

  cropOutputPixelType lowThreshold = -1300;
  thresholdConnected->SetConnectivity(neighborhood);
  thresholdConnected->SetLower( lowThreshold );
  thresholdConnected->SetUpper( upThreshold );
  thresholdConnected->SetReplaceValue( labelvalue );


  if( seed.size() > 0 )
    {
    cropOutputImageType::PointType lpsPoint;
    cropOutputImageType::IndexType index;
    for( ::size_t i = 0; i < seed.size(); ++i )
      {
      // seeds come in ras, convert to lps
      lpsPoint[0] = -seed[i][0];
      lpsPoint[1] = -seed[i][1];
      lpsPoint[2] = seed[i][2];

      tracheaReader->GetOutput()->TransformPhysicalPointToIndex(lpsPoint,
index);
      thresholdConnected->AddSeed(index);

//       std::cout << "LPS: " << lpsPoint << std::endl;
//       std::cout << "IJK: " << index << std::endl;
      }
    }
  else
    {
    std::cerr << "No seeds specified." << std::endl;
    return -1;
    }

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

  return 0;
}


Thank you again,

Pietro

2012/5/28 Luis Ibanez <luis.ibanez at kitware.com>:
> Hi Pietro,
>
> Could you please post to the list the line of code
> (and even better, include the 10 lines above and
> below this line) so we can see it in context ?
>
> There is something suspicious in your message,
> since, it should be trivial to cast a non-const pointer
> to a const pointer.  Please make sure that you copy
> paste the error message in your email exactly as
> you get it from the compiler.
>
>
>   Thanks
>
>
>      Luis
>
>
> -------------------------------------
>
> On Thu, May 24, 2012 at 7:53 AM, Pietro Nardelli <pie.nardelli at gmail.com>
> wrote:
>>
>> Sorry Luis,
>>
>> I have another quetion: I'm trying to create a pipeline which uses a
>> cropped volume (obtained by using the RegionOfInterestImageFilter) as
>> input of a connected threshold filter. Building it, though, I get an
>> error: it says that it cannot convert parameter 1 from
>> 'itk::Image<TPixel,VImageDimension> *' to 'const
>> itk::Image<TPixel,VImageDimension> *'.
>>
>> What might it be? Any Idea?
>>
>> Thank you again
>>
>> Pietro
>>
>> 2012/5/24 Luis Ibanez <luis.ibanez at kitware.com>:
>> > HI Pietro,
>> >
>> > It depends where in your application you are extracting the region.
>> >
>> > If you are doing it in the code that deals with visualization and the
>> > resulting region is only going to be used for display, then the VTK
>> > filter will do just fine.
>> >
>> > On the other hand, if you are doing this in the middle of a pipeline
>> > and the extracted region is going to be used for subsequent
>> > processing, you probably should prefer the ITK Region of Interest
>> > filter:
>> >
>> >
>> > http://www.itk.org/Doxygen/html/classitk_1_1RegionOfInterestImageFilter.html
>> >
>> > since it preserves the orientation (direction cosines) of the image,
>> > and properly sets the origin in physical coordinates of the output
>> > image region.  This are very important attributes if you are using
>> > the resulting image in registration, or if you plan to correlate it
>> > spatially to the input image.
>> >
>> >     Hope this helps,
>> >
>> >          Luis
>> >
>> >
>> > ------------------------------------
>> > On Wed, May 23, 2012 at 12:08 PM, Pietro Nardelli
>> > <pie.nardelli at gmail.com>
>> > wrote:
>> >>
>> >> Sorry guys,
>> >>
>> >> I have another question: I'm using the class Region of interest to
>> >> crop a volume of a DICOM dataset, but I found out that in VTK there is
>> >> the class Extract ROI which seems very similar. Does anyone know what
>> >> the difference between the 2 classes? Which is the best one?
>> >>
>> >> Thank you,
>> >>
>> >> Pietro
>> >> _____________________________________
>> >> 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://www.itk.org/mailman/listinfo/insight-users
>> >
>> >
>
>


More information about the Insight-users mailing list