Dear authors,

Thanks for implementing these filters in ITK. However, I agree with Fethallah Benmansour that there was a loss of precision in the structure tensor code and the outer product components of the structure tensor needed to be smoothed using an additional parameter rho.

I found another error in the function "UpdateDiffusionTensorImage()" of all three filters: itkAnisotropicCoherenceEnhancingDiffusionImageFilter.txx, itkAnisotropicEdgeEnhancementDiffusionImageFilter.txx, itkAnisotropicHybridDiffusionImageFilter.txx.

The eigenvalues of the structure tensor are neatly ordered in decreasing order, however the corresponding eigenvectors of the structure tensor are not! This will cause all diffusion filters (EED, CED, HDCS) to work incorrectly...

The diffusion tensor is based on the orientation of the structure tensor. The first eigenvector of the structure tensor is assumed to point into the direction with the highest contrast (largest eigenvalue), however, in the current version of the code, the first eigenvector points into the direction of the lowest contrast (still sorted in ascending order). Because the eigenvalues are sorted, the largest eigenvalue is 'coupled' to the eigenvector with originally the lowest eigenvalue, therefore filtering across edges will take place instead of along edges...

The eigenvectors can be ordered after ordering the eigenvalues in the function UpdateDiffusionTensorImage(), by for example adding the following code:

    unsigned int vectorLength = 3; // Eigenvector length
    EigenVectorMatrixType eigenVectorMatrix;
    eigenVectorMatrix = eigenVectorImageIterator.Get();

    itk::VariableLengthVector<double> firstEigenVector( vectorLength );
    itk::VariableLengthVector<double> secondEigenVector( vectorLength );
    itk::VariableLengthVector<double> thirdEigenVector( vectorLength );
    for ( unsigned int i=0; i < vectorLength; i++ ) {
      // Get eigenvectors belonging to eigenvalue order
      firstEigenVector[i] = eigenVectorMatrix[largestEigenValueIndex][i];
      secondEigenVector[i] = eigenVectorMatrix[middleEigenValueIndex][i];
      thirdEigenVector[i] = eigenVectorMatrix[smallestEigenValueIndex][i];
     
      // Set eigenVectorMatrix in correct order
      eigenVectorMatrix[0][i] = firstEigenVector[i];
      eigenVectorMatrix[1][i] = secondEigenVector[i];
      eigenVectorMatrix[2][i] = thirdEigenVector[i];
    }

Dear authors,

I found a couple of errors in the provided implementation of the structure tensor.

-- First, in line 260 of itkStructureTensorRecursiveGaussianImageFilter.txx.

std::vector<PixelType> tmp( numberTensorElements );

PixelType might be uchar or ushort, and therefore yields a loss of precision. I replaced PixelType by InternalRealType.

-- Second, the smoothing of the outer product' components was not performed.

This additional filtering have been added.

Attached the modified *.h and *.txx files.

Regards,