itkmribiasfieldcorrectionfiltertest.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMRIBiasFieldCorrectionFilterTest.cxx,v $
  Language:  C++
  Date:      $Date: 2007-12-29 13:36:00 $
  Version:   $Revision: 1.13 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include <iostream>
#include "vnl/vnl_vector.h"

#include "itkMRIBiasFieldCorrectionFilter.h"
#include "itkImage.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkGaussianImageSource.h"
#include "itkMultivariateLegendrePolynomial.h"
#include "itkCompositeValleyFunction.h"
#include "itkNormalVariateGenerator.h"
#include "itkArray.h"
#include "itkImageFileWriter.h"
#include "itkSphereSpatialFunction.h"

int itkMRIBiasFieldCorrectionFilterTest ( int , char* [] )
{
  const unsigned int ImageDimension = 3;

  typedef itk::Image< float, ImageDimension >            ImageType;
  typedef itk::Image< float, ImageDimension >            MaskType;
  typedef itk::ImageRegionIteratorWithIndex< ImageType > ImageIteratorType;

  bool SaveImages = false;
  ImageType::SizeType   imageSize;
  ImageType::IndexType  imageIndex;
  ImageType::RegionType imageRegion;
  imageSize[0] = 30;
  imageSize[1] = 30;
  imageSize[2] = 15;
  std::cout << "Random Test image size: " << imageSize[0] << "x" 
       << imageSize[1] << "x" << imageSize[2] << std::endl; 

  imageIndex.Fill( 0 );
  float spacing[ImageDimension] = {1.0, 1.0, 1.0};
  float origin[ImageDimension] = {0, 0, 0};

  imageRegion.SetSize( imageSize );
  imageRegion.SetIndex( imageIndex );

  // creates an image that will stores the Gaussian pixel * bias values
  ImageType::Pointer imageWithBias = ImageType::New();
  imageWithBias->SetBufferedRegion( imageRegion );
  imageWithBias->SetLargestPossibleRegion( imageRegion );
  imageWithBias->SetSpacing( spacing );
  imageWithBias->SetOrigin( origin );
  imageWithBias->Allocate();

  // creates the image source with a sphere.
  ImageType::Pointer image = ImageType::New();
  image->SetBufferedRegion( imageRegion );
  image->SetLargestPossibleRegion( imageRegion );
  image->SetSpacing( spacing );
  image->SetOrigin( origin );
  image->Allocate();

  // creates an image for bias
  ImageType::Pointer biasImage = ImageType::New();
  biasImage->SetBufferedRegion( imageRegion );
  biasImage->SetLargestPossibleRegion( imageRegion );
  biasImage->SetSpacing( spacing );
  biasImage->SetOrigin( origin );
  biasImage->Allocate();

  // class statistics for two classes: a bright sphere and background 
  itk::Array<double> classMeans(2);
  itk::Array<double> classSigmas(2);

  classMeans[0] = 10.0;
  classMeans[1] = 200.0;

  classSigmas[0] = 10.0;
  classSigmas[1] = 20.0;

  // creats a normal random variate generator
  itk::Statistics::NormalVariateGenerator::Pointer randomGenerator =
    itk::Statistics::NormalVariateGenerator::New();

  // fills the image with a sphere filled with intensity values from a
  // normal distribution. 
  typedef itk::SphereSpatialFunction<ImageDimension> SphereType;
  SphereType::Pointer sphere = SphereType::New();

  SphereType::InputType center;
  center[0] = imageSize[0]/2;
  center[1] = imageSize[1]/2;
  center[2] = imageSize[2]/2;
  sphere->SetCenter( center );
  sphere->SetRadius( 5.0 );

  randomGenerator->Initialize( 2003 );
  ImageIteratorType i_iter( image , imageRegion );
  SphereType::InputType point;
  while ( !i_iter.IsAtEnd() )
    {
    image->TransformIndexToPhysicalPoint( i_iter.GetIndex() , point );
    if ( sphere->Evaluate( point ) == 1 ) // inside or on surface
      {
      i_iter.Set( randomGenerator->GetVariate() * classSigmas[1] 
                  + classMeans[1] );
      }
    else
      {
      i_iter.Set( classMeans[0] );
      }
    ++i_iter;
    }
  
  // creates a bias field
  typedef itk::MultivariateLegendrePolynomial BiasFieldType;
  BiasFieldType::DomainSizeType biasSize(3);
  int biasDegree = 3;
  biasSize[0] = imageSize[0];
  biasSize[1] = imageSize[1];
  biasSize[2] = imageSize[2];
  BiasFieldType bias(biasSize.size(), 
                     biasDegree, // bias field degree 
                     biasSize);

  // generates the coefficients using the normal random variate generator.
  BiasFieldType::CoefficientArrayType 
    coefficients(bias.GetNumberOfCoefficients());
  BiasFieldType::CoefficientArrayType 
    initCoefficients(bias.GetNumberOfCoefficients());

  randomGenerator->Initialize( (int) 2003 );
  for ( unsigned int i = 0; i < bias.GetNumberOfCoefficients(); ++i )
    {
    coefficients[i] = ( randomGenerator->GetVariate() + 1 ) * 0.1;
    initCoefficients[i] = 0;
    }
  bias.SetCoefficients(coefficients);

  // set the imageWithBias pixel values with imageSource pixel value +
  // bias.
  ImageIteratorType ib_iter( imageWithBias, 
                             imageWithBias->GetLargestPossibleRegion() );

  BiasFieldType::SimpleForwardIterator b_iter( &bias );
  i_iter.GoToBegin();
  ib_iter.GoToBegin();
  float temp;
  while ( !i_iter.IsAtEnd() )
    {
    temp = i_iter.Get() * (2 + b_iter.Get()); // this is a multiplicative bias field
    ib_iter.Set( temp );
    ++i_iter;
    ++ib_iter;
    ++b_iter;
    }

  // creates a bias correction filter and run it.
  typedef itk::MRIBiasFieldCorrectionFilter<ImageType, ImageType, MaskType> 
    FilterType;

  FilterType::Pointer filter = FilterType::New();


  // To see the debug output for each iteration, uncomment the
  // following line. 
  // filter->DebugOn();

  double sumOfError = 0.0;
  i_iter.GoToBegin();
  ib_iter.GoToBegin();
  while ( !i_iter.IsAtEnd() )
    {
    sumOfError += vnl_math_abs( ib_iter.Get() - i_iter.Get() );
    ++i_iter;
    ++ib_iter;
    }
  std::cout << "Absolute Avg. error before correction = " 
            << sumOfError / (imageSize[0] * imageSize[1] * imageSize[2]) 
            << std::endl;
  double origSumError = sumOfError;

  std::cout << "Computing bias correction without mask, 2 classes 10,10 - 200,20" << std::endl; 
  filter->SetInput( imageWithBias.GetPointer() );
  filter->IsBiasFieldMultiplicative( true ); // correct with multiplicative bias 
  filter->SetBiasFieldDegree( biasDegree ); // default value = 3
  filter->SetTissueClassStatistics( classMeans, classSigmas );
  //filter->SetOptimizerGrowthFactor( 1.01 ); // default value
  //filter->SetOptimizerInitialRadius( 0.02 ); // default value
  filter->SetUsingInterSliceIntensityCorrection( true ); // default value
  filter->SetVolumeCorrectionMaximumIteration( 200 ); // default value = 100
  filter->SetInterSliceCorrectionMaximumIteration( 100 ); // default value = 100
  filter->SetUsingSlabIdentification( true ); // default value = false
  filter->SetSlabBackgroundMinimumThreshold( 0 ); // default value
  filter->SetSlabNumberOfSamples( 10 ); // default value 
  filter->SetSlabTolerance(0.0); // default value
  filter->SetSlicingDirection( 2 ); // default value
  filter->SetUsingBiasFieldCorrection( true ); // default value
  filter->SetGeneratingOutput( true ); // default value

  filter->SetInitialBiasFieldCoefficients(initCoefficients); //default value is all zero

  //timing
  long int t1 = time(NULL);
  filter->Update();