/******************************************************************************
* CagdOslo.c - an implementation of the oslo algorithm (1).		      *
*******************************************************************************
* Written by Gershon Elber, Aug. 92.					      *
******************************************************************************/

#include <cagd_loc.h>

#define KNOT_INFINITY 1e20

#ifdef DEBUG
static void PrintAlphaMat(BspKnotAlphaCoeffType *A);
#endif /* DEBUG */

/******************************************************************************
* Computes the values of the alpha coefficients, Ai,k(j) of order k:	      *
*                                                                             *
*         n              n    m                   m    n		      *
*         _              _    _                   _    _		      *
*  C(t) = \ Pi Bi,k(t) = \ Pi \ Ai,k(j) Nj,k(t) = \  ( \ Pi Ai,k(j) ) Nj,k(t) *
*         /              /    /                   /    /		      *
*         -              -    -                   -    -		      *
*        i=0            i=0  j=0                 j=0  i=0		      *
*                                                                             *
* Let T be the original knot vector and let t be the refined one, i.e. T is   *
* a subset of t. 							      *
* The Ai,k(j) are computed from the following recursive	definition:           *
*                                                                             *
*             1, T(i) <= t(i) < T(i+1)                                        *
*           /                                                                 *
* Ai,1(j) =                                                                   *
*           \                                                                 *
*             0, otherwise.                                                   *
*                                                                             *
*                                                                             *
*           T(j+k-1) - T(i)             T(i+k) - T(j+k-1)                     *
* Ai,k(j) = --------------- Ai,k-1(j) + --------------- Ai+1,k-1(j)           *
*           T(i+k-1) - T(i)             T(i+k) - T(i+1)                       *
*                                                                             *
* LengthKVT + k is the length of KVT and similarly LengthKVt + k is the       *
* length of KVt. In other words, LengthKVT and LengthKVt are the control      *
* points len...								      *
*                                                                             *
* The output matrix has LengthKVT rows and LengthKVt columns (#cols > #rows)  *
* ColIndex/Length hold LengthKVt pairs of first non zero scalar and length of *
* non zero values in that column, so not all LengthKVT scalars are blended.   *
******************************************************************************/
BspKnotAlphaCoeffType *BspKnotEvalAlphaCoef(int k, CagdRType *KVT,
				int LengthKVT, CagdRType *KVt, int LengthKVt)
{
    int i, j, o;
    CagdRType *m, **Rows;
    BspKnotAlphaCoeffType
	*A = (BspKnotAlphaCoeffType *)
				     IritMalloc(sizeof(BspKnotAlphaCoeffType));

    A -> Order = k;
    A -> Length = LengthKVT;
    A -> RefLength = LengthKVt;
    A -> Matrix = (CagdRType *) IritMalloc(sizeof(CagdRType) * (LengthKVT + 1)
							     * LengthKVt);
    Rows = A -> Rows = (CagdRType **) IritMalloc(sizeof(CagdRType *) *
							      (LengthKVT + 1));
    A -> ColIndex = (int *) IritMalloc(sizeof(int) * LengthKVt);
    A -> ColLength = (int *) IritMalloc(sizeof(int) * LengthKVt);

    /* Update the row pointers to point onto the matrix rows. */
    for (i = 0, j = 0; i <= LengthKVT; i++, j += LengthKVt)
	Rows[i] = &A -> Matrix[j];

    /* Initialize the matrix with according to order 1: */
    for (m = A -> Matrix, i = (LengthKVT + 1) * LengthKVt; i > 0; i--)
	*m++ = 0.0;
    for (i = 0; i < LengthKVT; i++) {
	CagdRType
	    *RowI = Rows[i],
	    KVTI = KVT[i],
	    KVTI1 = KVT[i + 1];

	for (j = 0; j < LengthKVt; j++, RowI++) {
	    if (KVTI <= KVt[j] && KVt[j] < KVTI1)
		*RowI = 1.0;
	}
    }

    /* Iterate upto order k: */
    for (o = 2; o <= k; o++) {
	for (i = 0; i < LengthKVT; i++) {
	    CagdRType
		*RowI = Rows[i],
		*RowI1 = Rows[i + 1],
		KVTI = KVT[i],
		KVTIO = KVT[i + o],
		t1 = KVT[i + o - 1] - KVTI,
		t2 = KVTIO - KVT[i + 1];

	    /* If ti == 0, the whole term should be ignored. */
	    if (t1 < EPSILON)
		t1 = KNOT_INFINITY;
	    if (t2 < EPSILON)
		t2 = KNOT_INFINITY;

	    for (j = 0; j < LengthKVt; j++, RowI++, RowI1++) {
		int jo = j + o;

		*RowI = *RowI * (KVt[jo - 1] - KVTI) / t1 +
			*RowI1 * (KVTIO - KVt[jo - 1]) / t2;
	    }
	}
    }

    /* Update the row non zero indices. */
    for (i = LengthKVt - 1; i >= 0; i--) {
	for (j = 0; j < LengthKVT && APX_EQ(Rows[j][i], 0.0); j++);
	A -> ColIndex[i] = j;
	for (j = LengthKVT - 1; j >= 0 && APX_EQ(Rows[j][i], 0.0); j--);
	if ((A -> ColLength[i] = j + 1 - A -> ColIndex[i]) <= 0)
	    FATAL_ERROR(CAGD_ERR_DEGEN_ALPHA);
    }

    return A;
}

/******************************************************************************
* Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to      *
* form the new knot vector KVt.						      *
******************************************************************************/
void BspKnotFreeAlphaCoef(BspKnotAlphaCoeffType *A)
{
    IritFree((VoidPtr) A -> Matrix);
    IritFree((VoidPtr) A -> Rows);
    IritFree((VoidPtr) A -> ColIndex);
    IritFree((VoidPtr) A -> ColLength);
    IritFree((VoidPtr) A);
}

/******************************************************************************
* Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to      *
* form the new knot vector KVt.						      *
******************************************************************************/
BspKnotAlphaCoeffType *BspKnotEvalAlphaCoefMerge(int k, CagdRType *KVT,
			       int LengthKVT, CagdRType *NewKV, int LengthNewKV)
{
    int LengthKVt;
    BspKnotAlphaCoeffType *A;

    if (NewKV == NULL || LengthNewKV == 0) {
	A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVT, LengthKVT);
    }
    else {
	CagdRType
	    *KVt = BspKnotMergeTwo(KVT, LengthKVT + k,
				   NewKV, LengthNewKV, 0, &LengthKVt);

	A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVt, LengthKVt - k);

	IritFree(KVt);
    }

    return A;
}

/******************************************************************************
* Prepare a refinement vector for the given knot vector domain with n	      *
* inserted knots equally spaced.					      *
******************************************************************************/
CagdRType *BspKnotPrepEquallySpaced(int n, CagdRType Tmin, CagdRType Tmax)
{
    int i;
    CagdRType dt, t, *RefKV;

    if (n <= 0) {
	FATAL_ERROR(CAGD_ERR_WRONG_INDEX);
	return NULL;
    }

    dt = (Tmax - Tmin) / (n + 1),
    t = Tmin + dt,
    RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * n);

    for (i = 0; i < n; i++, t += dt) {
	RefKV[i] = t;
    }
    return RefKV;
}

#ifdef DEBUG
/******************************************************************************
* Prints the content of the alpha matrix.				      *
******************************************************************************/
static void PrintAlphaMat(BspKnotAlphaCoeffType *A)
{
    int i, j;

    fprintf(stderr, "Order = %d, Length = %d\n", A -> Order, A -> Length);
    for (i = 0; i < A -> Length; i++) {
	for (j = 0; j < A -> RefLength; j++)
	    fprintf(stderr, " %9.5lg", A -> Rows[i][j]);
	fprintf(stderr, "\n");
    }

    fprintf(stderr, "    ");
    for (j = 0; j < A -> RefLength; j++)
	    fprintf(stderr, " %3d %3d |", A -> ColIndex[j], A -> ColLength[j]);
	fprintf(stderr, "\n");
}
#endif /* DEBUG */