還沒有親測，先放著，因為要用這個......

Eigenvalues of Hessian Matrix are used as a criteria for line-like structure such as pulmonary vessels and fungi vessels.ITK provides a class called ：itk::SymmetricEigenAnalysis to extract eigenvalues, eigenvalues vector as well as eigenvalues matrix from a given matrix. However, the ITK itself provided no direct implementation for the computation of Hessian Matrix. Instead, the ITK package provides a HessianRecursiveGaussianImageFilter

Iv´an Macıa explained[1] that the direct computation of the approximation of Hessian matrix for discrete signal would only enhance the noise component.  In order to overcome such disadvantage, the author introduced a new convolution operator--the Gaussian Derivative Operator.

Though the document provides detailed description of the new class, it is not a easy task for use. Based on a example provided at <http://itk-insight-users.2283740.n2.nabble.com/EigenVector-Computation-td6054230.html>, I implemented a class to extract the eigenvalues vector of each voxel of a CT image.

============================================================

#include <iostream>
#include <string>
#include "stdlib.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkHessianRecursiveGaussianImageFilter.h"
#include "itkCastImageFilter.h"

const unsigned short dimension = 3;
typedef short VoxelType;
typedef itk::Image<VoxelType, dimension> InputImageType;
typedef double HessianVoxelType;
typedef itk::Image<HessianVoxelType, dimension> HessianInnerImageType;
typedef itk::CastImageFilter<InputImageType, HessianInnerImageType> CastFilterType;
typedef itk::HessianRecursiveGaussianImageFilter<HessianInnerImageType> HFilterType;

typedef itk::Vector<HessianVoxelType, dimension> VecVoxType;
typedef itk::Matrix<HessianVoxelType, dimension, dimension> MatVoxType;
typedef itk::Image<VecVoxType, dimension> VecEigHImageType;
typedef itk::Image<MatVoxType, dimension> MatEigHImageType;
typedef itk::ImageRegionIterator<VecEigHImageType> OutputIterType;
typedef itk::SymmetricEigenAnalysis<MatVoxType, VecVoxType, MatVoxType> EigValAnalysisType;

typedef MatEigHImageType::Pointer MatEigHImagePointerType;
typedef MatEigHImageType::RegionType MatRegionType;
typedef MatEigHImageType::PointType MatPointType;
typedef MatEigHImageType::SpacingType MatSpacingType;
typedef VecEigHImageType::Pointer VecEigHImagePointerType;
typedef itk::ImageRegionIteratorWithIndex<HFilterType::OutputImageType> HIterType;

class EigValHessian {

private:

    MatRegionType region;
    MatSpacingType spacing;
    MatPointType origin;

    CastFilterType::Pointer CastFilter;
    HFilterType::Pointer HFilter;
    EigValAnalysisType EigValAnalysis;
    VecEigHImagePointerType VecEigHImagePointer;

public:

    EigValHessian(const InputImageType::Pointer& image) {
        VecVoxType EigVal;
        MatVoxType EigMat,tmpMat;
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                EigMat[i][j]=0;
        region=image->GetLargestPossibleRegion();
        spacing=image->GetSpacing();
        origin=image->GetOrigin();

        CastFilter = CastFilterType::New();
        HFilter = HFilterType::New();
        HFilter->SetSigma(3);//3 for sigma;read the definition of the Gaussian kernel;

        //EigValAnalysis = EigValAnalysisType::New();
        EigValAnalysis.SetDimension(3);
        CastFilter->SetInput(image);
        HFilter->SetInput(CastFilter->GetOutput());
        HFilter->Update();
        VecEigHImagePointer=VecEigHImageType::New();

        VecEigHImagePointer->SetRegions(region);
        VecEigHImagePointer->SetOrigin(origin);
        VecEigHImagePointer->SetSpacing(spacing);
        VecEigHImagePointer->Allocate();
        EigVal[0]=EigVal[1]=EigVal[2]=0;
        VecEigHImagePointer->FillBuffer(EigVal);
        HIterType HIter(HFilter->GetOutput(),region);
        OutputIterType OutputIter(VecEigHImagePointer,region);
        itk::SymmetricSecondRankTensor<double,3> Tensor;
        for(HIter.GoToBegin(),OutputIter.GoToBegin();!HIter.IsAtEnd()&&!OutputIter.IsAtEnd();++HIter,++OutputIter){
            Tensor=HIter.Get();
            tmpMat[0][0]=Tensor[0];
            tmpMat[0][1]=Tensor[1];
            tmpMat[1][0]=Tensor[1];
            tmpMat[0][2]=Tensor[2];
            tmpMat[2][0]=Tensor[2];
            tmpMat[1][1]=Tensor[3];
            tmpMat[2][1]=Tensor[4];
            tmpMat[1][2]=Tensor[4];
            tmpMat[2][2]=Tensor[5];
            EigValAnalysis.ComputeEigenValuesAndVectors(tmpMat,EigVal,EigMat);
            OutputIter.Set(EigVal);
        }
        try{
            VecEigHImagePointer->Update();
        }
        catch(itk::ExceptionObject &excp) {
            std::cerr << "Exception thrown while computing the eigenvalues vector!" << std::endl;
            std::cerr << excp << std::endl;
            exit(EXIT_FAILURE);
        }
    }
};


 

[1] Iv´an Macıa, Generalized Computation of Gaussian Derivatives Using ITK, October 5, 2007. 1VICOMTech, Paseo Mikeletegi 57, San Sebastián.