[Insight-users] making volume isotropic without rescaling

[iamrndm b]

hi ,
I am trying to make my stack of images isotropic without rescaling .
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;

http://dl.dropbox.com/u/2989703/Screen%20shot%202011-06-29%20at%2011.03.18%20AM.png

the image on the left is one slice of the original and one on the left is
the output.

// 25-6-2011
// wound micro enviornment :
// step1 : read and write an image

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkResampleImageFilter.h"


#include "itkRecursiveGaussianImageFilter.h"

using namespace std;
int main( int argc, char ** argv )
{
  if( argc < 3 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << " inputImageFile  outputImageFile " <<
std::endl;
    return EXIT_FAILURE;
    }

    typedef unsigned short PixelType;
    const unsigned int Dimension=3;
    typedef itk::Image< PixelType, Dimension> ImageType;
    typedef   unsigned char   OutputPixelType;

    typedef itk::Image< OutputPixelType,   Dimension >   OutputImageType;


    typedef itk::ImageFileReader<ImageType> ReaderType;
    typedef itk::ImageFileWriter<OutputImageType> WriterType;


    ReaderType::Pointer reader = ReaderType::New();
    WriterType::Pointer writer = WriterType::New();


    const char * inputFilename  = argv[1];
    const char * outputFilename = argv[2];

    reader->SetFileName( inputFilename );
      try
    {
    cerr << "reading file "<< inputFilename <<endl;
    reader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught!" << std::endl;
    std::cerr << excep << std::endl;
    }

    typedef itk::RecursiveGaussianImageFilter< ImageType, ImageType >
GaussianFilterType;
    GaussianFilterType::Pointer smootherX = GaussianFilterType::New();
    GaussianFilterType::Pointer smootherY = GaussianFilterType::New();

    smootherX->SetInput( reader->GetOutput() );
    smootherY->SetInput( smootherX->GetOutput() );

    ImageType::ConstPointer inputImage = reader->GetOutput();

    const ImageType::SpacingType & inputSpacing = inputImage->GetSpacing();
    const double isoSpacing = vcl_sqrt( inputSpacing[2] * inputSpacing[0] );

      //debug
      cerr <<" inputSpacing " << inputSpacing[0] << " + "<<
              inputSpacing[1] << " + " << inputSpacing[2]<<endl;

      cerr <<" isospacing "<< isoSpacing << endl;

    smootherX->SetSigma( isoSpacing );
    smootherY->SetSigma( isoSpacing );
    smootherX->SetDirection( 0 );
    smootherY->SetDirection( 1 );

    smootherX->SetNormalizeAcrossScale( true );
    smootherY->SetNormalizeAcrossScale( true );


    typedef   unsigned char   OutputPixelType;

    typedef itk::Image< OutputPixelType,   Dimension >   OutputImageType;
    typedef itk::ResampleImageFilter< ImageType, OutputImageType >
ResampleFilterType;

      ResampleFilterType::Pointer resampler = ResampleFilterType::New();

       typedef itk::IdentityTransform< double, Dimension >  TransformType;

     TransformType::Pointer transform = TransformType::New();
     transform->SetIdentity();

     resampler->SetTransform( transform );

     typedef itk::LinearInterpolateImageFunction <ImageType, double >
InterpolatorType;

     InterpolatorType::Pointer interpolator = InterpolatorType::New();

     resampler->SetInterpolator( interpolator );

     resampler->SetDefaultPixelValue( 4095 );

      OutputImageType::SpacingType spacing;

     spacing[0] = isoSpacing;
     spacing[1] = isoSpacing;
     spacing[2] = isoSpacing;

     resampler->SetOutputSpacing( spacing );
     resampler->SetOutputOrigin( inputImage->GetOrigin() );
     resampler->SetOutputDirection( inputImage->GetDirection() );

     ImageType::SizeType   inputSize =
inputImage->GetLargestPossibleRegion().GetSize();

     typedef ImageType::SizeType::SizeValueType SizeValueType;

  const double dx = inputSize[0] * inputSpacing[0] / isoSpacing;
  const double dy = inputSize[1] * inputSpacing[1] / isoSpacing;

  const double dz = (inputSize[2] - 1 ) * inputSpacing[2] / isoSpacing;

  ImageType::SizeType   size;

  size[0] = static_cast<SizeValueType>( dx );
  size[1] = static_cast<SizeValueType>( dy );
  size[2] = static_cast<SizeValueType>( dz );
  resampler->SetSize( size );
  resampler->SetInput( smootherY->GetOutput() );

  resampler->Update();
   try
    {
    cerr << " resampler update " << endl;
     resampler->Update();
    }
    catch (itk::ExceptionObject & err)
    {
           cerr << "Exception "<<endl;
           return EXIT_FAILURE;
       }


    writer->SetFileName( outputFilename );



  writer->SetInput( resampler->GetOutput() );

  try
    {
    cerr <<" the final write "<<endl;
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }
   return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110629/26dd4e85/attachment.htm>