itkmribiasfieldcorrectionfiltertest.cxx

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

  long int t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl;

  sumOfError = 0.0;
  ImageIteratorType o_iter( filter->GetOutput(), 
                            filter->GetOutput()->GetLargestPossibleRegion() );
  i_iter.GoToBegin();
  while ( !i_iter.IsAtEnd() )
    {
    sumOfError += vnl_math_abs( o_iter.Get() - i_iter.Get() );
    ++i_iter;
    ++o_iter;
    }

  if ( SaveImages )
    {
    typedef itk::ImageFileWriter< ImageType > WriterType;
    WriterType::Pointer writer = WriterType::New();
    
    writer->SetInput( image );
    writer->SetFileName( "MRISource.mhd" );
    writer->Update();
    
    WriterType::Pointer writer2 = WriterType::New();
    writer2->SetInput(imageWithBias);
    writer2->SetFileName( "MRISourceWithBias.mhd" );
    writer2->Update();
    
    WriterType::Pointer writer3 = WriterType::New();
    writer3->SetInput(filter->GetOutput());
    writer3->SetFileName( "MRICorrected.mhd" );
    writer3->Update();
  
    }

  std::cout << "Absolute Avg. error without input and output mask = " 
            << sumOfError / (imageSize[0] * imageSize[1] * imageSize[2]) 
            << std::endl;
  if (origSumError < sumOfError)
    {
    std::cout << "ERROR: sumOfError: " << sumOfError 
              << " is less than origSumError: " << origSumError
              << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Computing bias correction with mask" << std::endl;
  filter->SetInput( imageWithBias.GetPointer() );
  filter->SetInputMask( image.GetPointer() );
  filter->SetOutputMask( image.GetPointer() );
  filter->SetInitialBiasFieldCoefficients(initCoefficients);
  filter->SetVolumeCorrectionMaximumIteration( 200 ); // default value = 100
  filter->SetInterSliceCorrectionMaximumIteration( 100 ); // default value = 100
  //filter->SetOptimizerInitialRadius( 0.02 ); // default value 
  t1 = time(NULL);
  filter->Update();
  t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl;

  sumOfError = 0.0;
  ImageIteratorType o2_iter( filter->GetOutput(), 
                             filter->GetOutput()->GetLargestPossibleRegion() );
  i_iter.GoToBegin();
  while ( !i_iter.IsAtEnd() )
    {
    sumOfError += vnl_math_abs( o2_iter.Get() - i_iter.Get() );
    ++i_iter;
    ++o2_iter;
    }

  std::cout << "Absolute Avg. error with input and output mask = " 
            << sumOfError / (imageSize[0] * imageSize[1] * imageSize[2]) 
            << std::endl;
  if (origSumError < sumOfError)
    {
    std::cout << "ERROR: sumOfError: " << sumOfError 
              << " is less than origSumError: " << origSumError
              << std::endl;
    return EXIT_FAILURE;
    }

  // default schedule is 2 2 2 - 1 1 1, let's change this
  std::cout << "Computing bias correction only with 2,2,2 resolution & no interSlice/Slab" << std::endl;
  filter->SetUsingInterSliceIntensityCorrection( false ); // default value
  filter->SetUsingSlabIdentification( false ); // default value = false
  FilterType::ScheduleType schedule ( 1, ImageDimension);
  schedule.Fill( 2 );
  filter->SetNumberOfLevels( 1 ); // Important to set this first, otherwise the filter rejects the new schedule
  filter->SetSchedule( schedule );
  filter->SetInitialBiasFieldCoefficients(initCoefficients); 
  filter->SetVolumeCorrectionMaximumIteration( 200 ); // default value = 100
  filter->SetInterSliceCorrectionMaximumIteration( 100 ); // default value = 100
  //filter->SetOptimizerInitialRadius( 0.02 ); // default value
  t1 = time(NULL);
  filter->Update();
  t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl;

  sumOfError = 0.0;
  ImageIteratorType o3_iter( filter->GetOutput(), 
                             filter->GetOutput()->GetLargestPossibleRegion() );
  i_iter.GoToBegin();
  while ( !i_iter.IsAtEnd() )
    {
    sumOfError += vnl_math_abs( o3_iter.Get() - i_iter.Get() );
    ++i_iter;
    ++o3_iter;
    }

  std::cout << "Absolute Avg. error with input and output mask = " 
            << sumOfError / (imageSize[0] * imageSize[1] * imageSize[2]) 
            << std::endl;
  if (origSumError < sumOfError)
    {
    std::cout << "ERROR: sumOfError: " << sumOfError 
              << " is less than origSumError: " << origSumError
              << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Computing bias correction only with 4,4,4 resolution & no interSlice/Slab" << std::endl;
  filter->SetUsingInterSliceIntensityCorrection( false ); // default value
  filter->SetUsingSlabIdentification( false ); // default value = false
  schedule.Fill( 4 );
  filter->SetNumberOfLevels( 1 ); 
  //filter->SetOptimizerInitialRadius( 0.02 ); // default value
  filter->SetSchedule( schedule );
  filter->SetVolumeCorrectionMaximumIteration( 200 ); // default value = 100
  filter->SetInterSliceCorrectionMaximumIteration( 100 ); // default value = 100
  filter->SetInitialBiasFieldCoefficients(initCoefficients); 
  t1 = time(NULL);
  filter->Update();
  t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl;

  sumOfError = 0.0;
  ImageIteratorType o4_iter( filter->GetOutput(), 
                             filter->GetOutput()->GetLargestPossibleRegion() );
  i_iter.GoToBegin();
  while ( !i_iter.IsAtEnd() )
    {
    sumOfError += vnl_math_abs( o4_iter.Get() - i_iter.Get() );
    ++i_iter;
    ++o4_iter;
    }

  std::cout << "Absolute Avg. error with input and output mask = " 
            << sumOfError / (imageSize[0] * imageSize[1] * imageSize[2]) 
            << std::endl;
  if (origSumError < sumOfError)
    {
    std::cout << "ERROR: sumOfError: " << sumOfError 
              << " is less than origSumError: " << origSumError
              << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Computing bias correction only with 4,4,4 resolution & no interSlice/Slab & more iterations" << std::endl;
  initCoefficients = filter->GetEstimatedBiasFieldCoefficients();
  filter->SetInitialBiasFieldCoefficients(initCoefficients); 
  filter->SetVolumeCorrectionMaximumIteration( 2000 ); // default value = 100
  filter->SetInterSliceCorrectionMaximumIteration( 100 ); // default value = 100
  t1 = time(NULL);
  filter->Update();
  t2 = time(NULL);
  std::cout << "Run time (in s)" << t2-t1  << std::endl;

  double sumOfErrorFinal = 0.0;
  ImageIteratorType o5_iter( filter->GetOutput(), 
                             filter->GetOutput()->GetLargestPossibleRegion() );
  i_iter.GoToBegin();
  while ( !i_iter.IsAtEnd() )
    {
    sumOfErrorFinal += vnl_math_abs( o5_iter.Get() - i_iter.Get() );
    ++i_iter;
    ++o5_iter;
    }

  std::cout << "Absolute Avg. error with input and output mask = " 
            << sumOfErrorFinal / (imageSize[0] * imageSize[1] * imageSize[2]) 
            << std::endl;
  if (sumOfError < sumOfErrorFinal)
    {
    std::cout << "ERROR: sumOfError: " << sumOfError 
              << " is less than sumOfErrorFinal: " << sumOfErrorFinal
              << std::endl;
    return EXIT_FAILURE;
    }

  std::cout << "Using slab identification: " 
            << filter->GetUsingSlabIdentification() << std::endl;
  std::cout << "Slab identification background minimum intensity threshold: "
            << filter->GetSlabBackgroundMinimumThreshold() << std::endl;
  std::cout << "Slab number of samples per slice: " 
            << filter->GetSlabNumberOfSamples() << std::endl;
  std::cout << "Slab identification tolerance: "
            << filter->GetSlabTolerance() << std::endl;
  std::cout << "Using bias field correction: "
            << filter->GetUsingBiasFieldCorrection() << std::endl;
  std::cout << "Generating output: "
            << filter->GetGeneratingOutput() << std::endl;
  std::cout << "Bias field degree: "
            << filter->GetBiasFieldDegree() << std::endl;
  std::cout << "Volume bias field correction iterations: "
            << filter->GetVolumeCorrectionMaximumIteration() << std::endl;
  std::cout << "Interslice bias field correction iterations: "
            << filter->GetInterSliceCorrectionMaximumIteration() << std::endl;
  std::cout << "Optimizer initial radius: "
            << filter->GetOptimizerInitialRadius() << std::endl;
  std::cout << "Optimizer growth factor: "
            << filter->GetOptimizerGrowthFactor() << std::endl;
  std::cout << "Optimizer shrink factor: "
            << filter->GetOptimizerShrinkFactor() << std::endl;

  return EXIT_SUCCESS;
}