itkSymmetricSecondRankTensor.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 itkSymmetricSecondRankTensor_h
#define itkSymmetricSecondRankTensor_h
 
// Undefine an eventual SymmetricSecondRankTensor macro
#ifdef SymmetricSecondRankTensor
#  undef SymmetricSecondRankTensor
#endif
 
#include "itkFixedArray.h"
#include "itkMatrix.h"
#include "itkSymmetricEigenAnalysis.h"
 
namespace itk
{
 
template <typename TComponent, unsigned int VDimension = 3>
class ITK_TEMPLATE_EXPORT SymmetricSecondRankTensor : public FixedArray<TComponent, VDimension *(VDimension + 1) / 2>
{
public:
  using Self = SymmetricSecondRankTensor;
  using Superclass = FixedArray<TComponent, VDimension *(VDimension + 1) / 2>;

  static constexpr unsigned int Dimension = VDimension;
  static constexpr unsigned int InternalDimension = VDimension * (VDimension + 1) / 2;

  using BaseArray = FixedArray<TComponent, Self::InternalDimension>;

  using EigenValuesArrayType = FixedArray<TComponent, VDimension>;

  using MatrixType = Matrix<TComponent, VDimension, VDimension>;
  using EigenVectorsMatrixType = Matrix<TComponent, VDimension, VDimension>;

  using ComponentType = TComponent;
  using typename Superclass::ValueType;
  using AccumulateValueType = typename NumericTraits<ValueType>::RealType;
  using RealValueType = typename NumericTraits<ValueType>::RealType;
 
  using SymmetricEigenAnalysisType =
    SymmetricEigenAnalysisFixedDimension<Dimension, MatrixType, EigenValuesArrayType, EigenVectorsMatrixType>;

#ifdef ITK_FUTURE_LEGACY_REMOVE
  SymmetricSecondRankTensor() = default;
#else
  SymmetricSecondRankTensor() { this->Fill(0); }
#endif
 
  SymmetricSecondRankTensor(const ComponentType & r) { this->Fill(r); }

  template <typename TCoordinateB>
  SymmetricSecondRankTensor(const SymmetricSecondRankTensor<TCoordinateB, VDimension> & pa)
    : BaseArray(pa)
  {}
 
  using ComponentArrayType = ComponentType[Self::InternalDimension];

  SymmetricSecondRankTensor(const ComponentArrayType r)
    : BaseArray(r)
  {}

  template <typename TCoordinateB>
  Self &
  operator=(const SymmetricSecondRankTensor<TCoordinateB, VDimension> & pa)
  {
    BaseArray::operator=(pa);
    return *this;
  }

  Self &
  operator=(const ComponentType & r);
 
  Self &
  operator=(const ComponentArrayType r);

  Self
  operator+(const Self & r) const;
 
  Self
  operator-(const Self & r) const;
 
  const Self &
  operator+=(const Self & r);
 
  const Self &
  operator-=(const Self & r);

  Self
  operator*(const RealValueType & r) const;
 
  Self
  operator/(const RealValueType & r) const;
 
  const Self &
  operator*=(const RealValueType & r);
 
  const Self &
  operator/=(const RealValueType & r);

  static unsigned int
  GetNumberOfComponents()
  {
    return Self::InternalDimension;
  }

  [[nodiscard]] ComponentType
  GetNthComponent(int c) const
  {
    return this->operator[](c);
  }

  void
  SetNthComponent(int c, const ComponentType & v)
  {
    this->operator[](c) = v;
  }

  ValueType &
  operator()(unsigned int row, unsigned int col);
 
  const ValueType &
  operator()(unsigned int row, unsigned int col) const;

  void
  SetIdentity();

  [[nodiscard]] AccumulateValueType
  GetTrace() const;

  void
  ComputeEigenValues(EigenValuesArrayType & eigenValues) const;

  void
  ComputeEigenAnalysis(EigenValuesArrayType & eigenValues, EigenVectorsMatrixType & eigenVectors) const;

  template <typename TMatrixValueType>
  Self
  Rotate(const Matrix<TMatrixValueType, VDimension, VDimension> & m) const;
  template <typename TMatrixValueType>
  Self
  Rotate(const vnl_matrix_fixed<TMatrixValueType, VDimension, VDimension> & m) const
  {
    return this->Rotate(static_cast<Matrix<TMatrixValueType, VDimension, VDimension>>(m));
  }
  template <typename TMatrixValueType>
  Self
  Rotate(const vnl_matrix<TMatrixValueType> & m) const
  {
    return this->Rotate(static_cast<Matrix<TMatrixValueType>>(m));
  }


  [[nodiscard]] MatrixType
  PreMultiply(const MatrixType & m) const;

  [[nodiscard]] MatrixType
  PostMultiply(const MatrixType & m) const;
};

using OutputStreamType = std::ostream;
using InputStreamType = std::istream;
 
template <typename TComponent, unsigned int VDimension>
OutputStreamType &
operator<<(OutputStreamType & os, const SymmetricSecondRankTensor<TComponent, VDimension> & c);
 
template <typename TComponent, unsigned int VDimension>
InputStreamType &
operator>>(InputStreamType & is, SymmetricSecondRankTensor<TComponent, VDimension> & dt);
 
template <typename T>
inline void
swap(SymmetricSecondRankTensor<T> & a, SymmetricSecondRankTensor<T> & b) noexcept
{
  a.swap(b);
}
} // end namespace itk
 
#include "itkNumericTraitsTensorPixel.h"
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkSymmetricSecondRankTensor.hxx"
#endif
 
#endif