Hi All,<br><br>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.<br>the two images are avaliable at<br>http://www.duke.edu/~kurtzhao/temp.mha<br>http://www.duke.edu/~kurtzhao/temp1.zip<br><br>The output on my machine is <br>overlapping sum from statistical calculator: 4477<br>overlapping sum from volume calculator: 295072<br>sum of temp.mha from statistical calculator: 30425<br>sum of temp.mha from volume calculator: 292569<br>sum of temp1.mha from statistical calculator: 6980<br>sum of temp1.mha from volume calculator: 6980<br><br>Thanks<br><br>Kurt<br><br>#include <iostream><br>#include <itkImageFileReader.h><br>#include <itkImageMaskSpatialObject.h><br>#include <itkSpatialObjectToImageStatisticsCalculator.h><br><br>#include <itkAddImageFilter.h><br>#include <itkBinaryThresholdImageFilter.h><br><br>typedef itk::Image < unsigned char, 3 > ImageVType;<br><br><br>void GetSumV ( ImageVType::Pointer imageP, <br> ImageVType::Pointer maskImageP, double & sum )<br> {<br> const unsigned char PIXELONE = 1 ;<br><br> typedef itk::AddImageFilter< ImageVType, <br> ImageVType, ImageVType> AddType;<br> AddType::Pointer add = AddType::New ( ) ;<br> add -> SetInput1 ( imageP);<br> add -> SetInput2 ( maskImageP);<br> add -> Update ( ) ;<br> <br> typedef itk::BinaryThresholdImageFilter< ImageVType, <br> ImageVType> ThresholdType;<br> ThresholdType::Pointer thresholder <br> = ThresholdType::New ( ) ;<br> thresholder -> SetInput ( add -> GetOutput ( )) ;<br> thresholder -> SetOutsideValue ( 0 ) ;<br> thresholder -> SetInsideValue ( PIXELONE) ;<br> thresholder -> SetLowerThreshold ( PIXELONE ) ;<br> thresholder -> SetUpperThreshold ( PIXELONE * 2.0 ) ;<br> thresholder -> Update ( ) ;<br><br> ImageVType::Pointer image2P = thresholder -> GetOutput ( ) ;<br> sum = 0 ;<br> typedef itk::ImageRegionConstIterator< ImageVType > ConstIteratorType;<br> ConstIteratorType it( image2P , image2P -> GetLargestPossibleRegion ( ) );<br> for ( it.GoToBegin(); !it.IsAtEnd(); ++it)<br> {<br> if ( it.Get( ) == PIXELONE ) <br> {<br> sum ++ ;<br> }<br> }<br> }<br>void GetSumS ( ImageVType::Pointer imageP, <br> ImageVType::Pointer maskImageP, double & sum )<br> {<br> typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;<br> ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();<br> maskSO->SetImage( maskImageP );<br><br> typedef itk::SpatialObjectToImageStatisticsCalculator< <br> ImageVType, ImageMaskSpatialObject> CalculatorType;<br> CalculatorType::Pointer calculator = CalculatorType::New();<br> calculator->SetImage(imageP);<br> calculator->SetSpatialObject( maskSO );<br> calculator->Update();<br> /*<br> std::cout << "Sample mean = " << calculator->GetMean()<<" ";<br> std::cout << "Sample STD = " <br> << sqrt ( calculator->GetCovarianceMatrix()[0][0] ) <br> << std::endl;<br> std::cout << "Sample sum = " <br> << calculator -> GetSum ( ) <br> << std::endl;<br> */<br> sum = calculator -> GetSum ( ) ;<br> } <br>void GetOverlappingSum ( )<br> {<br><br> typedef itk::ImageFileReader < ImageVType > ImageReaderType;<br> ImageReaderType::Pointer reader = ImageReaderType::New();<br> reader -> SetFileName( "temp.mha" ) ;<br> reader -> Update();<br><br> ImageVType::Pointer image = reader -> GetOutput ( ) ;<br><br> ImageReaderType::Pointer reader1 = ImageReaderType::New();<br> reader1 -> SetFileName( "temp1.mha" ) ;<br> reader1 -> Update();<br> ImageVType::Pointer image1 = reader1 -> GetOutput ( ) ;<br><br> double sum;<br> GetSumS ( image, image1, sum ) ;<br> std::cout<<"overlapping sum from statistical calculator: "<br> <<sum<<std::endl;<br> GetSumV ( image, image1, sum ) ;<br> std::cout<<"overlapping sum from volume calculator: "<br> <<sum<<std::endl;<br> GetSumS ( image, image, sum ) ;<br> std::cout<<"sum of temp.mha from statistical calculator: "<br> <<sum<<std::endl;<br> GetSumV ( image, image, sum ) ;<br> std::cout<<"sum of temp.mha from volume calculator: "<br> <<sum<<std::endl;<br> GetSumS ( image1, image1, sum ) ;<br> std::cout<<"sum of temp1.mha from statistical calculator: "<br> <<sum<<std::endl;<br> GetSumV ( image1, image1, sum ) ;<br> std::cout<<"sum of temp1.mha from volume calculator: "<br> <<sum<<std::endl;<br> }<br><br><br>int main(int argc, char **argv)<br> {<br><br> GetOverlappingSum ( ) ;<br><br> return 1 ;<br> }<br><!-- footer --><br><br><br><br><br><div style="border-bottom:1px solid #999"></div><br>
<font color="black" style="font-size:14.8px">你 不 想 试 试 今 夏 最 “酷” 的 邮 箱 吗 ？</font>
<br>
<a href="http://www.126.com/" target="_blank" style="font-size:13px;line-height:160%;color:blue">蕴 涵 中 华 传 统 文 化 于 世 界 一 流 科 技 之 中，创 新 Ajax 技 术，126 “D 计 划”火 热 体 验 中 ！
</a>