[Insight-users] Problem with
SpatialObjectToImageStatisticsCalculator
Alain CORON
alain.coron at lip.bhdc.jussieu.fr
Wed Jun 14 17:32:45 EDT 2006
Hello,
I would like to extract the mean of a region of an image. So I try to use
SpatialObjectToImageStatisticsCalculator.
My spatial objects are an ellipse (to extract a circular region) or a plane
(to extract a square region). But sometimes the mean is not correctly
computed.
In order to illustrate the problem, I have written a program that creates a
constant image of 100x100 pixels. Then a plane spatial object and an
ellipse spatial object centered in the image are created. The spacing of
the image is chosen on the command line as well as the length of the edge of
the square (which is also the diameter of the circle).
Then the mean of the region is computed.
At the end the constant image and the two spatial objects (with a
SpatialObjectToImageFilter) are saved. They are written as VTK files and
their filenames are prefixed with the 2nd command line argument.
If I run the program on my PC
./ToyProgram Toy_ 0.1 2
with
Image spacing : 0.1
edge of the square, or diameter : 2
the results are
imageCenterPosition : [5, 5]
1st square corner : [3.95, 3.95]
2nd square corner : [6.05, 6.05]
Is the center of the image inside
- the PlaneSpatialObject : 1
- the EllipseSpatialObject : 1
Is (0,0) inside
- the PlaneSpatialObject : 0
- the EllipseSpatialObject : 0
Number of pixels in the region
- PlaneCalculator : 441
- EllipseCalculator : 305
Sample mean
- the PlaneCalculator : [1]
- the EllipseCalculator : [1]
ok no problem.
But if
Image spacing : 0.01
edge of the square, or diameter : 0.2
then we have,
imageCenterPosition : [0.5, 0.5]
1st square corner : [0.395, 0.395]
2nd square corner : [0.605, 0.605]
Is the center of the image inside
- the PlaneSpatialObject : 1
- the EllipseSpatialObject : 1
Is (0,0) inside
- the PlaneSpatialObject : 0
- the EllipseSpatialObject : 0
Number of pixels in the region
- PlaneCalculator : 0 <----
- EllipseCalculator : 0 <----
Sample mean
- the PlaneCalculator : [nan] <----
- the EllipseCalculator : [nan] <----
The images written on disk are correct but the mean is wrong.
If you try
Image spacing : 0.1
edge of the square, or diameter : 2
or higher values the results are correct but fails with
Image spacing : 0.01
edge of the square, or diameter : 0.2
The version of ITK I use is 2.6.
Is there something wrong in this program ?
Thanks.
//////// Start of ToyProgram.cxx
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
#include <iomanip>
#include "itkNumericTraits.h"
#include "itkImage.h"
#include "itkRandomImageSource.h"
#include "itkImageFileWriter.h"
#include "itkPlaneSpatialObject.h"
#include "itkEllipseSpatialObject.h"
#include "itkSpatialObjectToImageStatisticsCalculator.h"
#include "itkSpatialObjectToImageFilter.h"
typedef float PixelType; // Operations
const unsigned int Dimension = 2;
typedef itk::Image<PixelType, Dimension> ImageType;
typedef itk::ImageFileWriter<ImageType> WriterType;
typedef itk::EllipseSpatialObject<Dimension> EllipseSpatialObjectType;
typedef itk::PlaneSpatialObject<Dimension> PlaneSpatialObjectType;
typedef itk::SpatialObjectToImageStatisticsCalculator<ImageType, EllipseSpatialObjectType > EllipseCalculatorType;
typedef itk::SpatialObjectToImageFilter<EllipseSpatialObjectType, ImageType> EllipseSpatialObjectToImageFilterType;
typedef itk::SpatialObjectToImageStatisticsCalculator<ImageType, PlaneSpatialObjectType > PlaneCalculatorType;
typedef itk::SpatialObjectToImageFilter<PlaneSpatialObjectType, ImageType> PlaneSpatialObjectToImageFilterType;
typedef ImageType::IndexType ImageIndexType;
typedef ImageType::SpacingType ImageSpacingType;
typedef ImageType::RegionType ImageRegionType;
typedef ImageType::IndexType ImageIndexType;
typedef ImageType::SizeType ImageSizeType;
typedef ImageType::PointType PointType;
int main(int argc, char* argv[])
{
if( argc < 3)
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " spacing outputFileNamePrefix length"
<< std::endl;
return EXIT_FAILURE;
}
int ArgumentNumber = 0;
const char* outputFileNamePrefix = argv[++ArgumentNumber];
const float spacing = atof (argv[++ArgumentNumber]);
const float length = atof (argv[++ArgumentNumber]);
std::cout << "Image spacing : " << spacing << std::endl;
std::cout << "edge of the square, or diameter : " << length << std::endl;
ImageSizeType imageSize;
imageSize[0] = 100;
imageSize[1] = 100;
ImageSpacingType imageSpacing;
imageSpacing[0] = spacing;
imageSpacing[1] = spacing;
ImageType::IndexType start;
start[0] = 0;
start[1] = 0;
ImageType::Pointer image = ImageType::New();
image->SetSpacing (imageSpacing);
ImageType::RegionType region;
region.SetSize (imageSize);
region.SetIndex (start);
image->SetRegions (region);
image->Allocate ();
image->FillBuffer (itk::NumericTraits<PixelType>::One);
// Find the center of the image
ImageIndexType imageCenterIndex;
for (int i = 0 ; i < Dimension; ++i)
{
imageCenterIndex[i] = imageSize[i]/2;
}
PointType imageCenterPosition;
image->TransformIndexToPhysicalPoint (imageCenterIndex,
imageCenterPosition);
PointType LLCornerPosition, UUCornerPosition, zeroPosition;
for (int i = 0; i < Dimension; ++i)
{
LLCornerPosition[i] = imageCenterPosition[i] - length/2 - imageSpacing[i]/2;
UUCornerPosition[i] = imageCenterPosition[i] + length/2 + imageSpacing[i]/2;
zeroPosition[i] = 0.0;
}
std::cout << "imageCenterPosition : " << imageCenterPosition << std::endl;
std::cout << "1st square corner : " << LLCornerPosition << std::endl;
std::cout << "2nd square corner : " << UUCornerPosition << std::endl;
// Define a plane spatial object
PlaneSpatialObjectType::Pointer plane = PlaneSpatialObjectType::New();
plane->SetLowerPoint (LLCornerPosition);
plane->SetUpperPoint (UUCornerPosition);
PlaneCalculatorType::Pointer planeCalculator = PlaneCalculatorType::New();
planeCalculator->SetSpatialObject (plane);
planeCalculator->SetImage (image);
planeCalculator->Update ();
// Define an ellipse spatial object
EllipseSpatialObjectType::Pointer ellipse = EllipseSpatialObjectType::New();
ellipse->SetRadius (length/2);
EllipseSpatialObjectType::VectorType offset;
offset[0] = imageCenterPosition[0];
offset[1] = imageCenterPosition[1];
ellipse->GetIndexToObjectTransform()->SetOffset (offset);
ellipse->ComputeObjectToParentTransform ();
EllipseCalculatorType::Pointer ellipseCalculator = EllipseCalculatorType::New();
ellipseCalculator->SetSpatialObject (ellipse);
ellipseCalculator->SetImage (image);
ellipseCalculator->Update ();
std::cout << "Is the center of the image inside " << std::endl
<< " - the PlaneSpatialObject : "
<< plane->IsInside (imageCenterPosition) << std::endl
<< " - the EllipseSpatialObject : "
<< ellipse->IsInside (imageCenterPosition) << std::endl;
std::cout << "Is (0,0) inside " << std::endl
<< " - the PlaneSpatialObject : "
<< plane->IsInside (zeroPosition) << std::endl
<< " - the EllipseSpatialObject : "
<< ellipse->IsInside (zeroPosition) << std::endl;
std::cout << "Number of pixels in the region " << std::endl
<< " - PlaneCalculator : " << planeCalculator->GetNumberOfPixels()
<< std::endl
<< " - EllipseCalculator : " << ellipseCalculator->GetNumberOfPixels()
<< std::endl;
std::cout << "Sample mean " << std::endl
<< " - the PlaneCalculator : "
<< planeCalculator->GetMean() << std::endl
<< " - the EllipseCalculator : "
<< ellipseCalculator->GetMean() << std::endl;
PlaneSpatialObjectToImageFilterType::Pointer planeSpatialObjectToImageFilter = PlaneSpatialObjectToImageFilterType::New();
planeSpatialObjectToImageFilter->SetInput (plane);
planeSpatialObjectToImageFilter->SetOrigin (image->GetOrigin());
planeSpatialObjectToImageFilter->SetSpacing (image->GetSpacing());
planeSpatialObjectToImageFilter->SetSize (image->GetLargestPossibleRegion().GetSize());
WriterType::Pointer writer1 = WriterType::New();
std::ostringstream outputFileName1;
outputFileName1 << outputFileNamePrefix
<< "PlaneSpatialImageFilter.vtk";
writer1->SetFileName (outputFileName1.str().c_str());
writer1->SetInput (planeSpatialObjectToImageFilter->GetOutput());
writer1->Update ();
EllipseSpatialObjectToImageFilterType::Pointer ellipseSpatialObjectToImageFilter = EllipseSpatialObjectToImageFilterType::New();
ellipseSpatialObjectToImageFilter->SetInput (ellipse);
ellipseSpatialObjectToImageFilter->SetOrigin (image->GetOrigin());
ellipseSpatialObjectToImageFilter->SetSpacing (image->GetSpacing());
ellipseSpatialObjectToImageFilter->SetSize (image->GetLargestPossibleRegion().GetSize());
WriterType::Pointer writer2 = WriterType::New();
std::ostringstream outputFileName2;
outputFileName2 << outputFileNamePrefix
<< "EllipseSpatialImageFilter.vtk";
writer2->SetFileName (outputFileName2.str().c_str());
writer2->SetInput (ellipseSpatialObjectToImageFilter->GetOutput());
writer2->Update ();
WriterType::Pointer writer3 = WriterType::New();
std::ostringstream outputFileName3;
outputFileName3 << outputFileNamePrefix
<< "Image.vtk";
writer3->SetFileName (outputFileName3.str().c_str());
writer3->SetInput (image);
writer3->Update ();
}
--
Alain CORON
More information about the Insight-users
mailing list