[Insight-users] itkSpatialObjectToImageStatisticsCalculator problem

[Benjamin King]

Hello all,

I'd like to compute the mean and standard deviation of a part of an image. 
That part is described by a sphere i.e. by the center's coordinates and a 
radius.
The following program is what I came up with, but it behaves strangely. 
Most of the time, the output is "-1.#IND 0", sometimes it looks like 
something sensible, but I can't be sure about that. If I have a working set 
of input parameters, little variations of x, y or z lead to desaster aka "- 
1.#IND 0". I can vary the radius without that effect.

Also, I think that statistics->Update() should Update() the imageReader and 
ellipse, but when I start the program with a filename that doesn't exist, 
there is no complaint.

Any help/ideas are much appreciated, here is my code:

---
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkSpatialObject.h"
#include "itkEllipseSpatialObject.h"
#include "itkSpatialObjectToImageStatisticsCalculator.h"
#include "itkAffineTransform.h"
#include <iostream>

using std::cout;
using std::endl;

int main(int argc, char** argv)
{
#define DIMENSION 3
  if (argc != 6) {
    cout << "\nVolume Statistics\n\n"
      "This program computes the mean and standard deviation of \n"
      "a ball shaped region in a " << DIMENSION << " dimensional 
image.\n\n"
      "USAGE:\n"
      "  volumestatistics <input file> <x> <y> <z> <radius>\n"
      "    <input file>: " << DIMENSION<< "d 16bit Analyze/Dicom file\n"
      "    <x>:          x coordinate of the ball's center\n"
      "    <y>:          y coordinate of the ball's center\n"
      "    <z>:          z coordinate of the ball's center\n"
      "    <radius>:     radius of the ball\n\n"
      "HINT: Define the center coordinate and the radius with respect to 
the\n"
      "  spacing and origin of your input image.\n\n"
      "OUTPUT:\n"
      "  A single line with two floating point numbers <mean> and <standard 
deviation>\n"
      "  separated by a space"
      << endl;

    return 1;
  }
  typedef itk::Image< unsigned short, DIMENSION > ImageType;
  typedef itk::ImageFileReader< ImageType > ReaderType;
  ReaderType::Pointer imageReader = ReaderType::New();
  imageReader->SetFileName(argv[1]);

  // generate a ball centered around (0, 0, 0)
  typedef itk::EllipseSpatialObject< DIMENSION> EllipseType;
  EllipseType::Pointer ellipse = EllipseType::New();
  ellipse->SetSpacing(imageReader->GetOutput()->GetSpacing());
  ellipse->SetRadius(atof(argv[5]));
  EllipseType::VectorType offset;
  offset[0] = atof(argv[2]);
  offset[1] = atof(argv[3]);
  offset[2] = atof(argv[4]);
  ellipse->GetIndexToObjectTransform()->SetOffset(offset);
  ellipse->ComputeObjectToParentTransform();

  // compute the statistics
  typedef itk::SpatialObjectToImageStatisticsCalculator< ImageType, 
EllipseType > StatisticsType;
  StatisticsType::Pointer statistic = StatisticsType::New();
  statistic->SetImage(imageReader->GetOutput());
  statistic->SetSpatialObject(ellipse);

  imageReader->Update();
  ellipse->Update();
  statistic->Update();

  cout << statistic->GetMean()[0] << " " << sqrt(statistic- 
>GetCovarianceMatrix()[0][0]) << endl;

  return 0;
}
---


Best regards,
  Benjamin


-- 
Benjamin King
Institut für Medizinische Informatik
Medizinische Hochschule Hannover
Tel.: +49  511  532-2663