/******************************************************************************
* Bsp_knot.c - Various bspline routines to handle knot vectors.		      *
*******************************************************************************
* Written by Gershon Elber, Aug. 90.					      *
******************************************************************************/

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

#define KNOT_IS_THE_SAME 1e-10

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given curve has no interior knot open end KV.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To check for KV that mimics Bezier polynomial curve.         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:  TRUE if same as Bezier curve, FALSE otherwise.               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvHasBezierKV, conversion                                            M
*****************************************************************************/
CagdBType BspCrvHasBezierKV(CagdCrvStruct *Crv)
{
    return BspKnotHasBezierKV(Crv -> KnotVector, Crv-> Length, Crv -> Order);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given surface has no interior knot open end KVs.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       To check for KVs that mimics Bezier polynomial surface.       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:  TRUE if same as Bezier curve, FALSE otherwise.               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfHasBezierKVs, conversion                                           M
*****************************************************************************/
CagdBType BspSrfHasBezierKVs(CagdSrfStruct *Srf)
{
    return
	BspKnotHasBezierKV(Srf -> UKnotVector, Srf-> ULength, Srf -> UOrder) &&
	BspKnotHasBezierKV(Srf -> VKnotVector, Srf-> VLength, Srf -> VOrder);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given knot vector has no interior knots and it has    M
* an open end cond.                                                          M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:   To check for open end and no interior knots conditions.    M
*   Len:          Of knot vector.                                            M
*   Order:        Of curve/surface the exploits this knot vector.            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:    TRUE if has open end conditions and no interior knots,     M
*                 FALSE otherwise.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotHasBezierKV, knot vectors, conversion                             M
*****************************************************************************/
CagdBType BspKnotHasBezierKV(CagdRType *KnotVector, int Len, int Order)
{
    return Len == Order * 2 && BspKnotHasOpenEC(KnotVector, Len, Order);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given Bspline curve has open end coditions.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      To check for open end conditions.                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:  TRUE, if curve has open end conditions, FALSE otherwise.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvHasOpenEC, open end conditions                                     M
*****************************************************************************/
CagdBType BspCrvHasOpenEC(CagdCrvStruct *Crv)
{
    return BspKnotHasOpenEC(Crv -> KnotVector, Crv -> Length, Crv -> Order);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given Bspline surface has open end coditions.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      To check for open end conditions.                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:  TRUE, if surface has open end conditions, FALSE otherwise.   M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfHasOpenEC, open end conditions                                     M
*****************************************************************************/
CagdBType BspSrfHasOpenEC(CagdSrfStruct *Srf)
{
    return
	BspKnotHasOpenEC(Srf -> UKnotVector, Srf -> ULength, Srf -> UOrder) &&
	BspKnotHasOpenEC(Srf -> VKnotVector, Srf -> VLength, Srf -> VOrder);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given Bspline surface has open end coditions in the   M
* specified direction.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      To check for open end conditions.                              M
*   Dir:      Either the U or the V parametric direction.                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:  TRUE, if surface has open end conditions, FALSE otherwise.   M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfHasOpenECDir, open end conditions                                  M
*****************************************************************************/
CagdBType BspSrfHasOpenECDir(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    return BspKnotHasOpenEC(Srf -> UKnotVector, Srf -> ULength,
				    Srf -> UOrder);
	case CAGD_CONST_V_DIR:
	    return BspKnotHasOpenEC(Srf -> VKnotVector, Srf -> VLength,
				    Srf -> VOrder);
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    return FALSE;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the given knot vector has open end conditions.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:   To check for open end condition.                           M
*   Len:          Of knot vector.                                            M
*   Order:        Of curve/surface the exploits this knot vector.            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:     TRUE if KV has open end conditions.                       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotHasOpenEC, knot vectors, open end conditions                      M
*****************************************************************************/
CagdBType BspKnotHasOpenEC(CagdRType *KnotVector, int Len, int Order)
{
    int i,
        LastIndex = Len + Order - 1;

    for (i = 1; i < Order; i++)
       if (!APX_EQ_EPS(KnotVector[0], KnotVector[i], IRIT_EPSILON))
	    return FALSE;

    for (i = LastIndex - 1; i >= Len; i--)
       if (!APX_EQ_EPS(KnotVector[LastIndex], KnotVector[i], IRIT_EPSILON))
	    return FALSE;

    return TRUE;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff t is in the parametric domain as define by the knot       M
* vector KnotVector, its length Len, and the order Order.	             M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To verify t is indeed in.                                 M
*   Len:           Length of curve/surface using KnotVector. This is NOT the M
*                  length of KnotVector which is equal to (Length + Order).  M
*   Order:         Order of curve/surface using KnotVector.                  M
*   Periodic:      TRUE if this KnotVector is periodic.                      M
*   t:             Parametric value to verify.                               M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:    TRUE, if t is contained in the parametric domain, FALSE    M
*                 otherwise.                                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotParamInDomain, parametric domain, knot vectors                    M
*****************************************************************************/
CagdBType BspKnotParamInDomain(CagdRType *KnotVector,
			       int Len,
			       int Order,
			       CagdBType Periodic,
			       CagdRType t)
{
    CagdRType
	r1 = KnotVector[Order - 1],
	r2 = KnotVector[Len + (Periodic ? Order - 1 : 0)];

    return (r1 < t || APX_EQ_EPS(r1, t, IRIT_EPSILON))
        && (r2 > t || APX_EQ_EPS(r2, t, IRIT_EPSILON));
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns the index of the last knot which is less than or equal to t in the M
* given knot vector KnotVector of length Len. APX_EQ_EPS is used in equality.M
*   Parameter t is assumed to be in the parametric domain for the knot       M
* vector.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To search for a knot with the LE relation to t.           M
*   Len:           Of knot vector. This is not the length of the             M
*                  curve/surface using this KnotVector.                      M
*   t:             The parameter value to search for the LE relation.        M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:           Index of last knot in KnotVector that is LE t, or -1 if   M
*                  t is below the first knot.				     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotLastIndexLE, knot vectors                                         M
*****************************************************************************/
int BspKnotLastIndexLE(CagdRType *KnotVector, int Len, CagdRType t)
{
    int i;

    for (i = 0;
	 i < Len && (KnotVector[i] <= t ||
		     APX_EQ_EPS(KnotVector[i], t, IRIT_EPSILON));
	 i++);

    return i - 1;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns the index of the last knot which is less t in the given knot       M
* vector KnotVector of length Len. APX_EQ_EPS is used for equality.          M
*   Parameter t is assumed to be in the parametric domain for the knot       M
* vector.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To search for a knot with the L relation to t.            M
*   Len:           Of knot vector. This is not the length of the             M
*                  curve/surface using this KnotVector.                      M
*   t:             The parameter value to search for the L relation.         M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:           Index of last knot in KnotVector that is less than t or   M
*                  -1 if t is below the first knot.			     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotLastIndexL, knot vectors                                          M
*****************************************************************************/
int BspKnotLastIndexL(CagdRType *KnotVector, int Len, CagdRType t)
{
    int i;

    for (i = 0;
	 i < Len &&
	 KnotVector[i] < t &&
	 !APX_EQ_EPS(KnotVector[i], t, IRIT_EPSILON);
	 i++);

    return i - 1;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns the index of the first knot which is greater than t in the given   M
* knot vector KnotVector of length Len. APX_EQ_EPS is used for equality.     M
*   Parameter t is assumed to be in the parametric domain for the knot       M
* vector.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To search for a knot with the G relation to t.            M
*   Len:           Of knot vector. This is not the length of the             M
*                  curve/surface using this KnotVector.                      M
*   t:             The parameter value to search for the G relation.         M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:           Index of first knot in KnotVector that is greater than t  M
*                  or Len if t is greater than last knot in KnotVector.      M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotFirstIndexG, knot vectors                                         M
*****************************************************************************/
int BspKnotFirstIndexG(CagdRType *KnotVector, int Len, CagdRType t)
{
    int i;

    for (i = Len - 1;
	 i >= 0 &&
	 KnotVector[i] > t &&
	 !APX_EQ_EPS(KnotVector[i], t, IRIT_EPSILON);
	 i--);

    return i + 1;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a uniform periodic knot vector for Len Control points and order    M
* Order. The actual length of the KV is Len + Order + Order - 1.	     M
* If KnotVector is NULL it is being allocated dynamically.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Len:          Of control polygon/mesh of curve/surface that are to use   M
*                 this knot vector.                                          M
*   Order:        Of the curve/surface that are to use this knot vector.     M
*   KnotVector:   If new knot vector is to be saved here, otherwise a new    M
*                 space is allocated.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  The created uniform periodic knot vector, either newly     M
*                 allocated or given in Knotvector and just initialized.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotUniformPeriodic, knot vectors, periodic end conditions, end       M
*   conditions 								     M
*****************************************************************************/
CagdRType *BspKnotUniformPeriodic(int Len, int Order, CagdRType *KnotVector)
{
    int i,
	KVLen = Len + Order + Order - 1;
    CagdRType *KV;

    if (KnotVector == NULL)
	KV = KnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * KVLen);
    else
	KV = KnotVector;

    for (i = 0; i < KVLen; i++)
	*KV++ = (CagdRType) i;

    return KnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a uniform floating knot vector for Len Control points and order    M
* Order. The actual length of the KV is Len + Order.			     M
* If KnotVector is NULL it is being allocated dynamically.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Len:          Of control polygon/mesh of curve/surface that are to use   M
*                 this knot vector.                                          M
*   Order:        Of the curve/surface that are to use this knot vector.     M
*   KnotVector:   If new knot vector is to be saved here, otherwise a new    M
*                 space is allocated.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  The created uniform floating knot vector, either newly     M
*                 allocated or given in Knotvector and just initialized.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotUniformFloat, knot vectors, floating end conditions, end          M
*   conditions								     M
*****************************************************************************/
CagdRType *BspKnotUniformFloat(int Len, int Order, CagdRType *KnotVector)
{
    int i;
    CagdRType *KV;

    if (KnotVector == NULL)
	KV = KnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) *
								(Len + Order));
    else
	KV = KnotVector;

    for (i = 0; i < Len + Order; i++)
	*KV++ = (CagdRType) i;

    return KnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a uniform open knot vector for Len Control points and order        M
* Order. The actual length of the KV is Len + Order.			     M
* If KnotVector is NULL it is being allocated dynamically.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Len:          Of control polygon/mesh of curve/surface that are to use   M
*                 this knot vector.                                          M
*   Order:        Of the curve/surface that are to use this knot vector.     M
*   KnotVector:   If new knot vector is to be saved here, otherwise a new    M
*                 space is allocated.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  The created uniform open knot vector, either newly         M
*                 allocated or given in Knotvector and just initialized.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotUniformOpen, knot vectors, open end conditions, end               M
*   conditions 								     M
*****************************************************************************/
CagdRType *BspKnotUniformOpen(int Len, int Order, CagdRType *KnotVector)
{
    int i, j;
    CagdRType *KV;

    if (KnotVector == NULL)
	KV = KnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) *
								(Len + Order));
    else
	KV = KnotVector;

    for (i = 0; i < Order; i++)
	*KV++ = 0.0;
    for (i = 1, j = Len - Order; i <= j; )
	*KV++ = (CagdRType) i++;
    for (j = 0; j < Order; j++)
	*KV++ = (CagdRType) i;

    return KnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Applies an affine transformation to the given knot vector. Note affine     M
* transformation on the knot vector does not change the Bspline curve.	     M
*   Knot vector is translated by Translate amount and scaled by Scale as     M
*  KV[i] = (KV[i] - KV[0]) * Scale + (KV0 + Translate).                      V
*   All transformation as taken place in place.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:     To affinely transform.                                   M
*   Len:            Of knot vector. This is not the length of the curve or   M
*                   surface using this knot vector.                          M
*   Translate:      Amount to translate the knot vector.                     M
*   Scale:          Amount to scale the knot vector.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotAffineTrans, knot vectors, affine transformation                  M
*****************************************************************************/
void BspKnotAffineTrans(CagdRType *KnotVector,
			int Len,
			CagdRType Translate,
			CagdRType Scale)
{
    int i;
    CagdRType
	KV0 = KnotVector[0];

    for (i = 0; i < Len; i++)
	KnotVector[i] = (KnotVector[i] - KV0) * Scale + KV0 + Translate;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Applies an affine transformation to the given knot vector. Note affine     M
* transformation on the knot vector does not change the Bspline curve.	     M
*   Knot vector is translated and scaled so as to span the domain from       M
* MinVal top MaxVal. This works for open end condition curves only.	     M
*  KV[i] = (KV[i] - KV[0]) * Scale + MinVal,                                 V
* where Scale = (MaxVal - MinVal) / (KV[Len - 1] - KV[0]).		     M
*   All transformation as taken place in place.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:     To affinely transform.                                   M
*   Len:            Of knot vector. This is not the length of the curve or   M
*                   surface using this knot vector.                          M
*   MinVal, MaxVal: New parametric domain of knot vector.                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotAffineTrans2, knot vectors, affine transformation                 M
*****************************************************************************/
void BspKnotAffineTrans2(CagdRType *KnotVector,
			 int Len,
			 CagdRType MinVal,
			 CagdRType MaxVal)
{
    int i;
    CagdRType
	KV0 = KnotVector[0],
	KVn = KnotVector[Len - 1],
	Scale = (MaxVal - MinVal) / (KVn - KV0);

    for (i = 0; i < Len; i++)
	KnotVector[i] = (KnotVector[i] - KV0) * Scale + MinVal;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Creates an identical copy of a given knot vector KnotVector of length Len. M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:   Knot vector to duplicate                                   M
*   Len:          Length of knot vector. This is not the length of the curve M
*                 or surface using this knot vector.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  The duplicated knot vector.                                M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotCopy, allocation, knot vectors                                    M
*****************************************************************************/
CagdRType *BspKnotCopy(CagdRType *KnotVector, int Len)
{
    CagdRType
	*NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len);

    CAGD_GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * Len);

    return NewKnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Reverse a knot vector of length Len.	Reversing of knot vector keeps the   M
* knots monotonically non-decreasing as well as the parametric domain. Only  M
* the spaces between the knots are being flipped. For example the knot       M
* vector:                                                                    M
*                        [0 0 0 0 1 2 2 6 6 6 6]                             V
* is reversed to be:                                                         M
*                        [0 0 0 0 4 4 5 6 6 6 6]                             V
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:   Knot vector to be reversed.                                M
*   Len:          Length of knot vector. This is not the length of the curve M
*                 or surface using this knot vector.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  The reversed knot vector.                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotReverse, knot vectors, reverse                                    M
*****************************************************************************/
CagdRType *BspKnotReverse(CagdRType *KnotVector, int Len)
{
    int i;
    CagdRType
	*NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len),
	t = KnotVector[Len - 1] + KnotVector[0];

    for (i = 0; i < Len; i++)
	NewKnotVector[i] = t - KnotVector[Len - i - 1];

    return NewKnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a knot vector that contains all the knots in KnotVector1 that are  M
* not in KnotVector2.							     M
*   NewLen is set to new KnotVector length.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector1:   First  knot vector.                                       M
*   Len1:          Length of KnotVector1. This is not the length of the      M
*                  curve or surface using this knot vector.                  M
*   KnotVector2:   Second knot vector.                                       M
*   Len1:          Length of KnotVector2. This is not the length of the      M
*                  curve or surface using this knot vector.                  M
*   NewLen:        To save the size of the knot vector that contains the     M
*                  computed subset of KnotVector1 / KnotVector2.             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   The subset of knot in KnotVector1 that are not in         M
*                  KnotVector2 (KnotVector1 / KnotVector2).                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotSubtrTwo, knot vectors, compatibility, refinement                 M
*****************************************************************************/
CagdRType *BspKnotSubtrTwo(CagdRType *KnotVector1,
			   int Len1,
			   CagdRType *KnotVector2,
			   int Len2,
			   int *NewLen)
{
    int i = 0,
	j = 0;
    CagdRType
	*NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len1),
	*t = NewKnotVector;

    *NewLen = 0;
    while (i < Len1 && j < Len2) {
	if (APX_EQ_EPS(KnotVector1[i], KnotVector2[j], IRIT_EPSILON)) {
	    i++;
	    j++;
	}
	else if (KnotVector1[i] > KnotVector2[j]) {
	    j++;
	}
	else {
	    /* Current KnotVector1 value is less than current KnotVector2: */
	    *t++ = KnotVector1[i++];
	    (*NewLen)++;
	}
    }

    return NewKnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Merges two knot vector KnotVector1 and KnotVector2 of length Len1 and Len2 M
* respectively into one. If Mult is not zero then knot multiplicity is       M
* tested not to be larger than Mult value.				     M
*   NewLen is set to new KnotVector length.                                  M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector1:   First  knot vector.                                       M
*   Len1:          Length of KnotVector1. This is not the length of the      M
*                  curve or surface using this knot vector.                  M
*   KnotVector2:   Second knot vector.                                       M
*   Len2:          Length of KnotVector2. This is not the length of the      M
*                  curve or surface using this knot vector.                  M
*   Mult:          Maximum multiplicity to allow in merged knot vector.      M
*   NewLen:        To save the size of the knot vector that contains the     M
*                  merged knot vectors.				             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   The merged knot verctor (KnotVector1 U KnotVector2).      M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotMergeTwo, knot vectors, compatibility, refinement                 M
*****************************************************************************/
CagdRType *BspKnotMergeTwo(CagdRType *KnotVector1,
			   int Len1,
			   CagdRType *KnotVector2,
			   int Len2,
			   int Mult,
			   int *NewLen)
{
    int i = 0,
	j = 0,
	k = 0,
        m = 0,
	Len = Len1 + Len2;
    CagdRType t,
	*NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len);

    if (Mult == 0)
	Mult = Len + 1;

    while (i < Len1 && j < Len2) {
	if (KnotVector1[i] < KnotVector2[j]) {
	    /* Use value from KnotVector1: */
	    t = KnotVector1[i++];
	}
	else {
	    /* Use value from KnotVector2: */
	    t = KnotVector2[j++];
	}

	if (k == 0 || (k > 0 &&
		       !APX_EQ_EPS(NewKnotVector[k - 1], t, IRIT_EPSILON))) {
	    NewKnotVector[k++] = t;
	    m = 1;
	}
	else if (m < Mult) {
	    NewKnotVector[k++] = t;
	    m++;
	}
    }

    while (i < Len1)
	NewKnotVector[k++] = KnotVector1[i++];
    while (j < Len2)
	NewKnotVector[k++] = KnotVector1[j++];

    /* It should be noted that k <= Len so there is a chance some of the new */
    /* KnotVector space will not be used (it is not memory leak!). If you    */
    /* really care that much - copy it to a new allocated vector of size k   */
    /* and return it while freeing the original of size Len.		     */
    *NewLen = k;

    return NewKnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Merges two knot vector KnotVector1 and KnotVector2 of length Len1 and Len2 M
* respectively into one, from geometries of orders Order1 and Order2.	     M
*   Merged knot vector is for order ResOrder so that the resulting curve can M
* represent the discontinuities in both geometries.		             M
*   Assumes both knot vectors are open end spanning the same domain.         M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector1:   First  knot vector.                                       M
*   Len1:          Length of KnotVector1. This is not the length of the      M
*                  curve or surface using this knot vector.                  M
*   Order1:        Order of first knot vector's geometry.                    M
*   KnotVector2:   Second knot vector.                                       M
*   Len2:          Length of KnotVector2. This is not the length of the      M
*                  curve or surface using this knot vector.                  M
*   Order2:        Order of secon knot vector's geometry.                    M
*   ResOrder:      Expected order of geometry that will use the merged       M
*                  knot vector.                                              M
*   NewLen:        To save the size of the knot vector that contains the     M
*                  merged knot vectors.				             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   The merged knot verctor (KnotVector1 U KnotVector2).      M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotContinuityMergeTwo, knot vectors, compatibility, refinement       M
*****************************************************************************/
CagdRType *BspKnotContinuityMergeTwo(CagdRType *KnotVector1,
				     int Len1,
				     int Order1,
				     CagdRType *KnotVector2,
				     int Len2,
				     int Order2,
				     int ResOrder,
				     int *NewLen)
{
    int l1, l2, m, cont,
	i = 0,
	j = 0,
	k = 0,
	Len = (Len1 + Len2) * (1 + ResOrder - MAX(Order1, Order2));
    CagdRType t,
	*NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len);

    while (i < Len1 && j < Len2) {
	if (APX_EQ_EPS(KnotVector1[i], KnotVector2[j], IRIT_EPSILON)) {
	    /* Compute multiplicity and continuity. */
	    for (l1 = 1;
		 i + l1 < Len1 &&
		 APX_EQ_EPS(KnotVector1[i], KnotVector1[i + l1], IRIT_EPSILON);
		 l1++);
	    for (l2 = 1;
		 j + l2 < Len2 &&
		 APX_EQ_EPS(KnotVector2[j], KnotVector2[j + l2], IRIT_EPSILON);
		 l2++);
	    cont = MIN(Order1 - l1, Order2 - l2) - 1;

	    /* Use value from KnotVector1: */
	    t = KnotVector1[i];
	    i += l1;
	    j += l2;
	}
	else if (KnotVector1[i] < KnotVector2[j]) {
	    /* Compute multiplicity and continuity. */
	    for (l1 = 1;
		 i + l1 < Len1 &&
		 APX_EQ_EPS(KnotVector1[i], KnotVector1[i + l1], IRIT_EPSILON);
		 l1++);
	    cont = Order1 - l1 - 1;

	    /* Use value from KnotVector1: */
	    t = KnotVector1[i];
	    i += l1;
	}
	else {
	    /* Compute multiplicity and continuity. */
	    for (l2 = 1;
		 j + l2 < Len2 &&
		 APX_EQ_EPS(KnotVector2[j], KnotVector2[j + l2], IRIT_EPSILON);
		 l2++);
	    cont = Order2 - l2 - 1;

	    /* Use value from KnotVector2: */
	    t = KnotVector2[j++];
	    j += l2;
	}

	for (m = 0; m < ResOrder - cont - 1; m++)
	    NewKnotVector[k++] = t;
    }

    /* It should be noted that k <= Len so there is a chance some of the new */
    /* KnotVector space will not be used (it is not memory leak!). If you    */
    /* really care that much - copy it to a new allocated vector of size k   */
    /* and return it while freeing the original of size Len.		     */
    *NewLen = k;

    return NewKnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Creates a new knot vector from the given KnotVector that averages Ave	     M
* consequetive knots.							     M
*  Resulting vector will have (Len - Ave + 1) elements.			     M
*  See also BspKnotNodes.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To average out.                                           M
*   Len:           Length of KnotVector. This is not the length of the       M
*                  curve or surface using this knot vector.                  M
*   Ave:           How many knots to average each time.                      M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   The averaged knot vector of length (Len - Ave + 1).       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotAverage, knot vectors, node values                                M
*****************************************************************************/
CagdRType *BspKnotAverage(CagdRType *KnotVector, int Len, int Ave)
{
    int i,
	AveLen = Len - Ave + 1;
    CagdRType
	*AveVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * AveLen);

    if (Ave > Len || Ave < 1)
	CAGD_FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);

    /* Compute the first average value. */
    for (AveVector[0] = 0.0, i = 0; i < Ave; i++)
	AveVector[0] += KnotVector[i];

    for (i = 1; i < AveLen; i++)
	AveVector[i] = AveVector[i - 1] + KnotVector[i + Ave - 1]
					- KnotVector[i - 1];

    for (i = 0; i < AveLen; i++)
	AveVector[i] /= Ave;

    return AveVector;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Creates a new vector with the given KnotVector Node values.		     M
*   The nodes are the approximated parametric value associated with the each M
* control point. Therefore for a knot vector of length Len and order Order   M
* there are Len - Order control points and therefore nodes.		     M
*   Nodes are defined as (k = Order - 1 or degree):			     M
*	    								     M
*	  i+k								     V
*	   -			 	First Node N(i = 0)		     V
*	   \								     V
*          / KnotVector(j)		Last Node  N(i = Len - k - 2)	     V
*	   -								     V
*        j=i+1								     V
* N(i) = -----------------						     V
*	        k							     V
*	    								     M
*  See also BspKnotAverage.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To average out as nodes.                                  M
*   Len:           Length of KnotVector. This is not the length of the       M
*                  curve or surface using this knot vector.                  M
*   Order:         Of curve or surface that exploits this knot vector.       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   The nodes computed for the given knot vector.             M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotNodes, knot vectors, node values                                  M
*****************************************************************************/
CagdRType *BspKnotNodes(CagdRType *KnotVector, int Len, int Order)
{
    int i, j,
	k = MAX(Order - 1, 1),
	NodeLen = Len - Order;
    CagdRType
	TMin = KnotVector[k],
	TMax = KnotVector[NodeLen],
	*NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * NodeLen);

    for (i = 0; i < NodeLen; i++) {
	NodeVector[i] = 0.0;
	for (j = i + 1; j <= i + k; j++)
	    NodeVector[i] += BOUND(KnotVector[j], TMin, TMax);
	NodeVector[i] /= k;
    }

    return NodeVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns the nodes of a freeform curve.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      To compute node values for.                                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *: Node values of the given curve.                             M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdCrvNodes, node values                                                M
*****************************************************************************/
CagdRType *CagdCrvNodes(CagdCrvStruct *Crv)
{
    int i,
        Order = Crv -> Order,
	Length = CAGD_CRV_PT_LST_LEN(Crv);
    CagdRType *NodeVector;

    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	    NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Order);
	    for (i = 0; i < Order; i++)
	        NodeVector[i] = i / ((RealType) (Order - 1));
	    break;
	case CAGD_CBSPLINE_TYPE:
	    NodeVector = BspKnotNodes(Crv -> KnotVector, Length + Order, Order);
	    break;
	default:
	    NodeVector = NULL;
	    break;
    }

    return NodeVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns the nodes of a freeform surface.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      To compute node values for.                                    M
*   Dir:      Either the U or the V parametric direction.                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *: Node values of the given surface and given parametric       M
*                direction.                                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdSrfNodes, node values                                                M
*****************************************************************************/
CagdRType *CagdSrfNodes(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    int i,
        Order = Dir == CAGD_CONST_U_DIR ? Srf -> UOrder : Srf -> VOrder,
	Length = Dir == CAGD_CONST_U_DIR ? Srf -> ULength : Srf -> VLength;
    CagdRType *NodeVector;

    if (Dir != CAGD_CONST_U_DIR && Dir != CAGD_CONST_V_DIR) {
        CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	return NULL;
    }

    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	    NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Order);
	    for (i = 0; i < Order; i++)
	        NodeVector[i] = i / ((RealType) (Order - 1));
	    break;
	case CAGD_SBSPLINE_TYPE:
	    NodeVector = BspKnotNodes(Dir == CAGD_CONST_U_DIR ?
				          Srf -> UKnotVector :
				          Srf -> VKnotVector,
				      Length + Order, Order);
	    break;
	default:
	    NodeVector = NULL;
	    break;
    }

    return NodeVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Finds the parameter value with the largest coefficient of the curve using  M
* nodes values to estimate the coefficients' parameters.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:           To compute the parameter node value of the largest        M
*                  coefficient.                                              M
*   Axis:          Which axis should we search for maximal coefficient?      M
*                  1 for X, 2 for Y, etc.                                    M
*   MaxVal:        The coefficient itself will be place herein.              M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType:     The node parameter value of the detected maximal          M
*                  coefficient.						     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvMaxCoefParam, extremum                                             M
*****************************************************************************/
CagdRType BspCrvMaxCoefParam(CagdCrvStruct *Crv, int Axis, CagdRType *MaxVal)
{
    int i,
	Index = 0,
	Length = Crv -> Length,
	Order = Crv -> Order;
    CagdRType
	*Points = Crv -> Points[Axis],
	R = *Points,
        *Nodes = BspKnotNodes(Crv -> KnotVector, Length + Order, Order);

    for (i = 0; i < Length; i++, Points++) {
	if (*Points > R) {
	    R = *Points;
	    Index = i;
	}
    }
    *MaxVal = R;

    R = Nodes[Index];
    IritFree(Nodes);

    return R;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Finds the parameter value with the largest coefficient of the surface      M
* using nodes values to estimate the coefficients' parameters.		     M
* Returns a pointer to a static array of two elements holding U and V.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:           To compute the parameter node value of the largest        M
*                  coefficient.                                              M
*   Axis:          Which axis should we search for maximal coefficient?      M
*                  1 for X, 2 for Y, etc.                                    M
*   MaxVal:        The coefficient itself will be place herein.              M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   The node UV parameter values of the detected maximal      M
*                  coefficient.						     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfMaxCoefParam, extremum                                             M
*****************************************************************************/
CagdRType *BspSrfMaxCoefParam(CagdSrfStruct *Srf, int Axis, CagdRType *MaxVal)
{
    static CagdRType UV[2];
    int i,
	UIndex = 0,
	VIndex = 0,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder;
    CagdRType
	*Points = Srf -> Points[Axis],
	R = *Points,
        *UNodes = BspKnotNodes(Srf -> UKnotVector, ULength + UOrder, UOrder),
        *VNodes = BspKnotNodes(Srf -> VKnotVector, VLength + VOrder, VOrder);

    for (i = 0; i < ULength * VLength; i++, Points++) {
	if (*Points > R) {
	    R = *Points;
	    UIndex = i % ULength;
	    VIndex = i / ULength;
	}
    }
    *MaxVal = R;

    UV[0] = UNodes[UIndex];
    UV[1] = VNodes[VIndex];
    IritFree(UNodes);
    IritFree(VNodes);
    return UV;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Creates a new vector with t inserted as a new knot. Attempt is made to   M
* make sure t is in the knot vector domain.				     M
*   No test is made for the current multiplicity of knot t in KnotVector.    M
*   This function only constructs a refined knot vector and does not         M
* compute the actual refined coefficients.                                   M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To insert new knot t in.                                  M
*   Order:         Of geometry that exploits KnotVector.                     M
*   Len:           Length of curve/surface using KnotVector. This is NOT the M
*                  length of KnotVector which is equal to (Length + Order).  M
*   t:             The new knot t to insert.                                 M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   A new knot vector larger by one than KnotVector that      M
*                  contains t.                                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotInsertOne, knot vectors, knot insertion, refinement               M
*****************************************************************************/
CagdRType *BspKnotInsertOne(CagdRType *KnotVector,
			    int Order,
			    int Len,
			    CagdRType t)
{
    int Mult = BspKnotFindMult(KnotVector, Order, Len, t) + 1;

    return BspKnotInsertMult(KnotVector, Order, &Len, t, Mult);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Inserts Mult knots with value t into the knot vector KnotVector.	     M
*   Attempt is made to make sure t in knot vector domain.		     M
*   If a knot equal to t (up to APX_EQ) already exists with multiplicity i   M
* only (Mult - i) knot are being inserted into the new knot vector.	     M
*   Len is updated to the resulting knot vector.			     M
*   It is possible to DELETE a knot using this routine by specifying         M
* multiplicity less then current multiplicity!				     M
*   This function only constructs a refined knot vector and does not         M
* compute the actual refined coefficients.                                   M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To insert new knot t in.                                  M
*   Order:         Of geometry that exploits KnotVector.                     M
*   Len:           Length of curve/surface using KnotVector. This is NOT the M
*                  length of KnotVector which is equal to (Length + Order).  M
*   t:             The new knot t to insert.                                 M
*   Mult:          The multiplicity that this knot should have in resulting  M
*                  knot vector.						     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   A new knot vector derived from KnotVector that has        M
*                  a mltiplicity of exacly Mult at the knot t.               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotInsertMult, knot vectors, knot insertion, refinement              M
*****************************************************************************/
CagdRType *BspKnotInsertMult(CagdRType *KnotVector,
			     int Order,
			     int *Len,
			     CagdRType t,
			     int Mult)
{
    int i, CurrentMult, NewLen, FirstIndex;
    CagdRType *NewKnotVector;

    if (!BspKnotParamInDomain(KnotVector, *Len, Order, FALSE, t))
	CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

    CurrentMult = BspKnotFindMult(KnotVector, Order, *Len, t);
    NewLen = *Len + Mult - CurrentMult;
    NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) *
    							(NewLen + Order));
    FirstIndex = BspKnotLastIndexL(KnotVector, *Len + Order, t) + 1;

    /* Copy all the knot before the knot t. */
    CAGD_GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * FirstIndex);

    /* Insert Mult knot of value t. */
    for (i = FirstIndex; i < FirstIndex + Mult; i++)
	NewKnotVector[i] = t;

    /* And copy the second part. */
    CAGD_GEN_COPY(&NewKnotVector[FirstIndex + Mult],
		  &KnotVector[FirstIndex + CurrentMult],
		  sizeof(CagdRType) * (*Len + Order - FirstIndex - Mult + 1));

    *Len = NewLen;
    return NewKnotVector;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns the multiplicity of knot t in knot vector KnotVector, zero if      M
* none.                                                                      M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To test multiplicity of knot value t at.                  M
*   Order:         Of geometry that exploits KnotVector.                     M
*   Len:           Length of curve/surface using KnotVector. This is NOT the M
*                  length of KnotVector which is equal to (Length + Order).  M
*   t:             The knot to verify the multiplicity of.                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:           Multiplicity of t in KnotVector.                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotFindMult, knot vectors                                            M
*****************************************************************************/
int BspKnotFindMult(CagdRType *KnotVector, int Order, int Len, CagdRType t)
{
    int LastIndex, FirstIndex;

    if (!BspKnotParamInDomain(KnotVector, Len, Order, FALSE, t))
	CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

    FirstIndex = BspKnotLastIndexL(KnotVector, Len, t) + 1;

    for (LastIndex = FirstIndex;
	 LastIndex < Len && APX_EQ_EPS(KnotVector[LastIndex], t, IRIT_EPSILON);
	 LastIndex++);

    return LastIndex - FirstIndex;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Scans the given knot vector to a potential C1 discontinuity.		     M
*   Returns TRUE if found one and set t to its parameter value.		     M
*   Assumes knot vector has open end condition.				     M
*   A knot vector with multiplicity of a knot of (Order - 1) can be C1       M
* discontinuous at that knot. However, this is only a necessary condition    M
* for C1 discontinuity in the geometry.                                      M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To test for potential C1 dicontinuities.                  M
*   Order:         Of geometry that exploits KnotVector.                     M
*   Length:        Length of curve/surface using KnotVector. This is NOT the M
*                  length of KnotVector which is equal to (Length + Order).  M
*   t:             Where to put the parameter value (knot) that can be C2    M
*                  discontinuous.                                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:     TRUE if found a potential C1 discontinuity, FALSE         M
*		   otherwise.						     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotC1Discont, knot vectors, continuity, discontinuity                M
*****************************************************************************/
CagdBType BspKnotC1Discont(CagdRType *KnotVector,
			   int Order,
			   int Length,
			   CagdRType *t)
{
    int i, Count;
    CagdRType
	LastT = KnotVector[0] - 1.0;

    for (i = Order, Count = 0; i < Length; i++) {
	if (APX_EQ_EPS(LastT, KnotVector[i], IRIT_EPSILON))
	    Count++;
	else {
	    Count = 1;
	    LastT = KnotVector[i];
	}

	if (Count >= Order - 1) {
	    *t = LastT;
	    return TRUE;
	}
    }

    return FALSE;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Scans the given knot vector for all potential C1 discontinuity. Returns    M
* a vector holding the parameter values of the potential C1 discontinuities, M
* NULL of none found.							     M
*   Sets n to length of returned vector.				     M
*   Assumes knot vector has open end condition.				     M
*   A knot vector with multiplicity of a knot of (Order - 1) can be C1       M
* discontinuous at that knot. However, this is only a necessary condition    M
* for C1 discontinuity in the geometry.                                      M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To test for potential C1 dicontinuities.                  M
*   Order:         Of geometry that exploits KnotVector.                     M
*   Length:        Length of curve/surface using KnotVector. This is NOT the M
*                  length of KnotVector which is equal to (Length + Order).  M
*   n:             Length of returned vector - number of potential C1        M
*                  discontinuities found.                                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   Vector holding all parametr values with potential C1      M
*                  discontinuities.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotAllC1Discont, knot vectors, continuity, discontinuity             M
*****************************************************************************/
CagdRType *BspKnotAllC1Discont(CagdRType *KnotVector,
			       int Order,
			       int Length,
			       int *n)
{
    int i, Count,
	C1DiscontCount = 0;
    CagdRType
	LastT = KnotVector[0] - 1.0,
	*C1Disconts = (CagdRType *) IritMalloc(sizeof(CagdRType) * Length);

    for (i = Order, Count = 0; i < Length; i++) {
	if (APX_EQ_EPS(LastT, KnotVector[i], IRIT_EPSILON))
	    Count++;
	else {
	    Count = 1;
	    LastT = KnotVector[i];
	}

	if (Count >= Order - 1) {
	    C1Disconts[C1DiscontCount++] = LastT;
	    Count = 0;
	}
    }

    if ((*n = C1DiscontCount) > 0) {
	return C1Disconts;
    }
    else {
	IritFree((VoidPtr) C1Disconts);
	return NULL;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Routine to determine where to sample along the provided paramtric domain, M
* given the C1 discontinuities along it.				     M
*   Returns a vector of length NumSamples.				     M
*   If C1Disconts != NULL (NumC1Disconts > 0), C1Discont is being freed.     M
*                                                                            *
* PARAMETERS:                                                                M
*   PMin:            Minimum of parametric domain.                           M
*   PMax:            Maximum of parametric domain.                           M
*   NumSamples:      To allocate for the vector of samples.                  M
*   C1Disconts:      A vector of potential C1 discontinuities in the         M
*                    (PMin, PMax) domain. This vector is freed by this       M
*                    routine, if it is not NULL.                             M
*   NumC1Disconts:   Length of C1Discont. if zero then C1Discont == NULL.    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:    A vector of the suggested set of sampling locations.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotParamValues, piecewise linear approximation, knot vectors         M
*****************************************************************************/
CagdRType *BspKnotParamValues(CagdRType PMin,
			      CagdRType PMax,
			      int NumSamples,
			      CagdRType *C1Disconts,
			      int NumC1Disconts)
{
    int i, CrntIndex, NextIndex, CrntDiscont;
    CagdRType Step,
	*Samples = (CagdRType *) IritMalloc(sizeof(CagdRType) * NumSamples);

    Samples[0] = PMin;
    if (NumSamples <= 1)
	return Samples;
    Samples[NumSamples - 1] = PMax;
    if (NumSamples <= 2)
	return Samples;

    if (NumC1Disconts > NumSamples - 2) {
    	/* More C1 discontinuities than we are sampling. Grab a subset of    */
	/* The discontinuities as the sampling set.			     */
    	Step = ((CagdRType) (NumC1Disconts - 1)) / (NumSamples - 2)
								- IRIT_EPSILON;
	for (i = 0; i < NumSamples - 2; i++)
    	    Samples[i + 1] = C1Disconts[(int) (i * Step)];
    }
    else {
	/* More samples than C1 discontinuites. Uniformly distribute the C1  */
	/* discontinuites between the samples and linearly interpolate the   */
	/* sample values in between.					     */
    	Step = ((CagdRType) (NumC1Disconts + 1)) / (NumSamples - 2);
	CrntIndex = CrntDiscont = 0;

	while (CrntDiscont < NumC1Disconts) {
	    NextIndex = (int) ((CrntDiscont + 1) / Step);
	    Samples[NextIndex] = C1Disconts[CrntDiscont++];

	    for (i = CrntIndex + 1; i < NextIndex; i++) {
		CagdRType
		    t = (i - CrntIndex) / ((CagdRType) (NextIndex - CrntIndex)),
		    t1 = 1.0 - t;

		Samples[i] = Samples[CrntIndex] * t1 + Samples[NextIndex] * t;
	    }

	    CrntIndex = NextIndex;
	}

	/* Interpolate the last interval: */
    	for (i = CrntIndex + 1; i < NumSamples - 1; i++) {
    	    CagdRType
    		t = (i - CrntIndex) / ((CagdRType) (NumSamples - 1 - CrntIndex)),
    		t1 = 1.0 - t;

    	    Samples[i] = Samples[CrntIndex] * t1 + Samples[NumSamples - 1] * t;
    	}
    }

    if (C1Disconts != NULL)
	IritFree((VoidPtr) C1Disconts);

    return Samples;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a knot vector, make sure adjacent knots that are close "enough" are  M
* actually identical. Important for robustness of subdiv/refinement algs.    M
*                                                                            M
*                                                                            M
*                                                                            *
* PARAMETERS:                                                                M
*   KV:        Knot vector to make robust, in place.                         M
*   Len:       Length of knot vector KV.                                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspKnotMakeRobustKVm knot vectors                                        M
*****************************************************************************/
void BspKnotMakeRobustKV(CagdRType *KV, int Len)
{
    int i;

    for (i = 1; i < Len; i++)
	if (KV[i] - KV[i - 1] < KNOT_IS_THE_SAME && KV[i] != KV[i - 1])
	    KV[i] = KV[i - 1];
}
