[Insight-users] Filtering a tensor component accessed with NthElementImageAdaptor

[Vicente Grau]

Dear all,

I am writing a filter to calculate the gradient structure tensor. This
involves calculating an NxN matrix for each voxel of the image, and then
applying a gaussian filter to each one of the components of the matrix (i.e.,
taking one component of the matrix at every voxel to create a new image, and
applying a gaussian filter on this new image). Based on
HessianRecursiveGausianImageFilter, I am using NthElementImageAdaptor and
three instances of RecursiveGaussianImageFilter to filter in the x, y and z
directions. As I saw in previous submissions to the list, I defined the
first of the gaussian filters to accept the adaptor as an input.
I have been getting errors in the results (basically, the filters behave as
if the pixel spacing were always [1 1 1], though it's not), so I wrote this
very simple piece of code to check that image spacing has the expected
values. I have attached the code below. It works by reading a tensor image
from a file, and applying a gaussian filter to the first component of the
tensor (of course, I would have to define things as filter direction, sigma,
etc., but here I just wanted to check the spacing). The input image is all
zeros, though I don't think this matters in this case. With an input image
with spacing [2 2 2], the program output is:

m_ImageAdaptor->GetSpacing(): [2, 2, 2]
m_ComponentFilter->GetOutput()->GetSpacing(): [1, 1, 1]

which is exactly the same error I got in my structure tensor calculation
program. It seems that the output of the filter has a different spacing from
the input, so when I attach a second filter to the output of the first one,
the result is wrong. Could anybody tell me what I am doing wrong here?
Thanks a lot,

Vicente

/**************************** Code follows ***********************/

#include "itkSymmetricSecondRankTensor.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkRecursiveGaussianImageFilter.h"
#include "itkNthElementImageAdaptor.h"



int main( int argc, char *argv[] )
{
    const     unsigned int        Dimension       = 3;
    typedef   double              InputPixelType;

    typedef itk::Image< itk::SymmetricSecondRankTensor< InputPixelType,
Dimension > , Dimension > InputImageType;
    typedef itk::Image< InputPixelType, Dimension >   OutputImageType;

    typedef   itk::ImageFileReader< InputImageType >  ReaderType;
    ReaderType::Pointer   reader = ReaderType::New();
   reader->SetFileName( argv[1] );

    typedef itk::NthElementImageAdaptor< InputImageType, InputPixelType >
OutputImageAdaptorType;
    OutputImageAdaptorType::Pointer  m_ImageAdaptor =
OutputImageAdaptorType::New();
    m_ImageAdaptor->SetImage( reader->GetOutput() );
    m_ImageAdaptor->Allocate();
    m_ImageAdaptor->SelectNthElement( 0 );

    typedef itk::RecursiveGaussianImageFilter< OutputImageAdaptorType,
OutputImageType >    ComponentFilterType;
    ComponentFilterType::Pointer m_ComponentFilter  =
ComponentFilterType::New();
    m_ComponentFilter->SetInput( m_ImageAdaptor );

    m_ComponentFilter->Update();

    std::cout << "m_ImageAdaptor->GetSpacing(): " <<
m_ImageAdaptor->GetSpacing() << std::endl;
    std::cout << "m_ComponentFilter->GetOutput()->GetSpacing(): " <<

m_ComponentFilter->GetOutput()->GetSpacing() << std::endl;

  return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20060505/35cca6e1/attachment.html