/******************************************************************************
* CBzrEval.c - Bezier curves handling routines - evaluation routines.	      *
*******************************************************************************
* Written by Gershon Elber, Mar. 90.					      *
******************************************************************************/

#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include "cagd_loc.h"

static int PowerCacheFineNess,
    BezierCacheEnabled = FALSE,
    CacheFineNess = 0;
static CagdRType *BezierCache[CAGD_MAX_BEZIER_CACHE_ORDER + 1]
			     [CAGD_MAX_BEZIER_CACHE_ORDER + 1];

static int IChooseK[CAGD_MAX_BEZIER_CACHE_ORDER + 1]
		   [CAGD_MAX_BEZIER_CACHE_ORDER + 1] =
{
    {    1,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0  },
    {    1,    1,    0,    0,    0,    0,    0,    0,    0,    0,    0  },
    {    1,    2,    1,    0,    0,    0,    0,    0,    0,    0,    0  },
    {    1,    3,    3,    1,    0,    0,    0,    0,    0,    0,    0  },
    {    1,    4,    6,    4,    1,    0,    0,    0,    0,    0,    0  },
    {    1,    5,   10,   10,    5,    1,    0,    0,    0,    0,    0  },
    {    1,    6,   15,   20,   15,    6,    1,    0,    0,    0,    0  },
    {    1,    7,   21,   35,   35,   21,    7,    1,    0,    0,    0  },
    {    1,    8,   28,   56,   70,   56,   28,    8,    1,    0,    0  },
    {    1,    9,   36,   84,  126,  126,   84,   36,    9,    1,    0  },
    {    1,   10,   45,  120,  210,  252,  210,  120,   45,   10,    1  }
};

static CagdRType BzrCrvEvalFuncAux(int i, int k, CagdRType t);
static CagdRType IntPow(CagdRType x, int i);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Sets the bezier sampling cache - if enabled, a Bezier can be evaluated     M
* directly from presampled basis function.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   FineNess:       Number of samples to support (power of 2).               M
*   EnableCache:    Are we really planning on using this thing?              M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvSetCache, evaluation, caching                                      M
*****************************************************************************/
void BzrCrvSetCache(int FineNess, CagdBType EnableCache)
{
    int i, j, k;

    if (BezierCacheEnabled == EnableCache && FineNess == CacheFineNess)
	return;

    if (BezierCacheEnabled) {
	/* Free old cache if was one. */
	for (k = 2; k <= CAGD_MAX_BEZIER_CACHE_ORDER; k++)
	    for (i = 0; i < k; i++)
		if (BezierCache[k][i] != NULL) {
		    IritFree((VoidPtr) BezierCache[k][i]);
		    BezierCache[k][i] = NULL;
		}
    }

    if ((BezierCacheEnabled = EnableCache) != FALSE) {
	/* Allocate the new cache and initalize it: */
	if (FineNess < 2)
	    FineNess = 2;
	if (FineNess > CAGD_MAX_BEZIER_CACHE_FINENESS)
	    FineNess = CAGD_MAX_BEZIER_CACHE_FINENESS;
	CacheFineNess = FineNess;
	PowerCacheFineNess = 1 << FineNess;

	for (k = 2; k <= CAGD_MAX_BEZIER_CACHE_ORDER; k++)/* For all orders. */
	    for (i = 0; i < k; i++) { 	 /* Allocate and set all basis func. */
		BezierCache[k][i] = (CagdRType *)
			IritMalloc(sizeof(CagdRType) * PowerCacheFineNess);
		for (j = 0; j < PowerCacheFineNess; j++)
		    BezierCache[k][i][j] = BzrCrvEvalFuncAux(i, k,
				((CagdRType) j) / (PowerCacheFineNess - 1));
	    }
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Assumes Vec holds control points for scalar bezier curve of order Order,   M
* and evaluates and returns that curve value at parameter value t.	     M
*   Vec is incremented by VecInc (usually by 1) after each iteration.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Vec:            Coefficents of a scalar Bspline univariate function.     M
*   VecInc:         Step to move along Vec.                                  M
*   Order:          Order of associated geometry.                            M
*   t:              Parameter value where to evaluate the curve.             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType:      Geometry's value at parameter value t.                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvEvalVecAtParam, evaluation                                         M
*****************************************************************************/
CagdRType BzrCrvEvalVecAtParam(CagdRType *Vec,
			       int VecInc,
			       int Order,
			       CagdRType t)
{
    int i;
    CagdRType
	R = 0.0;

    if (VecInc == 1)
	for (i = 0; i < Order; i++)
	    R += BzrCrvEvalFuncAux(i, Order, t) * *Vec++;
    else
	for (i = 0; i < Order; i++) {
	    R += BzrCrvEvalFuncAux(i, Order, t) * *Vec;
	    Vec += VecInc;
	}

    return R;
}
/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a pointer to a static data, holding the value of the curve at      M
* given parametric location t. The curve is assumed to be Bezier.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      To evaluate at the given parametric location t.                M
*   t:        The parameter value at which the curve Crv is to be evaluated. M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  A vector holding all the coefficients of all components    M
*                 of curve Crv's point type. If for example the curve's      M
*                 point type is P2, the W, X, and Y will be saved in the     M
*                 first three locations of the returned vector. The first    M
*                 location (index 0) of the returned vector is reserved for  M
*                 the rational coefficient W and XYZ always starts at second M
*                 location of the returned vector (index 1).                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvEvalAtParam, evaluation                                            M
*****************************************************************************/
CagdRType *BzrCrvEvalAtParam(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdRType Pt[CAGD_MAX_PT_COORD];
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j,
	k = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType B;

    for (j = IsNotRational; j <= MaxCoord; j++)
	Pt[j] = 0.0;

    for (i = 0; i < k; i++) {
	B = BzrCrvEvalFuncAux(i, k, t);
	for (j = IsNotRational; j <= MaxCoord; j++)
	    Pt[j] += B * Crv -> Points[j][i];
    }

    return Pt;
}
/*****************************************************************************
* DESCRIPTION:                                                               M
* Samples the curve at FineNess location equally spaced in the curve's       M
* parametric domain.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:         To approximate as a polyline.                               M
*   FineNess:    Control on number of samples.                               M
*   Points:      Where to put the resulting polyline.                        M
*   A:           Optional alpha matrix for refinement.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:         The actual number of samples placed in Points. Always       M
*                less than or eaul to FineNess.                              M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvEvalToPolyline, conversion, refinement, evaluation                 M
*****************************************************************************/
void BzrCrvEvalToPolyline(CagdCrvStruct *Crv,
			  int FineNess,
			  CagdRType *Points[])
{
    CagdBType
	UseCache = BezierCacheEnabled &&
		   FineNess == PowerCacheFineNess &&
		   Crv -> Order <= CAGD_MAX_BEZIER_CACHE_ORDER,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, Count,
	k = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType B;

    if (UseCache) {
	for (Count = 0; Count < PowerCacheFineNess; Count++) {
	    for (j = IsNotRational; j <= MaxCoord; j++)
		Points[j][Count] = 0.0;
	    for (i = 0; i < k; i++) {
		B = BezierCache[k][i][Count];
		for (j = IsNotRational; j <= MaxCoord; j++)
		    Points[j][Count] += B * Crv -> Points[j][i];
	    }
	}
    }
    else {
	for (Count = 0; Count < FineNess; Count++) {
	    for (j = IsNotRational; j <= MaxCoord; j++)
		Points[j][Count] = 0.0;
	    for (i = 0; i < k; i++) {
		B = BzrCrvEvalFuncAux(i, k,
				      ((CagdRType) Count) / (FineNess - 1));
		for (j = IsNotRational; j <= MaxCoord; j++)
		    Points[j][Count] += Crv -> Points[j][i] * B;
	    }
	}
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Evaluates the i'th Bezier basis function of order k, at parametric value t *
* (t in [0..1]).							     *
*   This function is:		   i	 k - i - 1          i		     *
*			Bi,k(t) = ( ) * t          * (1 - t)		     *
*				   k					     *
*                                                                            *
* PARAMETERS:                                                                *
*   i:   I'th basi function.                                                 *
*   k:   Degree of the curve.                                                *
*   t:   Parameter value at which to evaluate the Bezier basis function.     *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdRType:  Value of the i'th Bezier basis function of degree k at t.    *
*****************************************************************************/
static CagdRType BzrCrvEvalFuncAux(int i, int k, CagdRType t)
{
    if (APX_EQ(t, 0.0))
	return (CagdRType) (i == 0);
    else if (APX_EQ(t, 1.0))
	return (CagdRType) (i == k - 1);
    else
	return (k >= CAGD_MAX_BEZIER_CACHE_ORDER ?
		    CagdIChooseK(i, k - 1) : (CagdRType) IChooseK[k - 1][i]) *
	       IntPow(1.0 - t, k - i - 1) * IntPow(t, i);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Computes x to the power of i, i integer.				     *
*                                                                            *
*                                                                            *
* PARAMETERS:                                                                *
*   x, i: Description says it all.                                           *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdRType:  Integer power of x, computed using pwoer of 2's.             *
*****************************************************************************/
static CagdRType IntPow(CagdRType x, int i)
{
    CagdRType Power, RetVal;

    for (Power = x, RetVal = 1.0; i != 0; i >>= 1) {
	if (i & 0x01)
	   RetVal *= Power;
	
	Power = SQR(Power);
    }

    return RetVal;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluates the following (in floating point arithmetic):		     M
*			 k         k!					     V
*			( ) = -------------				     V
*			 i    i! * (k - i)!				     V
*                                                                            *
* PARAMETERS:                                                                M
*   i, k:   Coefficients of i choose k.                                      M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType:   Result of i choose k, in floating point, to prevent from    M
*                overflows.                                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdIChooseK, evaluation, combinatorics                                  M
*****************************************************************************/
CagdRType CagdIChooseK(int i, int k)
{
    int j;
    CagdRType
	c = 1.0;

    if ((k >> 1) > i) {				/* i is less than half of k: */
	for (j = k - i + 1; j <= k; j++)
	    c *= j;
	for (j = 2; j <= i; j++)
	    c /= j;
    }
    else {
	for (j = i + 1; j <= k; j++)
	    c *= j;
	for (j = 2; j <= k - i; j++)
	    c /= j;
    }

    return c;
}

#ifdef K_CHOOSE_I_GEN

/*****************************************************************************
* DESCRIPTION:                                                               *
* Evaluate the following (in integer arithmetic):			     *
*			 k         k!					     *
*			( ) = -------------				     *
*			 i    i! * (k - i)!				     *
*                                                                            *
* PARAMETERS:                                                                *
*   i, k:   Coefficients of i choose k.                                      *
*                                                                            *
* RETURN VALUE:                                                              *
*   int: Result of i choose k.                                               *
*****************************************************************************/
static int IChooseKGenOne(int i, int k)
{
    int j;
    long c = 1;

    if ((k >> 1) > i) {				/* i is less than half of k: */
	for (j = k - i + 1; j <= k; j++)
	    c *= j;
	for (j = 2; j <= i; j++)
	    c /= j;
    }
    else {
	for (j = i + 1; j <= k; j++)
	    c *= j;
	for (j = 2; j <= k - i; j++)
	    c /= j;
    }

    return (int) c;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Evaluates IChooseK for all possiblities and prints them for the	     *
* IChooseK table defined in the beginning of this file.			     *
*                                                                            *
* PARAMETERS:                                                                *
*   None                                                                     *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
void IChooseKGen(void)
{
    int i, j;

    for (i = 0; i <= MAX_BEZIER_CACHE_ORDER; i++) {
	printf(" },\n    {");
	for (j = 0; j <= MAX_BEZIER_CACHE_ORDER; j++)
	    printf(" %4d,", j <= i ? IChooseKGenOne(j ,i) : 0);
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Main routine to create the table in the beginning of this file.	     *
* PARAMETERS:                                                                *
*   None                                                                     *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
void main(void)
{
    IChooseKGen();
}

#endif /* K_CHOOSE_I_GEN */
