[Insight-users] ITK FastMarching example--help!!!--the resultion of output images are changed to [1 1 1]

[Baoyun Li]

Dear All:

When I run ITK example for fastmarching (see the attached code below), the output image resolution are changed to [1 1 1] begining at sigmoid filter output [writer3->SetFileName("../data/FastMarchingFilterOutput3.hdr")]. But my input resolution is [0.7 0.7 2.5]. The first two output files kept the origianl resolution.

I only changed some arguments of  the expmale code and the Dimension to 3. My command aruguments are: FastMarchingImageFilter SE1.hdr out.hdr 1000 1001.

ITK deals images on Physical space, now the spacial resolution of images are changed, so I do not how I can trust the results and cannot continue some works, eg mutiple resoltion processing

Can somebody help me to figure out this problem?  Thanks a lot.


Best regards

Baoyun


#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif


#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkImage.h"
#include "itkFastMarchingImageFilter.h"

// Software Guide : BeginCodeSnippet
#include "itkBinaryThresholdImageFilter.h"

// Software Guide : BeginCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
// 
#include "itkRescaleIntensityImageFilter.h"

int main( int argc, char *argv[] )
{
  if( argc < 3 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " inputImage  outputImage seedX seedY";
    std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingValue" << std::endl;
    return 1;
    }

  
  // Software Guide : BeginCodeSnippet
  typedef   float           InternalPixelType;
  const     unsigned int    Dimension = 3;
  typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
  // Software Guide : EndCodeSnippet
  typedef unsigned char OutputPixelType;
  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
  // Software Guide : EndCodeSnippet
  typedef itk::BinaryThresholdImageFilter< InternalImageType, 
                        OutputImageType    >    ThresholdingFilterType;
  ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
  
  const InternalPixelType  timeThreshold = atof( argv[3] );
  
  // Software Guide : BeginCodeSnippet
  thresholder->SetLowerThreshold(           0.0  );
  thresholder->SetUpperThreshold( timeThreshold  );
  thresholder->SetOutsideValue(  0  );
  thresholder->SetInsideValue(  255 );
  
  // Software Guide : BeginCodeSnippet
  typedef  itk::ImageFileReader< InternalImageType > ReaderType;
  typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;
  // Software Guide : EndCodeSnippet

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();
  reader->SetFileName( argv[1] );
  writer->SetFileName( argv[2] );

  //  The RescaleIntensityImageFilter type is declared below. This filter will
  //  renormalize image before sending them to writers.
  //
  typedef itk::RescaleIntensityImageFilter< 
                               InternalImageType, 
                               OutputImageType >   CastFilterType;
  typedef   itk::CurvatureAnisotropicDiffusionImageFilter< 
                               InternalImageType, 
                               InternalImageType >  SmoothingFilterType;
  
  // Software Guide : BeginCodeSnippet
  SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
 
  // Software Guide : BeginCodeSnippet
  typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter< 
                               InternalImageType, 
                               InternalImageType >  GradientFilterType;
  typedef   itk::SigmoidImageFilter<                               
                               InternalImageType, 
                               InternalImageType >  SigmoidFilterType;
  
  GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();
  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginCodeSnippet
  sigmoid->SetOutputMinimum(  0.0  );
  sigmoid->SetOutputMaximum(  1.0  );
  
  // Software Guide : BeginCodeSnippet
  typedef  itk::FastMarchingImageFilter< InternalImageType, 
                              InternalImageType >    FastMarchingFilterType;
  
  // Software Guide : BeginCodeSnippet
  FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
 
  // Software Guide : BeginCodeSnippet
  smoothing->SetInput( reader->GetOutput() );
  gradientMagnitude->SetInput( smoothing->GetOutput() );
  sigmoid->SetInput( gradientMagnitude->GetOutput() );
  fastMarching->SetInput( sigmoid->GetOutput() );
  thresholder->SetInput( fastMarching->GetOutput() );
  writer->SetInput( thresholder->GetOutput() );
  // Software Guide : EndCodeSnippet

  
  // Software Guide : BeginCodeSnippet
  smoothing->SetTimeStep( 0.05 );
  smoothing->SetNumberOfIterations(  10 );
  smoothing->SetConductanceParameter( 9..0 );
  // Software Guide : EndCodeSnippet

  
  const double sigma = (float) 2.5;//atof( argv[5] );
  // Software Guide : BeginCodeSnippet
  gradientMagnitude->SetSigma(  sigma  );
  // Software Guide : EndCodeSnippet

  const double alpha =  (float)-0.33;//atof( argv[6] );
  const double beta  =  (float) 0.5;//atof( argv[7] );

  // Software Guide : BeginCodeSnippet
  sigmoid->SetAlpha( alpha );
  sigmoid->SetBeta(  beta  );
 
 
  typedef FastMarchingFilterType::NodeContainer           NodeContainer;
  typedef FastMarchingFilterType::NodeType                NodeType;
  NodeContainer::Pointer seeds = NodeContainer::New();
  //  Software Guide : EndCodeSnippet 
  
  InternalImageType::IndexType  seedPosition;
  
  seedPosition[0] = 131;//atoi( argv[3] );
  seedPosition[1] = 273;//atoi( argv[4] );
  seedPosition[2] = 36;//atoi( argv[4] );

  
  NodeType node;
  const double seedValue = 0.0;
  
  node.SetValue( seedValue );
  node.SetIndex( seedPosition );
 
  seeds->Initialize();
  seeds->InsertElement( 0, node );
  //  Software Guide : EndCodeSnippet 
  // Software Guide : BeginCodeSnippet
  fastMarching->SetTrialPoints(  seeds  );
  // Software Guide : EndCodeSnippet

  //  Here we configure all the writers required to see the intermediate
  //  outputs of the pipeline. This is added here only for
  //  pedagogical/debugging purposes. These intermediate output are normaly not
  //  required. Only the output of the final thresholding filter should be
  //  relevant.  Observing intermediate output is helpful in the process of
  //  fine tuning the parameters of filters in the pipeline. 
  //
  CastFilterType::Pointer caster1 = CastFilterType::New();
  CastFilterType::Pointer caster2 = CastFilterType::New();
  CastFilterType::Pointer caster3 = CastFilterType::New();
  CastFilterType::Pointer caster4 = CastFilterType::New();
  WriterType::Pointer writer1 = WriterType::New();
  WriterType::Pointer writer2 = WriterType::New();
  WriterType::Pointer writer3 = WriterType::New();
  WriterType::Pointer writer4 = WriterType::New();
  caster1->SetInput( smoothing->GetOutput() );
  writer1->SetInput( caster1->GetOutput() );
  writer1->SetFileName("../data/FastMarchingFilterOutput1.hdr");
  caster1->SetOutputMinimum(   0 );
  caster1->SetOutputMaximum( 255 );
  writer1->Update();
  caster2->SetInput( gradientMagnitude->GetOutput() );
  writer2->SetInput( caster2->GetOutput() );
  writer2->SetFileName("../data/FastMarchingFilterOutput2.hdr");
  caster2->SetOutputMinimum(   0 );
  caster2->SetOutputMaximum( 255 );
  writer2->Update();
  caster3->SetInput( sigmoid->GetOutput() );
  writer3->SetInput( caster3->GetOutput() );
  writer3->SetFileName("../data/FastMarchingFilterOutput3.hdr");
  caster3->SetOutputMinimum(   0 );
  caster3->SetOutputMaximum( 255 );
  writer3->Update();
  caster4->SetInput( fastMarching->GetOutput() );
  writer4->SetInput( caster4->GetOutput() );
  writer4->SetFileName("../data/FastMarchingFilterOutput4.hdr");
  caster4->SetOutputMinimum(   0 );
  caster4->SetOutputMaximum( 255 );

  //  Software Guide : BeginLatex
 
  //
  //  Software Guide : EndLatex 
  // Software Guide : BeginCodeSnippet
  fastMarching->SetOutputSize( 
           reader->GetOutput()->GetBufferedRegion().GetSize() );
  // Software Guide : EndCodeSnippet

  
  //
  //  Software Guide : EndLatex 
  const double stoppingTime = atof( argv[4] );
  // Software Guide : BeginCodeSnippet
  fastMarching->SetStoppingValue(  stoppingTime  );
  // Software Guide : EndCodeSnippet

 
  // Software Guide : BeginCodeSnippet
  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }
  // Software Guide : EndCodeSnippet

  writer4->Update();

  // The following writer type is used to save the output of the
  // time-crossing map in a file with appropiate pixel representation. The
  // advantage of saving this image in native format is that it can be used
  // with a viewer to help determine an appropriate threshold to be used on
  // the output of the \code{fastmarching} filter.
  //
  typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
  InternalWriterType::Pointer mapWriter = InternalWriterType::New();
  mapWriter->SetInput( fastMarching->GetOutput() );
  mapWriter->SetFileName("../data/FastMarchingFilterOutput5.hdr");
  mapWriter->Update();
  InternalWriterType::Pointer speedWriter = InternalWriterType::New();
  speedWriter->SetInput( sigmoid->GetOutput() );
  speedWriter->SetFileName("../data/FastMarchingFilterOutput6.hdr");
  speedWriter->Update();
  InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
  gradientWriter->SetInput( gradientMagnitude->GetOutput() );
  gradientWriter->SetFileName("../data/FastMarchingFilterOutput7.hdr");
  gradientWriter->Update();

  return 0;
}



      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090225/9135bf99/attachment-0001.htm>