itkImageBase.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.
 *
 *=========================================================================*/
/*=========================================================================
 *
 *  Portions of this file are subject to the VTK Toolkit Version 3 copyright.
 *
 *  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 *
 *  For complete copyright, license and disclaimer of warranty information
 *  please refer to the NOTICE file at the top of the ITK source tree.
 *
 *=========================================================================*/
#ifndef itkImageBase_h
#define itkImageBase_h
 
#include "itkDataObject.h"
 
#include "itkImageRegion.h"
#include "itkMatrix.h"
#include "itkObjectFactory.h"
#include "itkOffset.h"
#include "itkFixedArray.h"
#include "itkImageHelper.h"
#include "itkFloatTypes.h"
 
#include <vxl_version.h>
#include "vnl/vnl_matrix_fixed.hxx" // Get the templates
 
namespace itk
{

inline constexpr double DefaultImageCoordinateTolerance = 1e-6;
inline constexpr double DefaultImageDirectionTolerance = 1e-6;
 

template <unsigned int VImageDimension = 2>
class ITK_TEMPLATE_EXPORT ImageBase : public DataObject
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(ImageBase);

  using Self = ImageBase;
  using Superclass = DataObject;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(ImageBase);

  using ImageDimensionType = unsigned int;

  static constexpr ImageDimensionType ImageDimension = VImageDimension;

  using IndexType = Index<VImageDimension>;
  using IndexValueType = typename IndexType::IndexValueType;

  using OffsetType = Offset<VImageDimension>;
  using OffsetValueType = typename OffsetType::OffsetValueType;

  using SizeType = Size<VImageDimension>;
  using SizeValueType = typename SizeType::SizeValueType;

  using RegionType = ImageRegion<VImageDimension>;

  using SpacingValueType = SpacePrecisionType;
  using SpacingType = Vector<SpacingValueType, VImageDimension>;

  using PointValueType = SpacePrecisionType;
  using PointType = Point<PointValueType, VImageDimension>;

  using DirectionType = Matrix<SpacePrecisionType, VImageDimension, VImageDimension>;

  void
  Initialize() override;

  static constexpr unsigned int
  GetImageDimension()
  {
    return VImageDimension;
  }

  itkSetMacro(Origin, PointType);
  virtual void
  SetOrigin(const double origin[VImageDimension]);
  virtual void
  SetOrigin(const float origin[VImageDimension]);
  virtual void
  SetDirection(const DirectionType & direction);

  itkGetConstReferenceMacro(Direction, DirectionType);

  itkGetConstReferenceMacro(InverseDirection, DirectionType);

  itkGetConstReferenceMacro(Spacing, SpacingType);

  itkGetConstReferenceMacro(Origin, PointType);

  virtual void
  Allocate(bool itkNotUsed(initialize) = false)
  {}

  void
  AllocateInitialized()
  {
    return this->Allocate(true);
  }

  virtual void
  SetLargestPossibleRegion(const RegionType & region);

  virtual const RegionType &
  GetLargestPossibleRegion() const
  {
    return m_LargestPossibleRegion;
  }

  virtual void
  SetBufferedRegion(const RegionType & region);

  virtual const RegionType &
  GetBufferedRegion() const
  {
    return m_BufferedRegion;
  }

  virtual void
  SetRequestedRegion(const RegionType & region);

  void
  SetRequestedRegion(const DataObject * data) override;

  virtual const RegionType &
  GetRequestedRegion() const
  {
    return m_RequestedRegion;
  }

  virtual void
  SetRegions(const RegionType & region)
  {
    this->SetLargestPossibleRegion(region);
    this->SetBufferedRegion(region);
    this->SetRequestedRegion(region);
  }
 
  virtual void
  SetRegions(const SizeType & size)
  {
    RegionType region;
    region.SetSize(size);
 
    this->Self::SetRegions(region);
  }


  const OffsetValueType *
  GetOffsetTable() const
  {
    return m_OffsetTable;
  }

  inline OffsetValueType
  ComputeOffset(const IndexType & ind) const
  {
    OffsetValueType offset = 0;
 
    ImageHelper<VImageDimension, VImageDimension>::ComputeOffset(
      this->GetBufferedRegion().GetIndex(), ind, m_OffsetTable, offset);
    return offset;
    /* NON TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING data version
     * Leaving here for documentation purposes
     * OffsetValueType ComputeOffset(const IndexType & ind) const
     * {
     *   // need to add bounds checking for the region/buffer?
     *   OffsetValueType   offset = 0;
     *   const IndexType & bufferedRegionIndex = this->GetBufferedRegion().GetIndex();
     *   // data is arranged as [][][][slice][row][col]
     *   // with Index[0] = col, Index[1] = row, Index[2] = slice
     *   for ( int i = VImageDimension - 1; i > 0; i-- )
     *     {
     *     offset += ( ind[i] - bufferedRegionIndex[i] ) * m_OffsetTable[i];
     *     }
     *   offset += ( ind[0] - bufferedRegionIndex[0] );
     *   return offset;
     * }
     */
  }

  inline IndexType
  ComputeIndex(OffsetValueType offset) const
  {
    IndexType         index;
    const IndexType & bufferedRegionIndex = this->GetBufferedRegion().GetIndex();
 
    ImageHelper<VImageDimension, VImageDimension>::ComputeIndex(bufferedRegionIndex, offset, m_OffsetTable, index);
    return index;
    /* NON TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING data version
     * Leaving here for documentation purposes
     * IndexType ComputeIndex(OffsetValueType offset) const
     * {
     *   IndexType         index;
     *   const IndexType & bufferedRegionIndex = this->GetBufferedRegion().GetIndex();
     *   for ( int i = VImageDimension - 1; i > 0; i-- )
     *     {
     *     index[i] = static_cast< IndexValueType >( offset / m_OffsetTable[i] );
     *     offset -= ( index[i] * m_OffsetTable[i] );
     *     index[i] += bufferedRegionIndex[i];
     *     }
     *   index[0] = bufferedRegionIndex[0] + static_cast< IndexValueType >( offset );
     *   return index;
     * }
     */
  }

  virtual void
  SetSpacing(const SpacingType & spacing);
  virtual void
  SetSpacing(const double spacing[VImageDimension]);
  virtual void
  SetSpacing(const float spacing[VImageDimension]);
  template <typename TCoordinate>
  [[nodiscard]] IndexType
  TransformPhysicalPointToIndex(const Point<TCoordinate, VImageDimension> & point) const
  {
    IndexType index;
 
    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      TCoordinate sum{};
      for (unsigned int j = 0; j < VImageDimension; ++j)
      {
        sum += this->m_PhysicalPointToIndex[i][j] * (point[j] - this->m_Origin[j]);
      }
      index[i] = Math::RoundHalfIntegerUp<IndexValueType>(sum);
    }
    return index;
  }

  template <typename TCoordinate>
  ITK_NODISCARD("Call the overload which has the point as the only parameter and returns the index")
  bool TransformPhysicalPointToIndex(const Point<TCoordinate, VImageDimension> & point, IndexType & index) const
  {
    index = TransformPhysicalPointToIndex(point);
 
    // Now, check to see if the index is within allowed bounds
    const bool isInside = this->GetLargestPossibleRegion().IsInside(index);
    return isInside;
  }
 

  template <typename TIndexRep, typename TCoordinate>
  [[nodiscard]] ContinuousIndex<TIndexRep, VImageDimension>
  TransformPhysicalPointToContinuousIndex(const Point<TCoordinate, VImageDimension> & point) const
  {
    ContinuousIndex<TIndexRep, VImageDimension> index;
    Vector<SpacePrecisionType, VImageDimension> cvector;
 
    for (unsigned int k = 0; k < VImageDimension; ++k)
    {
      cvector[k] = point[k] - this->m_Origin[k];
    }
    cvector = m_PhysicalPointToIndex * cvector;
    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      index[i] = static_cast<TIndexRep>(cvector[i]);
    }
    return index;
  }

  template <typename TCoordinate, typename TIndexRep>
  ITK_NODISCARD("Call the overload which has the point as the only parameter and returns the index")
  bool TransformPhysicalPointToContinuousIndex(const Point<TCoordinate, VImageDimension> &   point,
                                               ContinuousIndex<TIndexRep, VImageDimension> & index) const
  {
    index = TransformPhysicalPointToContinuousIndex<TIndexRep>(point);
 
    // Now, check to see if the index is within allowed bounds
    const bool isInside = this->GetLargestPossibleRegion().IsInside(index);
    return isInside;
  }

  template <typename TCoordinate, typename TIndexRep>
  void
  TransformContinuousIndexToPhysicalPoint(const ContinuousIndex<TIndexRep, VImageDimension> & index,
                                          Point<TCoordinate, VImageDimension> &               point) const
  {
    for (unsigned int r = 0; r < VImageDimension; ++r)
    {
      TCoordinate sum{};
      for (unsigned int c = 0; c < VImageDimension; ++c)
      {
        sum += this->m_IndexToPhysicalPoint(r, c) * index[c];
      }
      point[r] = sum + this->m_Origin[r];
    }
  }

  template <typename TCoordinate, typename TIndexRep>
  [[nodiscard]] Point<TCoordinate, VImageDimension>
  TransformContinuousIndexToPhysicalPoint(const ContinuousIndex<TIndexRep, VImageDimension> & index) const
  {
    Point<TCoordinate, VImageDimension> point;
    TransformContinuousIndexToPhysicalPoint(index, point);
    return point;
  }

  template <typename TCoordinate>
  void
  TransformIndexToPhysicalPoint(const IndexType & index, Point<TCoordinate, VImageDimension> & point) const
  {
    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      point[i] = this->m_Origin[i];
      for (unsigned int j = 0; j < VImageDimension; ++j)
      {
        point[i] += m_IndexToPhysicalPoint[i][j] * index[j];
      }
    }
  }

  template <typename TCoordinate>
  [[nodiscard]] Point<TCoordinate, VImageDimension>
  TransformIndexToPhysicalPoint(const IndexType & index) const
  {
    Point<TCoordinate, VImageDimension> point;
    TransformIndexToPhysicalPoint(index, point);
    return point;
  }

  template <typename TCoordinate>
  void
  TransformLocalVectorToPhysicalVector(const FixedArray<TCoordinate, VImageDimension> & inputGradient,
                                       FixedArray<TCoordinate, VImageDimension> &       outputGradient) const
  {
    const DirectionType & direction = this->GetDirection();
 
    itkAssertInDebugAndIgnoreInReleaseMacro(inputGradient.GetDataPointer() != outputGradient.GetDataPointer());
 
    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      using CoordSumType = typename NumericTraits<TCoordinate>::AccumulateType;
      CoordSumType sum{};
      for (unsigned int j = 0; j < VImageDimension; ++j)
      {
        sum += direction[i][j] * inputGradient[j];
      }
      outputGradient[i] = static_cast<TCoordinate>(sum);
    }
  }

  template <typename TVector>
  [[nodiscard]] TVector
  TransformLocalVectorToPhysicalVector(const TVector & inputGradient) const
  {
    TVector outputGradient;
    TransformLocalVectorToPhysicalVector(inputGradient, outputGradient);
    return outputGradient;
  }
 

  template <typename TCoordinate>
  void
  TransformPhysicalVectorToLocalVector(const FixedArray<TCoordinate, VImageDimension> & inputGradient,
                                       FixedArray<TCoordinate, VImageDimension> &       outputGradient) const
  {
    const DirectionType & inverseDirection = this->GetInverseDirection();
 
    itkAssertInDebugAndIgnoreInReleaseMacro(inputGradient.GetDataPointer() != outputGradient.GetDataPointer());
 
    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      using CoordSumType = typename NumericTraits<TCoordinate>::AccumulateType;
      CoordSumType sum{};
      for (unsigned int j = 0; j < VImageDimension; ++j)
      {
        sum += inverseDirection[i][j] * inputGradient[j];
      }
      outputGradient[i] = static_cast<TCoordinate>(sum);
    }
  }

  template <typename TVector>
  [[nodiscard]] TVector
  TransformPhysicalVectorToLocalVector(const TVector & inputGradient) const
  {
    TVector outputGradient;
    TransformPhysicalVectorToLocalVector(inputGradient, outputGradient);
    return outputGradient;
  }

  void
  CopyInformation(const DataObject * data) override;

  virtual void
  Graft(const Self * image);

  void
  UpdateOutputInformation() override;

  void
  UpdateOutputData() override;

  void
  SetRequestedRegionToLargestPossibleRegion() override;

  bool
  RequestedRegionIsOutsideOfTheBufferedRegion() override;

  bool
  VerifyRequestedRegion() override;

  bool
  IsCongruentImageGeometry(const ImageBase * otherImage, double coordinateTolerance, double directionTolerance) const;

  bool
  IsSameImageGeometryAs(const ImageBase * otherImage,
                        double            coordinateTolerance = DefaultImageCoordinateTolerance,
                        double            directionTolerance = DefaultImageDirectionTolerance) const;

  virtual unsigned int
  GetNumberOfComponentsPerPixel() const;
  virtual void
  SetNumberOfComponentsPerPixel(unsigned int);
protected:
  ImageBase() = default;
  ~ImageBase() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;

  void
  ComputeOffsetTable();

  virtual void
  ComputeIndexToPhysicalPointMatrices();
 
protected:
  SpacingType   m_Spacing{ MakeFilled<SpacingType>(1.0) };
  PointType     m_Origin{};
  DirectionType m_Direction{ DirectionType::GetIdentity() };
  DirectionType m_InverseDirection{ DirectionType::GetIdentity() };
  DirectionType m_IndexToPhysicalPoint{ DirectionType::GetIdentity() };
  DirectionType m_PhysicalPointToIndex{ DirectionType::GetIdentity() };
  virtual void
  InitializeBufferedRegion();

  OffsetValueType
  FastComputeOffset(const IndexType & ind) const
  {
    OffsetValueType offset = 0;
    ImageHelper<VImageDimension, VImageDimension>::ComputeOffset(
      Self::GetBufferedRegion().GetIndex(), ind, m_OffsetTable, offset);
    return offset;
  }

  IndexType
  FastComputeIndex(OffsetValueType offset) const
  {
    IndexType         index;
    const IndexType & bufferedRegionIndex = Self::GetBufferedRegion().GetIndex();
    ImageHelper<VImageDimension, VImageDimension>::ComputeIndex(bufferedRegionIndex, offset, m_OffsetTable, index);
    return index;
  }
 
  void
  Graft(const DataObject * data) override;
 
private:
  OffsetValueType m_OffsetTable[VImageDimension + 1]{};
 
  RegionType m_LargestPossibleRegion{};
  RegionType m_RequestedRegion{};
  RegionType m_BufferedRegion{};
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkImageBase.hxx"
#endif
 
#endif