[Insight-users] Re-sort a 3D Volume

Gaëtan Lehmann gaetan.lehmann at jouy.inra.fr
Wed May 16 17:19:08 EDT 2007


Hi all,

Just a comment for the future users (a with the hope to get some  
reviews): the SliceBySliceImageFilter should be able to do what is  
done here, but in a much simple and lisible way.
That code should be easier to modify to work with uncommon cases though.

Gaëtan



Le 16 mai 07 à 18:18, d f a écrit :

> Hi all,
>
> I recently needed something just like what Roger was nice enough to  
> share last month (http://public.kitware.com/pipermail/insight-users/ 
> 2007-April/021796.html ).  I modified it a bit for my needs, and  
> decided to send this along in case it helps anyone else (although  
> I'm not sure how often other people would find the additional  
> functionality useful).
>
> This just uses a basic Gaussian filter for the processing, but I  
> was also able to implement a full level set segmentation, operating  
> on 2D slices at a time, by using this as a template.  (In case you  
> are curious, I was comparing some tradeoffs of implementing 3D vs  
> 2D segmentation on 3D images)
>
> drf
>
>
> //  This programs reads in a 3D image volume, extracts 2D slices,
> //  and performs some 2D filtering on each. The results are inserted
> //  into a 3D output volume of identical dimensions to the input. This
> //  is used to perform 2D processing on 3D images without having to
> //  write large numbers of 2D images to disk,  simplifying the batch
> //  files used to perform processing.
> //
> //  Any 2D processing pipeline can be used here; the object that
> //  passes its output to the TileImageFilter must be updated inside
> //  the loop (lines 167,168).  In this case, a discrete gaussian image
> //  filter is applied to the 2D images before being passed to the  
> tiler.
> //
> //  Additionally, 2D slices along the x,y, or z dimensions can be
> //  operated on.  The input parameter is collapseDim (0,1,or 2).
> //  For example, a value of 2 "collapses" (or steps through) the
> //  z-dimension, yielding many 2D images from the x-y plane.   
> Likewise,
> //  collapseDim=0 collapses x, yielding images from the y-z plane.   
> The
> //  result is permuted so the output has the same dimensions as the  
> original.
> //
> //  Thanks to Luis Ibanez for the help and to Roger Bramon Feixas,
> //  who actually coded the loop that puts the images into the tiler:
> //  http://public.kitware.com/pipermail/insight-users/2007-April/ 
> 021796.html
>
>
> #if defined(_MSC_VER)
> #pragma warning ( disable : 4786 )
> #endif
>
> #ifdef __BORLANDC__
> #define ITK_LEAN_AND_MEAN
> #endif
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkExtractImageFilter.h"
> #include "itkTileImageFilter.h"
> #include "itkDiscreteGaussianImageFilter.h "
> #include "itkRescaleIntensityImageFilter.h"
> #include "itkImage.h"
> #include "itkPermuteAxesImageFilter.h"
>
>
> int main( int argc, char ** argv )
> {
>   //  Check input arguments
>   if( argc < 3 )
>     {
>     std::cerr << "Usage: " << std::endl;
>     std::cerr << argv[0] << " input3DImageFile  output3DImageFile  
> collapseDim" << std::endl;
>     return EXIT_FAILURE;
>     }
>
>
>   //  Typedefs
>   typedef float         InputPixelType;
>   typedef float         MiddlePixelType;
>   typedef float         OutputPixelType;
>   typedef unsigned char WritePixelType;
>
>   typedef itk::Image< InputPixelType,  3 >    InputImageType;
>   typedef itk::Image< MiddlePixelType, 2 >    MiddleImageType;
>   typedef itk::Image< OutputPixelType, 3 >    OutputImageType;
>   typedef itk::Image< WritePixelType,  3 >    WriteImageType;
>
>   typedef itk::ImageFileReader< InputImageType  >  ReaderType;
>   typedef itk::ImageFileWriter< WriteImageType  >  WriterType;
>
>   typedef itk::ExtractImageFilter<
>                InputImageType, MiddleImageType  >  ExtractFilterType;
>   typedef itk::DiscreteGaussianImageFilter<
>                MiddleImageType, MiddleImageType >  FilterType;
>   typedef itk::TileImageFilter<
>                MiddleImageType, OutputImageType >  TileFilterType;
>   typedef itk::RescaleIntensityImageFilter<
>                OutputImageType, WriteImageType  >  RescaleFilterType;
>   typedef itk::PermuteAxesImageFilter<
>                OutputImageType                  >  PermuteFilterType;
>
>
>   //  Input/Output filenames
>   const char * inputFilename  = argv[1];
>   const char * outputFilename = argv[2];
>
>
>   // Reader/Writer instantiations
>   ReaderType::Pointer reader = ReaderType::New();
>   WriterType::Pointer writer = WriterType::New();
>   reader->SetFileName( inputFilename  );
>   writer->SetFileName( outputFilename );
>
>
>   //  Read input file
>   try
>     {
>     reader->Update();
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cerr << "ExceptionObject caught !" << std::endl;
>     std::cerr << err << std::endl;
>     return EXIT_FAILURE;
>     }
>
>
>   //  Find input image size
>   InputImageType::RegionType inputRegion =
>            reader->GetOutput()->GetLargestPossibleRegion();
>   InputImageType::SizeType size = inputRegion.GetSize();
>   const unsigned int collapseDim  = atoi( argv[3] );
>   InputImageType::SizeType::SizeValueType numSlices = size 
> [collapseDim];
>   size[collapseDim] = 0;   // collapse dimension ( 3D->2D )
>
>
>   //  Set region for 2D extraction
>   InputImageType::IndexType start = inputRegion.GetIndex();
>   InputImageType::RegionType desiredRegion;
>   desiredRegion.SetSize(  size  );
>
>
>   //  Extract filter instantiation
>   ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
>
>
>   //  2D Gaussian filter instantiation
>   FilterType::Pointer filter = FilterType::New();
>   const double gaussianVariance = 1;
>   const unsigned int maxKernelWidth = 20;
>   filter->SetVariance( gaussianVariance );
>   filter->SetMaximumKernelWidth( maxKernelWidth );
>
>
>   //  Tiling filter for inserting 2D results into 3D output volume
>   //  layout=[1,1,0] is code to automatically create necessary  
> image size
>   //     for a given number of slices
>   TileFilterType::Pointer tileFilter = TileFilterType::New();
>   TileFilterType::LayoutArrayType layout;
>   layout[0] = 1;
>   layout[1] = 1;
>   layout[2] = 0;
>   tileFilter->SetLayout( layout );
>
>
>   //  Rescale used to cast float->char for image writing
>   RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
>   rescaler->SetOutputMinimum(   0 );
>   rescaler->SetOutputMaximum( 255 );
>
>
>   //  Set up pipeline
>   extractFilter->SetInput( reader->GetOutput() );
>   filter->SetInput( extractFilter->GetOutput() );
>   tileFilter->SetInput( filter->GetOutput() );
>
>
>   //  Vector of pointers to the many 2D extracted images
>   std::vector< MiddleImageType::Pointer > extracts;
>
>
>   //  Extract 2D slices, perform processing
>   for(unsigned int sliceNumber = 0; sliceNumber<numSlices;  
> sliceNumber++)
>   {
>       start[collapseDim] = sliceNumber;
>       desiredRegion.SetIndex( start );
>
>       extractFilter->SetExtractionRegion( desiredRegion );
>
>       filter->Update();
>       extracts.push_back( filter->GetOutput() );
>       extracts.back()->DisconnectPipeline();
>
>       if (sliceNumber != 0)
>       {
>           tileFilter->PushBackInput( extracts.back() );
>       }
>   }
>
>   //  Update tiling filter
>   tileFilter->Update();
>
>
>   //  Permute image axes to match input
>   PermuteFilterType::Pointer permute = PermuteFilterType::New();
>   unsigned int permuteAxes[3];
>   switch (collapseDim)
>   {
>     case 0:
>         {
>         permuteAxes[0] = 2;
>         permuteAxes[1] = 0;
>         permuteAxes[2] = 1;
>         break;
>         }
>     case 1:
>         {
>         permuteAxes[0] = 0;
>         permuteAxes[1] = 2;
>         permuteAxes[2] = 1;
>         break;
>         }
>     case 2:
>         {
>         permuteAxes[0] = 0;
>         permuteAxes[1] = 1;
>         permuteAxes[2] = 2;
>         break;
>         }
>     default:
>         {
>         permuteAxes[0] = 0;
>         permuteAxes[1] = 1;
>         permuteAxes[2] = 2;
>         break;
>         }
>   }
>   permute->SetOrder( permuteAxes );
>   permute->SetInput( tileFilter->GetOutput() );
>
>
>   //  Remainder of pipeline
>   //  Cast float->char for image writing
>   rescaler->SetInput( permute->GetOutput() );
>   writer->SetInput( rescaler->GetOutput() );
>
>
>   //  Update pipeline, write file
>   try
>     {
>     writer->Update();
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cerr << "ExceptionObject caught !" << std::endl;
>     std::cerr << err << std::endl;
>     return EXIT_FAILURE;
>     }
>
>   return EXIT_SUCCESS;
> }
>
>
>
>
> On 4/9/07, Roger Bramon Feixas <rogerbramon at gmail.com > wrote:
> > Hi Luis,
> >
> > Problem solved. Now, I would like share with us the most important
> > parts of  my algorithm, if somebody has a similar problem:
> >
> > Befor the loop:
> > ----------------------
> >
> > ExtractImageType::Pointer extractFilter = ExtractImageType::New();
> > extractFilter->SetInput( myItkImageType3D );
> >
> > TileFilterType::Pointer tileFilter = TileFilterType::New();
> >
> > TileFilterType::LayoutArrayType layout;
> > layout[0] = 1;
> > layout[1] = 1;
> > layout[2] = 0;
> >
> > tileFilter->SetLayout( layout );
> >
> > std::vector< ItkImageType2D::Pointer > extracts;
> >
> > Into the loop:
> > ------------------
> >
> > extractFilter->SetExtractionRegion( region );
> > extractFilter->Update();
> > extracts.push_back( extractFilter->GetOutput() );
> > extracts.back()->DisconnectPipeline();
> > tileFilter->PushBackInput( extracts.back() );
> >
> > After the loop:
> > -------------------
> > tileFilter->Update();
> >
> >
> >
> > Thanks for all,
> >
> > Roger
> >
> >
> > El 08/04/2007, a las 0:33, Luis Ibanez escribió:
> >
> > >
> > > Hi Roger,
> > >
> > > You will need to run the ExtractImageFilter first in a loop,
> > > and at each iteration extract a different slice.
> > >
> > > You should put the output of the ExtractImageFilter into a
> > > std::vector of SmartPointers, and after pushing it in the
> > > vector you *MUST* call DisconnectPipeline() on that image.
> > >
> > > Otherwise, all your "slice" images are using the *same  
> SmartPointer*.
> > > That is, the ExtractFilter output is reusing the same instance  
> of the
> > > image object every time.
> > >
> > > I'm guessing that what you see as output of the TileImageFilter
> > > is that all the tiles contain a replication of the "Last" slice
> > > that you extracted with the ExtractImageFilter.
> > >
> > >
> > > Adding the call to DisconnectPipeline() will force the Extract
> > > ImageFilter to create a new instance of the image object for
> > > storing its output at every iteration.
> > >
> > >
> > >    Regards,
> > >
> > >
> > >       Luis
> > >
> > >
> > > ---------------------------------
> > > Roger Bramon Feixas wrote:
> > >> Hi,
> > >> I need to re-order the slices of a 3D volume. I tried to extract
> > >> 2D  slices using a "itkExtractImageFilter" and create a new  
> volume
> > >> using  "itkTileImageFilter" but it doesn't work. The extraction
> > >> works good  because I write the output in an image and it's
> > >> correct. The problem  is in "itkTileImageFilter". I receive data
> > >> empty error when I want to  use tileImageFilter output. That  
> is my
> > >> algorism:
> > >> before the loop:
> > >>  - I instantiate the TileImageFilter
> > >>  - I set the layout in 1,1,0
> > >> into the loop over the images:
> > >>  - I allocated a different ExtractImageFilter for each image
> > >>  - I fed the TileImageFilter with the output of each different
> > >> ExtractImageFilter: tileFilter->PushBackInput( extractFilter-
> > >> >GetOutput() ). I also tried to do that using: tileFilter-
> > >> >SetInput ( i, extractFilter->GetOutput() )
> > >> after the loop:
> > >>  TileImageFilter->Update();
> > >> After the loop, the layout is : 1,1,X where X is the number of
> > >> slices  that I pushed. Also, I read in a insight-user mail that
> > >> they solved a  similar problem saving all of ExtractimageFilter
> > >> objects in a vector  but in my way it doesn't work.
> > >> I hope anybody know anything about my problem.
> > >> Thanks,
> > >> Roger
> > >> _______________________________________________
> > >> Insight-users mailing list
> > >> Insight-users at itk.org
> > >> http://www.itk.org/mailman/listinfo/insight-users
> >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users at itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
> >
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users

--
Gaëtan Lehmann
Biologie du Développement et de la Reproduction
INRA de Jouy-en-Josas (France)
tel: +33 1 34 65 29 66    fax: 01 34 65 29 09
http://voxel.jouy.inra.fr



-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 186 bytes
Desc: =?ISO-8859-1?Q?Ceci_est_une_signature_=E9lectronique_PGP?=
Url : http://public.kitware.com/pipermail/insight-users/attachments/20070516/97dd8a37/PGP-0001.pgp


More information about the Insight-users mailing list