convolve.h

Go to the documentation of this file.

#ifndef _DLR_NUMERIC_CONVOLVE_H_
#define _DLR_NUMERIC_CONVOLVE_H_

#include <dlrCommon/tags.h>
#include <dlrNumeric/array1D.h>
#include <dlrNumeric/convolutionStrategy.h>

namespace dlr {

  namespace numeric {

    template <class Type>
    inline Array1D<Type>
    convolve(const Array1D<Type>& kernel,
             const Array1D<Type>& signal,
             ConvolutionStrategy strategy = DLR_CONVOLVE_ZERO_PAD_RESULT);


    template <class Type0, class Type1>
    inline Array1D<Type1>
    convolve(const Array1D<Type0>& kernel,
             const Array1D<Type0>& signal,
             type_tag<Type1> resultTag,
             ConvolutionStrategy strategy = DLR_CONVOLVE_ZERO_PAD_RESULT);
  

    template <class Type>
    inline Array1D<Type>
    correlate(const Array1D<Type>& kernel,
              const Array1D<Type>& signal,
              ConvolutionStrategy strategy = DLR_CONVOLVE_ZERO_PAD_RESULT);


    template <class Type0, class Type1>
    Array1D<Type1>
    correlate(const Array1D<Type0>& kernel,
              const Array1D<Type0>& signal,
              type_tag<Type1> resultTag,
              ConvolutionStrategy strategy = DLR_CONVOLVE_ZERO_PAD_RESULT);
  
  
  } // namespace numeric

} // namespace dlr


/* ======= Declarations to maintain compatibility with legacy code. ======= */

namespace dlr {

  using numeric::convolve;
  using numeric::convolve;
  using numeric::correlate;
  using numeric::correlate;

} // namespace dlr


/* ================ Definitions follow ================ */


#include <algorithm> // For std::reverse_copy()
#include <dlrNumeric/convolve1D.h>

namespace dlr {

  namespace numeric {

    namespace privateCode {

      template <class Type0, class Type1>
      Array1D<Type1>
      oldCorrelateReflectSignal(const Array1D<Type0>& kernel,
                                const Array1D<Type0>& signal,
                                type_tag<Type1> resultTag)
      {
        if(kernel.size() % 2 != 1) {
          DLR_THROW(ValueException, "correlate()",
                    "Argument kernel must have an odd number of elements.");
        }
        if(kernel.size() > signal.size()) {
          DLR_THROW(ValueException, "correlate()",
                    "Argument kernel must not have more elements than "
                    "argument signal.");
        }
        
        Array1D<Type1> result(kernel.size() + signal.size() - 1);
        for(size_t borderIndex = 0; borderIndex < kernel.size() - 1;
            ++borderIndex) {
          Type1 dotProduct0 = static_cast<Type1>(0);
          Type1 dotProduct1 = static_cast<Type1>(0);
          size_t kernelIndex = 0;
          size_t signalIndex = (kernel.size() - borderIndex) - 1;
          while(kernelIndex < kernel.size() - borderIndex - 1) {          
            dotProduct0 += kernel[kernelIndex] * signal[signalIndex];
            dotProduct1 +=
              kernel[(kernel.size() - kernelIndex) - 1]
              * signal[(signal.size() - signalIndex) - 1];
            ++kernelIndex;
            --signalIndex;
          }
          while(kernelIndex < kernel.size()) {          
            dotProduct0 += kernel[kernelIndex] * signal[signalIndex];
            dotProduct1 +=
              kernel[(kernel.size() - kernelIndex) - 1]
              * signal[(signal.size() - signalIndex) - 1];
            ++kernelIndex;
            ++signalIndex;
          }
          result[borderIndex] = dotProduct0;
          result[(result.size() - borderIndex) - 1] = dotProduct1;
        }

        correlate1DCommon<Type1, Type0, Type0>(
          kernel,
          signal.begin(),
          signal.end() - 2 * kernel.size() + 1,
          result.begin() + kernel.size() - 1);
        return result;
      }

    } // namespace privateCode
    
    template <class Type>
    inline Array1D<Type>
    convolve(const Array1D<Type>& kernel,
             const Array1D<Type>& signal,
             ConvolutionStrategy strategy)
    {
      return convolve(kernel, signal, type_tag<Type>(), strategy);
    }  

  
    template <class Type0, class Type1>
    inline Array1D<Type1>
    convolve(const Array1D<Type0>& kernel,
             const Array1D<Type0>& signal,
             type_tag<Type1> resultTag,
             ConvolutionStrategy strategy)
    {
      Array1D<Type0> reversedKernel(kernel.size());
      std::reverse_copy(kernel.begin(), kernel.end(), reversedKernel.begin());
      return correlate(reversedKernel, signal, resultTag, strategy);
    }
    

  
    template <class Type>
    inline Array1D<Type>
    correlate(const Array1D<Type>& kernel,
              const Array1D<Type>& signal,
              ConvolutionStrategy strategy)
    {
      return correlate(kernel, signal, type_tag<Type>(), strategy);
    }


    template <class Type0, class Type1>
    Array1D<Type1>
    correlate(const Array1D<Type0>& kernel,
              const Array1D<Type0>& signal,
              type_tag<Type1> resultTag,
              ConvolutionStrategy strategy)
    {
      switch(strategy) {
      case DLR_CONVOLVE_TRUNCATE_RESULT:
        return correlate1D<Type1>(
          kernel, signal, strategy, DLR_CONVOLVE_ROI_VALID);
        break;
      case DLR_CONVOLVE_ZERO_PAD_RESULT:
        return correlate1D<Type1>(
          kernel, signal, DLR_CONVOLVE_PAD_RESULT, DLR_CONVOLVE_ROI_SAME,
          static_cast<Type1>(0));
        break;
      case DLR_CONVOLVE_ZERO_PAD_SIGNAL:
      case DLR_CONVOLVE_WRAP_SIGNAL:
        return correlate1D<Type1>(
          kernel, signal, strategy, DLR_CONVOLVE_ROI_FULL);
        break;
      case DLR_CONVOLVE_REFLECT_SIGNAL:
        return privateCode::oldCorrelateReflectSignal(
          kernel, signal, type_tag<Type1>());
        break;
      default:
        DLR_THROW(LogicException, "correlate()",
                  "Illegal value for strategy argument.  "
                  "This function is deprecated.  "
                  "Please consider using convolve1D() or correlate1D() "
                  "instead.");
        break;
      }
      return Array1D<Type1>();
    }
  
  } // namespace numeric

} // namespace dlr

#endif // #ifndef _DLR_NUMERIC_CONVOLVE_H_