Comparing ITK 4.13 implementation of GaussianDerivativeImageFunction::EvaluateAtIndex to http://review.source.kitware.com/#/c/22926/ (WIP)

GaussianDerivativeImageFunction_old_versus_WIP22926.cpp

// Author: Niels Dekker, LKEB, Leiden University Medical Center, The Netherlands

// Large parts of this code are from ITK, which has the following license:

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


#define _SCL_SECURE_NO_WARNINGS

#include <itkGaussianDerivativeImageFunction.h>
#include <itkImage.h>
#include <stdlib.h>  // For rand

namespace
{
  template< typename TInputImage, typename TOutput = double >
  class ImageFunctionWithoutBlur : public itk::GaussianDerivativeImageFunction<TInputImage, TOutput>
  {
  private:
    const double * const m_Sigma;
    mutable OperatorArrayType         m_OperatorArray;
    OperatorImageFunctionPointer m_OperatorImageFunction;
    const double * const m_Extent;
    GaussianDerivativeFunctionPointer m_GaussianDerivativeFunction;

    ImageFunctionWithoutBlur()
      :
      m_Sigma(GetSigma()),
      m_OperatorImageFunction(OperatorImageFunctionType::New()),
      m_Extent(GetExtent()),
      m_GaussianDerivativeFunction(GaussianDerivativeFunctionType::New())
    {
      m_GaussianDerivativeFunction->SetNormalized(false); // faster
    }


    void RecomputeGaussianDerivativeKernel()
    {
      unsigned int direction = 0;

      for (unsigned int op = 0; op < itkGetStaticConstMacro(ImageDimension2); ++op)
      {
        // Set the derivative of the Gaussian first
        OperatorNeighborhoodType dogNeighborhood;
        typename GaussianDerivativeFunctionType::InputType pt;
        typename NeighborhoodType::SizeType size;
        size.Fill(0);
        size[direction] = static_cast<itk::SizeValueType>(m_Sigma[direction] * m_Extent[direction]);
        dogNeighborhood.SetRadius(size);

        typename GaussianDerivativeFunctionType::ArrayType s;
        s[0] = 1;
        m_GaussianDerivativeFunction->SetSigma(s);

        typename OperatorNeighborhoodType::Iterator it = dogNeighborhood.Begin();

        unsigned int i = 0;
        while (it != dogNeighborhood.End())
        {
          pt[0] = dogNeighborhood.GetOffset(i)[direction];

          (*it) = m_GaussianDerivativeFunction->Evaluate(pt);
          ++i;
          ++it;
        }

        m_OperatorArray[op * 2] = dogNeighborhood;
        ++direction;
      }
    }


    void SetInputImage(const InputImageType *ptr) ITK_OVERRIDE
    {
      this->itk::GaussianDerivativeImageFunction<TInputImage, TOutput>::SetInputImage(ptr);
      m_OperatorImageFunction->SetInputImage(GetInputImage());
      RecomputeGaussianDerivativeKernel();
    }


    OutputType EvaluateAtIndex(const IndexType & index) const
    {
      OutputType gradient;

      for (unsigned int ii = 0; ii < itkGetStaticConstMacro(ImageDimension2); ++ii)
      {
        // Apply each Gaussian kernel to a subset of the image
        typedef typename OutputType::RealValueType OutputRealValueType;

        // Estimate derivative in the direction 'ii'
        const unsigned int idx = 2 * ii;
        const signed int center = (unsigned int)((m_OperatorArray[idx].GetSize()[ii] - 1) / 2);
        m_OperatorArray[idx][center] = 0;
        m_OperatorImageFunction->SetOperator(m_OperatorArray[idx]);
        OutputRealValueType value = m_OperatorImageFunction->EvaluateAtIndex(index);

        gradient[ii] = static_cast<typename OutputType::ComponentType>(value);
      }

      return gradient;
    }

    typedef ImageFunctionWithoutBlur Self;
  public:
    itkNewMacro(Self);
  };


  // Returns zero when EvaluateAt returns the same for both image functions.
  template <typename TImage, typename TOutput>
  int CompareEvaluateAtIndex(
    itk::ImageFunction<TImage, itk::Vector< TOutput, TImage::ImageDimension >, TOutput >& imageFunction1,
    itk::ImageFunction<TImage, itk::Vector< TOutput, TImage::ImageDimension >, TOutput >& imageFunction2)
  {
    int result = 0;
    const typename TImage::Pointer image = TImage::New();
    const unsigned imageSizeX = 65;
    const unsigned imageSizeY = 65;
    const typename TImage::SizeType imageSize = { { imageSizeX, imageSizeY } };
    image->SetRegions(imageSize);
    image->Allocate();
    
    typename TImage::PixelType* const ptr = image->GetBufferPointer();

    const unsigned numberOfPixels = imageSizeX * imageSizeY;

    for (unsigned i = 0; i < numberOfPixels; ++i)
    {
      ptr[i] = static_cast<typename TImage::PixelType>(rand());
    }

    imageFunction1.SetInputImage(image);
    imageFunction2.SetInputImage(image);

    for (unsigned y = 0; y < imageSizeY; ++y)
    {
      for (unsigned x = 0; x < imageSizeX; ++x)
      {
        const typename TImage::IndexType index = { { x, y } };
        
        const itk::Vector< TOutput, TImage::ImageDimension > vector1 = imageFunction1.EvaluateAtIndex(index);
        const itk::Vector< TOutput, TImage::ImageDimension > vector2 = imageFunction2.EvaluateAtIndex(index);

        if (vector1 != vector2)
        {
          ++result;
          std::cerr << "imageFunction1 and imageFunction1 produce different vectors!!!" << std::endl;
        }
        std::cout << vector1 << std::endl;
      }
    }
    return result;
  }
}


int main()
{
  typedef itk::Image<unsigned short> ImageType;

  return CompareEvaluateAtIndex<ImageType, double>(
    *ImageFunctionWithoutBlur<ImageType>::New(),
    *itk::GaussianDerivativeImageFunction<ImageType>::New()
    );
}