itk::ERROR: MeanSampleFilter: Total frequency was too close to zero

Beginner Questions
1

Hello,

I’m quite new to C++ and ITK (I have mainly used Python and SimpleITK up until now).
I’m trying to read an image and get some statistics from this image, but only in an ellipse-shaped region.
For this, I am using the SpatialObjectToImageStatisticsCalculator (the code is added below).
However, I always get the following error:

.../Modules/Numerics/Statistics/include/itkMeanSampleFilter.hxx:182:
itk::ERROR: MeanSampleFilter(0x7fd1bed06c40): Total frequency was too close to zero: 0
Abort trap: 6

I have taken a look at ‘itkMeanSampleFilter.hxx’ and it seems that the frequency of the iterator stays 0.
But, with my limited knowledge of C++ and ITK, I am not sure what this means or what causes it.

Other things I’ve tried:

  • I’ve saved the ellipse to an image and checked that it is actually in the correct position (which it is).

  • I have tried the example in the ‘SpatialObjectToImageStatisticsCalculator.cxx’ file, where a random test image is created (not read from a file), and this works.

  • I have also tried saving this random test image and then reading it to use with the statistics calculator, but then I get the same error. So, it seems that something goes wrong in reading the file.

  • Searching online didn’t get me any further.

Can anyone help me understand what I’m doing wrong? Thanks in advance!
This is the code I used:

// Read image
typedef float PixelType;
const unsigned int Dimension = 3;
typedef itk::Image< PixelType, Dimension > ImageType;
typedef itk::ImageFileReader< ImageType > ReaderType;

ReaderType::Pointer reader = ReaderType::New();
const char * InputFilename = argv[1];
reader->SetFileName( InputFilename );
reader->UpdateOutputInformation();
ImageType::Pointer image = reader->GetOutput();


// Create ellipse-shape spatial object
itk::FixedArray< int , 3> radius_vx;
radius_vx[0] = 21; 
radius_vx[1] = 21; 
radius_vx[2] = 3; 

typedef itk::EllipseSpatialObject<3> EllipseType;
EllipseType::Pointer ellipse = EllipseType::New();
ellipse->SetRadius(radius_vx);
EllipseType::VectorType offset;
ellipse->ComputeObjectToParentTransform();


// Calculate statistics
typedef itk::SpatialObjectToImageStatisticsCalculator<ImageType, EllipseType> CalculatorType;
CalculatorType::Pointer calculator = CalculatorType::New();
calculator->SetImage(image);
calculator->SetSpatialObject(ellipse);
calculator->Update();

std::cout << "Sample mean = " << calculator->GetMean() << std::endl;
std::cout << "Sample covariance = " << calculator->GetCovarianceMatrix() << std::endl;
2

reader->UpdateOutputInformation(); only updates meta-information (size, spacing etc), it does not read pixel values. Try replacing this by reader->Update();.

3

Thanks for the reply!
I’ve changed it, but still get the same error.

4

Can you turn this into a runnable example? That will allow reproducing your problem easily.

5

After trying some more things, I noticed the problem is caused by the origin of the image.

The code below generates the same error, but when commenting out the change of the origin, it works fine.
Also, when I read in my actual image (instead of creating a random test image), and change the origin to [0,0,0], it works.

Any idea why an origin different from [0,0,0] would cause this error?

// Create image
typedef itk::Image< float, 3 >      ImageType;
typedef itk::RandomImageSource< ImageType > RandomImageSourceType;
RandomImageSourceType::Pointer randomImageSource
= RandomImageSourceType::New();
ImageType::SizeValueType size[3];
size[0] = 512;
size[1] = 512;
size[2] = 75;
randomImageSource->SetSize(size);
randomImageSource->Update();
ImageType::Pointer image = randomImageSource->GetOutput();

// This part causes the error /////////////////////////
ImageType::PointType origin;
origin[0] = -172;
origin[1] = 179;
origin[2] = -368;
image->SetOrigin( origin );
///////////////////////////////////////////////////////

// Create ellipse
itk::FixedArray<int, 3> radius_vx;
radius_vx[0] = 20;
radius_vx[1] = 20;
radius_vx[2] = 5;

// Create ellipse with radius in voxels
typedef itk::EllipseSpatialObject<3> EllipseType;
EllipseType::Pointer ellipse = EllipseType::New();
ellipse->SetRadius(radius_vx);
EllipseType::VectorType offset;
ellipse->ComputeObjectToParentTransform();

// Calculate coeffiecient of variation
typedef itk::SpatialObjectToImageStatisticsCalculator<ImageType, EllipseType> CalculatorType;
CalculatorType::Pointer calculator = CalculatorType::New();
calculator->SetImage(image);
calculator->SetSpatialObject(ellipse);
calculator->Update();

std::cout << "Sample mean = " << calculator->GetMean() << std::endl;
std::cout << "Sample covariance = " << calculator->GetCovarianceMatrix();
6

After a lot of trial and error I noticed I was confusing the world space and object space.

In case anyone else has the same issue (please correct me if I’m wrong):
Changing the origin of the image causes a misalignment between the ellipse (defined in object space) and the image (defined in world space).
Because they no longer overlap, there are no voxels from the image contained in the ellipse.

The ellipse can be moved using

ellipse->GetIndexToObjectTransform()->SetOffset(offset)

where ‘offset’ contains the index of the desired position in world coordinates, multiplied with the spacing and increased with the origin. This is the same as using TransformIndexToPhysicalPoint.

I haven’t found yet how to correctly compensate for a direction matrix that is not the identity matrix. Using TransformIndexToPhysicalPoint doesn’t seem to solve this either.

For now, I can fix this by setting the direction matrix of the image to the identity matrix.
However, it would be better to not have to do this, but be able to adapt the ellipse positioning coordinates correctly.

If anyone knows how to do this, please let me know.
Thanks!

7

You could construct image’s index to physical transform, and set it to your ellipse spatial object.

8

Thanks a lot!

I’ve tried many different things regarding the index to physical point transform, but whenever the image direction is different from the identity matrix, I can’t get the image and ellipse to align properly.

I ended up using a neighborhooditerator to extract the values in an elliptical neighborhood and calculate the statistics I need.