itkLinearInterpolateImageFunction.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 itkLinearInterpolateImageFunction_h
#define itkLinearInterpolateImageFunction_h
 
#include "itkInterpolateImageFunction.h"
#include "itkVariableLengthVector.h"
#include <algorithm> // For max.
 
namespace itk
{
template <typename TInputImage, typename TCoordinate = double>
class ITK_TEMPLATE_EXPORT LinearInterpolateImageFunction : public InterpolateImageFunction<TInputImage, TCoordinate>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(LinearInterpolateImageFunction);

  using Self = LinearInterpolateImageFunction;
  using Superclass = InterpolateImageFunction<TInputImage, TCoordinate>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkOverrideGetNameOfClassMacro(LinearInterpolateImageFunction);

  itkNewMacro(Self);

  using typename Superclass::OutputType;

  using typename Superclass::InputImageType;

  using typename Superclass::InputPixelType;

  using typename Superclass::RealType;

  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;

  using typename Superclass::IndexType;

  using typename Superclass::SizeType;

  using typename Superclass::ContinuousIndexType;
  using InternalComputationType = typename ContinuousIndexType::ValueType;

  OutputType
  EvaluateAtContinuousIndex(const ContinuousIndexType & index) const override
  {
    return this->EvaluateOptimized(Dispatch<ImageDimension>(), index);
  }
 
  SizeType
  GetRadius() const override
  {
    return SizeType::Filled(1);
  }
 
protected:
  LinearInterpolateImageFunction() = default;
  ~LinearInterpolateImageFunction() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;
 
private:
  struct DispatchBase
  {};
  template <unsigned int>
  struct Dispatch : public DispatchBase
  {};
 
  inline OutputType
  EvaluateOptimized(const Dispatch<0> &, const ContinuousIndexType &) const
  {
    return 0;
  }
 
  inline OutputType
  EvaluateOptimized(const Dispatch<1> &, const ContinuousIndexType & index) const
  {
    IndexType basei;
    basei[0] = std::max(Math::Floor<IndexValueType>(index[0]), this->m_StartIndex[0]);
 
    const InternalComputationType distance = index[0] - static_cast<InternalComputationType>(basei[0]);
 
    const TInputImage * const inputImagePtr = this->GetInputImage();
    const RealType &          val0 = inputImagePtr->GetPixel(basei);
    if (distance <= 0.)
    {
      return (static_cast<OutputType>(val0));
    }
 
    ++basei[0];
    if (basei[0] > this->m_EndIndex[0])
    {
      return (static_cast<OutputType>(val0));
    }
    const RealType & val1 = inputImagePtr->GetPixel(basei);
 
    return (static_cast<OutputType>(val0 + (val1 - val0) * distance));
  }
 
  inline OutputType
  EvaluateOptimized(const Dispatch<2> &, const ContinuousIndexType & index) const
  {
    IndexType basei;
 
    basei[0] = std::max(Math::Floor<IndexValueType>(index[0]), this->m_StartIndex[0]);
    const InternalComputationType distance0 = index[0] - static_cast<InternalComputationType>(basei[0]);
 
    basei[1] = std::max(Math::Floor<IndexValueType>(index[1]), this->m_StartIndex[1]);
    const InternalComputationType distance1 = index[1] - static_cast<InternalComputationType>(basei[1]);
 
    const TInputImage * const inputImagePtr = this->GetInputImage();
    const RealType &          val00 = inputImagePtr->GetPixel(basei);
    if (distance0 <= 0. && distance1 <= 0.)
    {
      return (static_cast<OutputType>(val00));
    }
    if (distance1 <= 0.) // if they have the same "y"
    {
      ++basei[0]; // then interpolate across "x"
      if (basei[0] > this->m_EndIndex[0])
      {
        return (static_cast<OutputType>(val00));
      }
      const RealType & val10 = inputImagePtr->GetPixel(basei);
      return (static_cast<OutputType>(val00 + (val10 - val00) * distance0));
    }
    else if (distance0 <= 0.) // if they have the same "x"
    {
      ++basei[1]; // then interpolate across "y"
      if (basei[1] > this->m_EndIndex[1])
      {
        return (static_cast<OutputType>(val00));
      }
      const RealType & val01 = inputImagePtr->GetPixel(basei);
      return (static_cast<OutputType>(val00 + (val01 - val00) * distance1));
    }
    // fall-through case:
    // interpolate across "xy"
    ++basei[0];
    if (basei[0] > this->m_EndIndex[0]) // interpolate across "y"
    {
      --basei[0];
      ++basei[1];
      if (basei[1] > this->m_EndIndex[1])
      {
        return (static_cast<OutputType>(val00));
      }
      const RealType & val01 = inputImagePtr->GetPixel(basei);
      return (static_cast<OutputType>(val00 + (val01 - val00) * distance1));
    }
    const RealType & val10 = inputImagePtr->GetPixel(basei);
 
    const RealType & valx0 = val00 + (val10 - val00) * distance0;
 
    ++basei[1];
    if (basei[1] > this->m_EndIndex[1]) // interpolate across "x"
    {
      return (static_cast<OutputType>(valx0));
    }
    const RealType & val11 = inputImagePtr->GetPixel(basei);
    --basei[0];
    const RealType & val01 = inputImagePtr->GetPixel(basei);
 
    const RealType & valx1 = val01 + (val11 - val01) * distance0;
 
    return (static_cast<OutputType>(valx0 + (valx1 - valx0) * distance1));
  }
 
  inline OutputType
  EvaluateOptimized(const Dispatch<3> &, const ContinuousIndexType & index) const
  {
    IndexType basei;
    basei[0] = std::max(Math::Floor<IndexValueType>(index[0]), this->m_StartIndex[0]);
    const InternalComputationType distance0 = index[0] - static_cast<InternalComputationType>(basei[0]);
 
    basei[1] = std::max(Math::Floor<IndexValueType>(index[1]), this->m_StartIndex[1]);
    const InternalComputationType distance1 = index[1] - static_cast<InternalComputationType>(basei[1]);
 
    basei[2] = std::max(Math::Floor<IndexValueType>(index[2]), this->m_StartIndex[2]);
    const InternalComputationType distance2 = index[2] - static_cast<InternalComputationType>(basei[2]);
 
    const TInputImage * const inputImagePtr = this->GetInputImage();
    const RealType &          val000 = inputImagePtr->GetPixel(basei);
    if (distance0 <= 0. && distance1 <= 0. && distance2 <= 0.)
    {
      return (static_cast<OutputType>(val000));
    }
 
    if (distance2 <= 0.)
    {
      if (distance1 <= 0.) // interpolate across "x"
      {
        ++basei[0];
        if (basei[0] > this->m_EndIndex[0])
        {
          return (static_cast<OutputType>(val000));
        }
        const RealType & val100 = inputImagePtr->GetPixel(basei);
 
        return static_cast<OutputType>(val000 + (val100 - val000) * distance0);
      }
      if (distance0 <= 0.) // interpolate across "y"
      {
        ++basei[1];
        if (basei[1] > this->m_EndIndex[1])
        {
          return (static_cast<OutputType>(val000));
        }
        const RealType & val010 = inputImagePtr->GetPixel(basei);
 
        return static_cast<OutputType>(val000 + (val010 - val000) * distance1);
      }
      else // interpolate across "xy"
      {
        ++basei[0];
        if (basei[0] > this->m_EndIndex[0]) // interpolate across "y"
        {
          --basei[0];
          ++basei[1];
          if (basei[1] > this->m_EndIndex[1])
          {
            return (static_cast<OutputType>(val000));
          }
          const RealType & val010 = inputImagePtr->GetPixel(basei);
          return static_cast<OutputType>(val000 + (val010 - val000) * distance1);
        }
        const RealType & val100 = inputImagePtr->GetPixel(basei);
        const RealType & valx00 = val000 + (val100 - val000) * distance0;
 
        ++basei[1];
        if (basei[1] > this->m_EndIndex[1]) // interpolate across "x"
        {
          return (static_cast<OutputType>(valx00));
        }
        const RealType & val110 = inputImagePtr->GetPixel(basei);
 
        --basei[0];
        const RealType & val010 = inputImagePtr->GetPixel(basei);
        const RealType & valx10 = val010 + (val110 - val010) * distance0;
 
        return static_cast<OutputType>(valx00 + (valx10 - valx00) * distance1);
      }
    }
    else
    {
      if (distance1 <= 0.)
      {
        if (distance0 <= 0.) // interpolate across "z"
        {
          ++basei[2];
          if (basei[2] > this->m_EndIndex[2])
          {
            return (static_cast<OutputType>(val000));
          }
          const RealType & val001 = inputImagePtr->GetPixel(basei);
 
          return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
        }
        // interpolate across "xz"
        ++basei[0];
        if (basei[0] > this->m_EndIndex[0]) // interpolate across "z"
        {
          --basei[0];
          ++basei[2];
          if (basei[2] > this->m_EndIndex[2])
          {
            return (static_cast<OutputType>(val000));
          }
          const RealType & val001 = inputImagePtr->GetPixel(basei);
 
          return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
        }
        const RealType & val100 = inputImagePtr->GetPixel(basei);
 
        const RealType & valx00 = val000 + (val100 - val000) * distance0;
 
        ++basei[2];
        if (basei[2] > this->m_EndIndex[2]) // interpolate across "x"
        {
          return (static_cast<OutputType>(valx00));
        }
        const RealType & val101 = inputImagePtr->GetPixel(basei);
 
        --basei[0];
        const RealType & val001 = inputImagePtr->GetPixel(basei);
 
        const RealType & valx01 = val001 + (val101 - val001) * distance0;
 
        return static_cast<OutputType>(valx00 + (valx01 - valx00) * distance2);
      }
      else if (distance0 <= 0.) // interpolate across "yz"
      {
        ++basei[1];
        if (basei[1] > this->m_EndIndex[1]) // interpolate across "z"
        {
          --basei[1];
          ++basei[2];
          if (basei[2] > this->m_EndIndex[2])
          {
            return (static_cast<OutputType>(val000));
          }
          const RealType & val001 = inputImagePtr->GetPixel(basei);
 
          return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
        }
        const RealType & val010 = inputImagePtr->GetPixel(basei);
 
        const RealType & val0x0 = val000 + (val010 - val000) * distance1;
 
        ++basei[2];
        if (basei[2] > this->m_EndIndex[2]) // interpolate across "y"
        {
          return (static_cast<OutputType>(val0x0));
        }
        const RealType & val011 = inputImagePtr->GetPixel(basei);
 
        --basei[1];
        const RealType & val001 = inputImagePtr->GetPixel(basei);
 
        const RealType & val0x1 = val001 + (val011 - val001) * distance1;
 
        return static_cast<OutputType>(val0x0 + (val0x1 - val0x0) * distance2);
      }
      else // interpolate across "xyz"
      {
        ++basei[0];
        if (basei[0] > this->m_EndIndex[0]) // interpolate across "yz"
        {
          --basei[0];
          ++basei[1];
          if (basei[1] > this->m_EndIndex[1]) // interpolate across "z"
          {
            --basei[1];
            ++basei[2];
            if (basei[2] > this->m_EndIndex[2])
            {
              return (static_cast<OutputType>(val000));
            }
            const RealType & val001 = inputImagePtr->GetPixel(basei);
 
            return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
          }
          const RealType & val010 = inputImagePtr->GetPixel(basei);
          const RealType & val0x0 = val000 + (val010 - val000) * distance1;
 
          ++basei[2];
          if (basei[2] > this->m_EndIndex[2]) // interpolate across "y"
          {
            return (static_cast<OutputType>(val0x0));
          }
          const RealType & val011 = inputImagePtr->GetPixel(basei);
 
          --basei[1];
          const RealType & val001 = inputImagePtr->GetPixel(basei);
 
          const RealType & val0x1 = val001 + (val011 - val001) * distance1;
 
          return static_cast<OutputType>(val0x0 + (val0x1 - val0x0) * distance2);
        }
        const RealType & val100 = inputImagePtr->GetPixel(basei);
 
        const RealType & valx00 = val000 + (val100 - val000) * distance0;
 
        ++basei[1];
        if (basei[1] > this->m_EndIndex[1]) // interpolate across "xz"
        {
          --basei[1];
          ++basei[2];
          if (basei[2] > this->m_EndIndex[2]) // interpolate across "x"
          {
            return (static_cast<OutputType>(valx00));
          }
          const RealType & val101 = inputImagePtr->GetPixel(basei);
 
          --basei[0];
          const RealType & val001 = inputImagePtr->GetPixel(basei);
 
          const RealType & valx01 = val001 + (val101 - val001) * distance0;
 
          return static_cast<OutputType>(valx00 + (valx01 - valx00) * distance2);
        }
        const RealType & val110 = inputImagePtr->GetPixel(basei);
 
        --basei[0];
        const RealType & val010 = inputImagePtr->GetPixel(basei);
 
        const RealType & valx10 = val010 + (val110 - val010) * distance0;
 
        const RealType & valxx0 = valx00 + (valx10 - valx00) * distance1;
 
        ++basei[2];
        if (basei[2] > this->m_EndIndex[2]) // interpolate across "xy"
        {
          return (static_cast<OutputType>(valxx0));
        }
        const RealType & val011 = inputImagePtr->GetPixel(basei);
 
        ++basei[0];
        const RealType & val111 = inputImagePtr->GetPixel(basei);
 
        --basei[1];
        const RealType & val101 = inputImagePtr->GetPixel(basei);
 
        --basei[0];
        const RealType & val001 = inputImagePtr->GetPixel(basei);
 
        const RealType & valx01 = val001 + (val101 - val001) * distance0;
        const RealType & valx11 = val011 + (val111 - val011) * distance0;
        const RealType & valxx1 = valx01 + (valx11 - valx01) * distance1;
 
        return (static_cast<OutputType>(valxx0 + (valxx1 - valxx0) * distance2));
      }
    }
  }
 
  inline OutputType
  EvaluateOptimized(const DispatchBase &, const ContinuousIndexType & index) const
  {
    return this->EvaluateUnoptimized(index);
  }

  virtual inline OutputType
  EvaluateUnoptimized(const ContinuousIndexType & index) const;

  template <typename RealTypeScalarRealType>
  void
  MakeZeroInitializer(const TInputImage * const                      inputImagePtr,
                      VariableLengthVector<RealTypeScalarRealType> & tempZeros) const
  {
    // Variable length vector version to get the size of the pixel correct.
    constexpr typename TInputImage::IndexType idx = { { 0 } };
    const typename TInputImage::PixelType &   tempPixel = inputImagePtr->GetPixel(idx);
    const unsigned int                        sizeOfVarLengthVector = tempPixel.GetSize();
    tempZeros.SetSize(sizeOfVarLengthVector);
    tempZeros.Fill(RealTypeScalarRealType{});
  }
 
  template <typename RealTypeScalarRealType>
  void
  MakeZeroInitializer(const TInputImage * const itkNotUsed(inputImagePtr), RealTypeScalarRealType & tempZeros) const
  {
    // All other cases
    tempZeros = RealTypeScalarRealType{};
  }
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkLinearInterpolateImageFunction.hxx"
#endif
 
#endif