[Insight-users] An issue in the sum calculation of Spatia lObjectToImageStatisticsCalculator

[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. Furthermore, I compared the results with the sum calculation of ume calculator: 292569LabelStatisticsImageFilter. 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 label calculator: 4477
overlapping sum from statistical calculator: 4477
overlapping sum from volume calculator: 4477
sum of temp.mha from label calculator: 292569
sum of temp.mha from statistical calculator: 30425
sum of temp.mha from volume calculator: 292569
sum of temp1.mha from label calculator: 6980
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 * 2.0 ) ;
  thresholder -> SetUpperThreshold ( PIXELONE * 2.0 ) ;
  thresholder -> Update ( ) ;

  /*
  typedef itk::ImageFileWriter< ImageVType > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput( thresholder -> GetOutput ( ) );
  writer -> UseCompressionOn ( ) ;
  writer->SetFileName("temp4.mha");
  writer->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 GetSumL ( ImageVType::Pointer imageP, 
  ImageVType::Pointer maskImageP, double & sum )
  {
  typedef itk::LabelStatisticsImageFilter< 
    ImageVType, ImageVType > LabelFilterType;
  LabelFilterType::Pointer calculator = LabelFilterType::New();
  calculator->SetInput(imageP);
  calculator->SetLabelInput( maskImageP );
  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 ( 1.0 ) ;
  } 
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;
  GetSumL ( image, image1, sum ) ;
  std::cout<<"overlapping sum from label calculator: "
    <<sum<<std::endl;
  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;
  GetSumL ( image, image, sum ) ;
  std::cout<<"sum of temp.mha from label 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;
  GetSumL ( image1, image1, sum ) ;
  std::cout<<"sum of temp1.mha from label 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/6abb0663/attachment.htm