/************************************************************************/
/*                                                                      */
/*               Copyright 1998-2004 by Ullrich Koethe                  */
/*       Cognitive Systems Group, University of Hamburg, Germany        */
/*                                                                      */
/*    This file is part of the VIGRA computer vision library.           */
/*    The VIGRA Website is                                              */
/*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
/*    Please direct questions, bug reports, and contributions to        */
/*        koethe@informatik.uni-hamburg.de          or                  */
/*        vigra@kogs1.informatik.uni-hamburg.de                         */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/

#ifndef VIGRA_RESAMPLING_CONVOLUTION_HXX
#define VIGRA_RESAMPLING_CONVOLUTION_HXX

#include <cmath>
#include "vigra/stdimage.hxx"
#include "vigra/array_vector.hxx"
#include "vigra/rational.hxx"
#include "vigra/functortraits.hxx"

namespace vigra {

namespace resampling_detail
{

struct MapTargetToSourceCoordinate
{
    MapTargetToSourceCoordinate(Rational<int> const & samplingRatio,
                                Rational<int> const & offset)
    : a(samplingRatio.denominator()*offset.denominator()),
      b(samplingRatio.numerator()*offset.numerator()),
      c(samplingRatio.numerator()*offset.denominator())
    {}

//        the following funcions are more efficient realizations of:
//             rational_cast<T>(i / samplingRatio + offset);
//        we need efficiency because this may be called in the inner loop

    int operator()(int i) const
    {
        return (i * a + b) / c;
    }

    double toDouble(int i) const
    {
        return double(i * a + b) / c;
    }

    Rational<int> toRational(int i) const
    {
        return Rational<int>(i * a + b, c);
    }

    int a, b, c;
};

} // namespace resampling_detail

template <>
class FunctorTraits<resampling_detail::MapTargetToSourceCoordinate>
: public FunctorTraitsBase<resampling_detail::MapTargetToSourceCoordinate>
{
  public:
    typedef VigraTrueType isUnaryFunctor;
};

template <class SrcIter, class SrcAcc,
          class DestIter, class DestAcc,
          class KernelArray,
          class Functor>
void
resamplingConvolveLine(SrcIter s, SrcIter send, SrcAcc src,
                       DestIter d, DestIter dend, DestAcc dest,
                       KernelArray const & kernels,
                       Functor mapTargetToSourceCoordinate)
{
    typedef typename
        NumericTraits<typename SrcAcc::value_type>::RealPromote
        TmpType;
    typedef typename KernelArray::value_type Kernel;
    typedef typename Kernel::const_iterator KernelIter;

    int wo = send - s;
    int wn = dend - d;
    int wo2 = 2*wo - 2;

    int i;
    typename KernelArray::const_iterator kernel = kernels.begin();
    for(i=0; i<wn; ++i, ++d, ++kernel)
    {
        // use the kernels periodically
        if(kernel == kernels.end())
            kernel = kernels.begin();

        // calculate current target point into source location
        int is = mapTargetToSourceCoordinate(i);

        TmpType sum = NumericTraits<TmpType>::zero();

        int lbound = is - kernel->right(),
            hbound = is - kernel->left();

        KernelIter k = kernel->center() + kernel->right();
        if(lbound < 0 || hbound >= wo)
        {
            vigra_precondition(-lbound < wo && wo2 - hbound >= 0,
                "resamplingConvolveLine(): kernel or offset larger than image.");
            for(int m=lbound; m <= hbound; ++m, --k)
            {
                int mm = (m < 0) ?
                            -m :
                            (m >= wo) ?
                                wo2 - m :
                                m;
                sum += *k * src(s, mm);
            }
        }
        else
        {
            SrcIter ss = s + lbound;
            SrcIter ssend = s + hbound;

            for(; ss <= ssend; ++ss, --k)
            {
                sum += *k * src(ss);
            }
        }

        dest.set(sum, d);
    }
}

template <class Kernel, class MapCoordinate, class KernelArray>
void
createResamplingKernels(Kernel const & kernel,
             MapCoordinate const & mapCoordinate, KernelArray & kernels)
{
    for(unsigned int idest = 0; idest < kernels.size(); ++idest)
    {
        int isrc = mapCoordinate(idest);
        double idsrc = mapCoordinate.toDouble(idest);
        double offset = idsrc - isrc;
        double radius = kernel.radius();
        int left = int(ceil(-radius - offset));
        int right = int(floor(radius - offset));
        kernels[idest].initExplicitly(left, right);

        double x = left + offset;
        for(int i = left; i <= right; ++i, ++x)
            kernels[idest][i] = kernel(x);
        kernels[idest].normalize(1.0, kernel.derivativeOrder(), offset);
    }
}

/** \addtogroup ResamplingConvolutionFilters Resampling Convolution Filters

    These functions implement the convolution operation when the source and target images
    have different sizes. This is realized by accessing a continous kernel at the
    appropriate non-integer positions. The technique is, for example, described in
    D. Schumacher: <i>General Filtered Image Rescaling</i>, in: Graphics Gems III,
    Academic Press, 1992.
*/
//@{

/********************************************************/
/*                                                      */
/*                  resamplingConvolveX                 */
/*                                                      */
/********************************************************/

/** \brief Apply a resampling filter in the x-direction.

    This function implements a convolution operation in x-direction
    (i.e. applies a 1D filter to every row) where the width of the source
    and destination images differ. This is typically used to avoid aliasing if
    the image is scaled down, or to interpolate smoothly if the image is scaled up.
    The target coordinates are transformed into source coordinates by

    \code
    xsource = (xtarget - offset) / samplingRatio
    \endcode

    The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
    in order to avoid rounding errors in this transformation. It is required that for all
    pixels of the target image, <tt>xsource</tt> remains within the range of the source
    image (i.e. <tt>0 <= xsource <= sourceWidth-1</tt>. Since <tt>xsource</tt> is
    in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
    arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
    which specifies the support (non-zero interval) of the kernel. VIGRA already
    provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
    \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
    \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and
    resamplingConvolveY().

    <b> Declarations:</b>

    pass arguments explicitly:
    \code
    namespace vigra {
        template <class SrcIter, class SrcAcc,
                  class DestIter, class DestAcc,
                  class Kernel>
        void
        resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
                            DestIter dul, DestIter dlr, DestAcc dest,
                            Kernel const & kernel,
                            Rational<int> const & samplingRatio, Rational<int> const & offset);
    }
    \endcode


    use argument objects in conjunction with \ref ArgumentObjectFactories:
    \code
    namespace vigra {
        template <class SrcIter, class SrcAcc,
                  class DestIter, class DestAcc,
                  class Kernel>
        void
        resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
                            triple<DestIter, DestIter, DestAcc> dest,
                            Kernel const & kernel,
                            Rational<int> const & samplingRatio, Rational<int> const & offset);
    }
    \endcode

    <b> Usage:</b>

    <b>\#include</b> "<a href="resampling__convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"


    \code
    Rational<int> ratio(2), offset(0);

    FImage src(w,h),
           dest(rational_cast<int>(ratio*w), h);

    float sigma = 2.0;
    Gaussian<float> smooth(sigma);
    ...

    // simpultaneously enlarge and smooth source image
    resamplingConvolveX(srcImageRange(src), destImageRange(dest),
                        smooth, ratio, offset);
    \endcode

    <b> Required Interface:</b>

    \code
    Kernel kernel;
    int kernelRadius = kernel.radius();
    double x = ...;  // must be <= radius()
    double value = kernel(x);
    \endcode
*/
template <class SrcIter, class SrcAcc,
          class DestIter, class DestAcc,
          class Kernel>
void
resamplingConvolveX(SrcIter sul, SrcIter slr, SrcAcc src,
                    DestIter dul, DestIter dlr, DestAcc dest,
                    Kernel const & kernel,
                    Rational<int> const & samplingRatio, Rational<int> const & offset)
{
    int wold = slr.x - sul.x;
    int wnew = dlr.x - dul.x;

    vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
                "resamplingConvolveX(): sampling ratio must be > 0 and < infinity");
    vigra_precondition(!offset.is_inf(),
                "resamplingConvolveX(): offset must be < infinity");

    int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());
    resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);

    ArrayVector<Kernel1D<double> > kernels(period);

    createResamplingKernels(kernel, mapCoordinate, kernels);

    for(; sul.y < slr.y; ++sul.y, ++dul.y)
    {
        typename SrcIter::row_iterator sr = sul.rowIterator();
        typename DestIter::row_iterator dr = dul.rowIterator();
        resamplingConvolveLine(sr, sr+wold, src, dr, dr+wnew, dest,
                               kernels, mapCoordinate);
    }
}

template <class SrcIter, class SrcAcc,
          class DestIter, class DestAcc,
          class Kernel>
inline void
resamplingConvolveX(triple<SrcIter, SrcIter, SrcAcc> src,
                    triple<DestIter, DestIter, DestAcc> dest,
                    Kernel const & kernel,
                    Rational<int> const & samplingRatio, Rational<int> const & offset)
{
    resamplingConvolveX(src.first, src.second, src.third,
                        dest.first, dest.second, dest.third,
                        kernel, samplingRatio, offset);
}

/********************************************************/
/*                                                      */
/*                  resamplingConvolveY                 */
/*                                                      */
/********************************************************/

/** \brief Apply a resampling filter in the y-direction.

    This function implements a convolution operation in y-direction
    (i.e. applies a 1D filter to every column) where the height of the source
    and destination images differ. This is typically used to avoid aliasing if
    the image is scaled down, or to interpolate smoothly if the image is scaled up.
    The target coordinates are transformed into source coordinates by

    \code
    ysource = (ytarget - offset) / samplingRatio
    \endcode

    The <tt>samplingRatio</tt> and <tt>offset</tt> must be given as \ref vigra::Rational
    in order to avoid rounding errors in this transformation. It is required that for all
    pixels of the target image, <tt>ysource</tt> remains within the range of the source
    image (i.e. <tt>0 <= ysource <= sourceHeight-1</tt>. Since <tt>ysource</tt> is
    in general not an integer, the <tt>kernel</tt> must be a functor that can be accessed at
    arbitrary (<tt>double</tt>) coordinates. It must also provide a member function <tt>radius()</tt>
    which specifies the support (non-zero interval) of the kernel. VIGRA already
    provides a number of suitable functors, e.g. \ref vigra::Gaussian, \ref vigra::BSpline
    \ref vigra::CatmullRomSpline, and \ref vigra::CoscotFunction. The function
    \ref resizeImageSplineInterpolation() is implemented by means resamplingConvolveX() and
    resamplingConvolveY().

    <b> Declarations:</b>

    pass arguments explicitly:
    \code
    namespace vigra {
        template <class SrcIter, class SrcAcc,
                  class DestIter, class DestAcc,
                  class Kernel>
        void
        resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
                            DestIter dul, DestIter dlr, DestAcc dest,
                            Kernel const & kernel,
                            Rational<int> const & samplingRatio, Rational<int> const & offset);
    }
    \endcode


    use argument objects in conjunction with \ref ArgumentObjectFactories:
    \code
    namespace vigra {
        template <class SrcIter, class SrcAcc,
                  class DestIter, class DestAcc,
                  class Kernel>
        void
        resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
                            triple<DestIter, DestIter, DestAcc> dest,
                            Kernel const & kernel,
                            Rational<int> const & samplingRatio, Rational<int> const & offset);
    }
    \endcode

    <b> Usage:</b>

    <b>\#include</b> "<a href="resampling__convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"


    \code
    Rational<int> ratio(2), offset(0);

    FImage src(w,h),
           dest(w, rational_cast<int>(ratio*h));

    float sigma = 2.0;
    Gaussian<float> smooth(sigma);
    ...

    // simpultaneously enlarge and smooth source image
    resamplingConvolveY(srcImageRange(src), destImageRange(dest),
                        smooth, ratio, offset);
    \endcode

    <b> Required Interface:</b>

    \code
    Kernel kernel;
    int kernelRadius = kernel.radius();
    double y = ...;  // must be <= radius()
    double value = kernel(y);
    \endcode
*/
template <class SrcIter, class SrcAcc,
          class DestIter, class DestAcc,
          class Kernel>
void
resamplingConvolveY(SrcIter sul, SrcIter slr, SrcAcc src,
                    DestIter dul, DestIter dlr, DestAcc dest,
                    Kernel const & kernel,
                    Rational<int> const & samplingRatio, Rational<int> const & offset)
{
    int hold = slr.y - sul.y;
    int hnew = dlr.y - dul.y;

    vigra_precondition(!samplingRatio.is_inf() && samplingRatio > 0,
                "resamplingConvolveY(): sampling ratio must be > 0 and < infinity");
    vigra_precondition(!offset.is_inf(),
                "resamplingConvolveY(): offset must be < infinity");

    int period = lcm(samplingRatio.numerator(), samplingRatio.denominator());

    resampling_detail::MapTargetToSourceCoordinate mapCoordinate(samplingRatio, offset);

    ArrayVector<Kernel1D<double> > kernels(period);

    createResamplingKernels(kernel, mapCoordinate, kernels);

    for(; sul.x < slr.x; ++sul.x, ++dul.x)
    {
        typename SrcIter::column_iterator sc = sul.columnIterator();
        typename DestIter::column_iterator dc = dul.columnIterator();
        resamplingConvolveLine(sc, sc+hold, src, dc, dc+hnew, dest,
                               kernels, mapCoordinate);
    }
}

template <class SrcIter, class SrcAcc,
          class DestIter, class DestAcc,
          class Kernel>
inline void
resamplingConvolveY(triple<SrcIter, SrcIter, SrcAcc> src,
                    triple<DestIter, DestIter, DestAcc> dest,
                    Kernel const & kernel,
                    Rational<int> const & samplingRatio, Rational<int> const & offset)
{
    resamplingConvolveY(src.first, src.second, src.third,
                        dest.first, dest.second, dest.third,
                        kernel, samplingRatio, offset);
}

/********************************************************/
/*                                                      */
/*               resamplingConvolveImage                */
/*                                                      */
/********************************************************/

/** \brief Apply two separable resampling filters successively, the first in x-direction,
           the second in y-direction.

    This function is a shorthand for the concatenation of a call to
    \link ResamplingConvolutionFilters#resamplingConvolveX resamplingConvolveX\endlink()
    and \link ResamplingConvolutionFilters#resamplingConvolveY resamplingConvolveY\endlink()
    with the given kernels. See there for detailed documentation.

    <b> Declarations:</b>

    pass arguments explicitly:
    \code
    namespace vigra {
        template <class SrcIterator, class SrcAccessor,
                  class DestIterator, class DestAccessor,
                  class KernelX, class KernelY>
        void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
                           DestIterator dul, DestIterator dlr, DestAccessor dest,
                           KernelX const & kx,
                           Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
                           KernelY const & ky,
                           Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
    }
    \endcode


    use argument objects in conjunction with \ref ArgumentObjectFactories:
    \code
    namespace vigra {
        template <class SrcIterator, class SrcAccessor,
                  class DestIterator, class DestAccessor,
                  class KernelX, class KernelY>
        void
        resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
                           triple<DestIterator, DestIterator, DestAccessor> dest,
                           KernelX const & kx,
                           Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
                           KernelY const & ky,
                           Rational<int> const & samplingRatioY, Rational<int> const & offsetY);
    }
    \endcode

    <b> Usage:</b>

    <b>\#include</b> "<a href="resampling__convolution_8hxx-source.html">vigra/resampling_convolution.hxx</a>"


    \code
    Rational<int> xratio(2), yratio(3), offset(0);

    FImage src(w,h),
           dest(rational_cast<int>(xratio*w), rational_cast<int>(yratio*h));

    float sigma = 2.0;
    Gaussian<float> smooth(sigma);
    ...

    // simpultaneously enlarge and smooth source image
    resamplingConvolveImage(srcImageRange(src), destImageRange(dest),
                            smooth, xratio, offset,
                            smooth, yratio, offset);

    \endcode

*/
template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class KernelX, class KernelY>
void resamplingConvolveImage(SrcIterator sul,SrcIterator slr, SrcAccessor src,
                   DestIterator dul, DestIterator dlr, DestAccessor dest,
                   KernelX const & kx,
                   Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
                   KernelY const & ky,
                   Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
{
    typedef typename
        NumericTraits<typename SrcAccessor::value_type>::RealPromote
        TmpType;

    BasicImage<TmpType> tmp(dlr.x - dul.x, slr.y - sul.y);

    resamplingConvolveX(srcIterRange(sul, slr, src),
                        destImageRange(tmp),
                        kx, samplingRatioX, offsetX);
    resamplingConvolveY(srcImageRange(tmp),
                        destIterRange(dul, dlr, dest),
                        ky, samplingRatioY, offsetY);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class KernelX, class KernelY>
inline void
resamplingConvolveImage(triple<SrcIterator, SrcIterator, SrcAccessor> src,
                   triple<DestIterator, DestIterator, DestAccessor> dest,
                   KernelX const & kx,
                   Rational<int> const & samplingRatioX, Rational<int> const & offsetX,
                   KernelY const & ky,
                   Rational<int> const & samplingRatioY, Rational<int> const & offsetY)
{
    resamplingConvolveImage(src.first, src.second, src.third,
                            dest.first, dest.second, dest.third,
                            kx, samplingRatioX, offsetX,
                            ky, samplingRatioY, offsetY);
}

} // namespace vigra


#endif /* VIGRA_RESAMPLING_CONVOLUTION_HXX */