// Geometric Tools, Inc.
// http://www.geometrictools.com
// Copyright (c) 1998-2005.  All Rights Reserved
//
// The Wild Magic Library (WM3) source code is supplied under the terms of
// the license agreement
//     http://www.geometrictools.com/License/WildMagic3License.pdf
// and may not be copied or disclosed except in accordance with the terms
// of that agreement.

#include "Wm3BSplineReduction.h"
#include "Wm3Integrate1.h"
#include "Wm3Intersector1.h"
#include "Wm3LinearSystem.h"

namespace Wm3
{
//----------------------------------------------------------------------------
template <class Real, class TVector>
BSplineReduction<Real,TVector>::BSplineReduction (int iQuantity,
    const TVector* akCtrl, int iDegree, Real fFraction, int& riQuantityOut,
    TVector*& rakCtrlOut)
{
    // Check for valid input.  If invalid, return a "null" curve.
    assert( iQuantity >= 2 );
    assert( akCtrl );
    assert( 1 <= iDegree && iDegree < iQuantity );
    if ( iQuantity < 2 || !akCtrl || iDegree < 1 || iDegree >= iQuantity )
    {
        riQuantityOut = 0;
        rakCtrlOut = 0;
        return;
    }

    // clamp the number of control points to [degree+1,quantity-1]
    riQuantityOut = (int)(fFraction*iQuantity);
    if ( riQuantityOut > iQuantity )
    {
        riQuantityOut = iQuantity;
        rakCtrlOut = new TVector[riQuantityOut];
        memcpy(rakCtrlOut,akCtrl,iQuantity*sizeof(TVector));
        return;
    }
    if ( riQuantityOut < iDegree+1 )
        riQuantityOut = iDegree+1;

    // allocate output control points and set all to the zero vector
    rakCtrlOut = new TVector[riQuantityOut];

    // Set up basis function parameters.  Function 0 corresponds to the
    // output curve.  Function 1 corresponds to the input curve.
    m_iDegree = iDegree;
    m_aiQuantity[0] = riQuantityOut;
    m_aiQuantity[1] = iQuantity;

    for (int j = 0; j <= 1; j++)
    {
        m_aiNumKnots[j] = m_aiQuantity[j] + m_iDegree + 1;
        m_aafKnot[j] = new Real[m_aiNumKnots[j]];

        int i;
        for (i = 0; i <= m_iDegree; i++)
            m_aafKnot[j][i] = (Real)0.0;

        Real fFactor = ((Real)1.0)/(Real)(m_aiQuantity[j]-m_iDegree);
        for (/**/; i < m_aiQuantity[j]; i++)
            m_aafKnot[j][i] = (i-m_iDegree)*fFactor;

        for (/**/; i < m_aiNumKnots[j]; i++)
            m_aafKnot[j][i] = (Real)1.0;
    }

    // construct matrix A (depends only on the output basis function)
    Real fValue, fTMin, fTMax;
    int i0, i1;

    m_aiBasis[0] = 0;
    m_aiBasis[1] = 0;

    BandedMatrix<Real> kA(m_aiQuantity[0],m_iDegree,m_iDegree);
    for (i0 = 0; i0 < m_aiQuantity[0]; i0++)
    {
        m_aiIndex[0] = i0;
        fTMax = MaxSupport(0,i0);

        for (i1 = i0; i1 <= i0 + m_iDegree && i1 < m_aiQuantity[0]; i1++)
        {
            m_aiIndex[1] = i1;
            fTMin = MinSupport(0,i1);

            fValue = Integrate1<Real>::RombergIntegral(8,fTMin,fTMax,
                Integrand,this);
            kA(i0,i1) = fValue;
            kA(i1,i0) = fValue;
        }
    }

    // construct A^{-1}
    GMatrix<Real> kInvA(m_aiQuantity[0],m_aiQuantity[0]);
    bool bResult = LinearSystem<Real>().Invert(kA,kInvA);
    assert( bResult );
    (void)bResult;  // avoid compiler warning in release build

    // construct B (depends on both input and output basis functions)
    m_aiBasis[1] = 1;
    GMatrix<Real> kB(m_aiQuantity[0],m_aiQuantity[1]);
    for (i0 = 0; i0 < m_aiQuantity[0]; i0++)
    {
        m_aiIndex[0] = i0;
        Real fTMin0 = MinSupport(0,i0);
        Real fTMax0 = MaxSupport(0,i0);

        for (i1 = 0; i1 < m_aiQuantity[1]; i1++)
        {
            m_aiIndex[1] = i1;
            Real fTMin1 = MinSupport(1,i1);
            Real fTMax1 = MaxSupport(1,i1);

            Intersector1<Real> kIntr(fTMin0,fTMax0,fTMin1,fTMax1);
            kIntr.Find();
            int iIQuantity = kIntr.GetQuantity();

            if ( iIQuantity == 2 )
            {
                fValue = Integrate1<Real>::RombergIntegral(8,
                    kIntr.GetOverlap(0),kIntr.GetOverlap(1),Integrand,this);
                kB[i0][i1] = fValue;
            }
            else
            {
                kB[i0][i1] = (Real)0.0;
            }
        }
    }

    // construct A^{-1}*B
    GMatrix<Real> kProd = kInvA*kB;

    // construct the control points for the least-squares curve
    memset(rakCtrlOut,0,riQuantityOut*sizeof(TVector));
    for (i0 = 0; i0 < m_aiQuantity[0]; i0++)
    {
        for (i1 = 0; i1 < m_aiQuantity[1]; i1++)
            rakCtrlOut[i0] += akCtrl[i1]*kProd[i0][i1];
    }
}
//----------------------------------------------------------------------------
template <class Real, class TVector>
BSplineReduction<Real,TVector>::~BSplineReduction ()
{
    delete[] m_aafKnot[0];
    delete[] m_aafKnot[1];
}
//----------------------------------------------------------------------------
template <class Real, class TVector>
Real BSplineReduction<Real,TVector>::MinSupport (int iBasis, int i) const
{
    return m_aafKnot[iBasis][i];
}
//----------------------------------------------------------------------------
template <class Real, class TVector>
Real BSplineReduction<Real,TVector>::MaxSupport (int iBasis, int i) const
{
    return m_aafKnot[iBasis][i+1+m_iDegree];
}
//----------------------------------------------------------------------------
template <class Real, class TVector>
Real BSplineReduction<Real,TVector>::F (int iBasis, int i, int j, Real fTime)
{
    if ( j > 0 )
    {
        Real fResult = (Real)0.0;

        Real fDenom = m_aafKnot[iBasis][i+j] - m_aafKnot[iBasis][i];
        if ( fDenom > (Real)0.0 )
        {
            fResult += (fTime - m_aafKnot[iBasis][i]) *
                F(iBasis,i,j-1,fTime)/fDenom;
        }

        fDenom = m_aafKnot[iBasis][i+j+1] - m_aafKnot[iBasis][i+1];
        if ( fDenom > (Real)0.0 )
        {
            fResult += (m_aafKnot[iBasis][i+j+1] - fTime) *
                F(iBasis,i+1,j-1,fTime)/fDenom;
        }

        return fResult;
    }

    if ( m_aafKnot[iBasis][i] <= fTime && fTime < m_aafKnot[iBasis][i+1] )
        return (Real)1.0;
    else
        return (Real)0.0;
}
//----------------------------------------------------------------------------
template <class Real, class TVector>
Real BSplineReduction<Real,TVector>::Integrand (Real fTime, void* pvThis)
{
    BSplineReduction<Real,TVector>& rkSelf =
        *(BSplineReduction<Real,TVector>*)pvThis;

    Real fValue0 = rkSelf.F(rkSelf.m_aiBasis[0],rkSelf.m_aiIndex[0],
        rkSelf.m_iDegree,fTime);

    Real fValue1 = rkSelf.F(rkSelf.m_aiBasis[1],rkSelf.m_aiIndex[1],
        rkSelf.m_iDegree,fTime);

    return fValue0*fValue1;
}
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// explicit instantiation
//----------------------------------------------------------------------------
// float
template WM3_ITEM
class BSplineReduction<float,Vector2f>;

template WM3_ITEM
class BSplineReduction<float,Vector3f>;

// double
template WM3_ITEM
class BSplineReduction<double,Vector2d>;

template WM3_ITEM
class BSplineReduction<double,Vector3d>;
//----------------------------------------------------------------------------
}



