GteIntpBicubic2.h

Go to the documentation of this file.

// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2017
// Distributed under the Boost Software License, Version 1.0.
// http://www.boost.org/LICENSE_1_0.txt
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// File Version: 3.0.0 (2016/06/19)

#pragma once

#include <LowLevel/GteLogger.h>

// The interpolator is for uniformly spaced (x,y)-values.  The input samples
// F must be stored in row-major order to represent f(x,y); that is,
// F[c + xBound*r] corresponds to f(x,y), where c is the index corresponding
// to x and r is the index corresponding to y.  Exact interpolation is
// achieved by setting catmullRom to 'true', giving you the Catmull-Rom
// blending matrix.  If a smooth interpolation is desired, set catmullRom to
// 'false' to obtain B-spline blending.

namespace gte
{

template <typename Real>
class IntpBicubic2
{
public:
    // Construction.
    IntpBicubic2(int xBound, int yBound, Real xMin, Real xSpacing,
        Real yMin, Real ySpacing, Real const* F, bool catmullRom);

    // Member access.
    inline int GetXBound() const;
    inline int GetYBound() const;
    inline int GetQuantity() const;
    inline Real const* GetF() const;
    inline Real GetXMin() const;
    inline Real GetXMax() const;
    inline Real GetXSpacing() const;
    inline Real GetYMin() const;
    inline Real GetYMax() const;
    inline Real GetYSpacing() const;

    // Evaluate the function and its derivatives.  The functions clamp the
    // inputs to xmin <= x <= xmax and ymin <= y <= ymax.  The first operator
    // is for function evaluation.  The second operator is for function or
    // derivative evaluations.  The xOrder argument is the order of the
    // x-derivative and the yOrder argument is the order of the y-derivative.
    // Both orders are zero to get the function value itself.
    Real operator()(Real x, Real y) const;
    Real operator()(int xOrder, int yOrder, Real x, Real y) const;

private:
    int mXBound, mYBound, mQuantity;
    Real mXMin, mXMax, mXSpacing, mInvXSpacing;
    Real mYMin, mYMax, mYSpacing, mInvYSpacing;
    Real const* mF;
    Real mBlend[4][4];
};


template <typename Real>
IntpBicubic2<Real>::IntpBicubic2(int xBound, int yBound, Real xMin,
    Real xSpacing, Real yMin, Real ySpacing, Real const* F, bool catmullRom)
    :
    mXBound(xBound),
    mYBound(yBound),
    mQuantity(xBound * yBound),
    mXMin(xMin),
    mXSpacing(xSpacing),
    mYMin(yMin),
    mYSpacing(ySpacing),
    mF(F)
{
    // At least a 3x3 block of data points are needed to construct the
    // estimates of the boundary derivatives.
    LogAssert(mXBound >= 3 && mYBound >= 3 && mF, "Invalid input.");
    LogAssert(mXSpacing > (Real)0 && mYSpacing > (Real)0, "Invalid input.");

    mXMax = mXMin + mXSpacing * static_cast<Real>(mXBound - 1);
    mInvXSpacing = ((Real)1) / mXSpacing;
    mYMax = mYMin + mYSpacing * static_cast<Real>(mYBound - 1);
    mInvYSpacing = ((Real)1) / mYSpacing;

    if (catmullRom)
    {
        mBlend[0][0] = (Real)0;
        mBlend[0][1] = (Real)-0.5;
        mBlend[0][2] = (Real)1;
        mBlend[0][3] = (Real)-0.5;
        mBlend[1][0] = (Real)1;
        mBlend[1][1] = (Real)0;
        mBlend[1][2] = (Real)-2.5;
        mBlend[1][3] = (Real)1.5;
        mBlend[2][0] = (Real)0;
        mBlend[2][1] = (Real)0.5;
        mBlend[2][2] = (Real)2;
        mBlend[2][3] = (Real)-1.5;
        mBlend[3][0] = (Real)0;
        mBlend[3][1] = (Real)0;
        mBlend[3][2] = (Real)-0.5;
        mBlend[3][3] = (Real)0.5;
    }
    else
    {
        mBlend[0][0] = (Real)1 / (Real)6;
        mBlend[0][1] = (Real)-3 / (Real)6;
        mBlend[0][2] = (Real)3 / (Real)6;
        mBlend[0][3] = (Real)-1 / (Real)6;;
        mBlend[1][0] = (Real)4 / (Real)6;
        mBlend[1][1] = (Real)0 / (Real)6;
        mBlend[1][2] = (Real)-6 / (Real)6;
        mBlend[1][3] = (Real)3 / (Real)6;
        mBlend[2][0] = (Real)1 / (Real)6;
        mBlend[2][1] = (Real)3 / (Real)6;
        mBlend[2][2] = (Real)3 / (Real)6;
        mBlend[2][3] = (Real)-3 / (Real)6;
        mBlend[3][0] = (Real)0 / (Real)6;
        mBlend[3][1] = (Real)0 / (Real)6;
        mBlend[3][2] = (Real)0 / (Real)6;
        mBlend[3][3] = (Real)1 / (Real)6;
    }
}

template <typename Real> inline
int IntpBicubic2<Real>::GetXBound() const
{
    return mXBound;
}

template <typename Real> inline
int IntpBicubic2<Real>::GetYBound() const
{
    return mYBound;
}

template <typename Real> inline
int IntpBicubic2<Real>::GetQuantity() const
{
    return mQuantity;
}

template <typename Real> inline
Real const* IntpBicubic2<Real>::GetF() const
{
    return mF;
}

template <typename Real> inline
Real IntpBicubic2<Real>::GetXMin() const
{
    return mXMin;
}

template <typename Real> inline
Real IntpBicubic2<Real>::GetXMax() const
{
    return mXMax;
}

template <typename Real> inline
Real IntpBicubic2<Real>::GetXSpacing() const
{
    return mXSpacing;
}

template <typename Real> inline
Real IntpBicubic2<Real>::GetYMin() const
{
    return mYMin;
}

template <typename Real> inline
Real IntpBicubic2<Real>::GetYMax() const
{
    return mYMax;
}

template <typename Real> inline
Real IntpBicubic2<Real>::GetYSpacing() const
{
    return mYSpacing;
}

template <typename Real>
Real IntpBicubic2<Real>::operator()(Real x, Real y) const
{
    // Compute x-index and clamp to image.
    Real xIndex = (x - mXMin) * mInvXSpacing;
    int ix = static_cast<int>(xIndex);
    if (ix < 0)
    {
        ix = 0;
    }
    else if (ix >= mXBound)
    {
        ix = mXBound - 1;
    }

    // Compute y-index and clamp to image.
    Real yIndex = (y - mYMin) * mInvYSpacing;
    int iy = static_cast<int>(yIndex);
    if (iy < 0)
    {
        iy = 0;
    }
    else if (iy >= mYBound)
    {
        iy = mYBound - 1;
    }

    Real U[4];
    U[0] = (Real)1;
    U[1] = xIndex - ix;
    U[2] = U[1] * U[1];
    U[3] = U[1] * U[2];

    Real V[4];
    V[0] = (Real)1;
    V[1] = yIndex - iy;
    V[2] = V[1] * V[1];
    V[3] = V[1] * V[2];

    // Compute P = M*U and Q = M*V.
    Real P[4], Q[4];
    for (int row = 0; row < 4; ++row)
    {
        P[row] = (Real)0;
        Q[row] = (Real)0;
        for (int col = 0; col < 4; ++col)
        {
            P[row] += mBlend[row][col] * U[col];
            Q[row] += mBlend[row][col] * V[col];
        }
    }

    // Compute (M*U)^t D (M*V) where D is the 4x4 subimage containing (x,y).
    --ix;
    --iy;
    Real result = (Real)0;
    for (int row = 0; row < 4; ++row)
    {
        int yClamp = iy + row;
        if (yClamp < 0)
        {
            yClamp = 0;
        }
        else if (yClamp > mYBound - 1)
        {
            yClamp = mYBound - 1;
        }

        for (int col = 0; col < 4; ++col)
        {
            int xClamp = ix + col;
            if (xClamp < 0)
            {
                xClamp = 0;
            }
            else if (xClamp > mXBound - 1)
            {
                xClamp = mXBound - 1;
            }

            result += P[col] * Q[row] * mF[xClamp + mXBound * yClamp];
        }
    }

    return result;
}

template <typename Real>
Real IntpBicubic2<Real>::operator()(int xOrder, int yOrder, Real x, Real y)
const
{
    // Compute x-index and clamp to image.
    Real xIndex = (x - mXMin) * mInvXSpacing;
    int ix = static_cast<int>(xIndex);
    if (ix < 0)
    {
        ix = 0;
    }
    else if (ix >= mXBound)
    {
        ix = mXBound - 1;
    }

    // Compute y-index and clamp to image.
    Real yIndex = (y - mYMin) * mInvYSpacing;
    int iy = static_cast<int>(yIndex);
    if (iy < 0)
    {
        iy = 0;
    }
    else if (iy >= mYBound)
    {
        iy = mYBound - 1;
    }

    Real U[4], dx, xMult;
    switch (xOrder)
    {
    case 0:
        dx = xIndex - ix;
        U[0] = (Real)1;
        U[1] = dx;
        U[2] = dx * U[1];
        U[3] = dx * U[2];
        xMult = (Real)1;
        break;
    case 1:
        dx = xIndex - ix;
        U[0] = (Real)0;
        U[1] = (Real)1;
        U[2] = ((Real)2) * dx;
        U[3] = ((Real)3) * dx * dx;
        xMult = mInvXSpacing;
        break;
    case 2:
        dx = xIndex - ix;
        U[0] = (Real)0;
        U[1] = (Real)0;
        U[2] = (Real)2;
        U[3] = (Real)6 * dx;
        xMult = mInvXSpacing * mInvXSpacing;
        break;
    case 3:
        U[0] = (Real)0;
        U[1] = (Real)0;
        U[2] = (Real)0;
        U[3] = (Real)6;
        xMult = mInvXSpacing * mInvXSpacing * mInvXSpacing;
        break;
    default:
        return (Real)0;
    }

    Real V[4], dy, yMult;
    switch (yOrder)
    {
    case 0:
        dy = yIndex - iy;
        V[0] = (Real)1;
        V[1] = dy;
        V[2] = dy * V[1];
        V[3] = dy * V[2];
        yMult = (Real)1;
        break;
    case 1:
        dy = yIndex - iy;
        V[0] = (Real)0;
        V[1] = (Real)1;
        V[2] = ((Real)2) * dy;
        V[3] = ((Real)3) * dy * dy;
        yMult = mInvYSpacing;
        break;
    case 2:
        dy = yIndex - iy;
        V[0] = (Real)0;
        V[1] = (Real)0;
        V[2] = (Real)2;
        V[3] = ((Real)6) * dy;
        yMult = mInvYSpacing * mInvYSpacing;
        break;
    case 3:
        V[0] = (Real)0;
        V[1] = (Real)0;
        V[2] = (Real)0;
        V[3] = (Real)6;
        yMult = mInvYSpacing * mInvYSpacing * mInvYSpacing;
        break;
    default:
        return (Real)0;
    }

    // Compute P = M*U and Q = M*V.
    Real P[4], Q[4];
    for (int row = 0; row < 4; ++row)
    {
        P[row] = (Real)0;
        Q[row] = (Real)0;
        for (int col = 0; col < 4; ++col)
        {
            P[row] += mBlend[row][col] * U[col];
            Q[row] += mBlend[row][col] * V[col];
        }
    }

    // Compute (M*U)^t D (M*V) where D is the 4x4 subimage containing (x,y).
    --ix;
    --iy;
    Real result = (Real)0;
    for (int row = 0; row < 4; ++row)
    {
        int yClamp = iy + row;
        if (yClamp < 0)
        {
            yClamp = 0;
        }
        else if (yClamp > mYBound - 1)
        {
            yClamp = mYBound - 1;
        }

        for (int col = 0; col < 4; ++col)
        {
            int xClamp = ix + col;
            if (xClamp < 0)
            {
                xClamp = 0;
            }
            else if (xClamp > mXBound - 1)
            {
                xClamp = mXBound - 1;
            }

            result += P[col] * Q[row] * mF[xClamp + mXBound * yClamp];
        }
    }
    result *= xMult * yMult;

    return result;
}


}