/*=========================================================================
 *
 *  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
{

/** The default tolerance when comparing the geometry of two images.
 *
 * The value was chosen based on precisions of file formats such as DICOM,
 * to enable interoperability with images saved to file formats with higher precision.
 */
inline constexpr double DefaultImageCoordinateTolerance = 1e-6;
inline constexpr double DefaultImageDirectionTolerance = 1e-6;


/** \class ImageBase
 * \brief Base class for templated image classes.
 *
 * ImageBase is the base class for the templated Image
 * classes. ImageBase is templated over the dimension of the image. It
 * provides the API and ivars that depend solely on the dimension of
 * the image.  ImageBase does not store any of the image (pixel) data.
 * Storage for the pixels and the pixel access methods are defined in
 * subclasses of ImageBase, namely Image and ImageAdaptor.
 *
 * ImageBase manages the geometry of an image. The geometry of an
 * image is defined by its position, orientation, spacing, and extent.
 *
 * The position and orientation of an image is defined by its "Origin"
 * and its "Directions".  The "Origin" is the physical position of the
 * pixel whose "Index" is all zeros. The "Direction" of an image is a
 * matrix whose columns indicate the direction in physical space that
 * each dimension of the image traverses. The first column defines the
 * direction that the fastest moving index in the image traverses in
 * physical space while the last column defines the direction that the
 * slowest moving index in the image traverses in physical space.
 *
 * The extent of an image is defined by the pixel spacing and a set of
 * regions. The "Spacing" is the size of a pixel in physical space
 * along each dimension.  Regions describe a portion of an image grid
 * via a starting index for the image array and a size (or number of
 * pixels) in each dimension. The ivar LargestPossibleRegion defines
 * the size and starting index of the image dataset. The entire image
 * dataset, however, may not be resident in memory. The region of the
 * image that is resident in memory is defined by the
 * "BufferedRegion". The Buffer is a contiguous block of memory.  The
 * third set of meta-data defines a region of interest, called the
 * "RequestedRegion". The RequestedRegion is used by the pipeline
 * execution model to define what a filter is requested to produce.
 *
 * [RegionIndex, RegionSize] C [BufferIndex, BufferSize]
 *                           C [ImageIndex, ImageSize]
 *
 * ImageBase provides all the methods for converting between the
 * physical space and index coordinate
 * frames. TransformIndexToPhysicalPoint() converts an Index in the
 * pixel array into its coordinates in physical space.
 * TransformPhysicalPointToIndex() converts a position in physical
 * space into an Index into the pixel array (using
 * rounding). Subpixel locations are supported by methods that
 * convert to and from ContinuousIndex types.
 *
 * ImageBase also provides helper routines for the ImageIterators
 * which convert an Index to an offset in memory from the first pixel
 * address as well as convert an offset in memory from the first pixel
 * address to an Index.
 *
 * \ingroup ImageObjects
 * \ingroup ITKSystemObjects
 *
 * \ingroup ITKCommon
 */
template <unsigned int VImageDimension = 2>
class ITK_TEMPLATE_EXPORT ImageBase : public DataObject
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(ImageBase);

  /** Standard type alias. */
  using Self = ImageBase;
  using Superclass = DataObject;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** \see LightObject::GetNameOfClass() */
  itkOverrideGetNameOfClassMacro(ImageBase);

  /** Type of image dimension */
  using ImageDimensionType = unsigned int;

  /** Dimension of the image.  This constant is used by functions that are
   * templated over image type (as opposed to being templated over pixel
   * type and dimension) when they need compile time access to the dimension
   * of the image. */
  static constexpr ImageDimensionType ImageDimension = VImageDimension;

  /** Index type alias support. An index is used to access pixel values. */
  using IndexType = Index<VImageDimension>;
  using IndexValueType = typename IndexType::IndexValueType;

  /** Offset type alias support. An offset represent relative position
   * between indices. */
  using OffsetType = Offset<VImageDimension>;
  using OffsetValueType = typename OffsetType::OffsetValueType;

  /** Size type alias support. A size is used to define region bounds. */
  using SizeType = Size<VImageDimension>;
  using SizeValueType = typename SizeType::SizeValueType;

  /** Region type alias support. A region is used to specify a subset of an image. */
  using RegionType = ImageRegion<VImageDimension>;

  /** Spacing type alias support.  Spacing holds the size of a pixel.
   * The spacing is the geometric distance between image samples along
   * each dimension. ITK only supports positive spacing values:
   * negative values may cause undesirable results.  */
  using SpacingValueType = SpacePrecisionType;
  using SpacingType = Vector<SpacingValueType, VImageDimension>;

  /** Origin type alias support.  The origin is the geometric coordinates
   * of the index (0,0). */
  using PointValueType = SpacePrecisionType;
  using PointType = Point<PointValueType, VImageDimension>;

  /** Direction type alias support.  The Direction is a matrix of
   * direction cosines that specify the direction in physical space
   * between samples along each dimension. */
  using DirectionType = Matrix<SpacePrecisionType, VImageDimension, VImageDimension>;

  /** Restore object to initialized state. */
  void
  Initialize() override;

  /** Image dimension. The dimension of an image is fixed at construction. */
  static unsigned int
  GetImageDimension()
  {
    return VImageDimension;
  }

  /** Set the origin of the image. The origin is the geometric
   * coordinates of the image origin (pixel [0,0]).  It is stored internally
   * as SpacePrecisionType but may be set from float or double.
   * \sa GetOrigin() */
  itkSetMacro(Origin, PointType);
  virtual void
  SetOrigin(const double origin[VImageDimension]);
  virtual void
  SetOrigin(const float origin[VImageDimension]);

  /** Set the direction cosines of the image. The direction cosines
   * are vectors that point from one pixel to the next.
   *
   * Each column of the matrix indicates the direction cosines of the unit vector
   * that is parallel to the lines of the image grid corresponding to that
   * dimension. For example, an image with Direction matrix
   *
   *    0.866   0.500
   *   -0.500   0.866
   *
   * has an image grid were the fastest changing index (dimension[0]) walks
   * over a line that in physical space is oriented parallel to the vector
   * (0.866, -0.5). The second fastest changing index (dimension[1]) walks along
   * a line that in Physical space is oriented parallel to the vector
   * (0.5, 0.866)
   *
   * The columns of the Direction matrix are expected to form an
   * orthogonal right handed coordinate system.  But this is not
   * checked nor enforced in itk::ImageBase.
   *
   * For details, please see:
   *
   * https://www.itk.org/Wiki/Proposals:Orientation#Some_notes_on_the_DICOM_convention_and_current_ITK_usage
   *
   * \sa GetDirection() */
  virtual void
  SetDirection(const DirectionType & direction);

  /** Get the direction cosines of the image. The direction cosines
   * are vectors that point from one pixel to the next.
   * For ImageBase and Image, the default direction is identity. */
  itkGetConstReferenceMacro(Direction, DirectionType);

  /** Get the inverse direction cosines of the image.
   * These are calculated automatically in SetDirection, thus there
   * is no Set accessor. */
  itkGetConstReferenceMacro(InverseDirection, DirectionType);

  /** Get the spacing (size of a pixel) `of the image. The
   * spacing is the geometric distance between image samples along
   * each dimension. The value returned is a Vector<double, VImageDimension>.
   * For ImageBase and Image, the default data spacing is unity. */
  itkGetConstReferenceMacro(Spacing, SpacingType);

  /** Get the origin of the image. The origin is the geometric
   * coordinates of the index (0,0).  The value returned is a
   * Point<double, VImageDimension>. For ImageBase and Image, the
   * default origin is 0. */
  itkGetConstReferenceMacro(Origin, PointType);

  /** Allocate the image memory. The size of the image must
   * already be set, e.g. by calling SetRegions() or SetBufferedRegion().
   *
   * This method should be pure virtual, if backwards compatibility
   * was not required.
   *
   */
  virtual void
  Allocate(bool initialize = false);

  /** Allocates the pixel buffer of the image, zero-initializing its pixels. `AllocateInitialized()` is equivalent to
   * `Allocate(true)`. It is just intended to make the code more readable. */
  void
  AllocateInitialized()
  {
    return this->Allocate(true);
  }

  /** Set the region object that defines the size and starting index
   * for the largest possible region this image could represent.  This
   * is used in determining how much memory would be needed to load an
   * entire dataset.  It is also used to determine boundary true
   * conditions.
   * \sa ImageRegion, SetBufferedRegion(), SetRequestedRegion() */
  virtual void
  SetLargestPossibleRegion(const RegionType & region);

  /** Get the region object that defines the size and starting index
   * for the largest possible region this image could represent.  This
   * is used in determining how much memory would be needed to load an
   * entire dataset.  It is also used to determine boundary true
   * conditions.
   * \sa ImageRegion, GetBufferedRegion(), GetRequestedRegion() */
  virtual const RegionType &
  GetLargestPossibleRegion() const
  {
    return m_LargestPossibleRegion;
  }

  /** Set the region object that defines the size and starting index
   * of the region of the image currently loaded in memory.
   * \sa ImageRegion, SetLargestPossibleRegion(), SetRequestedRegion() */
  virtual void
  SetBufferedRegion(const RegionType & region);

  /** Get the region object that defines the size and starting index
   * of the region of the image currently loaded in memory.
   * \sa ImageRegion, SetLargestPossibleRegion(), SetRequestedRegion() */
  virtual const RegionType &
  GetBufferedRegion() const
  {
    return m_BufferedRegion;
  }

  /** Set the region object that defines the size and starting index
   * for the region of the image requested (i.e., the region of the
   * image to be operated on by a filter). Setting the RequestedRegion
   * does not cause the object to be modified. This method is called
   * internally by the pipeline and therefore bypasses the modified
   * time calculation.
   * \sa ImageRegion, SetLargestPossibleRegion(), SetBufferedRegion() */
  virtual void
  SetRequestedRegion(const RegionType & region);

  /** Set the requested region from this data object to match the requested
   * region of the data object passed in as a parameter.  This method
   * implements the API from DataObject. The data object parameter must be
   * castable to an ImageBase. Setting the RequestedRegion does not cause
   * the object to be modified. This method is called internally by
   * the pipeline and therefore bypasses the modified time
   * calculation. */
  void
  SetRequestedRegion(const DataObject * data) override;

  /** Get the region object that defines the size and starting index
   * for the region of the image requested (i.e., the region of the
   * image to be operated on by a filter).
   * \sa ImageRegion, SetLargestPossibleRegion(), SetBufferedRegion() */
  virtual const RegionType &
  GetRequestedRegion() const
  {
    return m_RequestedRegion;
  }

  /** Convenience methods to set the LargestPossibleRegion,
   *  BufferedRegion and RequestedRegion. Allocate must still be called.
   */
  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);
  }

  /** Get the offset table.  The offset table gives increments for
   * moving from one pixel to next in the current row, column, slice,
   * etc..  This table if of size [VImageDimension+1], because its
   * values are computed progressively as: {1, N1, N1*N2,
   * N1*N2*N3,...,(N1*...*Nn)} Where the values {N1,...,Nn} are the
   * elements of the BufferedRegion::Size array.  The last element of
   * the OffsetTable is equivalent to the BufferSize.  Having a
   * [VImageDimension+1] size array, simplifies the implementation of
   * some data accessing algorithms. The entries in the offset table
   * are only valid after the BufferedRegion is set. */
  const OffsetValueType *
  GetOffsetTable() const
  {
    return m_OffsetTable;
  }

  /** Compute an offset from the beginning of the buffer for a pixel
   * at the specified index. The index is not checked as to whether it
   * is inside the current buffer, so the computed offset could
   * conceivably be outside the buffer. If bounds checking is needed,
   * one can call ImageRegion::IsInside(ind) on the BufferedRegion
   * prior to calling ComputeOffset. */
  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;
     * }
     */
  }

  /** Compute the index of the pixel at a specified offset from the
   * beginning of the buffered region. Bounds checking is not
   * performed. Thus, the computed index could be outside the
   * BufferedRegion. To ensure a valid index, the parameter "offset"
   * should be between 0 and the number of pixels in the
   * BufferedRegion (the latter can be found using
   * ImageRegion::GetNumberOfPixels()). */
  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;
     * }
     */
  }

  /** Set the spacing (size of a pixel) of the image. The spacing is
   * the geometric distance between image samples along each
   * dimension. It is stored internally as double, but may be set from
   * float. These methods also pre-compute the Index to Physical point
   * transforms of the image.
   * \sa GetSpacing() */
  virtual void
  SetSpacing(const SpacingType & spacing);
  virtual void
  SetSpacing(const double spacing[VImageDimension]);
  virtual void
  SetSpacing(const float spacing[VImageDimension]);

  /** Returns the index (discrete) of a voxel from a physical point.
   * Floating point index results are rounded to integers
   * \note This specific overload does not figure out whether or not
   * the returned index is inside the image. Of course, the user can
   * still test this afterwards by calling ImageRegion::IsInside(index):
     \code
      auto index = image->TransformPhysicalPointToIndex(point);
      if (image->GetLargestPossibleRegion().IsInside(index)) // Et cetera...
     \endcode
   * Which is equivalent to the following code, which calls the other overload:
     \code
      IndexType index;
      if (image->TransformPhysicalPointToIndex(point, index)) // Et cetera...
     \endcode
   * \sa Transform */
  template <typename TCoordRep>
  [[nodiscard]] IndexType
  TransformPhysicalPointToIndex(const Point<TCoordRep, VImageDimension> & point) const
  {
    IndexType index;

    for (unsigned int i = 0; i < VImageDimension; ++i)
    {
      TCoordRep 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;
  }

  /** Get the index (discrete) of a voxel from a physical point.
   * Floating point index results are rounded to integers
   * Returns true if the resulting index is within the image, false otherwise
   *
   * \note For performance reasons, if you do not need to use the `bool` return value, please call the corresponding
   * overload instead, which has only one parameter (the point), and returns the index.
   *
   * \sa Transform */
  template <typename TCoordRep>
  ITK_NODISCARD("Call the overload which has the point as the only parameter and returns the index")
  bool TransformPhysicalPointToIndex(const Point<TCoordRep, 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;
  }


  /** \brief Returns the continuous index from a physical point
   * \note This specific overload does not figure out whether or not
   *  the returned index is inside the image. Of course, the user can
   * still test this afterwards by calling ImageRegion::IsInside(index):
     \code
     auto index = image->TransformPhysicalPointToContinuousIndex<double>(point);
     if (image->GetLargestPossibleRegion().IsInside(index)) // Et cetera...
     \endcode
   * Which is equivalent to the following code, which calls the other overload:
     \code
     itk::ContinuousIndex<double, ImageDimension> index;
     if (image->TransformPhysicalPointToContinuousIndex(point, index)) // Et cetera...
     \endcode
   * \sa Transform */
  template <typename TIndexRep, typename TCoordRep>
  [[nodiscard]] ContinuousIndex<TIndexRep, VImageDimension>
  TransformPhysicalPointToContinuousIndex(const Point<TCoordRep, 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;
  }

  /** \brief Get the continuous index from a physical point
   *
   * Returns true if the resulting index is within the image, false otherwise.
   *
   * \note For performance reasons, if you do not need to use the `bool` return value, please call the corresponding
   * overload instead, which has only one parameter (the point), and returns the continuous index.
   *
   * \sa Transform */
  template <typename TCoordRep, typename TIndexRep>
  ITK_NODISCARD("Call the overload which has the point as the only parameter and returns the index")
  bool TransformPhysicalPointToContinuousIndex(const Point<TCoordRep, 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;
  }

  /** Get a physical point (in the space which
   * the origin and spacing information comes from)
   * from a continuous index (in the index space)
   * \sa Transform */
  template <typename TCoordRep, typename TIndexRep>
  void
  TransformContinuousIndexToPhysicalPoint(const ContinuousIndex<TIndexRep, VImageDimension> & index,
                                          Point<TCoordRep, VImageDimension> &                 point) const
  {
    for (unsigned int r = 0; r < VImageDimension; ++r)
    {
      TCoordRep sum{};
      for (unsigned int c = 0; c < VImageDimension; ++c)
      {
        sum += this->m_IndexToPhysicalPoint(r, c) * index[c];
      }
      point[r] = sum + this->m_Origin[r];
    }
  }

  /** Returns a physical point (in the space which
   * the origin and spacing information comes from)
   * from a continuous index (in the index space)
   * \sa Transform */
  template <typename TCoordRep, typename TIndexRep>
  [[nodiscard]] Point<TCoordRep, VImageDimension>
  TransformContinuousIndexToPhysicalPoint(const ContinuousIndex<TIndexRep, VImageDimension> & index) const
  {
    Point<TCoordRep, VImageDimension> point;
    TransformContinuousIndexToPhysicalPoint(index, point);
    return point;
  }

  /** Get a physical point (in the space which
   * the origin and spacing information comes from)
   * from a discrete index (in the index space)
   *
   * \sa Transform */
  template <typename TCoordRep>
  void
  TransformIndexToPhysicalPoint(const IndexType & index, Point<TCoordRep, 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];
      }
    }
  }

  /** Returns a physical point (in the space which
   * the origin and spacing information comes from)
   * from a discrete index (in the index space)
   *
   * \sa Transform */
  template <typename TCoordRep>
  [[nodiscard]] Point<TCoordRep, VImageDimension>
  TransformIndexToPhysicalPoint(const IndexType & index) const
  {
    Point<TCoordRep, VImageDimension> point;
    TransformIndexToPhysicalPoint(index, point);
    return point;
  }

  /** Take a vector or covariant vector that has been computed in the
   * coordinate system parallel to the image grid and rotate it by the
   * direction cosines in order to get it in terms of the coordinate system of
   * the image acquisition device.  This implementation in the Image
   * multiply the array (vector or covariant vector) by the matrix of Direction
   * Cosines. The arguments of the method are of type FixedArray to make
   * possible to use this method with both Vector and CovariantVector.
   * The Method is implemented differently in the itk::Image.
   *
   * The inputGradient and the outputGradient must not refer to the
   * same data.
   *
   * \sa Image
   */
  template <typename TCoordRep>
  void
  TransformLocalVectorToPhysicalVector(const FixedArray<TCoordRep, VImageDimension> & inputGradient,
                                       FixedArray<TCoordRep, 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<TCoordRep>::AccumulateType;
      CoordSumType sum{};
      for (unsigned int j = 0; j < VImageDimension; ++j)
      {
        sum += direction[i][j] * inputGradient[j];
      }
      outputGradient[i] = static_cast<TCoordRep>(sum);
    }
  }

  /** Take a vector or covariant vector that has been computed in the
   * coordinate system parallel to the image grid and rotate it by the
   * direction cosines in order to get it in terms of the coordinate system of
   * the image acquisition device. Returns the resulting gradient.
   * \sa Image
   */
  template <typename TVector>
  [[nodiscard]] TVector
  TransformLocalVectorToPhysicalVector(const TVector & inputGradient) const
  {
    TVector outputGradient;
    TransformLocalVectorToPhysicalVector(inputGradient, outputGradient);
    return outputGradient;
  }


  /** Take a vector or covariant vector that has been computed in terms of the
   * coordinate system of the image acquisition device, and rotate it by the
   * inverse direction cosines in order to get it in the coordinate system
   * parallel to the image grid. This implementation in the Image
   * multiply the array (vector or covariant vector) by the inverse of Direction
   * Cosines. The arguments of the method are of type FixedArray to make
   * possible to use this method with both Vector and CovariantVector.
   *
   * The inputGradient and the outputGradient must not refer to the
   * same data.
   *
   */
  template <typename TCoordRep>
  void
  TransformPhysicalVectorToLocalVector(const FixedArray<TCoordRep, VImageDimension> & inputGradient,
                                       FixedArray<TCoordRep, 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<TCoordRep>::AccumulateType;
      CoordSumType sum{};
      for (unsigned int j = 0; j < VImageDimension; ++j)
      {
        sum += inverseDirection[i][j] * inputGradient[j];
      }
      outputGradient[i] = static_cast<TCoordRep>(sum);
    }
  }

  /** Take a vector or covariant vector that has been computed in terms of the
   * coordinate system of the image acquisition device, and rotate it by the
   * inverse direction cosines in order to get it in the coordinate system
   * parallel to the image grid. Returns the result.
   *
   */
  template <typename TVector>
  [[nodiscard]] TVector
  TransformPhysicalVectorToLocalVector(const TVector & inputGradient) const
  {
    TVector outputGradient;
    TransformPhysicalVectorToLocalVector(inputGradient, outputGradient);
    return outputGradient;
  }

  /** Copy information from the specified data set.  This method is
   * part of the pipeline execution model. By default, a ProcessObject
   * will copy meta-data from the first input to all of its
   * outputs. See ProcessObject::GenerateOutputInformation().  Each
   * subclass of DataObject is responsible for being able to copy
   * whatever meta-data it needs from from another DataObject.
   * ImageBase has more meta-data than its DataObject.  Thus, it must
   * provide its own version of CopyInformation() in order to copy the
   * LargestPossibleRegion from the input parameter. */
  void
  CopyInformation(const DataObject * data) override;

  /** Graft the data and information from one image to another. This
   * is a convenience method to setup a second image with all the meta
   * information of another image and use the same pixel
   * container. Note that this method is different than just using two
   * SmartPointers to the same image since separate DataObjects are
   * still maintained. This method is similar to
   * ImageSource::GraftOutput(). The implementation in ImageBase
   * simply calls CopyInformation() and copies the region ivars.
   * Subclasses of ImageBase are responsible for copying the pixel
   * container. */
  virtual void
  Graft(const Self * image);

  /** Update the information for this DataObject so that it can be used
   * as an output of a ProcessObject.  This method is used the pipeline
   * mechanism to propagate information and initialize the meta data
   * associated with a DataObject. This method calls its source's
   * ProcessObject::UpdateOutputInformation() which determines modified
   * times, LargestPossibleRegions, and any extra meta data like spacing,
   * origin, etc. */
  void
  UpdateOutputInformation() override;

  /** UpdateOutputData() is part of the pipeline infrastructure to
   * communicate between ProcessObjects and DataObjects. The method of
   * the superclass is overridden to check if the requested image
   * region has zero pixels. This is needed so that filters can set an
   * input's requested region to zero, to indicate that it does not
   * need to be updated or executed.
   */
  void
  UpdateOutputData() override;

  /** Set the RequestedRegion to the LargestPossibleRegion.  This
   * forces a filter to produce all of the output in one execution
   * (i.e. not streaming) on the next call to Update(). */
  void
  SetRequestedRegionToLargestPossibleRegion() override;

  /** Determine whether the RequestedRegion is outside of the
   * BufferedRegion. This method returns true if the RequestedRegion
   * is outside the BufferedRegion (true if at least one pixel is
   * outside). This is used by the pipeline mechanism to determine
   * whether a filter needs to re-execute in order to satisfy the
   * current request.  If the current RequestedRegion is already
   * inside the BufferedRegion from the previous execution (and the
   * current filter is up to date), then a given filter does not need
   * to re-execute */
  bool
  RequestedRegionIsOutsideOfTheBufferedRegion() override;

  /** Verify that the RequestedRegion is within the
   * LargestPossibleRegion.  If the RequestedRegion is not within the
   * LargestPossibleRegion, then the filter cannot possible satisfy
   * the request. This method returns true if the request can be
   * satisfied and returns fails if the request cannot. This method is
   * used by PropagateRequestedRegion().  PropagateRequestedRegion()
   * throws a InvalidRequestedRegionError exception is the requested
   * region is not within the LargestPossibleRegion. */
  bool
  VerifyRequestedRegion() override;

  /** Checks whether the images' pixels with the same index occupy the same physical space.
   * Compares the origin, spacing, and direction for equality within provided tolerances.
   * There is no check for valid regions in between the images. */
  bool
  IsCongruentImageGeometry(const ImageBase * otherImage, double coordinateTolerance, double directionTolerance) const;

  /** Check whether this image and the other image have the same grid in physical space.
   * Compares largest possible regions for equality, and the origin, spacing,
   * and direction cosines for equality within provided tolerances.
   * See also: ImageToImageFilter, namely:
   * https://github.com/InsightSoftwareConsortium/ITK/blob/v5.3.0/Modules/Core/Common/include/itkImageToImageFilter.h#L78-L92
   * https://github.com/InsightSoftwareConsortium/ITK/blob/v5.3.0/Modules/Core/Common/src/itkImageToImageFilterCommon.cxx#L26-L27
   */
  bool
  IsSameImageGeometryAs(const ImageBase * otherImage,
                        double            coordinateTolerance = DefaultImageCoordinateTolerance,
                        double            directionTolerance = DefaultImageDirectionTolerance) const;

  /** INTERNAL This method is used internally by filters to copy meta-data from
   * the output to the input. Users should not have a need to use this method.
   *
   * Filters that override the ProcessObject's GenerateOutputInformation()
   * should generally have the following line if they want to propagate meta-
   * data for both Image and VectorImage
     \code
     outputImage->SetNumberOfComponentsPerPixel(
        inputImage->GetNumberOfComponentsPerPixel() )
     \endcode
   *
   * \sa ImageBase, VectorImage
   *
   * Returns/Sets the number of components in the image. Note that in the
   * ImageBase implementation, this always returns 1. Image returns the
   * # returned from NumericTraits for the pixel type, and VectorImage
   * returns the vector length set by the user.
   */
  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;

  /** Calculate the offsets needed to move from one pixel to the next
   * along a row, column, slice, volume, etc. These offsets are based
   * on the size of the BufferedRegion. This should be called after
   * the BufferedRegion is set. */
  void
  ComputeOffsetTable();

  /** Compute helper matrices used to transform Index coordinates to
   * PhysicalPoint coordinates and back. This method is virtual and will be
   * overloaded in derived classes in order to provide backward compatibility
   * behavior in classes that did not used to take image orientation into
   * account.  */
  virtual void
  ComputeIndexToPhysicalPointMatrices();

protected:
  /** Origin, spacing, and direction in physical coordinates. This variables are
   * protected for efficiency.  They are referenced frequently by
   * inner loop calculations. */
  SpacingType   m_Spacing{ MakeFilled<SpacingType>(1.0) };
  PointType     m_Origin{};
  DirectionType m_Direction{ DirectionType::GetIdentity() };
  DirectionType m_InverseDirection{ DirectionType::GetIdentity() };

  /** Matrices intended to help with the conversion of Index coordinates
   *  to PhysicalPoint coordinates */
  DirectionType m_IndexToPhysicalPoint{ DirectionType::GetIdentity() };
  DirectionType m_PhysicalPointToIndex{ DirectionType::GetIdentity() };

  /** Restores the buffered region to its default state
   *  This method does not call Modify because Initialization is
   *  called by ReleaseData and can not modify the MTime
   * \sa  ReleaseData, Initialize, SetBufferedRegion */
  virtual void
  InitializeBufferedRegion();

  /** Directly computes an offset from the beginning of the buffer for a pixel
   * at the specified index.
   * The index is not checked as to whether it is inside the current buffer, so
   * the computed offset could conceivably be outside the buffer. If bounds
   * checking is needed, one can call \c ImageRegion::IsInside(ind) on the
   * BufferedRegion prior to calling ComputeOffset.
   * \warning unlike \c ComputeOffset(), this version does not incur a
   * virtual call. It's meant to be used only for \c itk::Image<>, \c
   * itk::VectorImage<> and \c itk::SpecialCoordinatesImage<>.
   */
  OffsetValueType
  FastComputeOffset(const IndexType & ind) const
  {
    OffsetValueType offset = 0;
    ImageHelper<VImageDimension, VImageDimension>::ComputeOffset(
      Self::GetBufferedRegion().GetIndex(), ind, m_OffsetTable, offset);
    return offset;
  }

  /** Directly computes the index of the pixel at a specified offset from the
   * beginning of the buffered region.
   * Bounds checking is not performed. Thus, the computed index could be
   * outside the BufferedRegion. To ensure a valid index, the parameter
   * \c offset should be between 0 and the number of pixels in the
   * BufferedRegion (the latter can be found using
   * \c ImageRegion::GetNumberOfPixels()).
   * \warning unlike \c ComputeOffset(), this version does not incur a
   * virtual call. It's meant to be used only for \c itk::Image<>, \c
   * itk::VectorImage<> and \c itk::SpecialCoordinatesImage<>.
   */
  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