itkVectorGradientMagnitudeImageFilter.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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, RealVectorImageType);
 
  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 RealVectorImageType::ConstPointer m_RealValuedInputImage;
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkVectorGradientMagnitudeImageFilter.hxx"
#endif
 
#endif