[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