<commit_before>/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    itkMRIBiasFieldCorrectionFilterTest.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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] = 60 ;
  imageSize[1] = 60 ;
  imageSize[2] = 30 ;
  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->SetVolumeCorrectionMaximumIteration( 1000 ) ; // default value = 100
  filter->SetUsingInterSliceIntensityCorrection( true ) ; // default value
  filter->SetInterSliceCorrectionMaximumIteration( 200 ) ; // 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() ;
  long int t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl ;

  sumOfError = 0.0 ;<commit_msg>BUG: Testing time decreased by half<commit_after>/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    itkMRIBiasFieldCorrectionFilterTest.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  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] = 60 ;
  imageSize[1] = 60 ;
  imageSize[2] = 30 ;
  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->SetVolumeCorrectionMaximumIteration( 500 ) ; // default value = 100
  filter->SetUsingInterSliceIntensityCorrection( true ) ; // default value
  filter->SetInterSliceCorrectionMaximumIteration( 200 ) ; // 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() ;
  long int t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl ;

  sumOfError = 0.0 ;
{'fraction_non_alphanumeric': 0.13665749215899636, 'fraction_numerical': 0.01657812199961595, 'mean_word_length': 3.786764705882353, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 105, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5627986', 'n_tokens_mistral': 5214, 'n_tokens_neox': 4880, 'n_words': 1115}