itkVectorGradientMagnitudeImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkVectorGradientMagnitudeImageFilter_h
#define itkVectorGradientMagnitudeImageFilter_h

#include "itkNeighborhoodIterator.h"
#include "itkImageToImageFilter.h"
#include "itkImage.h"
#include "itkVector.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/algo/vnl_symmetric_eigensystem.h"
#include "itkMath.h"

namespace itk
{
template< typename TInputImage,
          typename TRealType = float,
          typename TOutputImage = Image< TRealType,
                                         TInputImage::ImageDimension >
          >
class ITK_TEMPLATE_EXPORT VectorGradientMagnitudeImageFilter:
  public ImageToImageFilter< TInputImage, TOutputImage >
{
public:
  ITK_DISALLOW_COPY_AND_ASSIGN(VectorGradientMagnitudeImageFilter);

  using Self = VectorGradientMagnitudeImageFilter;
  using Superclass = ImageToImageFilter< TInputImage, TOutputImage >;
  using Pointer = SmartPointer< Self >;
  using ConstPointer = SmartPointer< const Self >;

  itkNewMacro(Self);

  itkTypeMacro(VectorGradientMagnitudeImageFilter, ImageToImageFilter);

  using OutputPixelType = typename TOutputImage::PixelType;
  using InputPixelType = typename TInputImage::PixelType;

  using InputImageType = TInputImage;
  using OutputImageType = TOutputImage;
  using InputImagePointer = typename InputImageType::Pointer;
  using OutputImagePointer = typename OutputImageType::Pointer;

  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;

  static constexpr unsigned int VectorDimension = InputPixelType::Dimension;

  using RealType = TRealType;
  using RealVectorType =
      Vector< TRealType, InputPixelType::Dimension >;
  using RealVectorImageType =
      Image< RealVectorType, TInputImage::ImageDimension >;

  using ConstNeighborhoodIteratorType = ConstNeighborhoodIterator< RealVectorImageType >;
  using RadiusType = typename ConstNeighborhoodIteratorType::RadiusType;

  using OutputImageRegionType = typename Superclass::OutputImageRegionType;

  void GenerateInputRequestedRegion() override;

  void SetUseImageSpacingOn()
  { this->SetUseImageSpacing(true); }

  void SetUseImageSpacingOff()
  { this->SetUseImageSpacing(false); }

  void SetUseImageSpacing(bool);

  itkGetConstMacro(UseImageSpacing, bool);

  using WeightsType = FixedArray< TRealType, VectorDimension >;

  itkSetMacro(DerivativeWeights, WeightsType);
  itkGetConstReferenceMacro(DerivativeWeights, WeightsType);

  itkSetMacro(ComponentWeights, WeightsType);
  itkGetConstReferenceMacro(ComponentWeights, WeightsType);

  itkSetMacro(UsePrincipleComponents, bool);
  itkGetConstMacro(UsePrincipleComponents, bool);
  void SetUsePrincipleComponentsOn()
  {
    this->SetUsePrincipleComponents(true);
  }

  void SetUsePrincipleComponentsOff()
  {
    this->SetUsePrincipleComponents(false);
  }

  static int CubicSolver(double *, double *);

#ifdef ITK_USE_CONCEPT_CHECKING
  // Begin concept checking
  itkConceptMacro( InputHasNumericTraitsCheck,
                   ( Concept::HasNumericTraits< typename InputPixelType::ValueType > ) );
  itkConceptMacro( RealTypeHasNumericTraitsCheck,
                   ( Concept::HasNumericTraits< RealType > ) );
  // End concept checking
#endif

protected:
  VectorGradientMagnitudeImageFilter();
  ~VectorGradientMagnitudeImageFilter() override = default;

  void BeforeThreadedGenerateData() override;

  void DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;


  void PrintSelf(std::ostream & os, Indent indent) const override;

  using ImageBaseType = typename InputImageType::Superclass;

  itkGetConstObjectMacro(RealValuedInputImage, ImageBaseType);

  TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
  {
    unsigned  i, j;
    TRealType dx, sum, accum;

    accum = NumericTraits< TRealType >::ZeroValue();
    for ( i = 0; i < ImageDimension; ++i )
      {
      sum = NumericTraits< TRealType >::ZeroValue();
      for ( j = 0; j < VectorDimension; ++j )
        {
        dx =  m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
             * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
        sum += dx * dx;
        }
      accum += sum;
      }
    return std::sqrt(accum);
  }

  TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType & it) const
  {
    // WARNING:  ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
    unsigned int i, j;
    double       Lambda[3];
    double       CharEqn[3];
    double       ans;

    vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
    vnl_vector_fixed< TRealType, VectorDimension >
    d_phi_du[TInputImage::ImageDimension];

    // Calculate the directional derivatives for each vector component using
    // central differences.
    for ( i = 0; i < ImageDimension; i++ )
      {
      for ( j = 0; j < VectorDimension; j++ )
        {
        d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
                         * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
        }
      }

    // Calculate the symmetric metric tensor g
    for ( i = 0; i < ImageDimension; i++ )
      {
      for ( j = i; j < ImageDimension; j++ )
        {
        g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
        }
      }

    // Find the coefficients of the characteristic equation det(g - lambda I)=0
    //    CharEqn[3] = 1.0;

    CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );

    CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
                 - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );

    CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] -  g[1][1] * g[2][2]  )
                 + g[1][0] * (  g[2][2] * g[0][1] - g[0][2] * g[2][1] )
                 + g[2][0] * ( g[1][1] * g[0][2] - g[0][1] * g[1][2] );

    // Find the eigenvalues of g
    int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);

    // Define gradient magnitude as the difference between two largest
    // eigenvalues.  Other definitions may be appropriate here as well.
    if ( numberOfDistinctRoots == 3 ) // By far the most common case
      {
      if ( Lambda[0] > Lambda[1] )
        {
        if ( Lambda[1] > Lambda[2] )  // Most common, guaranteed?
          {
          ans = Lambda[0] - Lambda[1];
          }
        else
          {
          if ( Lambda[0] > Lambda[2] )
            {
            ans = Lambda[0] - Lambda[2];
            }
          else
            {
            ans = Lambda[2] - Lambda[0];
            }
          }
        }
      else
        {
        if ( Lambda[0] > Lambda[2] )
          {
          ans = Lambda[1] - Lambda[0];
          }
        else
          {
          if ( Lambda[1] > Lambda[2] )
            {
            ans = Lambda[1] - Lambda[2];
            }
          else
            {
            ans = Lambda[2] - Lambda[1];
            }
          }
        }
      }
    else if ( numberOfDistinctRoots == 2 )
      {
      if ( Lambda[0] > Lambda[1] )
        {
        ans = Lambda[0] - Lambda[1];
        }
      else
        {
        ans = Lambda[1] - Lambda[0];
        }
      }
    else if ( numberOfDistinctRoots == 1 )
      {
      ans = 0.0;
      }
    else
      {
      itkExceptionMacro(<< "Undefined condition. Cubic root solver returned "
                        << numberOfDistinctRoots << " distinct roots.");
      }

    return ans;
  }

  // Function is defined here because the templating confuses gcc 2.96 when
  // defined
  // in .hxx file.  jc 1/29/03
  TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
  {
    unsigned int i, j;

    vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
    vnl_vector_fixed< TRealType, VectorDimension >
    d_phi_du[TInputImage::ImageDimension];

    // Calculate the directional derivatives for each vector component using
    // central differences.
    for ( i = 0; i < ImageDimension; i++ )
      {
      for ( j = 0; j < VectorDimension; j++ )
        {
        d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
                         * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
        }
      }

    // Calculate the symmetric metric tensor g
    for ( i = 0; i < ImageDimension; i++ )
      {
      for ( j = i; j < ImageDimension; j++ )
        {
        g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
        }
      }

    // Find the eigenvalues of g
    vnl_symmetric_eigensystem< TRealType > E(g);

    // Return the difference in length between the first two principle axes.
    // Note that other edge strength metrics may be appropriate here instead..
    return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
  }

  WeightsType m_DerivativeWeights;

  WeightsType m_ComponentWeights;
  WeightsType m_SqrtComponentWeights;

private:
  bool m_UseImageSpacing;
  bool m_UsePrincipleComponents;

  ThreadIdType m_RequestedNumberOfThreads;

  typename ImageBaseType::ConstPointer m_RealValuedInputImage;

};
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkVectorGradientMagnitudeImageFilter.hxx"
#endif

#endif