itkvectorimagetest.cxx

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

  
  typedef itk::ImageRegionIteratorWithIndex< VectorImageType > IteratorType;
  IteratorType it( image, region );
  it.GoToBegin();

  while( !it.IsAtEnd() )
    {
    VectorPixelType f;
    f[0] = it.GetIndex()[0];
    f[1] = it.GetIndex()[1];
    f[2] = it.GetIndex()[2];
    f[3] = it.GetIndex()[0];
    f[4] = it.GetIndex()[1];
    f[5] = it.GetIndex()[2];
    it.Set( f );
    ++it;
    }

  typedef itk::ImageFileWriter< VectorImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput( image );
  writer->SetFileName( argv[2] );
  writer->Update();

  writer->SetFileName( argv[1] );
  writer->Update();
  }
  
  {
  // Now read it as a itk::VectorImage.
  typedef itk::VectorImage< PixelType, Dimension > VectorImageType;
  typedef itk::ImageFileReader< VectorImageType >  ReaderType; 
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( argv[1] );
  reader->Update();

  VectorImageType::Pointer vectorImage = reader->GetOutput();
  
  typedef itk::ImageRegionConstIteratorWithIndex< VectorImageType > IteratorType;
  IteratorType cit( vectorImage, vectorImage->GetBufferedRegion() );
  cit.GoToBegin();

  bool failed1 = false;
  while( !cit.IsAtEnd() )
    {
    if(   (cit.Get()[0] != cit.GetIndex()[0])
       || (cit.Get()[1] != cit.GetIndex()[1])
       || (cit.Get()[2] != cit.GetIndex()[2])
       || (cit.Get()[3] != cit.GetIndex()[0])
       || (cit.Get()[4] != cit.GetIndex()[1])
       || (cit.Get()[5] != cit.GetIndex()[2]))
      {
      failed1 = true;
      }
    ++cit;
    }

  if( failed1 )
    {
    std::cerr << "Read VectorImage [FAILED]" << std::endl;
    failed = true;
    }
  else
    {
    std::cout << "Read VectorImage [PASSED]" << std::endl;
    }

  
  // Now write this out this VectorImage and read it again
  typedef itk::ImageFileWriter< VectorImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput( vectorImage );
  writer->SetFileName(argv[1]);
  writer->Update();
  }
  
  
  {
  // Now read it as a itk::VectorImage.
  typedef itk::VectorImage< PixelType, Dimension > VectorImageType;
  typedef itk::ImageFileReader< VectorImageType >  ReaderType;  
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( argv[1] );
  reader->Update();

  VectorImageType::Pointer vectorImage = reader->GetOutput();
  
  typedef itk::ImageRegionConstIteratorWithIndex< VectorImageType > IteratorType;
  IteratorType cit( vectorImage, vectorImage->GetBufferedRegion() );
  cit.GoToBegin();

  bool failed1 = false;
  while( !cit.IsAtEnd() )
    {
    if(   (cit.Get()[0] != cit.GetIndex()[0])
       || (cit.Get()[1] != cit.GetIndex()[1])
       || (cit.Get()[2] != cit.GetIndex()[2])
       || (cit.Get()[3] != cit.GetIndex()[0])
       || (cit.Get()[4] != cit.GetIndex()[1])
       || (cit.Get()[5] != cit.GetIndex()[2]))
      {
      failed1 = true;
      }
    ++cit;
    }

  if( failed1 )
    {
    std::cerr << "Write VectorImage [FAILED]" << std::endl;
    failed = true;
    }
  else
    {
    std::cout << "Write VectorImage [PASSED]" << std::endl;
    }

  
    {
    // Check support for Neighborhood Iterators
    //
    // 1. Test ConstNeighborhoodIterator
    // 
    std::cout << "Testing ConstNeighborhoodIterator...." << std::endl;
    
    typedef itk::ConstNeighborhoodIterator< VectorImageType > 
                                         ConstNeighborhoodIteratorType;
    ConstNeighborhoodIteratorType::RadiusType radius;
    radius[0] = radius[1] = radius[2] = 1;

    ConstNeighborhoodIteratorType::RegionType region 
                            = vectorImage->GetBufferedRegion();
    ConstNeighborhoodIteratorType::SizeType size;
    size[0]= size[1] = size[2] = 4;
    ConstNeighborhoodIteratorType::IndexType index;
    index[0] = index[1] = index[2] = 1; 
    region.SetIndex(index); 
    region.SetSize( size );
    
    ConstNeighborhoodIteratorType cNit(radius, vectorImage, region);

    // Move Iterator to a point and see if it reads out the right value
    //
    const unsigned int centerIndex = static_cast< unsigned int >(
                ((radius[0]*2+1)*(radius[1]*2+1)*(radius[2]*2+1))/2);

    ConstNeighborhoodIteratorType::IndexType location;
    location[0] = 1; location[1] = 2;  location[2] = 3;
    cNit.SetLocation( location );
    
    if( cNit.GetPixel(centerIndex) != vectorImage->GetPixel( location ) )
      {
      std::cerr << "  SetLocation [FAILED]" << std::endl;
      failed=true;
      }

    // Test GoToBegin()
    cNit.GoToBegin();
    if( cNit.GetPixel(centerIndex) != vectorImage->GetPixel( index ) )
      {
      std::cerr << "  GoToBegin [FAILED]" << std::endl;
      failed=true;
      }

    // Test GoToEnd()
    cNit.GoToEnd(); 
    --cNit;
    ConstNeighborhoodIteratorType::IndexType endIndex;
    endIndex[0] = endIndex[1] = endIndex[2] = 4;
    if( cNit.GetPixel(centerIndex) != vectorImage->GetPixel( endIndex ) )
      {
      std::cerr << "  GoToEnd [FAILED]" << std::endl;
      failed=true;
      }

    // Test IsAtEnd()
    if( !((++cNit).IsAtEnd()) )
      {
      std::cerr << "  IsAtEnd() [FAILED]" << std::endl;  
      failed = true;
      }
    cNit.GoToBegin();
    unsigned int numPixelsTraversed = size[0]*size[1]*size[2];
    while (! cNit.IsAtEnd())
      {
      ++cNit;
      --numPixelsTraversed;
      }
    if( numPixelsTraversed ) { std::cerr << "  IsAtEnd() [FAILED]" << std::endl; }

    // Test operator-
    --cNit;
    ConstNeighborhoodIteratorType::OffsetType offset;
    offset[0] = offset[1] = offset[2] = 1;
    cNit -= offset;
    itk::VariableLengthVector< PixelType > pixel = cNit.GetCenterPixel();
    itk::VariableLengthVector< PixelType > correctAnswer( VectorLength );
    correctAnswer.Fill( 3 );
    if( pixel != correctAnswer ) 
      {
      std::cerr << "  operator- [FAILED]" << std::endl;
      failed = true;
      }

    // Test GetNeighborhood()
    cNit.SetLocation( location );
    ConstNeighborhoodIteratorType::NeighborhoodType 
                    neighborhood = cNit.GetNeighborhood();
    //const unsigned int neighborhoodSize = neighborhood.Size();
    //for( unsigned int i=0; i< neighborhoodSize; i++)
    //  { std::cout << neighborhood[i] << std::endl; }
    if( (neighborhood[0][0] != 0) || (neighborhood[0][5] != 2))
      {
      std::cerr << "  GetNeighborhood() on ConstNeighborhoodIterator [FAILED]" << std::endl;
      failed = true;
      }



    // 
    // 2. Test NeighborhoodIterator on VectorImage
    //
    
    std::cout << "Testing NeighborhoodIterator..." << std::endl;
    
    typedef itk::NeighborhoodIterator< VectorImageType > NeighborhoodIteratorType;
    NeighborhoodIteratorType nit(radius, vectorImage, region);
    nit.SetLocation( location );
    itk::VariableLengthVector< PixelType > p( VectorLength );
    p.Fill( 100.0 );
    nit.SetNext( 1, 1, p );
    
    // Test SetNext()
    NeighborhoodIteratorType::IndexType index1;
    index1 = location; index1[1] = location[1] + 1;
    nit.SetLocation( index1 );

    if( nit.GetCenterPixel() != p )
      {
      std::cerr << "  SetNext() [FAILED]" << std::endl;
      failed = true;
      }
    p[0] = p[3] = (float)index1[0];  
    p[1] = p[4] = (float)index1[1]; 
    p[2] = p[5] = (float)index1[2]; 
    nit.SetCenterPixel( p );
    if( nit.GetCenterPixel() != p )
      {
      std::cerr << "  SetCenterPixel() [FAILED]" << std::endl;
      failed = true;
      }
    
    // Test SetNeighborhood() and GetPrevious()
    nit.SetLocation( index1 );
    nit.SetNeighborhood( neighborhood );
    p[0] = p[3] = 1; p[1] = p[4] = 2; p[2] = p[5] = 2;
    if( nit.GetPrevious( 2, 1 ) != p )
      {
      std::cerr << "  SetNeighborhood() or GetPrevious() [FAILED]" << std::endl;
      failed = true;
      }
    

    // 
    // 3. Testing ConstShapedNeighborhoodIterator on VectorImage
    //

    // Go back to original image, where pixel values tell us the indices
    std::cout << "Testing ConstShapedNeighborhoodIterator on VectorImage..." 
                                                                  << std::endl;
    reader->SetFileName( "dummy.nrrd");
    reader->SetFileName( argv[1] );
    reader->Update();
    vectorImage = reader->GetOutput(); 

    typedef itk::ConstShapedNeighborhoodIterator< VectorImageType > 
                                   ConstShapedNeighborhoodIteratorType;
    ConstShapedNeighborhoodIteratorType cSnit( radius, vectorImage, region );
    cSnit.SetLocation( location );
    ConstShapedNeighborhoodIteratorType::OffsetType offset1;
    offset1[0] = 0; offset1[1] = 0; offset1[2] = 0;
    cSnit.ActivateOffset( offset1 ); //activate the center
    // activate the top plane
    offset1[0] = 0; offset1[1] = 0; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = 0; offset1[1] = 1; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = 0; offset1[1] = -1; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = 1; offset1[1] = 0; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = 1; offset1[1] = 1; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = 1; offset1[1] = -1; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = -1; offset1[1] = 0; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = -1; offset1[1] = 1; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
    offset1[0] = -1; offset1[1] = -1; offset1[2] = -1;
    cSnit.ActivateOffset( offset1 ); 
 
    ConstShapedNeighborhoodIteratorType::IndexListType l
                                  = cSnit.GetActiveIndexList();
    ConstShapedNeighborhoodIteratorType::IndexListType::const_iterator 
                                                        ali = l.begin();
    while (ali != l.end())
      {
      std::cout << *ali << " ";
      ++ali;
      }
    std::cout << std::endl;

    ConstShapedNeighborhoodIteratorType::ConstIterator ci = cSnit.Begin();
    while (! ci.IsAtEnd())
      {
      ConstShapedNeighborhoodIteratorType::OffsetType offset2 = ci.GetNeighborhoodOffset();
      if( (offset2[0] == -1) && (offset2[1]== -1) && (offset2[2]== -1) )
        { 
        if( ci.GetNeighborhoodIndex() != 0 )
          {
          failed = true;
          std::cerr << "GetNeighborhoodOffset() on ConstShapedNeighborhoodIterato [FAILED]" 
                                                                              << std::endl;
          }
        if( (ci.Get()[0]!=0) || (ci.Get()[1]!=1) || (ci.Get()[2]!=2) )
          {
          failed=true;
          std::cerr 
            << "ConstShapedNeighborhoodIterator returned incorrect index [FAILED]" 
                                                                        << std::endl;
          }
        }
      ci++;
      }

    //
    // 4. Test ShapedNeighborhoodIterator 
    //
    typedef itk::ShapedNeighborhoodIterator< VectorImageType > 
                                      ShapedNeighborhoodIteratorType;
    ShapedNeighborhoodIteratorType sNit( radius, vectorImage, region );
    
    offset1[0] = 0; offset1[1] = 0; offset1[2] = 0;
    sNit.ActivateOffset( offset1 ); //activate the center
    // activate the top plane
    offset1[0] = 0; offset1[1] = 0; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = 0; offset1[1] = 1; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = 0; offset1[1] = -1; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = 1; offset1[1] = 0; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = 1; offset1[1] = 1; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = 1; offset1[1] = -1; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = -1; offset1[1] = 0; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = -1; offset1[1] = 1; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    offset1[0] = -1; offset1[1] = -1; offset1[2] = -1;
    sNit.ActivateOffset( offset1 ); 
    
    sNit.SetLocation( location );
    ShapedNeighborhoodIteratorType::Iterator shit = sNit.Begin();
    shit = sNit.Begin(); 
    p[0] = p[3] = 10; p[1] = p[4] = 20; p[2] = p[5] = 30;
    shit.Set( p );
    index[0]=location[0]-1; index[1]=location[1]-1; index[2]=location[2]-1;
    cNit.SetLocation( index );
    if( cNit.GetCenterPixel() != p )
      {
      std::cerr << "ShapedNeighborhoodIterator Set() [FAILED]" << std::endl;
      failed=true;
      }

    
    }  // End Testing Neighborhood Iterators on VectorImage

  }

  if( failed )
    {
    return EXIT_FAILURE;
    }

  return EXIT_SUCCESS;
}