itkImageBase.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  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.
 *
 *=========================================================================*/
/*=========================================================================
 *
 *  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>
#if VXL_VERSION_DATE_FULL < 20160229
#include "vnl/vnl_matrix_fixed.txx" // Get the templates
#else
#include "vnl/vnl_matrix_fixed.hxx" // Get the templates
#endif

#include "itkImageTransformHelper.h"

namespace itk
{

/* Forward declaration (ImageTransformHelper include's ImageBase) */
template< unsigned int NImageDimension, unsigned int R, unsigned int C, typename TPointValue, typename TMatrixValue >
class ImageTransformHelper;

template< unsigned int VImageDimension = 2 >
class ImageBase:public DataObject
{
public:
  typedef ImageBase                  Self;
  typedef DataObject                 Superclass;
  typedef SmartPointer< Self >       Pointer;
  typedef SmartPointer< const Self > ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(ImageBase, DataObject);

  typedef unsigned int ImageDimensionType;

  itkStaticConstMacro(ImageDimension, ImageDimensionType, VImageDimension);

  typedef Index< VImageDimension >           IndexType;
  typedef typename IndexType::IndexValueType IndexValueType;

  typedef Offset< VImageDimension >            OffsetType;
  typedef typename OffsetType::OffsetValueType OffsetValueType;

  typedef Size< VImageDimension >          SizeType;
  typedef typename SizeType::SizeValueType SizeValueType;

  typedef ImageRegion< VImageDimension > RegionType;

  typedef SpacePrecisionType                          SpacingValueType;
  typedef Vector< SpacingValueType, VImageDimension > SpacingType;

  typedef SpacePrecisionType                       PointValueType;
  typedef Point< PointValueType, VImageDimension > PointType;

  typedef Matrix< SpacePrecisionType, VImageDimension, VImageDimension > DirectionType;

  virtual void Initialize() ITK_OVERRIDE;

  static 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 initialize=false);

  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);

  virtual void SetRequestedRegion( const DataObject *data ) ITK_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->SetLargestPossibleRegion(region);
    this->SetBufferedRegion(region);
    this->SetRequestedRegion(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 TCoordRep >
  bool TransformPhysicalPointToIndex(
    const Point< TCoordRep, VImageDimension > & point,
    IndexType & index) const
  {
    ImageTransformHelper< VImageDimension,VImageDimension - 1, VImageDimension - 1, TCoordRep, SpacePrecisionType >
      ::TransformPhysicalPointToIndex(this->m_PhysicalPointToIndex, this->m_Origin, point, index);

    // Now, check to see if the index is within allowed bounds
    const bool isInside = this->GetLargestPossibleRegion().IsInside(index);
    return isInside;
    /* NON TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING data version
     * Leaving here for documentation purposes
     * template< typename TCoordRep >
     * bool TransformPhysicalPointToIndex(
     *   const Point< TCoordRep, VImageDimension > & point,
     *   IndexType & index) const
     * {
     *   for ( unsigned int i = 0; i < VImageDimension; i++ )
     *     {
     *     TCoordRep sum = NumericTraits< TCoordRep >::ZeroValue();
     *     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);
     *     }
     *   // Now, check to see if the index is within allowed bounds
     *   const bool isInside = this->GetLargestPossibleRegion().IsInside(index);
     *   return isInside;
     * }
     */
  }

  template< typename TCoordRep, typename TIndexRep >
  bool TransformPhysicalPointToContinuousIndex(
    const Point< TCoordRep, VImageDimension > & point,
    ContinuousIndex< TIndexRep, VImageDimension > & index) const
  {
    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] );
      }

    // Now, check to see if the index is within allowed bounds
    const bool isInside = this->GetLargestPossibleRegion().IsInside(index);

    return isInside;
  }

  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 = NumericTraits< TCoordRep >::ZeroValue();
      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 TCoordRep >
  void TransformIndexToPhysicalPoint(
    const IndexType & index,
    Point< TCoordRep, VImageDimension > & point) const
  {
    ImageTransformHelper< VImageDimension, VImageDimension - 1, VImageDimension - 1,TCoordRep, SpacePrecisionType >::
      TransformIndexToPhysicalPoint(this->m_IndexToPhysicalPoint, this->m_Origin, index, point);
    /* NON TEMPLATE_META_PROGRAMMING_LOOP_UNROLLING data version
     * Leaving here for documentation purposes
     * 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];
     *       }
     *     }
     * }
     */
  }

  template< typename TCoordRep >
  void TransformLocalVectorToPhysicalVector(
    const FixedArray< TCoordRep, VImageDimension > & inputGradient,
    FixedArray< TCoordRep, VImageDimension > & outputGradient) const
  {
    //
    //TODO: This temporary implementation should be replaced with Template
    // MetaProgramming.
    //
    const DirectionType & direction = this->GetDirection();

    for ( unsigned int i = 0; i < VImageDimension; ++i )
      {
      typedef typename NumericTraits< TCoordRep >::AccumulateType CoordSumType;
      CoordSumType sum = NumericTraits< CoordSumType >::ZeroValue();
      for ( unsigned int j = 0; j < VImageDimension; ++j )
        {
        sum += direction[i][j] * inputGradient[j];
        }
      outputGradient[i] = static_cast< TCoordRep >( sum );
      }
  }

  template< typename TCoordRep >
  void TransformPhysicalVectorToLocalVector(
    const FixedArray< TCoordRep, VImageDimension > & inputGradient,
    FixedArray< TCoordRep, VImageDimension > & outputGradient) const
  {
    //
    //TODO: This temporary implementation should be replaced with Template
    // MetaProgramming.
    //
    const DirectionType & inverseDirection = this->GetInverseDirection();

    for ( unsigned int i = 0; i < VImageDimension; ++i )
      {
      typedef typename NumericTraits< TCoordRep >::AccumulateType CoordSumType;
      CoordSumType sum = NumericTraits< CoordSumType >::ZeroValue();
      for ( unsigned int j = 0; j < VImageDimension; ++j )
        {
        sum += inverseDirection[i][j] * inputGradient[j];
        }
      outputGradient[i] = static_cast< TCoordRep >( sum );
      }
  }

  virtual void CopyInformation(const DataObject *data) ITK_OVERRIDE;

  virtual void Graft(const DataObject *data) ITK_OVERRIDE;

  virtual void UpdateOutputInformation() ITK_OVERRIDE;

  virtual void UpdateOutputData() ITK_OVERRIDE;

  virtual void SetRequestedRegionToLargestPossibleRegion() ITK_OVERRIDE;

  virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() ITK_OVERRIDE;

  virtual bool VerifyRequestedRegion() ITK_OVERRIDE;

  virtual unsigned int GetNumberOfComponentsPerPixel() const;
  virtual void SetNumberOfComponentsPerPixel(unsigned int);

protected:
  ImageBase();
  ~ImageBase();
  virtual void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE;

  void ComputeOffsetTable();

  virtual void ComputeIndexToPhysicalPointMatrices();

protected:
  SpacingType   m_Spacing;
  PointType     m_Origin;
  DirectionType m_Direction;
  DirectionType m_InverseDirection;

  DirectionType m_IndexToPhysicalPoint;
  DirectionType m_PhysicalPointToIndex;

  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;
    }

private:
  ImageBase(const Self &) ITK_DELETE_FUNCTION;
  void operator=(const Self &) ITK_DELETE_FUNCTION;

  void InternalSetSpacing(const SpacingValueType spacing[VImageDimension])
    {
      SpacingType s(spacing);
      this->SetSpacing(s);
    }

  template <typename TSpacingValue>
  void InternalSetSpacing(const TSpacingValue spacing[VImageDimension])
    {
      Vector<TSpacingValue,VImageDimension> sf(spacing);
      SpacingType                           s;
      s.CastFrom(sf);
      this->SetSpacing(s);
    }

  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