[Insight-users] A SpatialObjectToImageStatisticsCalculato r issue

[kurt]

Hi All,

I tried to use the SpatialObjectToImageStatisticsCalculator to calculate the volume of the portion overlapping two binary images, temp.mha and temp1.mha, as well as the volume of the two images. I was used to use SpatialObjectToImageStatisticsCalculator to calculate the sum of the image, the volume can be then calculated. However, the sum seems not right. the sum is not equal to the value by iterating all pixels. My concern is that if the filter performs properly.
the two images are avaliable at
http://www.duke.edu/~kurtzhao/temp.mha
http://www.duke.edu/~kurtzhao/temp1.zip

The output on my machine is 
overlapping sum from statistical calculator: 4477
overlapping sum from volume calculator: 295072
sum of temp.mha from statistical calculator: 30425
sum of temp.mha from volume calculator: 292569
sum of temp1.mha from statistical calculator: 6980
sum of temp1.mha from volume calculator: 6980

Thanks

Kurt

#include <iostream>
#include <itkImageFileReader.h>
#include <itkImageMaskSpatialObject.h>
#include <itkSpatialObjectToImageStatisticsCalculator.h>

#include <itkAddImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

typedef itk::Image < unsigned char, 3 > ImageVType;


void GetSumV ( ImageVType::Pointer imageP, 
  ImageVType::Pointer maskImageP, double & sum )
  {
  const unsigned char PIXELONE = 1 ;

  typedef itk::AddImageFilter< ImageVType, 
          ImageVType, ImageVType>  AddType;
  AddType::Pointer add = AddType::New ( ) ;
  add -> SetInput1 ( imageP);
  add -> SetInput2 ( maskImageP);
  add -> Update ( ) ;
  
  typedef itk::BinaryThresholdImageFilter< ImageVType, 
          ImageVType>  ThresholdType;
  ThresholdType::Pointer thresholder 
    = ThresholdType::New ( ) ;
  thresholder -> SetInput ( add -> GetOutput ( )) ;
  thresholder -> SetOutsideValue ( 0 ) ;
  thresholder -> SetInsideValue ( PIXELONE) ;
  thresholder -> SetLowerThreshold ( PIXELONE  ) ;
  thresholder -> SetUpperThreshold ( PIXELONE * 2.0 ) ;
  thresholder -> Update ( ) ;

  ImageVType::Pointer image2P = thresholder -> GetOutput ( ) ;
  sum = 0 ;
  typedef itk::ImageRegionConstIterator< ImageVType > ConstIteratorType;
  ConstIteratorType it( image2P , image2P -> GetLargestPossibleRegion ( )  );
  for ( it.GoToBegin(); !it.IsAtEnd(); ++it)
    {
    if ( it.Get( ) == PIXELONE ) 
      {
      sum ++ ;
      }
    }
  }
void GetSumS ( ImageVType::Pointer imageP, 
  ImageVType::Pointer maskImageP, double & sum )
  {
  typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;
  ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();
  maskSO->SetImage( maskImageP );

  typedef itk::SpatialObjectToImageStatisticsCalculator< 
    ImageVType, ImageMaskSpatialObject> CalculatorType;
  CalculatorType::Pointer calculator = CalculatorType::New();
  calculator->SetImage(imageP);
  calculator->SetSpatialObject( maskSO );
  calculator->Update();
  /*
  std::cout << "Sample mean = " << calculator->GetMean()<<" ";
  std::cout << "Sample STD = " 
    << sqrt ( calculator->GetCovarianceMatrix()[0][0] ) 
    << std::endl;
  std::cout << "Sample sum = " 
    <<  calculator -> GetSum ( ) 
    << std::endl;
    */
  sum =  calculator -> GetSum ( ) ;
  } 
void GetOverlappingSum ( )
  {

  typedef itk::ImageFileReader < ImageVType > ImageReaderType;
  ImageReaderType::Pointer reader = ImageReaderType::New();
  reader -> SetFileName( "temp.mha" ) ;
  reader -> Update();

  ImageVType::Pointer image = reader -> GetOutput ( ) ;

  ImageReaderType::Pointer reader1 = ImageReaderType::New();
  reader1 -> SetFileName( "temp1.mha" ) ;
  reader1 -> Update();
  ImageVType::Pointer image1 = reader1 -> GetOutput ( ) ;

  double sum;
  GetSumS ( image, image1, sum ) ;
  std::cout<<"overlapping sum from statistical calculator: "
    <<sum<<std::endl;
  GetSumV ( image, image1, sum ) ;
  std::cout<<"overlapping sum from volume calculator: "
    <<sum<<std::endl;
  GetSumS ( image, image, sum ) ;
  std::cout<<"sum of temp.mha from statistical calculator: "
    <<sum<<std::endl;
  GetSumV ( image, image, sum ) ;
  std::cout<<"sum of temp.mha from volume calculator: "
    <<sum<<std::endl;
  GetSumS ( image1, image1, sum ) ;
  std::cout<<"sum of temp1.mha from statistical calculator: "
    <<sum<<std::endl;
  GetSumV ( image1, image1, sum ) ;
  std::cout<<"sum of temp1.mha from volume calculator: "
    <<sum<<std::endl;
  }


int main(int argc, char **argv)
  {

  GetOverlappingSum ( ) ;

  return 1 ;
  }
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20060913/1d32ce55/attachment.html