Sequence2.hpp

Go to the documentation of this file.

// Copyright (c) 2011 Michael D. Adams
// All rights reserved.

// __START_OF_LICENSE__
// 
// Copyright (c) 2015 Michael D. Adams
// All rights reserved.
// 
// This file is part of the Signal Processing Library (SPL).
// 
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 3,
// or (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public
// License along with this program; see the file LICENSE.  If not,
// see <http://www.gnu.org/licenses/>.
// 
// __END_OF_LICENSE__

#ifndef SPL_Sequence2_hpp
#define SPL_Sequence2_hpp

// Configuration

//#define SPL_SEQUENCE2_DEBUG

#define SPL_SEQUENCE2_USE_NEW_CONV

#if defined(SPL_SEQUENCE2_DEBUG)
#define SPL_SEQUENCE2_INLINE
#else
#define SPL_SEQUENCE2_INLINE inline
#endif

// Include files.

#include <SPL/config.hpp>
#include <iostream>
#include <vector>
#include <SPL/Array2.hpp>
#include <SPL/Sequence.hpp>
#include <SPL/Sequence1.hpp>
#include <SPL/Sequence.hpp>
#include <SPL/math.hpp>

namespace SPL {

// Types.

template <class T>
class Sequence2
{
public:

        // Types

        typedef T ElemType;

        typedef typename Array2<T>::ConstIterator ConstIterator;

        typedef typename Array2<T>::Iterator Iterator;

        typedef typename Array2<T>::ConstXIterator ConstXIterator;

        typedef typename Array2<T>::XIterator XIterator;

        typedef typename Array2<T>::ConstYIterator ConstYIterator;

        typedef typename Array2<T>::YIterator YIterator;

        // Constructors and destructor

        Sequence2();

        Sequence2(int startX, int startY, int width, int height);

        Sequence2(int startX, int startY, int width, int height, const T& data);

        template <class InputIterator>
        Sequence2(int startX, int startY, int width, int height,
          InputIterator data);

        Sequence2(const Sequence2& f);

        template <class OtherT>
        Sequence2(const Sequence2<OtherT>& f);

        Sequence2(const Array2<T>& data);

        Sequence2(int startX, int startY, const Array2<T>& data);

        ~Sequence2();

        // Assignment operators

        Sequence2& operator=(const Sequence2& f);

        template <class OtherT>
        Sequence2& operator=(const Sequence2<OtherT>& f);

        // Compound assignment operators

        Sequence2& operator+=(const Sequence2& f);

        Sequence2& operator-=(const Sequence2& f);

        Sequence2& operator*=(const Sequence2& f);

        Sequence2& operator/=(const Sequence2& f);

        Sequence2& operator+=(const T& value);

        Sequence2& operator-=(const T& value);

        Sequence2& operator*=(const T& value);

        Sequence2& operator/=(const T& value);

        // Basic queries

        int getStartX() const;

        int getStartY() const;

        int getEndX() const;

        int getEndY() const;

        int getWidth() const;

        int getHeight() const;

        int getSize() const;

        bool isShared() const;

        // Accessors

        T& operator()(int x, int y);

        const T& operator()(int x, int y) const;

        // Iterators

        ConstIterator begin() const;

        Iterator begin();

        ConstIterator end() const;

        Iterator end();

        ConstXIterator rowBegin(int y) const;

        XIterator rowBegin(int y);

        ConstXIterator rowEnd(int y) const;

        XIterator rowEnd(int y);

        ConstYIterator colBegin(int x) const;

        YIterator colBegin(int x);

        ConstYIterator colEnd(int x) const;

        YIterator colEnd(int x);

        // Resizing

        // Get various statistics

        T min() const;

        T max() const;

        T sum() const;

        // Input/ouput

        std::ostream& output(std::ostream& out, int fieldWidth) const;

        // Miscellaneous

        Array2<T> getArray() const;

        void swapArray(Array2<T>& data);

        void fill(const T& value);

        Sequence2& translate(int x, int y);

private:

        template <class> friend class Sequence2;

        template <class U>
        friend std::ostream& operator<<(std::ostream& out, const Sequence2<U>&);
        template <class U>
        friend Sequence2<U> translate(const Sequence2<U>&, int x, int y);

        typedef Array2<T> Array;

        int startX_;

        int startY_;

        Array data_;

};

//


void combineDomains(int firstStartX, int firstStartY, int firstEndX,
  int firstEndY, int secondStartX, int secondStartY, int secondEndX,
  int secondEndY, int& startX, int& startY, int& endX, int& endY);


// Constructors and destructor

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2() : startX_(0), startY_(0),
  data_()
{
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(int startX, int startY,
  int width, int height) : startX_(startX), startY_(startY),
  data_(width, height)
{
        assert(width >= 0 && height >= 0);
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(int startX, int startY,
  int width, int height, const T& value) : startX_(startX), startY_(startY),
  data_(width, height, value)
{
        assert(width >= 0 && height >= 0);
}

template <class T>
template <class InputIterator>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(int startX, int startY,
  int width, int height, InputIterator data) : startX_(startX),
  startY_(startY), data_(width, height, data)
{
        assert(width >= 0 && height >= 0);
        assert(getWidth() == width && getHeight() == height);
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(const Sequence2<T>& f) :
  startX_(f.startX_), startY_(f.startY_), data_(f.data_)
{
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(const Array2<T>& data) :
  startX_(0), startY_(0), data_(data)
{
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(int startX, int startY,
  const Array2<T>& data) :
  startX_(startX), startY_(startY), data_(data)
{
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>::~Sequence2()
{
}

template <class T>
template <class OtherT>
SPL_SEQUENCE2_INLINE Sequence2<T>::Sequence2(const Sequence2<OtherT>& f) :
  startX_(f.startX_), startY_(f.startY_), data_(f.data_)
{
}

// Assignment operators

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator=(
  const Sequence2<T>& f)
{
        startX_ = f.startX_;
        startY_ = f.startY_;
        data_ = f.data_;
        return *this;
}

template <class T>
template <class OtherT>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator=(
  const Sequence2<OtherT>& f)
{
        startX_ = f.startX_;
        startY_ = f.startY_;
        data_ = f.data_;
        return *this;
}

// Compound assignment operators

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator+=(
  const Sequence2<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((!getSize() && !f.getSize()) ||
        (getStartX() == f.getStartX() && getEndX() == f.getEndX() &&
        getStartY() == f.getStartY() && getEndY() == f.getEndY()));

        data_ += f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator-=(
  const Sequence2<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((!getSize() && !f.getSize()) ||
        (getStartX() == f.getStartX() && getEndX() == f.getEndX() &&
        getStartY() == f.getStartY() && getEndY() == f.getEndY()));

        data_ -= f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator*=(
  const Sequence2<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((!getSize() && !f.getSize()) ||
        (getStartX() == f.getStartX() && getEndX() == f.getEndX() &&
        getStartY() == f.getStartY() && getEndY() == f.getEndY()));

        data_ *= f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator/=(
  const Sequence2<T>& f)
{
        // Ensure that both sequences have the same domain.
        assert((!getSize() && !f.getSize()) ||
        (getStartX() == f.getStartX() && getEndX() == f.getEndX() &&
        getStartY() == f.getStartY() && getEndY() == f.getEndY()));

        data_ /= f.data_;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator+=(
  const T& value)
{
        data_ += value;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator-=(
  const T& value)
{
        data_ -= value;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator*=(
  const T& value)
{
        data_ *= value;
        return *this;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::operator/=(
  const T& value)
{
        data_ /= value;
        return *this;
}

// Basic queries

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getStartX() const
{
        return startX_;
}

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getStartY() const
{
        return startY_;
}

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getEndX() const
{
        return startX_ + getWidth();
}

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getEndY() const
{
        return startY_ + getHeight();
}

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getWidth() const
{
        return data_.getWidth();
}

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getHeight() const
{
        return data_.getHeight();
}

template <class T>
SPL_SEQUENCE2_INLINE int Sequence2<T>::getSize() const
{
        return data_.getWidth() * data_.getHeight();
}

template <class T>
SPL_SEQUENCE2_INLINE bool Sequence2<T>::isShared() const
{
        return data_.isShared();
}

// Accessors

template <class T>
SPL_SEQUENCE2_INLINE const T& Sequence2<T>::operator()(int x, int y) const
{
    assert(x >= startX_ && x < getEndX() && y >= startY_ && y < getEndY());
    return data_(x - startX_, y - startY_);
}

template <class T>
SPL_SEQUENCE2_INLINE T& Sequence2<T>::operator()(int x, int y)
{
    assert(x >= startX_ && x < getEndX() && y >= startY_ && y < getEndY());
    return data_(x - startX_, y - startY_);
}

// Iterators

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::ConstIterator Sequence2<T>::begin()
  const
{
        return data_.begin();
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::Iterator Sequence2<T>::begin()
{
        return data_.begin();
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::ConstIterator Sequence2<T>::end()
  const
{
        return data_.end();
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::Iterator Sequence2<T>::end()
{
        return data_.end();
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::ConstXIterator
  Sequence2<T>::rowBegin(int y) const
{
        assert(y >= startY_ && y < getEndY());
        return data_.rowBegin(y - startY_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::XIterator
  Sequence2<T>::rowBegin(int y)
{
        assert(y >= startY_ && y < getEndY());
        return data_.rowBegin(y - startY_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::ConstXIterator
  Sequence2<T>::rowEnd(int y) const
{
        assert(y >= startY_ && y < getEndY());
        return data_.rowEnd(y - startY_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::XIterator
  Sequence2<T>::rowEnd(int y)
{
        assert(y >= startY_ && y < getEndY());
        return data_.rowEnd(y - startY_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::ConstYIterator
  Sequence2<T>::colBegin(int x) const
{
        assert(x >= startX_ && x < getEndX());
        return data_.colBegin(x - startX_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::YIterator
  Sequence2<T>::colBegin(int x)
{
        assert(x >= startX_ && x < getEndX());
        return data_.colBegin(x - startX_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::ConstYIterator
  Sequence2<T>::colEnd(int x) const
{
        assert(x >= startX_ && x < getEndX());
        return data_.colEnd(x - startX_);
}

template <class T>
SPL_SEQUENCE2_INLINE typename Sequence2<T>::YIterator
  Sequence2<T>::colEnd(int x)
{
        assert(x >= startX_ && x < getEndX());
        return data_.colEnd(x - startX_);
}

// Resizing

// Get various statistics

template <class T>
SPL_SEQUENCE2_INLINE T Sequence2<T>::min() const
{
        assert(getSize() > 0);
        return data_.min();
}

template <class T>
SPL_SEQUENCE2_INLINE T Sequence2<T>::max() const
{
        assert(getSize() > 0);
        return data_.max();
}

template <class T>
SPL_SEQUENCE2_INLINE T Sequence2<T>::sum() const
{
        return data_.sum();
}

// Input/output

template <class T>
std::ostream& Sequence2<T>::output(std::ostream& out, int fieldWidth) const
{
        out << startX_ << " " << startY_ << " ";
        data_.output(out, fieldWidth);
        return out;
}

template <class T>
std::ostream& operator<<(std::ostream& out, const Sequence2<T>& f)
{
        out << f.getStartX() << " " << f.getStartY() << " ";
        out << f.data_;
        return out;
}

template <class T>
std::istream& operator>>(std::istream& in, Sequence2<T>& f)
{
        int startX;
        int startY;
        int width;
        int height;
        if (!(in >> startX)) {
                return in;
        }
        if (!(in >> startY)) {
                return in;
        }
        if (!(in >> width)) {
                return in;
        }
        if (!(in >> height)) {
                return in;
        }
        typename std::vector<T> data;
        data.reserve(width * height);
        for (int y = 0; y < height; ++y) {
                for (int x = 0; x < width; ++x) {
                        T value;
                        if (!(in >> value)) {
                                return in;
                        }
                        data.push_back(value);
                }
        }
        f = Sequence2<T>(startX, startY, width, height, data.begin());
        return in;
}

// Miscellaneous

template<class T>
SPL_SEQUENCE2_INLINE void Sequence2<T>::fill(const T& value)
{
        data_.fill(value);
}

template <class T>
SPL_SEQUENCE2_INLINE void Sequence2<T>::swapArray(Array2<T>& data)
{
        data_.swap(data);
}

template<class T>
SPL_SEQUENCE2_INLINE Array2<T> Sequence2<T>::getArray() const
{
        return data_;
}

template<class T>
SPL_SEQUENCE2_INLINE Sequence2<T>& Sequence2<T>::translate(int x, int y)
{
        startX_ += x;
        startY_ += y;
        return *this;
}

//

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator+(const Sequence2<T>& f,
  const Sequence2<T>& g)
{
        return Sequence2<T>(f) += g;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator-(const Sequence2<T>& f,
  const Sequence2<T>& g)
{
        return Sequence2<T>(f) -= g;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator*(const Sequence2<T>& f,
  const Sequence2<T>& g)
{
        return Sequence2<T>(f) *= g;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator/(const Sequence2<T>& f,
  const Sequence2<T>& g)
{
        return Sequence2<T>(f) /= g;
}

template <class T>
Sequence2<T> add(const Sequence2<T>& f, const Sequence2<T>& g)
{
        int startX;
        int startY;
        int endX;
        int endY;
        combineDomains(f.getStartX(), f.getStartY(), f.getEndX(), f.getEndY(),
          g.getStartX(), g.getStartY(), g.getEndX(), g.getEndY(), startX, startY,
          endX, endY);
        Sequence2<T> result(startX, startY, endX - startX, endY - startY, T(0));
        int xMin = std::max(f.getStartX(), startX);
        int xMax = std::min(f.getEndX(), endX) - 1;
        int yMin = std::max(f.getStartY(), startY);
        int yMax = std::min(f.getEndY(), endY) - 1;
        for (int y = yMin; y <= yMax; ++y) {
                for (int x = xMin; x <= xMax; ++x) {
                        result(x, y) += f(x, y);
                }
        }
        xMin = std::max(g.getStartX(), startX);
        xMax = std::min(g.getEndX(), endX) - 1;
        yMin = std::max(g.getStartY(), startY);
        yMax = std::min(g.getEndY(), endY) - 1;
        for (int y = yMin; y <= yMax; ++y) {
                for (int x = xMin; x <= xMax; ++x) {
                        result(x, y) += g(x, y);
                }
        }
        return result;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator+(const T& value,
  const Sequence2<T>& f)
{
        return Sequence2<T>(f) += value;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator+(const Sequence2<T>& f,
  const T& value)
{
        return Sequence2<T>(f) += value;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator-(const Sequence2<T>& f,
  const T& value)
{
        return Sequence2<T>(f) -= value;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator*(const T& value,
  const Sequence2<T>& f)
{
        return Sequence2<T>(f) *= value;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator*(const Sequence2<T>& f,
  const T& value)
{
        return Sequence2<T>(f) *= value;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> operator/(const Sequence2<T>& f,
  const T& value)
{
        return Sequence2<T>(f) /= value;
}

template <class T>
bool operator==(const Sequence2<T>& f, const Sequence2<T>& g)
{
        bool result;
        if (!f.getSize() && !g.getSize()) {
                // Both sequences are empty.
                result = true;
        } else if (f.getStartX() != g.getStartX() ||
                f.getEndX() != g.getEndX() || f.getStartY() != g.getStartY() ||
                  f.getEndY() != g.getEndY()) {
                // The sequences are not defined on the same domain.
                result = false;
        } else {
                // The sequences are defined on the same domain.
                result = std::equal(f.begin(), f.end(), g.begin());
        }
        return result;
}

template <class T>
SPL_SEQUENCE2_INLINE bool operator!=(const Sequence2<T>& f,
  const Sequence2<T>& g)
{
        return !(f == g);
}

template <class T>
SPL_SEQUENCE2_INLINE bool approxEqual(const Sequence2<T>& f,
  const Sequence2<T>& g, T threshold = 1e-9)
{
        bool result;
        if (!f.getSize() && !g.getSize()) {
                // Empty sequence.
                result = true;
        } else if (f.getStartX() != g.getStartX() ||
          f.getStartY() != g.getStartY() ||
          f.getEndX() != g.getEndX() ||
          f.getEndY() != g.getEndY()) {
                // Domains do not match.
                result = false;
        } else {
                assert(f.getWidth() == g.getWidth() && f.getHeight() == g.getHeight());
                result = true;
                typename Sequence2<T>::ConstIterator src = f.begin();
                typename Sequence2<T>::ConstIterator dst = g.begin();
                while (src != f.end()) {
                        if (absVal((*src) - (*dst)) > threshold) {
                                result = false;
                                break;
                        }
                        ++src;
                        ++dst;
                }
        }
        return result;
}

//

template <class T>
Sequence2<T> subsequence(const Sequence2<T>& f, int startX, int startY,
  int width, int height)
{
        assert(startX >= f.getStartX() && startX + width <= f.getEndX() &&
          startY >= f.getStartY() && startY + height <= f.getEndY());
        if (!width || !height) {
                width = 0;
                height = 0;
        }
        Sequence2<T> result(startX, startY, width, height);
        for (int i = 0; i < height; ++i) {
                assert(startX - f.getStartX() >= 0);
                typename Sequence2<T>::ConstXIterator src =
                  f.rowBegin(startY + i) + (startX - f.getStartX());
                typename Sequence2<T>::XIterator dst = result.rowBegin(startY + i);
                SPL::copy_n(src, width, dst);
        }
        return result;
}

template <class T>
SPL_SEQUENCE2_INLINE Sequence2<T> translate(const Sequence2<T>& f,
  int deltaX, int deltaY)
{
        return Sequence2<T>(f).translate(deltaX, deltaY);
}

// Convolution


template<class T>
SPL_SEQUENCE2_INLINE T getExtSeqValue(const Sequence2<T>& f, int x, int y,
  int mode)
{
        assert(mode == ConvolveMode::full ||
          mode == ConvolveMode::sameDomainConstExt ||
          mode == ConvolveMode::sameDomainZeroExt ||
          mode == ConvolveMode::sameDomainPerExt ||
          mode == ConvolveMode::sameDomainSymExt0);
        T value;
        switch (mode) {
        case ConvolveMode::full:
        case ConvolveMode::sameDomainZeroExt:
                if (x >= f.getStartX() && x < f.getEndX() &&
                  y >= f.getStartY() && y < f.getEndY()) {
                        value = f(x, y);
                } else {
                        value = T(0);
                }
                break;
        case ConvolveMode::sameDomainConstExt:
                if (x < f.getStartX()) {
                        x = f.getStartX();
                } else if (x >= f.getEndX()) {
                        x = f.getEndX() - 1;
                }
                if (y < f.getStartY()) {
                        y = f.getStartY();
                } else if (y >= f.getEndY()) {
                        y = f.getEndY() - 1;
                }
                value = f(x, y);
                break;
        case ConvolveMode::sameDomainPerExt:
                if (x < f.getStartX() || x >= f.getEndX()) {
                        x = f.getStartX() + mod(x - f.getStartX(), f.getWidth());
                }
                if (y < f.getStartY() || y >= f.getEndY()) {
                        y = f.getStartY() + mod(y - f.getStartY(), f.getHeight());
                }
                assert(x >= f.getStartX() && x < f.getEndX());
                assert(y >= f.getStartY() && y < f.getEndY());
                value = f(x, y);
                break;
        case ConvolveMode::sameDomainSymExt0:
                if (x < f.getStartX() || x >= f.getEndX()) {
                        int tmp = mod(x - f.getStartX(), 2 * f.getWidth() - 2);
                        x = f.getStartX() + std::min(tmp, 2 * f.getWidth() - 2 - tmp);
                }
                if (y < f.getStartY() || y >= f.getEndY()) {
                        int tmp = mod(y - f.getStartY(), 2 * f.getHeight() - 2);
                        y = f.getStartY() + std::min(tmp, 2 * f.getHeight() - 2 - tmp);
                }
                assert(x >= f.getStartX() && x < f.getEndX());
                assert(y >= f.getStartY() && y < f.getEndY());
                value = f(x, y);
                break;
        default:
                assert(0);
                abort();
                break;
        }
        return value;
}


template <class T>
Sequence2<T> convolveSeparableFull(const Sequence2<T>& f,
  const Sequence1<T>& horzFilt, const Sequence1<T>& vertFilt)
{
        if (!horzFilt.getSize() || !vertFilt.getSize()) {
                return Sequence2<T>();
        }

        // Perform the horizontal filtering.
        Sequence2<T> tmpResult(f.getStartX() + horzFilt.getStartInd(),
          f.getStartY(), f.getWidth() + horzFilt.getSize() - 1,
          f.getHeight());
        for (int y = f.getStartY(); y < f.getEndY(); ++y) {
                convolveHelper(f.rowBegin(y), f.rowEnd(y), horzFilt.begin(),
                  horzFilt.end(), tmpResult.rowBegin(y), T());
        }

        // Perform the vertical filtering.
        Sequence2<T> result(tmpResult.getStartX(), tmpResult.getStartY() +
          vertFilt.getStartInd(), tmpResult.getWidth(), tmpResult.getHeight() +
          vertFilt.getSize() - 1);
        for (int x = tmpResult.getStartX(); x < tmpResult.getEndX(); ++x) {
                convolveHelper(tmpResult.colBegin(x), tmpResult.colEnd(x),
                  vertFilt.begin(), vertFilt.end(), result.colBegin(x), T());
        }

        return result;
}


template <class T>
Sequence2<T> convolveSeparableSameDomainZeroExt(const Sequence2<T>& f,
  const Sequence1<T>& horzFilt, const Sequence1<T>& vertFilt)
{
        if (!horzFilt.getSize() || !vertFilt.getSize()) {
                return Sequence2<T>();
        }

        int skip;
        int shift;

        // Perform the horizontal filtering.
        Sequence2<T> tmpResult(f.getStartX(), f.getStartY(), f.getWidth(),
          f.getHeight());
        if (horzFilt.getStartInd() >= 0) {
                shift = horzFilt.getStartInd();
                skip = 0;
        } else {
                shift = 0;
                skip = -horzFilt.getStartInd();
        }
        for (int y = f.getStartY(); y < f.getEndY(); ++y) {
                convolveHelper2(f.rowBegin(y), f.rowEnd(y), horzFilt.begin(),
                  horzFilt.end(), tmpResult.rowBegin(y) + shift, skip,
                  tmpResult.getWidth() - shift, T());
        }

        // Perform the vertical filtering.
        if (vertFilt.getStartInd() >= 0) {
                shift = vertFilt.getStartInd();
                skip = 0;
        } else {
                shift = 0;
                skip = -vertFilt.getStartInd();
        }
        Sequence2<T> result(tmpResult.getStartX(), tmpResult.getStartY(),
          tmpResult.getWidth(), tmpResult.getHeight());
        for (int x = tmpResult.getStartX(); x < tmpResult.getEndX(); ++x) {
                convolveHelper2(tmpResult.colBegin(x), tmpResult.colEnd(x),
                  vertFilt.begin(), vertFilt.end(), result.colBegin(x) + shift,
                  skip, result.getHeight() - shift, T());
        }

        return result;
}

template<class T>
Sequence2<T> convolve(const Sequence2<T>& f, const Sequence2<T>& g, int mode)
{
        assert(mode == ConvolveMode::full ||
          mode == ConvolveMode::sameDomainConstExt ||
          mode == ConvolveMode::sameDomainZeroExt ||
          mode == ConvolveMode::sameDomainPerExt ||
          mode == ConvolveMode::sameDomainSymExt0);
        assert(f.getWidth() > 0 && f.getHeight() > 0);
        assert(g.getWidth() > 0 && g.getHeight() > 0);

        Sequence2<T> result;

        switch (mode) {
        case ConvolveMode::full:
                result = Sequence2<T>(f.getStartX() + g.getStartX(),
                  f.getStartY() + g.getStartY(), f.getWidth() + g.getWidth() - 1,
                  f.getHeight() + g.getHeight() - 1);
                break;
        case ConvolveMode::sameDomainConstExt:
        case ConvolveMode::sameDomainZeroExt:
        case ConvolveMode::sameDomainPerExt:
        case ConvolveMode::sameDomainSymExt0:
                result = Sequence2<T>(f.getStartX(), f.getStartY(), f.getWidth(),
                  f.getHeight());
                break;
        default:
                assert(0);
                abort();
                break;
        }

        for (int y = result.getStartY(); y < result.getEndY(); ++y) {
                for (int x = result.getStartX(); x < result.getEndX(); ++x) {
                        T sum = T(0);
                        for (int j = g.getStartY(); j < g.getEndY(); ++j) {
                                for (int i = g.getStartX(); i < g.getEndX(); ++i) {
                                        sum += getExtSeqValue(f, x - i, y - j, mode) * g(i, j);
                                }
                        }
                        result(x, y) = sum;
                }
        }

        return result;
}

template <class T>
Sequence2<T> convolveSeparable(const Sequence2<T>& f,
  const Sequence1<T>& horzFilt, const Sequence1<T>& vertFilt,
  int mode = ConvolveMode::full)
{
#if defined(SPL_SEQUENCE2_USE_NEW_CONV)
        switch (mode) {
        case ConvolveMode::full:
                return convolveSeparableFull(f, horzFilt, vertFilt);
                break;
        case ConvolveMode::sameDomainZeroExt:
                return convolveSeparableSameDomainZeroExt(f, horzFilt, vertFilt);
                break;
        }
#endif
        // The following code is quite inefficient.  Rewrite it someday.
        Sequence2<T> horzFilt0(horzFilt.getStartInd(), 0, horzFilt.getSize(), 1,
          horzFilt.begin());
        Sequence2<T> vertFilt0(0, vertFilt.getStartInd(), 1, vertFilt.getSize(),
          vertFilt.begin());
        Sequence2<T> result = convolve(f, horzFilt0, mode);
        result = convolve(result, vertFilt0, mode);
        return result;
}

// Downsampling and upsampling

template <class T>
Sequence2<T> downsample(const Sequence2<T>& f, int factorX, int factorY)
{
        assert(factorX >= 1 && factorY >= 1);
        int startX = ceilDiv(f.getStartX(), factorX) * factorX;
        int startY = ceilDiv(f.getStartY(), factorY) * factorY;
        int width = ceilDiv(f.getEndX() - startX, factorX);
        int height = ceilDiv(f.getEndY() - startY, factorY);
        Sequence2<T> result(roundTowardZeroDiv(startX, factorX),
          roundTowardZeroDiv(startY, factorY), width, height);

        for (int i = 0; i < height; ++i) {
                typename Sequence2<T>::ConstXIterator srcIter =
                  f.rowBegin(startY + factorY * i) + (startX - f.getStartX());
                typename Sequence2<T>::XIterator dstIter =
                  result.rowBegin(result.getStartY() + i);
                int n = result.getWidth();
                while (--n >= 0) {
                        *dstIter = *srcIter;
                        ++dstIter;
                        srcIter += factorX;
                }
        }

        return result;
}

template <class T>
Sequence2<T> upsample(const Sequence2<T>& f, int factorX, int factorY)
{
        return upsample(f, factorX, factorY, 0, 0);
}

template <class T>
Sequence2<T> upsample(const Sequence2<T>& f, int factorX, int factorY,
  int padX, int padY)
{
        assert(factorX >= 1 && factorY >= 1);
        assert(padX >= 0 && padX < factorX);
        assert(padY >= 0 && padY < factorY);

        Sequence2<T> result(factorX * f.getStartX(), factorY * f.getStartY(),
          f.getWidth() + std::max((f.getWidth() - 1), 0) * (factorX - 1) + padX,
          f.getHeight() + std::max((f.getHeight() - 1), 0) * (factorY - 1) + padY);
        int srcY = f.getStartY();
        int dstY = result.getStartY();
        int rowCount = f.getHeight();

        while (rowCount > 0) {
                typename Sequence2<T>::ConstXIterator srcIter = f.rowBegin(srcY);
                typename Sequence2<T>::XIterator dstIter = result.rowBegin(dstY);
                int count = f.getWidth();

                while (count > 0) {
                        *dstIter = *srcIter;
                        ++srcIter;
                        ++dstIter;
                        for (int k = (count > 1) * (factorX - 1); k > 0; --k) {
                                *dstIter = T(0);
                                ++dstIter;
                        }
                        --count;
                }
                for (int k = padX; k > 0; --k) {
                        *dstIter = T(0);
                        ++dstIter;
                }
                ++srcY;
                ++dstY;

                for (int k = (rowCount > 1) * (factorY - 1); k > 0; --k) {
                        std::fill_n(result.rowBegin(dstY), result.getWidth(), T(0));
                        ++dstY;
                }
                --rowCount;
        }

        for (int k = padY; k > 0; --k) {
                assert(dstY < result.getEndY());
                std::fill_n(result.rowBegin(dstY), result.getWidth(), T(0));
                ++dstY;
        }

        return result;
}

// Polyphase decompositions

template <class T>
Array2<Sequence2<T> > polyphaseSplit(const Sequence2<T>& seq, int typeX,
  int numPhasesX, int typeY, int numPhasesY)
{
        assert(numPhasesX >= 2);
        assert(numPhasesY >= 2);
        Array2<Sequence2<T> > result(numPhasesX, numPhasesY);
        for (int j = 0; j < numPhasesY; ++j) {
                for (int i = 0; i < numPhasesX; ++i) {
                        int offsetX = getCosetOffset(typeX, numPhasesX, i);
                        int offsetY = getCosetOffset(typeY, numPhasesY, j);
                        result(i, j) = downsample(translate(seq, offsetX, offsetY),
                          numPhasesX, numPhasesY);
                }
        }
        return result;
}

template <class T>
Sequence2<T> polyphaseJoin(const Array2<Sequence2<T> >& comps, int typeX,
  int typeY)
{
        const int numPhasesX = comps.getWidth();
        const int numPhasesY = comps.getHeight();
        assert(numPhasesX >= 2);
        assert(numPhasesY >= 2);
        Sequence2<T> seq;
        for (int j = 0; j < numPhasesY; ++j) {
                for (int i = 0; i < numPhasesX; ++i) {
                        int offsetX = getCosetOffset(typeX, numPhasesX, i);
                        int offsetY = getCosetOffset(typeY, numPhasesY, j);
//std::cerr << "XXX " << i << " " << j << " " << upsample(comps(i, j), numPhasesX, numPhasesY) << "\n";
                        seq = add(seq, translate(upsample(comps(i, j), numPhasesX,
                          numPhasesY), -offsetX, -offsetY));
                }
        }
        return seq;
}

// Types.

typedef Sequence2<double> RealSequence2;

typedef Sequence2<int> IntSequence2;

//

}

#endif