[Insight-users] A minor issue in SpatialObjectToImageStat isticsCalculator

[Julien Jomier]

Hi Kurt,

I just put a fix in the cvs repository. Please do a cvs update and let
us know if it works for you or not.

Thanks for the report,

Julien

kurt wrote:
> Hi guys,
> 
> I am using SpatialObjectToImageStatisticsCalculator to generate the sum 
> of pixel values. but I got "0". the SO is a ImageMaskSpatialObject.
> 
> This could be fixed real quick because the sum can be generated by 
> calculator->GetMean() * calculator -> GetNumberOfPixels ( ), which works 
> fine. The code and data are included.
> 
> Thanks !
> 
> Kurt
> 
> #if defined(_MSC_VER)
> #pragma warning ( disable : 4786 )
> #endif
> 
> #include <iostream>
> 
> 
> #include <itkImageFileWriter.h>
> #include <itkImageFileReader.h>
> #include <itkThresholdImageFilter.h>
> #include <itkImageMaskSpatialObject.h>
> #include <itkSpatialObjectToImageStatisticsCalculator.h>
> 
> #include <itkEllipseSpatialObject.h>
> #include "itkEllipseboundaryToImageFilter.h"
> 
> 
> 
> int itkEllipseBoundaryToImageFilterTest( )
> {
> typedef itk::EllipseSpatialObject<3> EllipseType;
> 
> EllipseType::Pointer ellipse = EllipseType::New();
> ellipse->SetRadius(10);
> 
> // Center the circle in the image
> EllipseType::TransformType::OffsetType offset;
> offset.Fill(25);
> ellipse->GetObjectToParentTransform()->SetOffset(offset);
> ellipse->ComputeObjectToWorldTransform();
> 
> EllipseType::Pointer ellipse1 = EllipseType::New();
> ellipse1->SetRadius(15);
> 
> // Center the circle in the image
> EllipseType::TransformType::OffsetType offset1;
> offset1.Fill(55);
> ellipse1->GetObjectToParentTransform()->SetOffset(offset1);
> ellipse1->ComputeObjectToWorldTransform();
> 
> typedef itk::Image<double,3> ImageType;
> 
> typedef itk::EllipseBoundaryToImageFilter<3,ImageType> 
> EllipseBoundaryToImageFilterType;
> EllipseBoundaryToImageFilterType::Pointer imageFilter = 
> EllipseBoundaryToImageFilterType::New();
> imageFilter->SetInput(0, ellipse);
> imageFilter->SetInput(1, ellipse1);
> imageFilter->SetInsideValue(2);
> imageFilter->GetInsideValue();
> imageFilter->SetOutsideValue(0);
> imageFilter->GetOutsideValue();
> imageFilter->SetChildrenDepth(1);
> imageFilter->GetChildrenDepth();
> ImageType::SizeType size;
> size[0]=80;
> size[1]=80;
> size[2]=80;
> imageFilter->SetSize(size);
> 
> // Testing spacing
> std::cout << "Testing Spacing: ";
> 
> float spacing_float[3];
> double spacing_double[3];
> 
> for(unsigned int i=0;i<3;i++)
> {
> spacing_float[i]=1.0;
> spacing_double[i]=1.0;
> }
> imageFilter->SetSpacing(spacing_float);
> imageFilter->SetSpacing(spacing_double);
> const double* spacing_result = imageFilter->GetSpacing();
> 
> for(unsigned int i=0;i<3;i++)
> {
> if(spacing_result[i]!=1.0)
> {
> std::cout << "[FAILURE]" << std::endl;
> return EXIT_FAILURE;
> }
> }
> 
> std::cout << "[PASSED]" << std::endl;
> 
> // Testing Origin
> std::cout << "Testing Origin: ";
> 
> float origin_float[3];
> double origin_double[3];
> 
> for(unsigned int i=0;i<3;i++)
> {
> origin_float[i]=0.0;
> origin_double[i]=0.0;
> }
> imageFilter->SetOrigin(origin_float);
> imageFilter->SetOrigin(origin_double);
> const double* origin_result = imageFilter->GetOrigin();
> 
> for(unsigned int i=0;i<3;i++)
> {
> if(origin_result[i]!=0.0)
> {
> std::cout << "[FAILURE]" << std::endl;
> return EXIT_FAILURE;
> }
> }
> 
> std::cout << "[PASSED]" << std::endl;
> 
> // Testing PrintSelf
> std::cout << imageFilter << std::endl;
> 
> //Update the filter
> imageFilter->Update();
> 
> ImageType::Pointer image = imageFilter->GetOutput();
> 
> std::cout << "Testing Output Image: ";
> 
> typedef itk::ImageFileWriter< ImageType > WriterType;
> WriterType::Pointer writer = WriterType::New();
> writer->SetInput( image );
> writer -> UseCompressionOn ( ) ;
> writer->SetFileName("temp3.mhd");
> writer->Update();
> 
> ImageType::IndexType index;
> //The center should be outsidevalue
> index [ 0 ] = 25;
> index [ 1 ] = 25;
> index [ 2 ] = 25 ;
> 
> if(image->GetPixel(index) == 2.0)
> {
> std::cout << "[FAILURE]" << std::endl;
> return EXIT_FAILURE;
> }
> index [ 0 ] = 27;
> index [ 1 ] = 19;
> index [ 2 ] = 18 ;
> 
> //Test one boundary point
> if(image->GetPixel(index) != 2.0)
> {
> std::cout << "[FAILURE]" << std::endl;
> return EXIT_FAILURE;
> }
> std::cout << "[PASSED]" << std::endl;
> 
> 
> std::cout << "[PASSED]" << std::endl;
> std::cout << "Test [DONE]" << std::endl;
> 
> return 1;
> }
> int TestSampler( )
> {
> typedef double PixelType ;
> typedef unsigned char OverlayPixelType ;
> typedef itk::Image < PixelType, 3 > ImageType;
> typedef itk::Image < OverlayPixelType, 3 > OverlayType;
> 
> typedef itk::ImageFileReader < ImageType > ImageReaderType;
> ImageReaderType::Pointer imageReader = ImageReaderType::New();
> imageReader -> SetFileName( "image.mhd" ) ;
> try
> {
> imageReader -> Update();
> }
> catch( ... )
> {
> std::cout << "Problems reading image file " <<
> imageReader-> GetFileName ( ) << std::endl;
> return EXIT_FAILURE ;
> }
> 
> typedef itk::ImageFileReader < OverlayType > OverlayReaderType;
> OverlayReaderType::Pointer overlayReader = OverlayReaderType::New();
> overlayReader -> SetFileName( "overlay.mhd" ) ;
> try
> {
> overlayReader -> Update();
> }
> catch( ... )
> {
> std::cout << "Problems reading overlay file " <<
> overlayReader-> GetFileName ( ) << std::endl;
> return EXIT_FAILURE;
> }
> 
> ImageType::Pointer image = imageReader -> GetOutput ( ) ;
> 
> OverlayType::Pointer overlay = overlayReader -> GetOutput ( ) ;
> 
> ImageType::IndexType indexInside ;
> indexInside [ 0 ] = 140 ;
> indexInside [ 1 ] = 122 ;
> indexInside [ 2 ] = 24 ;
> 
> typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;
> ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();
> maskSO->SetImage( overlay );
> //maskSO -> Print ( std::cout ) ;
> ImageMaskSpatialObject::PointType inside ;
> ImageMaskSpatialObject::PointType outside ;
> for ( int i = 0 ; i < 3 ; i ++)
> {
> inside [ i ] = indexInside [ i ] *
> image -> GetSpacing ( ) [ i ] ;
> }
> outside [ 0 ] = 197 ;
> outside [ 1 ] = 65 ;
> outside [ 2 ] = 14 ;
> std::cout << "Is my point " << indexInside << " inside my mask image? "
> << (int ) overlay -> GetPixel ( indexInside ) << std::endl;
> std::cout << "Is my point " << inside << " inside my mask? "
> << maskSO->IsInside(inside) << std::endl;
> std::cout << "Is my point " << outside << " outside my mask? "
> << !maskSO->IsInside(outside) << std::endl;
> 
> /*
> typedef itk::EllipseSpatialObject<3> EllipseType;
> EllipseType::Pointer ellipse = EllipseType::New();
> ellipse->SetRadius(10);
> EllipseType::VectorType offset;
> for ( int i = 0 ; i < 3 ; i ++)
> {
> offset[ i ] = indexInside [ i ] *
> image -> GetSpacing ( ) [ i ] ;
> }
> ellipse->GetIndexToObjectTransform()->SetOffset(offset);
> ellipse->ComputeObjectToParentTransform();
> */
> typedef itk::SpatialObjectToImageStatisticsCalculator<
> // ImageType, EllipseType> CalculatorType;
> ImageType, ImageMaskSpatialObject> CalculatorType;
> CalculatorType::Pointer calculator = CalculatorType::New();
> calculator->SetImage(image);
> calculator->SetSpatialObject( maskSO );
> //calculator->SetSpatialObject( ellipse );
> calculator->Update();
> //calculator->Print ( std::cout );
> std::cout << "Sample mean = " << calculator->GetMean() << std::endl ;
> std::cout << "Sample covariance = " << calculator->GetCovarianceMatrix();
> std::cout << "Sample sum = " << calculator->GetSum()<<std::endl;
> std::cout << "Sample num of pixels = " << calculator->GetNumberOfPixels()
> <<std::endl;
> std::cout << "Sample calculated sum = " <<
> calculator->GetMean() * calculator -> GetNumberOfPixels ( ) <<std::endl;
> return EXIT_SUCCESS ;
> }
> int main(int argc, char **argv)
> {
> //itkEllipseBoundaryToImageFilterTest ( ) ;
> TestSampler ( ) ;
> 
> return 1 ;
> }
> 
> 
> 
> 
> 
> 
> 
> 
> 你 不 想 试 试 今 夏 最 “酷” 的 邮 箱 吗 ？
> 蕴 涵 中 华 传 统 文 化 于 世 界 一 流 科 技 之 中，创 新 Ajax 技 术， 
> 126 “D 计 划”火 热 体 验 中 ！ <http://www.126.com/>
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users