itkbsplineresampleimagefiltertest.cxx

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

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkBSplineResampleImageFilterTest.cxx,v $
  Language:  C++
  Date:      $Date: 2007-08-10 14:34:01 $
  Version:   $Revision: 1.17 $

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

  Portions of this code are covered under the VTK copyright.
  See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.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 "itkImage.h"
#include "itkSize.h"
#include "itkBSplineDownsampleImageFilter.h"
#include "itkBSplineUpsampleImageFilter.h"
#include "itkImageLinearIteratorWithIndex.h"
#include "itkFilterWatcher.h"


#include "vnl/vnl_math.h"

typedef double InputPixelType;
typedef int    IntInputPixelType;

// Set up for 2D Images
enum { ImageDimension2D = 2 };

typedef itk::Image< InputPixelType, ImageDimension2D > ImageType2D;
typedef ImageType2D::Pointer ImageTypePtr2D;
typedef ImageType2D::SizeType SizeType2D;
typedef itk::ImageRegionIterator<ImageType2D>  InputIterator;

typedef itk::Image< IntInputPixelType, ImageDimension2D > IntImageType2D;
typedef IntImageType2D::Pointer IntImageTypePtr2D;
typedef IntImageType2D::SizeType IntSizeType2D;
typedef itk::ImageRegionIterator<IntImageType2D>  IntInputIterator;

void set2DData(ImageType2D::Pointer);

void PrintImageData(ImageTypePtr2D imgPtr)
{
  typedef itk::ImageLinearIteratorWithIndex<ImageType2D> Iterator;
 
  std::cout << "Size: " << imgPtr->GetLargestPossibleRegion().GetSize() << std::endl;
  int dim = ImageType2D::ImageDimension;

  std::cout << "Spacing: " << std::endl;
  for (int n = 0; n < dim ; n++)
    {
    std::cout << imgPtr->GetSpacing()[n] << ", ";
    }
  std::cout << std::endl;
  Iterator outIt = Iterator( imgPtr, imgPtr->GetLargestPossibleRegion() );
  outIt.SetDirection(0);

  SizeType2D size = imgPtr->GetLargestPossibleRegion().GetSize();

  std::cout << "Data: " <<std::endl;
  for (int n=0; n < dim - 1 ; n++)
    {
    for (unsigned int jj=0; jj < size[n + 1]; jj++)
      {

      while ( !outIt.IsAtEndOfLine() )
        {
        std::cout << outIt.Get() << ", " ;
        ++outIt;
        }
       outIt.NextLine();
       std::cout << std::endl;
      }
    }   
}

void set2DData(ImageType2D::Pointer imgPtr)
{
  SizeType2D size = { {4,4} };
  double mydata[ 49 ] = {  0, 1, 2, 3, 
    1, 2, 3, 4,
    2, 3, 4, 5,
    3, 4, 3, 2};

  ImageType2D::RegionType region;
  region.SetSize( size );

  imgPtr->SetLargestPossibleRegion( region );
  imgPtr->SetBufferedRegion( region );
  imgPtr->Allocate();

  // Set origin and spacing of physical coordinates 
  double origin [] = { 0.5, 1.0 };
  double spacing[] = { 0.1, 0.5  };
  imgPtr->SetOrigin(origin);
  imgPtr->SetSpacing(spacing);

  InputIterator inIter( imgPtr, region );

  int j = 0;
  while( !inIter.IsAtEnd() )
    {
    inIter.Set(mydata[j]);
    ++inIter;
    ++j;
    }  
}

void setInt2DData(IntImageType2D::Pointer imgPtr)
{
  IntSizeType2D size = { {4,4} };
  int mydata[ 49 ] = {  0, 1, 2, 3, 
    1, 2, 3, 4,
    2, 3, 4, 5,
    3, 4, 3, 2};

  IntImageType2D::RegionType region;
  region.SetSize( size );

  imgPtr->SetLargestPossibleRegion( region );
  imgPtr->SetBufferedRegion( region );
  imgPtr->Allocate();

  // Set origin and spacing of physical coordinates 
  double origin [] = { 0.5, 1.0 };
  double spacing[] = { 0.1, 0.5  };
  imgPtr->SetOrigin(origin);
  imgPtr->SetSpacing(spacing);

  IntInputIterator inIter( imgPtr, region );

  int j = 0;
  while( !inIter.IsAtEnd() )
    {
    inIter.Set(mydata[j]);
    ++inIter;
    ++j;
    }  
}

bool VerifyResultsHigherOrderSpline(ImageTypePtr2D ActualResults, double *ExpectedResults)
{
  double * ERptr;
  ERptr = ExpectedResults;

  InputIterator ActualResultsIter( ActualResults, ActualResults->GetLargestPossibleRegion() );
  double percentErr = 0;

  while (!ActualResultsIter.IsAtEnd() )
    {
    double val1 = ActualResultsIter.Get();

    percentErr += vnl_math_abs( ( val1 - * ERptr ) / val1 );

    ++ActualResultsIter;
    ++ERptr;
    }

  //Mean error is determined over the number of pixels, which in the test case is 16
  if( (percentErr/16) > 0.05 )
    {
    // std::cout << "*** Error: total error is more than 10%: " << percentErr << std::endl;
    return false;
    }
  return true;
}

bool VerifyResults3rdOrderSpline(ImageTypePtr2D ActualResults, double *ExpectedResults)
{
  double * ERptr;
  ERptr = ExpectedResults;

  InputIterator ActualResultsIter( ActualResults, ActualResults->GetLargestPossibleRegion() );

  while (!ActualResultsIter.IsAtEnd() )
    {
    double val1 = ActualResultsIter.Get();
    if( vnl_math_abs( val1 - * ERptr ) > 1e-6 )
      {
      // std::cout << "*** Error: value should be " << trueValue << std::endl;
      return false;
      }
    ++ActualResultsIter;
    ++ERptr;
    }
  return true;
}

bool VerifyResults2ndOrderSpline(ImageTypePtr2D ActualResults, double *ExpectedResults)
{
  double * ERptr;
  ERptr = ExpectedResults;

  InputIterator ActualResultsIter( ActualResults, ActualResults->GetLargestPossibleRegion() );
  double percentErr = 0;

  while (!ActualResultsIter.IsAtEnd() )
    {
    double val1 = ActualResultsIter.Get();

    percentErr += vnl_math_abs( ( val1 - * ERptr ) / val1 );

    ++ActualResultsIter;
    ++ERptr;
    }

  //Mean error is determined over the number of pixels, which in the test case is 16
  if( (percentErr/16) > 0.1 )
    {
    // std::cout << "*** Error: total error is more than 10%: " << percentErr << std::endl;
    return false;
    }
  return true;
}

//Spline order 1 or 0
bool VerifyResultsLowerOrderSpline(ImageTypePtr2D ActualResults, double *ExpectedResults)
{
  double * ERptr;
  ERptr = ExpectedResults;

  InputIterator ActualResultsIter( ActualResults, ActualResults->GetLargestPossibleRegion() );
  double percentErr = 0;

  while (!ActualResultsIter.IsAtEnd() )
    {
    double val1 = ActualResultsIter.Get();

    percentErr += vnl_math_abs( ( val1 - * ERptr ) / val1 );

    ++ActualResultsIter;
    ++ERptr;
    }

  //Mean error is determined over the number of pixels, which in the test case is 16
  //Threshold error is much higher since the order of the spline is lower.
  if( (percentErr)/16 > 0.6 )
    {
    return false;
    }
  return true;
}

   
int test2D_Standard_l2_NthOrderSpline_filter(unsigned int splineOrder)
{
  int flag = 0;

  /* Allocate a simple test image */
  ImageTypePtr2D image = ImageType2D::New();

  set2DData(image);
  double ExpectedResults[] = {0.517188, 1.575548, 2.634784, 2.847124,
                              1.547224, 2.323852, 3.101771, 3.257328,
                              2.578121, 3.073447, 3.570484, 3.669343,
                              2.784775, 3.223312, 3.663641, 3.751054};

  // l2 norm resampler.
  typedef itk::BSplineDownsampleImageFilter<ImageType2D,ImageType2D> DownsamplerType2D;
  typedef itk::BSplineUpsampleImageFilter<ImageType2D,ImageType2D> UpsamplerType2D;

  DownsamplerType2D::Pointer downSampler = DownsamplerType2D::New();
  FilterWatcher downWatcher(downSampler, "test2D_Standard_l2_filter");

  UpsamplerType2D::Pointer   upSampler =   UpsamplerType2D::New();
  FilterWatcher upWatcher(upSampler, "test2D_Standard_l2_filter");

  downSampler->SetSplineOrder(splineOrder);
  upSampler->SetSplineOrder(splineOrder);

  downSampler->SetInput(image);
  downSampler->Update();
  ImageTypePtr2D outImage1 = downSampler->GetOutput();
  PrintImageData(outImage1);
  upSampler->SetInput( outImage1 );
    upSampler->Update();
  ImageTypePtr2D outImage2 = upSampler->GetOutput();
  PrintImageData(outImage2);

  bool sameResults=false;
  if( splineOrder == 3 ) 
    {
    sameResults = VerifyResults3rdOrderSpline(outImage2, ExpectedResults);
    }
  else if( splineOrder == 2 ) 
    {
    sameResults = VerifyResultsLowerOrderSpline(outImage2, ExpectedResults);
    }
  else if( splineOrder == 1 ) 
    {
    sameResults = VerifyResultsLowerOrderSpline(outImage2, ExpectedResults);
    }
  else if( splineOrder == 0 ) 
    {
    sameResults = VerifyResultsLowerOrderSpline(outImage2, ExpectedResults);
    }

  if (!sameResults)
    {
    flag = 1;
    std::cout << "*** Error: unexpected value in Standard l2 - resampler with order " << splineOrder <<
      "  spline." << std::endl;
    std::cout << "" << std::endl;
    }
  else
    {
    std::cout << "Tests for Standard l2 - resampler with order " <<  splineOrder <<  "  spline PASSED " << std::endl;
    std::cout << "" << std::endl;
    }

  return flag;
}


int test2D_Standard_L2_NthOrderSpline_filter(unsigned int splineOrder)
{
  int flag = 0;

  // Allocate a simple test image 
  ImageTypePtr2D image = ImageType2D::New();

  set2DData(image);
  double ExpectedResults[] = {0.527796, 1.584850, 2.642784, 2.854860,
                              1.554069, 2.328161, 3.103548, 3.258594,
                              2.581206, 3.072767, 3.566038, 3.664140,
                              2.787102, 3.221627, 3.657944, 3.744552};

  // L2 norm resampler.
  typedef itk::BSplineL2ResampleImageFilterBase<ImageType2D, ImageType2D> ResamplerType;
  typedef itk::BSplineDownsampleImageFilter<ImageType2D,ImageType2D,ResamplerType> DownsamplerType2D;
  typedef itk::BSplineUpsampleImageFilter<ImageType2D,ImageType2D,ResamplerType> UpsamplerType2D;

  DownsamplerType2D::Pointer downSampler = DownsamplerType2D::New();
  FilterWatcher downWatcher(downSampler, "test2D_Standard_L2_filter");

  UpsamplerType2D::Pointer   upSampler =   UpsamplerType2D::New();
  FilterWatcher upWatcher(upSampler, "test2D_Standard_L2_filter");

  downSampler->SetSplineOrder(splineOrder);
  upSampler->SetSplineOrder(splineOrder);

  downSampler->SetInput(image);
  downSampler->Update();
  ImageTypePtr2D outImage1 = downSampler->GetOutput();
  PrintImageData(outImage1);
  upSampler->SetInput( outImage1 );