hi , <br>I am trying to make my stack of images isotropic without rescaling . <br>I used the example in the Examples folder called "Resamplevolumetobeisotropic.cxx" and and modified it a bit to remove the rescaling part, the resulting grey scales are visibly bad. I have added my code below;<br>
<br><a href="http://dl.dropbox.com/u/2989703/Screen%20shot%202011-06-29%20at%2011.03.18%20AM.png">http://dl.dropbox.com/u/2989703/Screen%20shot%202011-06-29%20at%2011.03.18%20AM.png</a><br><br>the image on the left is one slice of the original and one on the left is the output.<br>
<br>// 25-6-2011 <br>// wound micro enviornment :<br>// step1 : read and write an image <br><br>#include "itkImage.h"<br>#include "itkImageFileReader.h" <br>#include "itkImageFileWriter.h"<br>
#include "itkBinaryThresholdImageFilter.h"<br>#include "itkResampleImageFilter.h"<br><br><br>#include "itkRecursiveGaussianImageFilter.h"<br><br>using namespace std;<br>int main( int argc, char ** argv )<br>
{<br> if( argc < 3 )<br> {<br> std::cerr << "Usage: " << std::endl;<br> std::cerr << argv[0] << " inputImageFile outputImageFile " << std::endl;<br> return EXIT_FAILURE;<br>
}<br> <br> typedef unsigned short PixelType;<br> const unsigned int Dimension=3;<br> typedef itk::Image< PixelType, Dimension> ImageType;<br> typedef unsigned char OutputPixelType;<br><br> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;<br>
<br><br> typedef itk::ImageFileReader<ImageType> ReaderType;<br> typedef itk::ImageFileWriter<OutputImageType> WriterType;<br> <br> <br> ReaderType::Pointer reader = ReaderType::New();<br> WriterType::Pointer writer = WriterType::New();<br>
<br> <br> const char * inputFilename = argv[1];<br> const char * outputFilename = argv[2];<br> <br> reader->SetFileName( inputFilename );<br> try <br> {<br> cerr << "reading file "<< inputFilename <<endl;<br>
reader->Update();<br> }<br> catch( itk::ExceptionObject & excep )<br> {<br> std::cerr << "Exception caught!" << std::endl;<br> std::cerr << excep << std::endl;<br>
}<br> <br> typedef itk::RecursiveGaussianImageFilter< ImageType, ImageType > GaussianFilterType;<br> GaussianFilterType::Pointer smootherX = GaussianFilterType::New();<br> GaussianFilterType::Pointer smootherY = GaussianFilterType::New();<br>
<br> smootherX->SetInput( reader->GetOutput() );<br> smootherY->SetInput( smootherX->GetOutput() );<br> <br> ImageType::ConstPointer inputImage = reader->GetOutput();<br><br> const ImageType::SpacingType & inputSpacing = inputImage->GetSpacing();<br>
const double isoSpacing = vcl_sqrt( inputSpacing[2] * inputSpacing[0] );<br> <br> //debug <br> cerr <<" inputSpacing " << inputSpacing[0] << " + "<<<br> inputSpacing[1] << " + " << inputSpacing[2]<<endl;<br>
<br> cerr <<" isospacing "<< isoSpacing << endl;<br> <br> smootherX->SetSigma( isoSpacing );<br> smootherY->SetSigma( isoSpacing );<br> smootherX->SetDirection( 0 );<br>
smootherY->SetDirection( 1 );<br><br> smootherX->SetNormalizeAcrossScale( true );<br> smootherY->SetNormalizeAcrossScale( true );<br> <br> <br> typedef unsigned char OutputPixelType;<br><br>
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;<br> typedef itk::ResampleImageFilter< ImageType, OutputImageType > ResampleFilterType;<br><br> ResampleFilterType::Pointer resampler = ResampleFilterType::New();<br>
<br> typedef itk::IdentityTransform< double, Dimension > TransformType;<br><br> TransformType::Pointer transform = TransformType::New();<br> transform->SetIdentity();<br><br> resampler->SetTransform( transform );<br>
<br> typedef itk::LinearInterpolateImageFunction <ImageType, double > InterpolatorType;<br><br> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br><br> resampler->SetInterpolator( interpolator );<br>
<br> resampler->SetDefaultPixelValue( 4095 );<br> <br> OutputImageType::SpacingType spacing;<br><br> spacing[0] = isoSpacing;<br> spacing[1] = isoSpacing;<br> spacing[2] = isoSpacing;<br><br>
resampler->SetOutputSpacing( spacing );<br> resampler->SetOutputOrigin( inputImage->GetOrigin() );<br> resampler->SetOutputDirection( inputImage->GetDirection() );<br> <br> ImageType::SizeType inputSize = inputImage->GetLargestPossibleRegion().GetSize();<br>
<br> typedef ImageType::SizeType::SizeValueType SizeValueType;<br><br> const double dx = inputSize[0] * inputSpacing[0] / isoSpacing;<br> const double dy = inputSize[1] * inputSpacing[1] / isoSpacing;<br><br> const double dz = (inputSize[2] - 1 ) * inputSpacing[2] / isoSpacing;<br>
<br> ImageType::SizeType size;<br><br> size[0] = static_cast<SizeValueType>( dx );<br> size[1] = static_cast<SizeValueType>( dy );<br> size[2] = static_cast<SizeValueType>( dz );<br> resampler->SetSize( size );<br>
resampler->SetInput( smootherY->GetOutput() );<br><br> resampler->Update();<br> try <br> {<br> cerr << " resampler update " << endl;<br> resampler->Update();<br> }<br> catch (itk::ExceptionObject & err)<br>
{<br> cerr << "Exception "<<endl;<br> return EXIT_FAILURE; <br> }<br> <br> <br> writer->SetFileName( outputFilename );<br> <br> <br><br> writer->SetInput( resampler->GetOutput() );<br>
<br> try <br> {<br> cerr <<" the final write "<<endl;<br> writer->Update();<br> }<br> catch( itk::ExceptionObject & excep )<br> {<br> std::cerr << "Exception caught !" << std::endl;<br>
std::cerr << excep << std::endl;<br> }<br> return EXIT_SUCCESS;<br>}<br><br><br>