convolveND.h

Go to the documentation of this file.

#ifndef DLRNUMERIC_CONVOLVEND_H
#define DLRNUMERIC_CONVOLVEND_H

#include <iostream>
#include <string>
#include <dlrCommon/exception.h>
#include <dlrNumeric/arrayND.h>
#include <dlrNumeric/convolutionStrategy.h>


namespace dlr {

  namespace numeric {
    

    template <class OutputType, class AccumulatorType,
              class KernelType, class SignalType, size_t Dimension>
    ArrayND<Dimension, OutputType>
    convolve(const Array1D<KernelType>& kernel,
             const ArrayND<Dimension, SignalType>& signal,
             size_t axis,
             ConvolutionStrategy strategy = DLR_CONVOLVE_PAD_RESULT,
             ConvolutionROI roi = DLR_CONVOLVE_ROI_SAME);

    
  } // namespace numeric

} // namespace dlr



/* ============ Inline and template definitions follow. ============ */


namespace dlr {

  namespace numeric {

    namespace privateCode {

      inline bool
      convolveNDUpdatePosition(Array1D<size_t>& position,
                               Array1D<size_t> const& lowerBound,
                               Array1D<size_t> const& upperBound)
      {
        for(size_t ii = position.size() - 1; ii < position.size(); --ii) {
          if(++(position[ii]) < upperBound[ii]) {
            return true;
          }
          position[ii] = lowerBound[ii];
        }
        return false;
      }

    } // namespace privateCode
    
    
    template <class OutputType, class AccumulatorType,
              class KernelType, class SignalType, size_t Dimension>
    ArrayND<Dimension, OutputType>
    convolve(const Array1D<KernelType>& kernel,
             const ArrayND<Dimension, SignalType>& signal,
             size_t axis,
             ConvolutionStrategy strategy = DLR_CONVOLVE_PAD_RESULT,
             ConvolutionROI roi = DLR_CONVOLVE_ROI_SAME)
    {
      Array1D<AccumulatorType> reversedKernel(kernel.size());
      for(size_t ii = 0; ii < kernel.size(); ++ii) {
        reversedKernel[kernel.size() - ii - 1] =
          static_cast<AccumulatorType>(kernel[ii]);
      }
      
      ArrayND<Dimension, OutputType> result(signal.getShape());
      result = OutputType(0);
      
      if(strategy != DLR_CONVOLVE_PAD_RESULT) {
        DLR_THROW(NotImplementedException, "convolve()",
                  "Only ConvolutionStrategy DLR_CONVOLVE_PAD_RESULT is "
                  "currently implemented.");
      }
      if(roi != DLR_CONVOLVE_ROI_SAME) {
        DLR_THROW(NotImplementedException, "convolve()",
                  "Only ConvolutionROI DLR_CONVOLVE_ROI_SAME is "
                  "currently implemented.");
      }

      // We'll be using 1D indexing into the ArrayND instance.  This
      // stride tells us how much we need to increment the index in
      // order to move one element along the requested axis, and
      // resultOffset tells us the offset between the first element of
      // the convolution kernel and the center of the convolution
      // kernel.
      size_t stride = signal.getStride(axis);
      size_t resultOffset = (reversedKernel.size() / 2) * stride;

      // This array keeps track of where we are in the ND array.  It
      // will be stepped in raster order through each position at
      // which the convolution kernel is applied to signal.
      Array1D<size_t> elementPosition(Dimension);
      elementPosition = 0;

      // Establish upper and lower bounds on elementPosition for valid
      // indexing.
      Array1D<size_t> positionLowerBound = elementPosition.copy();
      Array1D<size_t> positionUpperBound = signal.getShape().copy();
      positionUpperBound[axis] -= (reversedKernel.size() - 1);

      do {
        size_t signalIndex = signal.flattenIndex(elementPosition);
        size_t resultIndex = signalIndex + resultOffset;
        
        AccumulatorType accumulator = AccumulatorType(0);
        for(size_t ii = 0; ii < reversedKernel.size(); ++ii) {
          accumulator += (reversedKernel[ii]
                          * static_cast<AccumulatorType>(signal(signalIndex)));
          signalIndex += stride;
        }
        result(resultIndex) = static_cast<OutputType>(accumulator);
      } while(privateCode::convolveNDUpdatePosition(
                elementPosition, positionLowerBound, positionUpperBound));
      return result;
    }
    

  } // namespace numeric

} // namespace dlr

#endif /* #ifdef DLRNUMERIC_CONVOLVEND_H */