[Insight-users] EigenVector Computation

[Sundaresan R]

Hi Kana,

Yes, you are right, the dimension not being set was causing the problem.

Thank you very much.

Sundar

On Wed, Feb 23, 2011 at 4:05 AM, Arunachalam Kana <
Kana.Arunachalam at fh-wels.at> wrote:

> 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/f74c5124/attachment-0001.htm>