[Insight-users] EigenVector Computation

[Arunachalam Kana]

Hi Sundar,

I use the eigenvalue analysis too. So i cross verified with my code whether you were missing something. My code is:

                typedef itk::SymmetricEigenAnalysis < HessianRecursiveGaussianFilterType::OutputImageType::PixelType, VectorPixelType, EV_PixelType > EigAnalysisType;
                //set eigen analysis
                EigAnalysisType eig;
                eig.SetDimension( 3 );
                eig.SetOrderEigenMagnitudes( true );

I found that the input matrix dimension was not set. Please check whether this is causing the error.
All the best.
Kana

From: insight-users-bounces at itk.org [mailto:insight-users-bounces at itk.org] On Behalf Of Sundaresan R
Sent: 22 February 2011 22:58
To: insight-users at itk.org
Subject: [Insight-users] EigenVector Computation

Hi

I'm computing the Hessians and then use SymmetricEigenAnalysis to compute the eigenVectors as follows:
I get a segmentation fault at the call to ComputeEigenValuesAndVectors method. Can anyone help me fix this issue?

Thanks,
Sundar

**************************************************************************************

#include "itkHessianRecursiveGaussianImageFilter.h"
#include "itkSymmetricEigenAnalysis.h"
#include "itkSymmetricEigenAnalysisImageFilter.h"
#include "itkSymmetricSecondRankTensor.h"
#include "itkCastImageFilter.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <iostream>


  const   unsigned int Dimension = 3;
  typedef double       InputPixelType;
  typedef double       OutputPixelType;
  typedef itk::Image< InputPixelType, Dimension >  InputImageType;
  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
  typedef itk::ImageFileReader< InputImageType >  ReaderType;
  typedef itk::ImageFileWriter< OutputImageType > WriterType;
  typedef itk::HessianRecursiveGaussianImageFilter< InputImageType > HessianFilterType;
  typedef itk::Vector< InputPixelType, Dimension > VectorType;
  typedef itk::Image< VectorType, Dimension > VectorImageType;
  typedef itk::Matrix< InputPixelType, Dimension, Dimension > MatrixType;
  typedef itk::Image< MatrixType, Dimension > MatrixImageType;
  typedef itk::ImageRegionConstIterator< MatrixImageType > MatrixIteratorType;
  typedef itk::ImageRegionIterator< MatrixImageType > OutputIteratorType;
  typedef itk::SymmetricEigenAnalysis< MatrixType, VectorType, MatrixType > EigenVectorAnalysisType;
  typedef MatrixImageType::Pointer        MatrixImagePointer;
  typedef MatrixImageType::RegionType     RegionType;
  typedef MatrixImageType::PointType      PointType;
  typedef MatrixImageType::SpacingType    SpacingType;
  typedef itk::ImageRegionIteratorWithIndex<HessianFilterType::OutputImageType>  HessianIteratorType;
  typedef itk::ImageFileWriter< MatrixImageType > MatrixWriterType;

int main( int argc, char *argv[] )
{
  if( argc < 3 )
    {
    std::cerr << "Usage: ./" << argv[0] << " inputImage sigma" << std::endl;
    }

  ReaderType::Pointer   reader = ReaderType::New();
  reader->SetFileName( argv[1] );

  HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
  hessianFilter->SetInput( reader->GetOutput() );
  if( argc >= 3 )
    hessianFilter->SetSigma( static_cast< double >( atof(argv[2])) );
  std::cout << "Calculating Hessian matrix" << std::endl;
  hessianFilter->Update();
  std::cout << "Done" << std::endl;

  EigenVectorAnalysisType  evap;
  HessianFilterType::OutputImageType::Pointer input = hessianFilter->GetOutput();
  RegionType region = input->GetLargestPossibleRegion();
  SpacingType spacing = input->GetSpacing();
  PointType origin = input->GetOrigin();

  MatrixType eigenMatrix, myMatrix;
  eigenMatrix[0][0] = 0; eigenMatrix[0][1]=0;eigenMatrix[0][2]=0;
  eigenMatrix[1][0] = 0; eigenMatrix[1][1]=0;eigenMatrix[1][2]=0;
  eigenMatrix[2][0] = 0; eigenMatrix[2][1]=0;eigenMatrix[2][2]=0;

  // Allocate the output image
  MatrixImagePointer m_Output = MatrixImageType::New();
  m_Output->SetRegions( region );
  m_Output->SetOrigin( origin );
  m_Output->SetSpacing( spacing );
  m_Output->Allocate();
  m_Output->FillBuffer( eigenMatrix );

  VectorType eigenVal;
  eigenVal[0]=eigenVal[1]=eigenVal[2]=0;
  HessianIteratorType mIt( input, region );
  OutputIteratorType oIt( m_Output, region );
  itk::SymmetricSecondRankTensor<double, 3> myTensor;
  mIt.GoToBegin();
  oIt.GoToBegin();
  while( !mIt.IsAtEnd() )
    {
      // Compute eigen values and eigen matrix
      evap.SetOrderEigenValues( 1 );
      myTensor = mIt.Get();
      myMatrix[0][0] = myTensor[0];
      myMatrix[0][1] = myTensor[1];
      myMatrix[0][2] = myTensor[2];
      myMatrix[1][0] = myTensor[1];
      myMatrix[1][1] = myTensor[3];
      myMatrix[1][2] = myTensor[4];
      myMatrix[2][0] = myTensor[2];
      myMatrix[2][1] = myTensor[4];
      myMatrix[2][2] = myTensor[5];
      std::cout << "myTensor " << myTensor << std::endl;
      std::cout << "myMatrix " << myMatrix << std::endl;
      std::cout << "eigenVal " << eigenVal << std::endl;
      std::cout << "eigenMatrix " << eigenMatrix << std::endl;
      evap.ComputeEigenValuesAndVectors( myMatrix, eigenVal, eigenMatrix );
      std::cout << " After ComputeEigenValuesAndVectors" << std::endl;
      oIt.Set( eigenMatrix );
      ++oIt;
      ++mIt;
    }
  MatrixWriterType::Pointer   writer = MatrixWriterType::New();
  writer->SetInput( m_Output );
  writer->SetFileName( "EigenMatrix.mhd" );
  try
  {
    writer->Update();
  }
  catch ( itk::ExceptionObject e )
  {
    std::cerr << "Error: " << e << std::endl;
  }
  return EXIT_SUCCESS;
}
***************************************************************************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110223/57f7ae89/attachment.htm>