[Insight-users] Re-sort a 3D Volume
d f
drf8888 at gmail.com
Wed May 16 12:18:51 EDT 2007
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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070516/8287a342/attachment.html
More information about the Insight-users
mailing list