Hi Kana,<br><br>Yes, you are right, the dimension not being set was causing the problem.<br><br>Thank you very much.<br><br>Sundar<br><br><div class="gmail_quote">On Wed, Feb 23, 2011 at 4:05 AM, Arunachalam Kana <span dir="ltr"><<a href="mailto:Kana.Arunachalam@fh-wels.at">Kana.Arunachalam@fh-wels.at</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div link="blue" vlink="purple" lang="DE-AT"><div><p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);">Hi Sundar,</span></p>
<p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);"> </span></p><p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB">I use the eigenvalue analysis too. So i cross verified with my code whether you were missing something. My code is:</span></p>
<p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> </span></p><p class="MsoNormal" style=""><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> typedef itk::SymmetricEigenAnalysis < HessianRecursiveGaussianFilterType::OutputImageType::PixelType, VectorPixelType, EV_PixelType > EigAnalysisType;</span></p>
<p class="MsoNormal" style=""><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> //set eigen analysis</span></p><p class="MsoNormal" style=""><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> EigAnalysisType eig;</span></p>
<p class="MsoNormal" style=""><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> <span style="background: red none repeat scroll 0% 0%; -moz-background-clip: -moz-initial; -moz-background-origin: -moz-initial; -moz-background-inline-policy: -moz-initial;">eig.SetDimension( 3 );</span></span></p>
<p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> eig.SetOrderEigenMagnitudes( true );</span></p><p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> </span></p>
<p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB">I found that the input matrix dimension was not set. Please check whether this is causing the error.</span></p><p class="MsoNormal">
<span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB">All the best.</span></p><p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB">Kana</span></p><p class="MsoNormal"><span style="font-size: 11pt; color: rgb(31, 73, 125);" lang="EN-GB"> </span></p>
<div style="border-style: solid none none; border-color: rgb(181, 196, 223) -moz-use-text-color -moz-use-text-color; border-width: 1pt medium medium; padding: 3pt 0cm 0cm;"><p class="MsoNormal"><b><span style="font-size: 10pt;" lang="EN-US">From:</span></b><span style="font-size: 10pt;" lang="EN-US"> <a href="mailto:insight-users-bounces@itk.org" target="_blank">insight-users-bounces@itk.org</a> [mailto:<a href="mailto:insight-users-bounces@itk.org" target="_blank">insight-users-bounces@itk.org</a>] <b>On Behalf Of </b>Sundaresan R<br>
<b>Sent:</b> 22 February 2011 22:58<br><b>To:</b> <a href="mailto:insight-users@itk.org" target="_blank">insight-users@itk.org</a><br><b>Subject:</b> [Insight-users] EigenVector Computation</span></p></div><div><div></div>
<div class="h5"><p class="MsoNormal"><span lang="EN-GB"> </span></p><p class="MsoNormal">Hi<br><br>I'm computing the Hessians and then use SymmetricEigenAnalysis to compute the eigenVectors as follows:<br>I get a segmentation fault at the call to ComputeEigenValuesAndVectors method. Can anyone help me fix this issue?<br>
<br>Thanks,<br>Sundar<br><br>**************************************************************************************<br> <br>#include "itkHessianRecursiveGaussianImageFilter.h" <br>#include "itkSymmetricEigenAnalysis.h"<br>
#include "itkSymmetricEigenAnalysisImageFilter.h" <br>#include "itkSymmetricSecondRankTensor.h" <br>#include "itkCastImageFilter.h" <br>#include "itkImage.h" <br>#include "itkImageFileReader.h" <br>
#include "itkImageFileWriter.h" <br>#include <iostream><br><br> <br> const unsigned int Dimension = 3; <br> typedef double InputPixelType; <br> typedef double OutputPixelType; <br> typedef itk::Image< InputPixelType, Dimension > InputImageType; <br>
typedef itk::Image< OutputPixelType, Dimension > OutputImageType; <br> typedef itk::ImageFileReader< InputImageType > ReaderType; <br> typedef itk::ImageFileWriter< OutputImageType > WriterType; <br>
typedef itk::HessianRecursiveGaussianImageFilter< InputImageType > HessianFilterType; <br> typedef itk::Vector< InputPixelType, Dimension > VectorType;<br> typedef itk::Image< VectorType, Dimension > VectorImageType; <br>
typedef itk::Matrix< InputPixelType, Dimension, Dimension > MatrixType;<br> typedef itk::Image< MatrixType, Dimension > MatrixImageType;<br> typedef itk::ImageRegionConstIterator< MatrixImageType > MatrixIteratorType;<br>
typedef itk::ImageRegionIterator< MatrixImageType > OutputIteratorType;<br> typedef itk::SymmetricEigenAnalysis< MatrixType, VectorType, MatrixType > EigenVectorAnalysisType;<br> typedef MatrixImageType::Pointer MatrixImagePointer;<br>
typedef MatrixImageType::RegionType RegionType;<br> typedef MatrixImageType::PointType PointType;<br> typedef MatrixImageType::SpacingType SpacingType;<br> typedef itk::ImageRegionIteratorWithIndex<HessianFilterType::OutputImageType> HessianIteratorType;<br>
typedef itk::ImageFileWriter< MatrixImageType > MatrixWriterType;<br> <br>int main( int argc, char *argv[] ) <br>{ <br> if( argc < 3 ) <br> { <br> std::cerr << "Usage: ./" << argv[0] << " inputImage sigma" << std::endl; <br>
} <br> <br> ReaderType::Pointer reader = ReaderType::New(); <br> reader->SetFileName( argv[1] ); <br> <br> HessianFilterType::Pointer hessianFilter = HessianFilterType::New(); <br> hessianFilter->SetInput( reader->GetOutput() ); <br>
if( argc >= 3 ) <br> hessianFilter->SetSigma( static_cast< double >( atof(argv[2])) ); <br> std::cout << "Calculating Hessian matrix" << std::endl; <br> hessianFilter->Update(); <br>
std::cout << "Done" << std::endl; <br> <br> EigenVectorAnalysisType evap; <br> HessianFilterType::OutputImageType::Pointer input = hessianFilter->GetOutput();<br> RegionType region = input->GetLargestPossibleRegion();<br>
SpacingType spacing = input->GetSpacing();<br> PointType origin = input->GetOrigin();<br><br> MatrixType eigenMatrix, myMatrix;<br> eigenMatrix[0][0] = 0; eigenMatrix[0][1]=0;eigenMatrix[0][2]=0;<br> eigenMatrix[1][0] = 0; eigenMatrix[1][1]=0;eigenMatrix[1][2]=0;<br>
eigenMatrix[2][0] = 0; eigenMatrix[2][1]=0;eigenMatrix[2][2]=0;<br><br> // Allocate the output image<br> MatrixImagePointer m_Output = MatrixImageType::New();<br> m_Output->SetRegions( region );<br> m_Output->SetOrigin( origin );<br>
m_Output->SetSpacing( spacing );<br> m_Output->Allocate();<br> m_Output->FillBuffer( eigenMatrix );<br><br> VectorType eigenVal;<br> eigenVal[0]=eigenVal[1]=eigenVal[2]=0;<br> HessianIteratorType mIt( input, region );<br>
OutputIteratorType oIt( m_Output, region );<br> itk::SymmetricSecondRankTensor<double, 3> myTensor;<br> mIt.GoToBegin();<br> oIt.GoToBegin();<br> while( !mIt.IsAtEnd() )<br> {<br> // Compute eigen values and eigen matrix<br>
evap.SetOrderEigenValues( 1 );<br> myTensor = mIt.Get();<br> myMatrix[0][0] = myTensor[0];<br> myMatrix[0][1] = myTensor[1];<br> myMatrix[0][2] = myTensor[2];<br> myMatrix[1][0] = myTensor[1];<br>
myMatrix[1][1] = myTensor[3];<br> myMatrix[1][2] = myTensor[4];<br> myMatrix[2][0] = myTensor[2];<br> myMatrix[2][1] = myTensor[4];<br> myMatrix[2][2] = myTensor[5];<br> std::cout << "myTensor " << myTensor << std::endl;<br>
std::cout << "myMatrix " << myMatrix << std::endl; <br> std::cout << "eigenVal " << eigenVal << std::endl;<br> std::cout << "eigenMatrix " << eigenMatrix << std::endl; <br>
evap.ComputeEigenValuesAndVectors( myMatrix, eigenVal, eigenMatrix );<br> std::cout << " After ComputeEigenValuesAndVectors" << std::endl; <br> oIt.Set( eigenMatrix );<br> ++oIt;<br>
++mIt;<br> }<br> MatrixWriterType::Pointer writer = MatrixWriterType::New(); <br> writer->SetInput( m_Output );<br> writer->SetFileName( "EigenMatrix.mhd" ); <br> try<br> {<br> writer->Update();<br>
}<br> catch ( itk::ExceptionObject e )<br> {<br> std::cerr << "Error: " << e << std::endl;<br> }<br> return EXIT_SUCCESS; <br>} <br>***************************************************************************************************</p>
</div></div></div></div></blockquote></div><br>