[Insight-developers] hard time with statistics framework (and also VectorImage and multicomponent pixels)

[Gaëtan Lehmann]

Hi,

I'm having a hard time with the statistics framework.
I don't have the feeling to try to do fancy things though: I want to  
look at the colocalization of the signal in two images, by looking at  
the 2D histogram generated from those two images. The histogram should  
be generated only in some part of the images given by a mask.

To do that, I put my two images in a VectorImage, and pass that  
VectorImage to an ImageToHistogramFilter. It doesn't build - seems to  
be a problem with the direct access to the PixelContainer in  
ImageToListSampleAdaptor which can't work with VectorImage.

Ok, so I switch to a standard image with a Vector pixel type. It  
builds, but fail at run time, because some parameter or not set:  
MarginalScale, HistogramSize, HistogramBinMinimum,  
HistogramBinMaximum. Come on! I just want my histogram! And I don't  
even know what is a marginal scale...

So I think I can help the future users by setting the default values  
in ImageToHistogramFilter. Quite easy for MarginalScale and  
AutoMinimumMaximum, but harder for HistogramSize, because the number  
of components may not be known at build time, in case of a  
VectorImage. I thought that this code can do it:

+  if( this->GetHistogramSizeInput() )
+    {
+    // user provided value
+    this->m_HistogramGenerator->SetHistogramSizeInput( this- 
 >GetHistogramSizeInput() );
+    }
+  else
+    {
+    // use a default value, which must be computed at run time for  
the VectorImage
+    HistogramSizeType numberOfBins( this->GetInput()- 
 >GetNumberOfComponentsPerPixel() );
+    numberOfBins.Fill(128);
+    this->m_HistogramGenerator->SetHistogramSize( numberOfBins );
+    }
+

No, it can't: GetNumberOfComponentsPerPixel() always return 1,  
excepted for VectorImage. I can't see any easy way to write that  
simple code without doing any kind of template specialization. OK,  
we'll probably have to do that anyway to fix the PixelContainer  
access. I'll provide the size I want by hand for now.
It builds and it run. Great! It looks quite slow though - something to  
investigate.

Then, to look at the produced histogram, I try to convert it an  
itk::Image with HistogramToIntensityImageFilter. The filter takes a  
pixel type and a dimension instead of the usual image type. Ok. Unlike  
all the other statistics classes, HistogramToIntensityImageFilter is  
not in the itk::Statistics namespace. Ok, I just remove  
the ::Statistics and it builds.
But doesn't run: segmentation error.

Ok, back in ccmake, switch to Debug build mode, rebuild, and run in  
gdb to get at stack trace.
The segmentation fault comes from the histogram->GetSize(i) call in  
HistogramToImageFilter::GenerateOutputInformation(). The histogram is  
not in the expected state in that method. Indeed, the histogram must  
be filled to be able to get it's size, so somehow,  
HistogramToImageFilter must update its input histogram to be able to  
generate its output information. Currently, those filters are not able  
to work properly in the pipeline.

Ok, so I'll update the histogram by hand for now, before running the  
HistogramToIntensityImageFilter. It works!

But I still have to figure how to constrain the histogram generation  
to a mask.

This is the short story of course - that was quite painful, and would  
probably be more painful for a new comer to ITK.
And I have no idea about how to do the same in python at this time,  
because the filters are not building with vector images.

So what I'd like to do (if I have some time):

* make GetNumberOfComponentsPerPixel() actually return the number of  
component per pixel, even if the image is not a vector image
* provide a default histogram size in ImageToHistogramFilter
* add the masking feature to ImageToHistogramFilter
* make ImageToListSampleAdaptor able to work with VectorImage
* maybe not use ImageToListSampleAdaptor and SampleToHistogramFilter  
internally in ImageToHistogramFilter, to try to speedup the execution
* make HistogramToImageFilter able to work with the pipeline
* check the dimension of the histogram before accessing it's size in  
HistogramToImageFilter
* make HistogramToImageFilter and subclasses take an Image as template  
parameter instead of the dimension and the type
* suppress the Decorated object layer in HistogramToImageFilter, now  
that Histogram is a subclass of DataObject
* remove the Statistics name space

Quite a big list of things to do!

Of course any comment or help is welcome.

Regards,

Gaëtan

-- 
Gaëtan Lehmann
Biologie du Développement et de la Reproduction
INRA de Jouy-en-Josas (France)
tel: +33 1 34 65 29 66    fax: 01 34 65 29 09
http://voxel.jouy.inra.fr  http://www.itk.org
http://www.mandriva.org  http://www.bepo.fr

-------------- next part --------------
A non-text attachment was scrubbed...
Name: PGP.sig
Type: application/pgp-signature
Size: 203 bytes
Desc: Ceci est une signature ?lectronique PGP
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20110429/8398bc37/attachment.pgp>