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

#include <cagd_loc.h>

#define KNOT_INFINITY		1e20
#define OSLO_EPSILON		1e-20
#define OSLO_APX_EQ(x, y)	(ABS((x) - (y)) < OSLO_EPSILON)

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

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes the values of the alpha coefficients, Ai,k(j) of order k:	     M
*                                                                            M
*         n              n    m                   m    n		     V
*         _              _    _                   _    _		     V
*  C(t) = \ Pi Bi,k(t) = \ Pi \ Ai,k(j) Nj,k(t) = \  ( \ Pi Ai,k(j) ) Nj,k(t)V
*         /              /    /                   /    /		     V
*         -              -    -                   -    -		     V
*        i=0            i=0  j=0                 j=0  i=0		     V
*                                                                            M
* Let T be the original knot vector and let t be the refined one, i.e. T is  M
* a subset of t. 							     M
*   The Ai,k(j) are computed from the following recursive definition:        M
*                                                                            M
*             1, T(i) <= t(i) < T(i+1)                                       V
*           /                                                                V
* Ai,1(j) =                                                                  V
*           \                                                                V
*             0, otherwise.                                                  V
*                                                                            V
*                                                                            V
*           T(j+k-1) - T(i)             T(i+k) - T(j+k-1)                    V
* Ai,k(j) = --------------- Ai,k-1(j) + --------------- Ai+1,k-1(j)          V
*           T(i+k-1) - T(i)             T(i+k) - T(i+1)                      V
*                                                                            M
* LengthKVT + k is the length of KVT and similarly LengthKVt + k is the      M
* length of KVt. In other words, LengthKVT and LengthKVt are the control     M
* points len...								     M
*                                                                            M
* The output matrix has LengthKVT rows and LengthKVt columns (#cols > #rows) M
* ColIndex/Length hold LengthKVt pairs of first non zero scalar and length ofM
* non zero values in that column, so not all LengthKVT scalars are blended.  M
*                                                                            *
* PARAMETERS:                                                                M
*   k:           Order of geometry.                                          M
*   KVT:         Original knot vector.                                       M
*   LengthKVT:   Length of original knot vector.                             M
*   KVt:         Refined knot vector. Must contain all knots of KVT.         M
*   LengthKVt:   Length of refined knot vector.                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   BspKnotAlphaCoeffType *:   A matrix to multiply the coefficients of the  M
*                              geometry using KVT, in order to get the       M
*                              coefficients under the space defined using    M
*                              KVt that represent the same geometry.         M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotEvalAlphaCoef, alpha matrix, refinement                           M
*****************************************************************************/
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 < OSLO_EPSILON)
		t1 = KNOT_INFINITY;
	    if (t2 < OSLO_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 && OSLO_APX_EQ(Rows[j][i], 0.0); j++);
	A -> ColIndex[i] = j;
	for (j = LengthKVT - 1; j >= 0 && OSLO_APX_EQ(Rows[j][i], 0.0); j--);
	if ((A -> ColLength[i] = j + 1 - A -> ColIndex[i]) <= 0)
	    CAGD_FATAL_ERROR(CAGD_ERR_DEGEN_ALPHA);
    }

    return A;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Frees the BspKnotAlphaCoeffType data structrure.                           M
*                                                                            *
* PARAMETERS:                                                                M
*   A:      Alpha matrix to free.                                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotFreeAlphaCoef, alpha matrix, refinement                           M
*****************************************************************************/
void BspKnotFreeAlphaCoef(BspKnotAlphaCoeffType *A)
{
    IritFree((VoidPtr) A -> Matrix);
    IritFree((VoidPtr) A -> Rows);
    IritFree((VoidPtr) A -> ColIndex);
    IritFree((VoidPtr) A -> ColLength);
    IritFree((VoidPtr) A);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to     M
* form the new knot vector KVt.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   k:           Order of geometry.                                          M
*   KVT:         Original knot vector.                                       M
*   LengthKVT:   Length of original knot vector.                             M
*   NewKV:       A sequence of new knots to introduce into KVT.              M
*   LengthNewKV: Length of new knot sequence.                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   BspKnotAlphaCoeffType *:   A matrix to multiply the coefficients of the  M
*                              geometry using KVT, in order to get the       M
*                              coefficients under the space defined using    M
*                              KVt that represent the same geometry.         M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotEvalAlphaCoefMerge, alpha matrix, refinement                      M
*****************************************************************************/
BspKnotAlphaCoeffType *BspKnotEvalAlphaCoefMerge(int k,
						 CagdRType *KVT,
						 int LengthKVT,
						 CagdRType *NewKV,
						 int LengthNewKV)
{
    BspKnotAlphaCoeffType *A;

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

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

	IritFree(KVt);
    }

    return A;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Prepares a refinement vector for the given knot vector domain with n	     M
* inserted knots equally spaced.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   n:       Number of knots to introduce.                                   M
*   Tmin:    Minimum domain to introduce knots.                              M
*   Tmax:    Maximum domain to introduce knots.                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *: A vector of n knots uniformly spaced between TMin and TMax. M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotPrepEquallySpaced, refinement                                     M
*****************************************************************************/
CagdRType *BspKnotPrepEquallySpaced(int n, CagdRType Tmin, CagdRType Tmax)
{
    int i;
    CagdRType dt, t, *RefKV;

    if (n <= 0) {
	CAGD_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_OSLO
/*****************************************************************************
* DESCRIPTION:                                                               *
* Prints the content of the alpha matrix.				     *
*                                                                            *
* PARAMETERS:                                                                *
*   A:        Alpha matrix to print.                                         *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
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_OSLO */
