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
 *
 *         https://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
{
 
 * \par Constraints
 * The filter requires an image with at least two dimensions and a vector
 * length of at least 2.  The theory supports extension to scalar images, but
 * the implementation of the itk vector classes do not
 *
 * The template parameter TRealType must be floating point (float or double) or
 * a user-defined "real" numerical type with arithmetic operations defined
 * sufficient to compute derivatives.
 *
 * \par Performance
 * This filter will automatically multithread if run with
 * SetUsePrincipleComponents=Off or on 3D data in UsePrincipleComponents=On
 * mode.   Unfortunately the ND eigen solver used is not thread safe (a special
 * 3D solver is), so it cannot multithread for data other than 3D in
 * UsePrincipleComponents=On mode.
 *
 * For algorithmic details see \cite sapiro1996.
 *
 * \ingroup GradientFilters
 *
 * \sa Image
 * \sa Neighborhood
 * \sa NeighborhoodOperator
 * \sa NeighborhoodIterator
 * \ingroup ITKImageGradient
 */
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_MOVE(VectorGradientMagnitudeImageFilter);

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

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(VectorGradientMagnitudeImageFilter);

  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 InputImageDimension = TInputImage::ImageDimension;
  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 typename Superclass::OutputImageRegionType;

  void
  GenerateInputRequestedRegion() override;

  void
  SetUseImageSpacing(bool);
  itkGetConstMacro(UseImageSpacing, bool);
  itkBooleanMacro(UseImageSpacing);
 
#if !defined(ITK_FUTURE_LEGACY_REMOVE)
  void
  SetUseImageSpacingOn()
  {
    this->SetUseImageSpacing(true);
  }

  void
  SetUseImageSpacingOff()
  {
    this->SetUseImageSpacing(false);
  }
#endif
 
  using ComponentWeightsType = FixedArray<TRealType, VectorDimension>;
  using DerivativeWeightsType = FixedArray<TRealType, ImageDimension>;
#if !defined(ITK_LEGACY_REMOVE)
  using WeightsType [[deprecated("Use DerivativeWeightsType or ComponentWeightsType instead.")]] = ComponentWeightsType;
#endif

  itkSetMacro(DerivativeWeights, DerivativeWeightsType);
  itkGetConstReferenceMacro(DerivativeWeights, DerivativeWeightsType);

  itkSetMacro(ComponentWeights, ComponentWeightsType);
  itkGetConstReferenceMacro(ComponentWeights, ComponentWeightsType);

  itkSetMacro(UsePrincipleComponents, bool);
  itkGetConstMacro(UsePrincipleComponents, bool);
  itkBooleanMacro(UsePrincipleComponents);
 
#if !defined(ITK_FUTURE_LEGACY_REMOVE)
  void
  SetUsePrincipleComponentsOn()
  {
    this->SetUsePrincipleComponents(true);
  }

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

  static int
  CubicSolver(const double *, double *);
 
  itkConceptMacro(SameDimensionCheck, (Concept::SameDimension<InputImageDimension, ImageDimension>));
  itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits<typename InputPixelType::ValueType>));
  itkConceptMacro(RealTypeHasNumericTraitsCheck, (Concept::HasNumericTraits<RealType>));
 
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 int i;
    unsigned int j;
    TRealType    dx;
    TRealType    sum;
    auto         accum = TRealType{};
    for (i = 0; i < ImageDimension; ++i)
    {
      sum = TRealType{};
      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;
    unsigned int 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
    const 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;
    unsigned int 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
    const 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));
  }

  DerivativeWeightsType m_DerivativeWeights = DerivativeWeightsType::Filled(1);

  ComponentWeightsType m_ComponentWeights = ComponentWeightsType::Filled(1);
  ComponentWeightsType m_SqrtComponentWeights = ComponentWeightsType::Filled(1);
 
private:
  bool m_UseImageSpacing{ true };
  bool m_UsePrincipleComponents{};
 
  ThreadIdType m_RequestedNumberOfWorkUnits{};
 
  typename RealVectorImageType::ConstPointer m_RealValuedInputImage{};
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkVectorGradientMagnitudeImageFilter.hxx"
#endif
 
#endif