/******************************************************************************
* 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 BezierCacheEnabled = FALSE;
static int CacheFineNess = 0, PowerCacheFineNess;
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);

/******************************************************************************
* Sets the bezier sampling cache - if enabled, a bezier can be evaluated      *
* directly from pre sampled basis function.				      *
******************************************************************************/
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));
	    }
    }
}

/******************************************************************************
* Assumes Vec holds control points for scalar bezier curve of order Order,    *
* and evaluates and returns that curve value at parameter value t.	      *
* Vec is incremented by VecInc (usually by 1) after each iteration.	      *
******************************************************************************/
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;
}

/******************************************************************************
* Returns a pointer to a static data, holding the value of the curve at given *
* parametric location t. The curve is assumed to be Bezier.		      *
******************************************************************************/
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;
}

/******************************************************************************
* Samples the curves at FineNess location equally spaced in the Bezier        *
* parametric domain [0..1]. If Cache is enabled, and FineNess is power of     *
* two, upto or equal to CacheFineNess, the cache is used, otherwise the       *
* points are evaluated manually for each of the samples.		      *
* Data is saved at the Points array of vectors (according to Curve PType),    *
* each vector is assumed to be allocated for FineNess CagdRType points.	      *
******************************************************************************/
void BzrCrvEvalToPolyline(CagdCrvStruct *Crv, int FineNess,
							CagdRType *Points[])
{
    CagdBType
	UseCache = BezierCacheEnabled &&
		   FineNess == CacheFineNess &&
		   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 {
	FineNess = 1 << FineNess;
	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;
	    }
	}
    }
}

/******************************************************************************
* Evaluate the i bezier 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					      *
******************************************************************************/
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);
}

/******************************************************************************
* Computes x to the power of i, i integer.				      *
******************************************************************************/
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;
}

/******************************************************************************
* Evaluate the following (in floating point arithmetic):		      *
*			 k         k!					      *
*			( ) = -------------				      *
*			 i    i! * (k - i)!				      *
******************************************************************************/
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

/******************************************************************************
* Evaluate the following (in integer arithmetic):			      *
*			 k         k!					      *
*			( ) = -------------				      *
*			 i    i! * (k - i)!				      *
******************************************************************************/
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;
}

/******************************************************************************
* Evaluate IChooseK for all possiblities and print them for the		      *
* IChooseK table defined in the beginning of this file.			      *
******************************************************************************/
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);
    }
}

/******************************************************************************
* main to create the table in the befinning of this file.		      *
******************************************************************************/
void main(void)
{
    IChooseKGen();
}

#endif /* K_CHOOSE_I_GEN */