// Geometric Tools, LLC
// Copyright (c) 1998-2012
// 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: 5.0.2 (2010/12/11)

#include "BSplineBasis.h"

namespace NURBS
{
//----------------------------------------------------------------------------
template <typename Real>
BSplineBasis<Real>::BSplineBasis ()
{
}
//----------------------------------------------------------------------------
template <typename Real>
BSplineBasis<Real>::BSplineBasis (int numCtrlPoints, int degree, bool open)
{
    Create(numCtrlPoints, degree, open);
}
//----------------------------------------------------------------------------
template <typename Real>
void BSplineBasis<Real>::Create (int numCtrlPoints, int degree, bool open)
{
    mUniform = true;

    int i, numKnots = Initialize(numCtrlPoints, degree, open);
    Real factor = ((Real)1)/(mNumCtrlPoints - mDegree);
    if (mOpen)
    {
        for (i = 0; i <= mDegree; ++i)
        {
            mKnot[i] = (Real)0;
        }

        for (/**/; i < mNumCtrlPoints; ++i)
        {
            mKnot[i] = (i - mDegree)*factor;
        }

        for (/**/; i < numKnots; ++i)
        {
            mKnot[i] = (Real)1;
        }
    }
    else
    {
        for (i = 0; i < numKnots; ++i)
        {
            mKnot[i] = (i - mDegree)*factor;
        }
    }
}

template <typename Real>
Array1D_Real BSplineBasis<Real>::uniformKnotVector(int num_ctrl_points, int degree)
{
	Real factor = ((Real)1)/(num_ctrl_points - degree);
	int numKnots = num_ctrl_points + degree + 1;

	Array1D_Real knot(numKnots,0);

	int i=0;

	for (i = 0; i <= degree; ++i)	knot[i] = 0.0;

	for (/**/; i < num_ctrl_points; ++i)
	{
		knot[i] = (i - degree) * factor;
	}

	for (/**/; i < numKnots; ++i)	knot[i] = 1.0;

	return knot;
}

//----------------------------------------------------------------------------
template <typename Real>
BSplineBasis<Real>::BSplineBasis (int numCtrlPoints, int degree,
    const Real* interiorKnot)
{
    Create(numCtrlPoints, degree, interiorKnot);
}
//----------------------------------------------------------------------------
template <typename Real>
void BSplineBasis<Real>::Create (int numCtrlPoints, int degree,
    const Real* interiorKnot)
{
    mUniform = false;

    int i, numKnots = Initialize(numCtrlPoints, degree, true);
    for (i = 0; i <= mDegree; ++i)
    {
        mKnot[i] = (Real)0;
    }

    for (int j = 0; i < mNumCtrlPoints; ++i, ++j)
    {
        mKnot[i] = interiorKnot[j];
    }

    for (/**/; i < numKnots; ++i)
    {
        mKnot[i] = (Real)1;
    }
}
//----------------------------------------------------------------------------
template <typename Real>
int BSplineBasis<Real>::GetNumCtrlPoints () const
{
    return mNumCtrlPoints;
}
//----------------------------------------------------------------------------
template <typename Real>
int BSplineBasis<Real>::GetDegree () const
{
    return mDegree;
}
//----------------------------------------------------------------------------
template <typename Real>
bool BSplineBasis<Real>::IsOpen () const
{
    return mOpen;
}
//----------------------------------------------------------------------------
template <typename Real>
bool BSplineBasis<Real>::IsUniform () const
{
    return mUniform;
}
//----------------------------------------------------------------------------
template <typename Real>
Real BSplineBasis<Real>::GetD0 (int i) const
{
    return mBD0[mDegree][i];
}
//----------------------------------------------------------------------------
template <typename Real>
Real BSplineBasis<Real>::GetD1 (int i) const
{
    return mBD1[mDegree][i];
}
//----------------------------------------------------------------------------
template <typename Real>
Real BSplineBasis<Real>::GetD2 (int i) const
{
    return mBD2[mDegree][i];
}
//----------------------------------------------------------------------------
template <typename Real>
Real BSplineBasis<Real>::GetD3 (int i) const
{
    return mBD3[mDegree][i];
}
//----------------------------------------------------------------------------
template <typename Real>
int BSplineBasis<Real>::Initialize (int numCtrlPoints, int degree, bool open)
{
    assertion(numCtrlPoints >= 2, "Invalid input\n");
    assertion(1 <= degree && degree <= numCtrlPoints-1, "Invalid input\n");

    mNumCtrlPoints = numCtrlPoints;
    mDegree = degree;
    mOpen = open;

    int numKnots = mNumCtrlPoints + mDegree + 1;
    mKnot = new1<Real>(numKnots);

	mBD0 = Allocate();
	mBD1.clear();
	mBD2.clear();
	mBD3.clear();

    return numKnots;
}
//----------------------------------------------------------------------------
template <typename Real>
void BSplineBasis<Real>::SetKnot (int j, Real value)
{
    if (!mUniform)
    {
        // Access only allowed to elements d+1 <= i <= n.
        int i = j + mDegree + 1;
        if (mDegree + 1 <= i && i <= mNumCtrlPoints)
        {
            mKnot[i] = value;
        }
        else
        {
            assertion(false, "Knot index out of range.\n");
        }
    }
    else
    {
        assertion(false, "Knots cannot be set for uniform splines.\n");
    }
}
//----------------------------------------------------------------------------
template <typename Real>
Real BSplineBasis<Real>::GetKnot (int j) const
{
    // Access only allowed to elements d+1 <= i <= n.
    int i = j + mDegree + 1;
    if (mDegree + 1 <= i && i <= mNumCtrlPoints)
    {
        return mKnot[i];
    }

    assertion(false, "Knot index out of range.\n");
    return std::numeric_limits<Real>::max();
}
//----------------------------------------------------------------------------
template <typename Real>
int BSplineBasis<Real>::GetKey (Real& t) const
{
    if (mOpen)
    {
        // Open splines clamp to [0,1].
        if (t <= (Real)0)
        {
            t = (Real)0;
            return mDegree;
        }
        else if (t >= (Real)1)
        {
            t = (Real)1;
            return mNumCtrlPoints - 1;
        }
    }
    else
    {
        // Periodic splines wrap to [0,1).
        if (t < (Real)0 || t >= (Real)1)
        {
            t -= floor(t);
        }
    }


    int i;

    if (mUniform)
    {
        i = mDegree + (int)((mNumCtrlPoints - mDegree)*t);
    }
    else
    {
        for (i = mDegree + 1; i <= mNumCtrlPoints; ++i)
        {
            if (t < mKnot[i])
            {
                break;
            }
        }
        --i;
    }

    return i;
}

//----------------------------------------------------------------------------
template <typename Real>
Array2D_Real BSplineBasis<Real>::Allocate()
{
	int numRows = mDegree + 1;
	int numCols = mNumCtrlPoints + mDegree;

	return Array2D_Real( numRows, Array1D_Real( numCols, 0.0 ) );
}

//----------------------------------------------------------------------------
template <typename Real>
void BSplineBasis<Real>::Compute (Real t, unsigned int order, int& minIndex, int& maxIndex)
{
    assertion(order <= 3, "Only derivatives to third order supported\n");

	if (order >= 1)
	{
		if (!mBD1.size())
		{
			mBD1 = Allocate();
		}

		if (order >= 2)
		{
			if (!mBD2.size())
			{
				mBD2 = Allocate();
			}

			if (order >= 3)
			{
				if (!mBD3.size())
				{
					mBD3 = Allocate();
				}
			}
		}
	}

    int i = GetKey(t);
    mBD0[0][i] = (Real)1;

    if (order >= 1)
    {
        mBD1[0][i] = (Real)0;
        if (order >= 2)
        {
            mBD2[0][i] = (Real)0;
            if (order >= 3)
            {
                mBD3[0][i] = (Real)0;
            }
        }
    }

    Real n0 = t - mKnot[i], n1 = mKnot[i+1] - t;
    Real invD0, invD1;
    int j;
    for (j = 1; j <= mDegree; j++)
    {
        invD0 = ((Real)1)/(mKnot[i+j] - mKnot[i]);
        invD1 = ((Real)1)/(mKnot[i+1] - mKnot[i-j+1]);

        mBD0[j][i] = n0*mBD0[j-1][i]*invD0;
        mBD0[j][i-j] = n1*mBD0[j-1][i-j+1]*invD1;

        if (order >= 1)
        {
            mBD1[j][i] = (n0*mBD1[j-1][i] + mBD0[j-1][i])*invD0;
            mBD1[j][i-j] = (n1*mBD1[j-1][i-j+1] - mBD0[j-1][i-j+1])*invD1;

            if (order >= 2)
            {
                mBD2[j][i] = (n0*mBD2[j-1][i] + ((Real)2)*mBD1[j-1][i])*invD0;
                mBD2[j][i-j] = (n1*mBD2[j-1][i-j+1] -
                    ((Real)2)*mBD1[j-1][i-j+1])*invD1;

                if (order >= 3)
                {
                    mBD3[j][i] = (n0*mBD3[j-1][i] +
                        ((Real)3)*mBD2[j-1][i])*invD0;
                    mBD3[j][i-j] = (n1*mBD3[j-1][i-j+1] -
                        ((Real)3)*mBD2[j-1][i-j+1])*invD1;
                }
            }
        }
    }

    for (j = 2; j <= mDegree; ++j)
    {
        for (int k = i-j+1; k < i; ++k)
        {
            n0 = t - mKnot[k];
            n1 = mKnot[k+j+1] - t;
            invD0 = ((Real)1)/(mKnot[k+j] - mKnot[k]);
            invD1 = ((Real)1)/(mKnot[k+j+1] - mKnot[k+1]);

            mBD0[j][k] = n0*mBD0[j-1][k]*invD0 + n1*mBD0[j-1][k+1]*invD1;

            if (order >= 1)
            {
                mBD1[j][k] = (n0*mBD1[j-1][k]+mBD0[j-1][k])*invD0 +
                    (n1*mBD1[j-1][k+1]-mBD0[j-1][k+1])*invD1;

                if (order >= 2)
                {
                    mBD2[j][k] = (n0*mBD2[j-1][k] +
                        ((Real)2)*mBD1[j-1][k])*invD0 +
                        (n1*mBD2[j-1][k+1] - ((Real)2)*mBD1[j-1][k+1])*invD1;

                    if (order >= 3)
                    {
                        mBD3[j][k] = (n0*mBD3[j-1][k] +
                            ((Real)3)*mBD2[j-1][k])*invD0 +
                            (n1*mBD3[j-1][k+1] - ((Real)3)*
                            mBD2[j-1][k+1])*invD1;
                    }
                }
            }
        }
    }

    minIndex = i - mDegree;
    maxIndex = i;
}
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// Explicit instantiation.
//----------------------------------------------------------------------------
//template
//class BSplineBasis<float>;

template
class BSplineBasis<double>;
//----------------------------------------------------------------------------
}