Hi<br><br>I&#39;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 &quot;itkHessianRecursiveGaussianImageFilter.h&quot;
<br>#include &quot;itkSymmetricEigenAnalysis.h&quot;<br>#include &quot;itkSymmetricEigenAnalysisImageFilter.h&quot;
<br>#include &quot;itkSymmetricSecondRankTensor.h&quot;
<br>#include &quot;itkCastImageFilter.h&quot;
<br>#include &quot;itkImage.h&quot;
<br>#include &quot;itkImageFileReader.h&quot;
<br>#include &quot;itkImageFileWriter.h&quot;
<br>#include &lt;iostream&gt;<br><br> <br>  const   unsigned int Dimension = 3;
<br>  typedef double       InputPixelType;
<br>  typedef double       OutputPixelType;
<br>  typedef itk::Image&lt; InputPixelType, Dimension &gt;  InputImageType;
<br>  typedef itk::Image&lt; OutputPixelType, Dimension &gt; OutputImageType;
<br>  typedef itk::ImageFileReader&lt; InputImageType &gt;  ReaderType;
<br>  typedef itk::ImageFileWriter&lt; OutputImageType &gt; WriterType;
<br>  typedef itk::HessianRecursiveGaussianImageFilter&lt; InputImageType &gt; HessianFilterType;
<br>  typedef itk::Vector&lt; InputPixelType, Dimension &gt; VectorType;<br>  typedef itk::Image&lt; VectorType, Dimension &gt; VectorImageType; <br>  typedef itk::Matrix&lt; InputPixelType, Dimension, Dimension &gt; MatrixType;<br>
  typedef itk::Image&lt; MatrixType, Dimension &gt; MatrixImageType;<br>  typedef itk::ImageRegionConstIterator&lt; MatrixImageType &gt; MatrixIteratorType;<br>  typedef itk::ImageRegionIterator&lt; MatrixImageType &gt; OutputIteratorType;<br>
  typedef itk::SymmetricEigenAnalysis&lt; MatrixType, VectorType, MatrixType &gt; 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&lt;HessianFilterType::OutputImageType&gt;  HessianIteratorType;<br>
  typedef itk::ImageFileWriter&lt; MatrixImageType &gt; MatrixWriterType;<br> <br>int main( int argc, char *argv[] )
<br>{
<br>  if( argc &lt; 3 )
<br>    {
<br>    std::cerr &lt;&lt; &quot;Usage: ./&quot; &lt;&lt; argv[0] &lt;&lt; &quot; inputImage sigma&quot; &lt;&lt; std::endl;
<br>    }
<br>  
<br>  ReaderType::Pointer   reader = ReaderType::New();
<br>  reader-&gt;SetFileName( argv[1] );
<br>    
<br>  HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
<br>  hessianFilter-&gt;SetInput( reader-&gt;GetOutput() );
<br>  if( argc &gt;= 3 )
<br>    hessianFilter-&gt;SetSigma( static_cast&lt; double &gt;( atof(argv[2])) );
<br>  std::cout &lt;&lt; &quot;Calculating Hessian matrix&quot; &lt;&lt; std::endl;
<br>  hessianFilter-&gt;Update();
<br>  std::cout &lt;&lt; &quot;Done&quot; &lt;&lt; std::endl;
<br>   <br>  EigenVectorAnalysisType  evap; <br>  HessianFilterType::OutputImageType::Pointer input = hessianFilter-&gt;GetOutput();<br>  RegionType region = input-&gt;GetLargestPossibleRegion();<br>  SpacingType spacing = input-&gt;GetSpacing();<br>
  PointType origin = input-&gt;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-&gt;SetRegions( region );<br>  m_Output-&gt;SetOrigin( origin );<br>
  m_Output-&gt;SetSpacing( spacing );<br>  m_Output-&gt;Allocate();<br>  m_Output-&gt;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&lt;double, 3&gt; 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 &lt;&lt; &quot;myTensor &quot; &lt;&lt; myTensor &lt;&lt; std::endl;<br>
      std::cout &lt;&lt; &quot;myMatrix &quot; &lt;&lt; myMatrix &lt;&lt; std::endl;  <br>      std::cout &lt;&lt; &quot;eigenVal &quot; &lt;&lt; eigenVal &lt;&lt; std::endl;<br>      std::cout &lt;&lt; &quot;eigenMatrix &quot; &lt;&lt; eigenMatrix &lt;&lt; std::endl;     <br>
      evap.ComputeEigenValuesAndVectors( myMatrix, eigenVal, eigenMatrix );<br>      std::cout &lt;&lt; &quot; After ComputeEigenValuesAndVectors&quot; &lt;&lt; std::endl; <br>      oIt.Set( eigenMatrix );<br>      ++oIt;<br>
      ++mIt;<br>    }<br>  MatrixWriterType::Pointer   writer = MatrixWriterType::New(); <br>  writer-&gt;SetInput( m_Output );<br>  writer-&gt;SetFileName( &quot;EigenMatrix.mhd&quot; );  <br>  try<br>  {<br>    writer-&gt;Update();<br>
  }<br>  catch ( itk::ExceptionObject e )<br>  {<br>    std::cerr &lt;&lt; &quot;Error: &quot; &lt;&lt; e &lt;&lt; std::endl;<br>  }<br>  return EXIT_SUCCESS;
<br>}
<br>***************************************************************************************************<br>