/******************************************************************************
* SBsp-Aux.c - Bspline surface auxilary routines.			      *
*******************************************************************************
* Written by Gershon Elber, July. 90.					      *
******************************************************************************/

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

/* Define some marcos to make some of the routines below look better. They  */
/* calculate the index of the U, V point of the control mesh in Points.	    */
#define DERIVED_SRF(U, V)	CAGD_MESH_UV(DerivedSrf, U, V)
#define RAISED_SRF(U, V)	CAGD_MESH_UV(RaisedSrf, U, V)
#define SRF(U, V)		CAGD_MESH_UV(Srf, U, V)

#define NORMAL_EPSILON		1e-4

static CagdBType EvalTangentVector(CagdVType Vec,
				   CagdCrvStruct *Crv,
				   CagdRType t,
				   CagdRType TMin,
				   CagdRType TMax);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bspline surface - subdivides it into two sub-surfaces at the given M
* parametric value.                                                          M
*   Returns pointer to first surface in a list of two subdivided surfaces.   M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      To subdivide at parameter value t.                             M
*   t:        Parameter value to subdivide Srf at.                           M
*   Dir:      Direction of subdivision. Either U or V.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A list of the two subdivided surfaces.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfSubdivAtParam, subdivision, refinement                             M
*****************************************************************************/
CagdSrfStruct *BspSrfSubdivAtParam(CagdSrfStruct *Srf,
				   CagdRType t,
				   CagdSrfDirType Dir)
{
    CagdBType
	NewSrf = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf);
    int i, j, Row, Col, KVLen, Index1, Index2, Mult,
	LULength, RULength, LVLength, RVLength,	ULength, VLength,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdRType *RefKV, **Pts, **LPts, **RPts, *LOnePts, *ROnePts, *OnePts,
    	UMin, UMax, VMin, VMax;
    CagdSrfStruct *RSrf, *LSrf;
    BspKnotAlphaCoeffType *A;

    if (CAGD_IS_PERIODIC_SRF(Srf)) {
	NewSrf = TRUE;
	Srf = CnvrtPeriodic2FloatSrf(Srf);
    }

    ULength = Srf -> ULength;
    VLength = Srf -> VLength;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    RefKV = Srf -> UKnotVector;
	    KVLen = UOrder + ULength;
	    Index1 = BspKnotLastIndexL(RefKV, KVLen, t);
	    if (Index1 + 1 < UOrder)
		Index1 = UOrder - 1;
	    Index2 = BspKnotFirstIndexG(RefKV, KVLen, t);
	    if (Index2 > ULength)
		Index2 = ULength;
	    LSrf = BspSrfNew(Index1 + 1, VLength,
			     UOrder, VOrder, Srf -> PType);
	    RSrf = BspSrfNew(ULength - Index2 + UOrder, VLength,
			     UOrder, VOrder, Srf -> PType);
	    Mult = UOrder - 1 - (Index2 - Index1 - 1);

	    /* Update the new knot vectors. */
	    CAGD_GEN_COPY(LSrf -> UKnotVector,
			  Srf -> UKnotVector,
			  sizeof(CagdRType) * (Index1 + 1));
	    /* Close the knot vector with multiplicity Order: */
	    for (j = Index1 + 1; j <= Index1 + UOrder; j++)
		LSrf -> UKnotVector[j] = t;
	    CAGD_GEN_COPY(&RSrf -> UKnotVector[UOrder],
			  &Srf -> UKnotVector[Index2],
			  sizeof(CagdRType) * (ULength + UOrder - Index2));
	    /* Make sure knot vector starts with multiplicity Order: */
	    for (j = 0; j < UOrder; j++)
		RSrf -> UKnotVector[j] = t;
	    
	    /* And copy the other direction knot vectors. */
	    CAGD_GEN_COPY(LSrf -> VKnotVector,
			  Srf -> VKnotVector,
			  sizeof(CagdRType) * (VOrder + VLength));
	    CAGD_GEN_COPY(RSrf -> VKnotVector,
			  Srf -> VKnotVector,
			  sizeof(CagdRType) * (VOrder + VLength));
	    break;
	case CAGD_CONST_V_DIR:
	    RefKV = Srf -> VKnotVector;
	    KVLen = VOrder + VLength;
	    Index1 = BspKnotLastIndexL(RefKV, KVLen, t);
	    if (Index1 + 1 < VOrder)
		Index1 = VOrder - 1;
	    Index2 = BspKnotFirstIndexG(RefKV, KVLen, t);
	    if (Index2 > VLength)
		Index2 = VLength;
	    LSrf = BspSrfNew(ULength, Index1 + 1,
			     UOrder, VOrder, Srf -> PType);
	    RSrf = BspSrfNew(ULength, VLength - Index2 + VOrder,
			     UOrder, VOrder, Srf -> PType);
	    Mult = VOrder - 1 - (Index2 - Index1 - 1);

	    /* Update the new knot vectors. */
	    CAGD_GEN_COPY(LSrf -> VKnotVector,
			  Srf -> VKnotVector,
			  sizeof(CagdRType) * (Index1 + 1));
	    /* Close the knot vector with multiplicity Order: */
	    for (j = Index1 + 1; j <= Index1 + VOrder; j++)
		LSrf -> VKnotVector[j] = t;
	    CAGD_GEN_COPY(&RSrf -> VKnotVector[VOrder],
			  &Srf -> VKnotVector[Index2],
			  sizeof(CagdRType) * (VLength + VOrder - Index2));
	    /* Make sure knot vector starts with multiplicity Order: */
	    for (j = 0; j < VOrder; j++)
		RSrf -> VKnotVector[j] = t;
	    
	    /* And copy the other direction knot vectors. */
	    CAGD_GEN_COPY(LSrf -> UKnotVector,
			  Srf -> UKnotVector,
			  sizeof(CagdRType) * (UOrder + ULength));
	    CAGD_GEN_COPY(RSrf -> UKnotVector,
			  Srf -> UKnotVector,
			  sizeof(CagdRType) * (UOrder + ULength));
	    break;
	default:
	    Mult = 1;
	    RefKV = NULL;
	    LSrf = RSrf = NULL;
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    Pts = Srf -> Points;
    LPts = LSrf -> Points;
    RPts = RSrf -> Points;
    LULength = LSrf -> ULength,
    RULength = RSrf -> ULength;
    LVLength = LSrf -> VLength,
    RVLength = RSrf -> VLength;

    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    /* Compute the Alpha refinement matrix. */
	    if (Mult > 0) {
		CagdRType
		    *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);

		CAGD_DOMAIN_T_VERIFY(t, UMin, UMax);
		if (t == UMax)
		    t -= CAGD_DOMAIN_EPSILON;
		for (i = 0; i < Mult; i++)
		    NewKV[i] = t;
		A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength,
								NewKV, Mult);
		IritFree((VoidPtr) NewKV);
	    }
	    else
		A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength, NULL, 0);

	    /* Now work on the two new surfaces meshes. */

	    /* Note that Mult can be negative in cases where original       */
	    /* multiplicity was order or more and we need to compensate     */
	    /* here, since Alpha matrix will be just a unit matrix then.    */
	    Mult = Mult >= 0 ? 0 : -Mult;

	    for (Row = 0; Row < VLength; Row++) {
	        /* Blend Srf into LSrf. */
	        for (j = IsNotRational; j <= MaxCoord; j++) {
		    LOnePts = &LPts[j][Row * LULength];
		    OnePts = &Pts[j][Row * ULength];
		    for (i = 0; i < LULength; i++, LOnePts++)
		        CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, LOnePts);
		}

	        /* Blend Srf into LSrf. */
		for (j = IsNotRational; j <= MaxCoord; j++) {
		    ROnePts = &RPts[j][Row * RULength];
		    OnePts = &Pts[j][Row * ULength];
		    for (i = LSrf -> ULength - 1 + Mult;
			 i < LSrf -> ULength + RSrf -> ULength - 1 + Mult;
			 i++, ROnePts++)
			CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, ROnePts);
		}
	    }

	    BspKnotFreeAlphaCoef(A);
	    break;
	case CAGD_CONST_V_DIR:
	    /* Compute the Alpha refinement matrix. */
	    if (Mult > 0) {
		CagdRType
		    *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);

		CAGD_DOMAIN_T_VERIFY(t, VMin, VMax);
		if (t == VMax)
		    t -= CAGD_DOMAIN_EPSILON;
		for (i = 0; i < Mult; i++)
		    NewKV[i] = t;
		A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength,
						  NewKV, Mult);
		IritFree((VoidPtr) NewKV);
	    }
	    else
		A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength, NULL, 0);

	    /* Now work on the two new surfaces meshes. */

	    /* Note that Mult can be negative in cases where original       */
	    /* multiplicity was order or more and we need to compensate     */
	    /* here, since Alpha matrix will be just a unit matrix then.    */
	    Mult = Mult >= 0 ? 0 : -Mult;

	    for (Col = 0; Col < ULength; Col++) {
		LULength = LSrf -> ULength;
		RULength = RSrf -> ULength;

		/* Blend Srf into LSrf. */
		for (j = IsNotRational; j <= MaxCoord; j++) {
		    LOnePts = &LPts[j][Col];
		    OnePts = &Pts[j][Col];
		    for (i = 0;
		         i < LVLength;
		         i++, LOnePts += LULength)
		        CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength,
					      LOnePts, ULength);
		}

		/* Blend Srf into RSrf. */
		for (j = IsNotRational; j <= MaxCoord; j++) {
		    ROnePts = &RPts[j][Col];
		    OnePts = &Pts[j][Col];
		    for (i = LVLength - 1 + Mult;
			 i < LVLength + RVLength - 1 + Mult;
			 i++, ROnePts += RULength)
			CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength,
					      ROnePts, ULength);
		}
	    }

	    BspKnotFreeAlphaCoef(A);
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    BspKnotMakeRobustKV(RSrf -> UKnotVector,
			RSrf -> UOrder + RSrf -> ULength);
    BspKnotMakeRobustKV(RSrf -> VKnotVector,
			RSrf -> VOrder + RSrf -> VLength);

    BspKnotMakeRobustKV(LSrf -> UKnotVector,
			LSrf -> UOrder + LSrf -> ULength);
    BspKnotMakeRobustKV(LSrf -> VKnotVector,
			LSrf -> VOrder + LSrf -> VLength);

    LSrf -> Pnext = RSrf;

    if (NewSrf)
	CagdSrfFree(Srf);

    return LSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Inserts n knot, all with the value t in direction Dir. In no case will    M
* the multiplicity of a knot be greater or equal to the curve order.         M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To refine by insertion (upto) n knot of value t.             M
*   Dir:        Direction of refinement. Either U or V.                      M
*   t:          Parameter value of new knot to insert.                       M
*   n:          Maximum number of times t should be inserted.                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  Refined Srf with n knots of value t in direction Dir.  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfKnotInsertNSame, refinement, subdivision                           M
*****************************************************************************/
CagdSrfStruct *BspSrfKnotInsertNSame(CagdSrfStruct *Srf,
				     CagdSrfDirType Dir,
				     CagdRType t,
				     int n)
{
    CagdBType
	NewSrf = FALSE;
    int i, CrntMult, Mult;
    CagdSrfStruct *RefinedSrf;

    if (CAGD_IS_PERIODIC_SRF(Srf)) {
        NewSrf = TRUE;
        Srf = CnvrtPeriodic2FloatSrf(Srf);
    }

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    CrntMult = BspKnotFindMult(Srf -> UKnotVector, Srf -> UOrder,
							 Srf -> ULength, t),
	    Mult = MIN(n, Srf -> UOrder - CrntMult - 1);
	    break;
	case CAGD_CONST_V_DIR:
	    CrntMult = BspKnotFindMult(Srf -> VKnotVector, Srf -> VOrder,
							 Srf -> VLength, t),
	    Mult = MIN(n, Srf -> VOrder - CrntMult - 1);
	    break;
	default:
	    Mult = 0;
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    if (Mult > 0) {
	CagdRType
	    *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);

	for (i = 0; i < Mult; i++)
	    NewKV[i] = t;

	RefinedSrf = BspSrfKnotInsertNDiff(Srf, Dir, FALSE, NewKV, Mult);

	IritFree((VoidPtr) NewKV);
    }
    else {
	RefinedSrf = CagdSrfCopy(Srf);
    }

    if (NewSrf)
	CagdSrfFree(Srf);

    return RefinedSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Inserts n knot with different values as defined by the vector t. If,      M
* however, Replace is TRUE, the knot are simply replacing the current knot   M
* vector.                                                                    M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To refine by insertion (upto) n knot of value t.             M
*   Dir:        Direction of refinement. Either U or V.                      M
*   Replace:    if TRUE, the n knots in t should replace the knot vector     M
*               of size n of Srf. Sizes must match. If False, n new knots    M
*               as defined by t will be introduced into Srf.                 M
*   t:          New knots to introduce/replace knot vector of Srf.           M
*   n:          Size of t.                                                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:    Refined Srf with n new knots in direction Dir.       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfKnotInsertNDiff, refinement, subdivision                           M
*****************************************************************************/
CagdSrfStruct *BspSrfKnotInsertNDiff(CagdSrfStruct *Srf,
				     CagdSrfDirType Dir,
				     int Replace,
				     CagdRType *t,
				     int n)
{
    CagdBType
	NewSrf = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, Row, Col, ULength, VLength,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct
	*RefSrf = NULL;

    if (CAGD_IS_PERIODIC_SRF(Srf)) {
        NewSrf = TRUE;
        Srf = CnvrtPeriodic2FloatSrf(Srf);
    }

    ULength = Srf -> ULength;
    VLength = Srf -> VLength;

    if (Replace) {
	for (i = 1; i < n; i++)
	    if (t[i] < t[i - 1])
		CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);

    	switch (Dir) {
	    case CAGD_CONST_U_DIR:
		if (Srf -> UOrder + Srf -> ULength != n)
		    CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);

		RefSrf = CagdSrfCopy(Srf);
		for (i = 0; i < n; i++)
		    RefSrf -> UKnotVector[i] = *t++;
		break;
	    case CAGD_CONST_V_DIR:
		if (Srf -> VOrder + Srf -> VLength != n)
		    CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);

		RefSrf = CagdSrfCopy(Srf);
		for (i = 0; i < n; i++)
		    RefSrf -> VKnotVector[i] = *t++;
		break;
	    default:
		CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
		break;
	}
    }
    else if (n == 0) {
	RefSrf = CagdSrfCopy(Srf);
    }
    else {
	int j, LengthKVt, RULength, RVLength;
	BspKnotAlphaCoeffType *A;
	CagdRType *MergedKVt, UMin, UMax, VMin, VMax,
	    *UKnotVector = Srf -> UKnotVector,
	    *VKnotVector = Srf -> VKnotVector;

	CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);

	for (i = 1; i < n; i++)
	    if (t[i] < t[i - 1])
		CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);

	switch (Dir) {
	    case CAGD_CONST_U_DIR:
	        for (i = 0; i < n; i++) {
		    CAGD_DOMAIN_T_VERIFY(t[i], UMin, UMax);
		    if (t[i] == UMax)
			t[i] -= CAGD_DOMAIN_EPSILON;
		}

		/* Compute the Alpha refinement matrix. */
		MergedKVt = BspKnotMergeTwo(UKnotVector, ULength + UOrder,
					    t, n, 0, &LengthKVt);
		A = BspKnotEvalAlphaCoef(UOrder, UKnotVector, ULength,
					 MergedKVt, LengthKVt - UOrder);

	        RefSrf = BspSrfNew(ULength + n, VLength, UOrder, VOrder,
							      Srf -> PType);
		IritFree((VoidPtr) RefSrf -> UKnotVector);
		IritFree((VoidPtr) RefSrf -> VKnotVector);
		RefSrf -> UKnotVector = MergedKVt;
		RefSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
					Srf -> VLength + Srf -> VOrder);

		RULength = RefSrf -> ULength;

		/* Update the control mesh */
		for (Row = 0; Row < VLength; Row++) {
		    for (j = IsNotRational; j <= MaxCoord; j++) {
			CagdRType
			    *ROnePts = &RefSrf -> Points[j][Row * RULength],
			    *OnePts = &Srf -> Points[j][Row * ULength];

			for (i = 0; i < RULength; i++, ROnePts++)
			    CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength,
					     ROnePts);
		    }
		}

		BspKnotFreeAlphaCoef(A);
		break;
	    case CAGD_CONST_V_DIR:
		for (i = 0; i < n; i++) {
		    CAGD_DOMAIN_T_VERIFY(t[i], VMin, VMax);
		    if (t[i] == VMax)
			t[i] -= CAGD_DOMAIN_EPSILON;
		}

		/* Compute the Alpha refinement matrix. */
		MergedKVt = BspKnotMergeTwo(VKnotVector, VLength + VOrder,
					    t, n, 0, &LengthKVt);
		A = BspKnotEvalAlphaCoef(VOrder, VKnotVector, VLength,
					 MergedKVt, LengthKVt - VOrder);

	        RefSrf = BspSrfNew(ULength, VLength + n, UOrder, VOrder,
							      Srf -> PType);
		IritFree((VoidPtr) RefSrf -> UKnotVector);
		IritFree((VoidPtr) RefSrf -> VKnotVector);
		RefSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
					Srf -> ULength + Srf -> UOrder);
		RefSrf -> VKnotVector = MergedKVt;

		RULength = RefSrf -> ULength;
		RVLength = RefSrf -> VLength;

		/* Update the control mesh */
		for (Col = 0; Col < ULength; Col++) {
		    for (j = IsNotRational; j <= MaxCoord; j++) {
			CagdRType
			    *ROnePts = &RefSrf -> Points[j][Col],
			    *OnePts = &Srf -> Points[j][Col];

			for (i = 0; i < RVLength; i++, ROnePts += RULength)
			    CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength,
						  ROnePts, ULength);
		    }
		}

		BspKnotFreeAlphaCoef(A);
		break;
	    default:
		CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
		break;
	}
    }

    BspKnotMakeRobustKV(RefSrf -> UKnotVector,
			RefSrf -> UOrder + RefSrf -> ULength);
    BspKnotMakeRobustKV(RefSrf -> VKnotVector,
			RefSrf -> VOrder + RefSrf -> VLength);

    if (NewSrf)
	CagdSrfFree(Srf);

    return RefSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new Bspline surface, identical to the original but with one      M
* degree higher, in the requested direction Dir.                             M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To raise it degree by one.                                   M
*   Dir:        Direction of degree raising. Either U or V.                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A surface with one degree higher in direction Dir,     M
*                     representing the same geometry as Srf.		     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfDegreeRaise, degree raising                                        M
*****************************************************************************/
CagdSrfStruct *BspSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    CagdBType
	NewSrf = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, i2, j, RaisedLen, Row, Col, Length,
	Order = Dir == CAGD_CONST_V_DIR ? Srf -> UOrder : Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct
	*RaisedSrf = NULL;

    if (CAGD_IS_PERIODIC_SRF(Srf)) {
        NewSrf = TRUE;
        Srf = CnvrtPeriodic2FloatSrf(Srf);
    }

    Length = Dir == CAGD_CONST_V_DIR ? Srf -> ULength : Srf -> VLength;

    if (Order > 2) {
	CagdSrfStruct *UnitSrf;
	int UKvLen1 = Srf -> UOrder + Srf -> ULength - 1,
	    VKvLen1 = Srf -> VOrder + Srf -> VLength - 1;
	CagdRType
	    *UKv = Srf -> UKnotVector,
	    *VKv = Srf -> VKnotVector;

	/* Degree raise by multiplying by a constant 1 linear surface in the */
	/* raised direction and constant 1 constant surface in the other.    */

	switch (Dir) {
	    case CAGD_CONST_U_DIR:
		UnitSrf = BspSrfNew(1, 2, 1, 2,
				    CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
		for (i = 0; i < 2; i++)
		    UnitSrf -> UKnotVector[i] = i > 0 ? UKv[UKvLen1] : UKv[0];
		for (i = 0; i < 4; i++)
		    UnitSrf -> VKnotVector[i] = i > 1 ? VKv[VKvLen1] : VKv[0];
		break;
	    case CAGD_CONST_V_DIR:
		UnitSrf = BspSrfNew(2, 1, 2, 1,
				    CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
		for (i = 0; i < 4; i++)
		    UnitSrf -> UKnotVector[i] = i > 1 ? UKv[UKvLen1] : UKv[0];
		for (i = 0; i < 2; i++)
		    UnitSrf -> VKnotVector[i] = i > 0 ? VKv[VKvLen1] : VKv[0];
		break;
	    default:
		UnitSrf = NULL;
		CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
		break;
	}
	for (i = 1; i <= MaxCoord; i++)
	    UnitSrf -> Points[i][0] = UnitSrf -> Points[i][1] = 1.0;

	RaisedSrf = BspSrfMult(Srf, UnitSrf);

	CagdSrfFree(UnitSrf);

	if (NewSrf)
	    CagdSrfFree(Srf);

	return RaisedSrf;
    }

    /* If surface is linear, degree raising means basically to increase the  */
    /* knot multiplicity of each segment by one and add a middle point for   */
    /* each such segment.						     */
    RaisedLen = Length * 2 - 1;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    RaisedSrf = BspSrfNew(Srf -> ULength, RaisedLen,
				  Srf -> UOrder, Order + 1, Srf -> PType);

	    /* Update the knot vectors. */
	    CAGD_GEN_COPY(RaisedSrf -> UKnotVector,
			  Srf -> UKnotVector,
			  sizeof(CagdRType) * (Srf -> ULength + Srf -> UOrder));
	    for (i = 0; i < 3; i++)
		RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[0];
	    for (i = 2, j = 3; i < Length; i++, j += 2)
		RaisedSrf -> VKnotVector[j] = RaisedSrf -> VKnotVector[j + 1] = 
		    Srf -> VKnotVector[i];
	    for (i = j; i < j + 3; i++)
		RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[Length];

	    /* Update the mesh. */
       	    for (Col = 0; Col < Srf -> ULength; Col++) {
		for (j = IsNotRational; j <= MaxCoord; j++)  /* First point. */
		    RaisedSrf -> Points[j][RAISED_SRF(Col, 0)] =
			Srf -> Points[j][SRF(Col, 0)];

		for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
		    for (j = IsNotRational; j <= MaxCoord; j++) {
			RaisedSrf -> Points[j][RAISED_SRF(Col, i2)] =
				Srf -> Points[j][SRF(Col, i - 1)] * 0.5 +
				Srf -> Points[j][SRF(Col, i)] * 0.5;
			RaisedSrf -> Points[j][RAISED_SRF(Col, i2 + 1)] =
				Srf -> Points[j][SRF(Col, i)];
		    }
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    RaisedSrf = BspSrfNew(RaisedLen, Srf -> VLength,
				  Order + 1, Srf -> VOrder, Srf -> PType);

	    /* Update the knot vectors. */
	    CAGD_GEN_COPY(RaisedSrf -> VKnotVector,
			  Srf -> VKnotVector,
			  sizeof(CagdRType) * (Srf -> VLength + Srf -> VOrder));
	    for (i = 0; i < 3; i++)
		RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[0];
	    for (i = 2, j = 3; i < Length; i++, j += 2)
		RaisedSrf -> UKnotVector[j] = RaisedSrf -> UKnotVector[j + 1] = 
		    Srf -> UKnotVector[i];
	    for (i = j; i < j + 3; i++)
		RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[Length];

	    /* Update the mesh. */
       	    for (Row = 0; Row < Srf -> VLength; Row++) {
		for (j = IsNotRational; j <= MaxCoord; j++)  /* First point. */
		    RaisedSrf -> Points[j][RAISED_SRF(0, Row)] =
			Srf -> Points[j][SRF(0, Row)];

		for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
		    for (j = IsNotRational; j <= MaxCoord; j++) {
			RaisedSrf -> Points[j][RAISED_SRF(i2, Row)] =
				Srf -> Points[j][SRF(i - 1, Row)] * 0.5 +
				Srf -> Points[j][SRF(i, Row)] * 0.5;
			RaisedSrf -> Points[j][RAISED_SRF(i2 + 1, Row)] =
				Srf -> Points[j][SRF(i, Row)];
		    }
	    }
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    if (NewSrf)
	CagdSrfFree(Srf);

    return RaisedSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new surface equal to the given surface, differentiated once in   M
* the direction Dir.                                                         M
*   Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M
* Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2.  V
* This is applied to all rows/cols of the surface.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To differentiate.                                            M
*   Dir:        Direction of differentiation. Either U or V.                 M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   Differentiated surface.                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfDerive, derivatives, partial derivatives                           M
*****************************************************************************/
CagdSrfStruct *BspSrfDerive(CagdSrfStruct *Srf,
			    CagdSrfDirType Dir)
{
    CagdBType
	NewSrf = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, Row, Col, NewUOrder, NewVOrder, NewULength, NewVLength,
	ULength, VLength,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdRType **DPoints, *UKv, *VKv, **Points;
    CagdSrfStruct
        *DerivedSrf = NULL;

    if (CAGD_IS_PERIODIC_SRF(Srf)) {
        NewSrf = TRUE;
        Srf = CnvrtPeriodic2FloatSrf(Srf);
    }

    if (!IsNotRational)
	return BspSrfDeriveRational(Srf, Dir);

    ULength = Srf -> ULength;
    VLength = Srf -> VLength;
    UKv = Srf -> UKnotVector;
    VKv = Srf -> VKnotVector;
    Points = Srf -> Points;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    NewULength = UOrder < 2 ? ULength : ULength - 1;
	    NewUOrder = MAX(UOrder - 1, 1);
	    DerivedSrf = BspSrfNew(NewULength, VLength,
				   NewUOrder, VOrder, Srf -> PType);
	    CAGD_GEN_COPY(DerivedSrf -> UKnotVector, &UKv[UOrder < 2 ? 0 : 1],
			  sizeof(CagdRType) * (NewULength + NewUOrder));
	    CAGD_GEN_COPY(DerivedSrf -> VKnotVector, VKv,
			  sizeof(CagdRType) * (VLength + VOrder));
	    DPoints = DerivedSrf -> Points;

       	    for (Row = 0; Row < VLength; Row++)
		for (i = 0; i < NewULength; i++) {
		    CagdRType
			Denom = UKv[i + UOrder] - UKv[i + 1];

		    for (j = IsNotRational; j <= MaxCoord; j++) {
			if (APX_EQ(Denom, 0.0))
			    Denom = INFINITY;

			DPoints[j][DERIVED_SRF(i, Row)] =
			    UOrder < 2 ? 0.0
				       : (UOrder - 1) *
					    (Points[j][SRF(i + 1, Row)] -
					     Points[j][SRF(i, Row)]) / Denom;
		    }
		}
	    break;
	case CAGD_CONST_V_DIR:
	    NewVLength = VOrder < 2 ? VLength : VLength - 1;
	    NewVOrder = MAX(VOrder - 1, 1);
	    DerivedSrf = BspSrfNew(ULength, NewVLength,
	    			   UOrder, NewVOrder, Srf -> PType);
	    CAGD_GEN_COPY(DerivedSrf -> UKnotVector, UKv,
			  sizeof(CagdRType) * (ULength + UOrder));
	    CAGD_GEN_COPY(DerivedSrf -> VKnotVector, &VKv[VOrder < 2 ? 0 : 1],
			  sizeof(CagdRType) * (NewVLength + NewVOrder));
	    DPoints = DerivedSrf -> Points;

	    for (Col = 0; Col < ULength; Col++)
		for (i = 0; i < NewVLength; i++) {
		    CagdRType
			Denom = VKv[i + VOrder] - VKv[i + 1];

		    for (j = IsNotRational; j <= MaxCoord; j++) {
			if (APX_EQ(Denom, 0.0))
			    Denom = INFINITY;

			DPoints[j][DERIVED_SRF(Col, i)] =
			    VOrder < 2 ? 0.0
				       : (VOrder - 1) *
					    (Points[j][SRF(Col, i + 1)] -
					     Points[j][SRF(Col, i)]) / Denom;
		    }
		}
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    if (NewSrf)
	CagdSrfFree(Srf);

    return DerivedSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluates the unit tangent to a surface at a given parametric location     M
( u, v) and given direction Dir.                                             M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Bspline surface to evaluate tangent vector for.              M
*   u, v:       Parametric location of required unit tangent.                M
*   Dir:        Direction of tangent vector. Either U or V.                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the unit tangent  M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfTangent, tangent                                                   M
*****************************************************************************/
CagdVecStruct *BspSrfTangent(CagdSrfStruct *Srf,
			     CagdRType u,
			     CagdRType v,
			     CagdSrfDirType Dir)
{
    CagdVecStruct
	*Tangent = NULL;
    CagdCrvStruct *Crv;

    switch (Dir) {
	case CAGD_CONST_V_DIR:
	    Crv = BspSrfCrvFromSrf(Srf, v, Dir);
	    Tangent = BspCrvTangent(Crv, u);
	    CagdCrvFree(Crv);
	    break;
	case CAGD_CONST_U_DIR:
	    Crv = BspSrfCrvFromSrf(Srf, u, Dir);
	    Tangent = BspCrvTangent(Crv, v);
	    CagdCrvFree(Crv);
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return Tangent;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluate the unit normal of a surface at a given parametric location.      M
*   If we fail to compute the normal at given location we retry by moving a  M
* tad.                                                                       M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Bspline surface to evaluate normal vector for.               M
*   u, v:       Parametric location of required unit normal.                 M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the unit normal   M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfNormal, normal                                                     M
*****************************************************************************/
CagdVecStruct *BspSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
{
    static CagdVecStruct Normal;
    CagdVecStruct *V, T1, T2;
    CagdRType UMin, UMax, VMin, VMax;

    CAGD_DOMAIN_GET_AND_VERIFY_SRF(u, v, Srf, UMin, UMax, VMin, VMax);

    V = BspSrfTangent(Srf, u, v, CAGD_CONST_U_DIR);
    if (CAGD_LEN_VECTOR(*V) < IRIT_EPSILON)
	V = BspSrfTangent(Srf,
			  u > UMin + EPSILON ? u - EPSILON : u + EPSILON,
			  v > VMin + EPSILON ? v - EPSILON : v + EPSILON,
			  CAGD_CONST_U_DIR);
    CAGD_COPY_VECTOR(T1, *V);

    V = BspSrfTangent(Srf, u, v, CAGD_CONST_V_DIR);
    if (CAGD_LEN_VECTOR(*V) < IRIT_EPSILON)
	V = BspSrfTangent(Srf,
			  u > UMin + EPSILON ? u - EPSILON : u + EPSILON,
			  v > VMin + EPSILON ? v - EPSILON : v + EPSILON,
			  CAGD_CONST_V_DIR);
    CAGD_COPY_VECTOR(T2, *V);

    /* The normal is the cross product of T1 and T2: */
    Normal.Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1];
    Normal.Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2];
    Normal.Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0];

    CAGD_NORMALIZE_VECTOR(Normal);		    /* Normalize the vector. */

    return &Normal;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluates the unit normals of a surface at a mesh defined by subdividing   M
* the parametric space into a grid of size UFineNess by VFineNess.	     M
*   The normals are saved in a linear CagdVecStruct vector which is          M
* allocated dynamically. Data is saved u inc. first.			     M
*   This routine is much faster than evaluating normal for each point,       M
* individually.                                                              M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:          To compute normals on a grid of its parametric domain.     M
*   UFineNess:    U Fineness of imposed grid on Srf's parametric domain.     M
*   VFineNess:    V Fineness of imposed grid on Srf's parametric domain.     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:   An vector of unit normals (u increments first).       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfMeshNormals, normal                                                M
*****************************************************************************/
CagdVecStruct *BspSrfMeshNormals(CagdSrfStruct *Srf,
				 int UFineNess,
				 int VFineNess)
{
    int i, j;
    CagdRType UMin, UMax, VMin, VMax;
    CagdVecStruct *Normals, *NPtr;
    CagdSrfStruct
	*DuSrf = CagdSrfDerive(Srf, CAGD_CONST_U_DIR),
	*DvSrf = CagdSrfDerive(Srf, CAGD_CONST_V_DIR);

    UFineNess = MAX(2, UFineNess);
    VFineNess = MAX(2, VFineNess);

    Normals = CagdVecArrayNew(UFineNess * VFineNess);

    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);

    NPtr = Normals;
    for (i = 0; i < UFineNess; i++) {
	CagdRType
	   U = UMin + (UMax - UMin) * i / (UFineNess - 1);
	CagdCrvStruct
	   *DuCrv2 = NULL,
	   *DvCrv2 = NULL,
	   *DuCrv = CagdCrvFromSrf(DuSrf, U, CAGD_CONST_U_DIR),
	   *DvCrv = CagdCrvFromSrf(DvSrf, U, CAGD_CONST_U_DIR);

	for (j = 0; j < VFineNess; j++) {
	    CagdVType Du, Dv;
	    CagdRType
		V = VMin + (VMax - VMin) * j / (VFineNess - 1);

	    if (!EvalTangentVector(Du, DuCrv, V, VMin, VMax)) {
		if (DuCrv2 == NULL) {
		    DuCrv2 = CagdCrvFromSrf(DuSrf,
					    U > UMin + NORMAL_EPSILON ?
					        U - NORMAL_EPSILON :
					        U + NORMAL_EPSILON,
					    CAGD_CONST_U_DIR);
		}
		EvalTangentVector(Du, DuCrv2, V, VMin, VMax);
	    }

	    if (!EvalTangentVector(Dv, DvCrv, V, VMin, VMax)) {
		if (DvCrv2 == NULL) {
		    DvCrv2 = CagdCrvFromSrf(DvSrf,
					    U > UMin + NORMAL_EPSILON ?
					        U - NORMAL_EPSILON :
					        U + NORMAL_EPSILON,
					    CAGD_CONST_U_DIR);
		}
		EvalTangentVector(Dv, DvCrv2, V, VMin, VMax);
	    }

	    CROSS_PROD(NPtr -> Vec, Dv, Du);
	    NPtr++;
	}

	if (DuCrv2 != NULL)
	    CagdCrvFree(DuCrv2);
	if (DvCrv2 != NULL)
	    CagdCrvFree(DvCrv2);
	CagdCrvFree(DuCrv);
	CagdCrvFree(DvCrv);
    }

    /* Normalize the results. */
    NPtr = Normals;
    for (i = 0; i < UFineNess; i++) {
	for (j = 0; j < VFineNess; j++) {
	    CagdRType
		Len = PT_LENGTH(NPtr -> Vec);

	    if (Len < PT_NORMALIZE_ZERO) {
		CagdVecStruct *NPtrTmp;

		/* Try its neighbors. */
		if (i > 0) {
		    NPtrTmp = NPtr - VFineNess;
		    PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
		    Len = PT_LENGTH(NPtr -> Vec);
		}
		if (Len < PT_NORMALIZE_ZERO && i < UFineNess - 1) {
		    NPtrTmp = NPtr + VFineNess;
		    PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
		    Len = PT_LENGTH(NPtr -> Vec);
		}
		if (Len < PT_NORMALIZE_ZERO && j > 0) {
		    NPtrTmp = NPtr - 1;
		    PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
		    Len = PT_LENGTH(NPtr -> Vec);
		}
		if (Len < PT_NORMALIZE_ZERO && j < VFineNess - 1) {
		    NPtrTmp = NPtr + 1;
		    PT_COPY(NPtr -> Vec, NPtrTmp -> Vec);
		    Len = PT_LENGTH(NPtr -> Vec);
		}
		if (Len > PT_NORMALIZE_ZERO) {
		    Len = 1 / Len;
		    PT_SCALE(NPtr -> Vec, Len);
		}
		else {
		    /* Do something. */
		    NPtr -> Vec[0] = NPtr -> Vec[1] = 0.0;
		    NPtr -> Vec[2] = 1.0;
		}
	    }
	    else {
		Len = 1 / Len;
		PT_SCALE(NPtr -> Vec, Len);
	    }
	    NPtr++;
	}
    }


    CagdSrfFree(DuSrf);
    CagdSrfFree(DvSrf);

    return Normals;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Evaluates a vector field Crv at param t. If zero length, tries to move a *
* tad.                                                                       *
*                                                                            *
* PARAMETERS:                                                                *
*   Vec:     Where to place the result.                                      *
*   Crv:     Vector field to evaluate.                                       *
*   t:       Parameter value where to evaluate the vector field.             *
*   TMin:    Minimum of vector field domain.                                 *
*   TMax:    Maximum of vector field domain.                                 *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdBType:  TRUE, if successful, FALSE otherwise.                        *
*****************************************************************************/
static CagdBType EvalTangentVector(CagdVType Vec,
				   CagdCrvStruct *Crv,
				   CagdRType t,
				   CagdRType TMin,
				   CagdRType TMax)
{
    CagdRType Len,
	*R = CagdCrvEval(Crv, t);
    
    CagdCoerceToE3(Vec, &R, -1, Crv -> PType);
    Len = PT_LENGTH(Vec);

    if (Len < PT_NORMALIZE_ZERO && t - NORMAL_EPSILON > TMin) {
	CagdCrvEval(Crv, t - NORMAL_EPSILON);
	CagdCoerceToE3(Vec, &R, -1, Crv -> PType);
	Len = PT_LENGTH(Vec);
    }
    if (Len < PT_NORMALIZE_ZERO && t + NORMAL_EPSILON < TMax) {
	CagdCrvEval(Crv, t + NORMAL_EPSILON);
	CagdCoerceToE3(Vec, &R, -1, Crv -> PType);
	Len = PT_LENGTH(Vec);
    }

    if (Len > PT_NORMALIZE_ZERO) {
	Len = 1.0 / Len;
	PT_SCALE(Vec, Len);
	return TRUE;
    }
    else
        return FALSE;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluates the unit normals of a surface at a mesh defined by subdividing   M
* the parametric space into a grid of size UFineNess by VFineNess.	     M
*   The normals are saved in a linear CagdVecStruct vector which is          M
* allocated dynamically. Data is saved u inc. first.			     M
*   This routine is much faster than evaluating normal for each point,       M
* individually.                                                              M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:          To compute normals on a grid of its parametric domain.     M
*   UFineNess:    U Fineness of imposed grid on Srf's parametric domain.     M
*   VFineNess:    V Fineness of imposed grid on Srf's parametric domain.     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:   An vector of unit normals (u increments first).       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfMeshNormalsSymb, normal                                            M
*****************************************************************************/
CagdVecStruct *BspSrfMeshNormalsSymb(CagdSrfStruct *Srf,
				     int UFineNess,
				     int VFineNess)
{
    int i, j;
    CagdRType U, V, UMin, UMax, VMin, VMax, **Points;
    CagdVecStruct *Normals, *NPtr;
    CagdSrfStruct *TstNormalSrf,
	*NormalSrf = SymbSrfNormalSrf(Srf);

    /* Verify we have a valid regular surface. */
    TstNormalSrf = CagdCoerceSrfTo(NormalSrf, CAGD_PT_E3_TYPE);
    for (Points = TstNormalSrf -> Points, i = 0;
	 i < TstNormalSrf -> ULength * TstNormalSrf -> VLength;
	 i++)
	if (!APX_EQ(Points[1][i], 0.0 ) ||
	    !APX_EQ(Points[2][i], 0.0 ) ||
	    !APX_EQ(Points[3][i], 0.0 ))
	    break;
    CagdSrfFree(TstNormalSrf);
    if (i >= TstNormalSrf -> ULength * TstNormalSrf -> VLength) {
	/* Not a regular surface - ignore it, */
	return NULL;
    }

    UFineNess = MAX(2, UFineNess);
    VFineNess = MAX(2, VFineNess);

    Normals = CagdVecArrayNew(UFineNess * VFineNess);

    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);

    NPtr = Normals;
    for (i = 0; i < UFineNess; i++) {
	for (j = 0; j < VFineNess; j++) {
	    U = UMin + (UMax - UMin) * i / (UFineNess - 1);
	    V = VMin + (VMax - VMin) * j / (VFineNess - 1);
	    CagdEvaluateSurfaceVecField(NPtr -> Vec, NormalSrf, U, V);

	    NPtr++;
	}
    }

    CagdSrfFree(NormalSrf);

    return Normals;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a Bspline surface into a Bspline surface with floating end        M
* conditions.                                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       Bspline surface to convert to floating end conditions. Assume M
*              Crv is either periodic or has floating end condition.         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A Bspline surface with floating end conditions,        M
*                     representing the same geometry as Srf.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtPeriodic2FloatSrf, conversion                                       M
*****************************************************************************/
CagdSrfStruct *CnvrtPeriodic2FloatSrf(CagdSrfStruct *Srf)
{
    int i, j,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	MaxAxis = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct *NewSrf;

    if (!CAGD_IS_UPERIODIC_SRF(Srf) && !CAGD_IS_VPERIODIC_SRF(Srf)) {
	CAGD_FATAL_ERROR(CAGD_ERR_PERIODIC_EXPECTED);
	return NULL;
    }

    NewSrf = BspSrfNew(CAGD_SRF_UPT_LST_LEN(Srf), CAGD_SRF_VPT_LST_LEN(Srf),
		       UOrder, VOrder, Srf -> PType);

    NewSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
					CAGD_SRF_UPT_LST_LEN(Srf) + UOrder);
    NewSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
					CAGD_SRF_VPT_LST_LEN(Srf) + UOrder);

    for (i = !CAGD_IS_RATIONAL_PT(Srf -> PType); i <= MaxAxis; i++) {
	CagdRType *NewPts,
	    *Pts = Srf -> Points[i];

	NewSrf -> Points[i] = (CagdRType *) IritMalloc(sizeof(CagdRType) *
						   CAGD_SRF_UPT_LST_LEN(Srf) *
						   CAGD_SRF_VPT_LST_LEN(Srf));
	NewPts = NewSrf -> Points[i];

	for (j = 0; j < VLength; j++, Pts += ULength) {
	    CAGD_GEN_COPY(NewPts, Pts, sizeof(CagdRType) * ULength);
	    NewPts += ULength;
	    if (Srf -> UPeriodic) {
		CAGD_GEN_COPY(NewPts, Pts,
			      sizeof(CagdRType) * (UOrder - 1));
		NewPts += UOrder - 1;
	    }
	}
	if (Srf -> VPeriodic) {
	    CAGD_GEN_COPY(NewPts, NewSrf -> Points[i],
			  sizeof(CagdRType) * (VOrder - 1) *
			  CAGD_SRF_UPT_LST_LEN(Srf));
	}
    }

    for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
	NewSrf -> Points[i] = NULL;
    
    return NewSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a Bspline surface to a Bspline surface with open end conditions.  M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       Bspline surface to convert to open end conditions.            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A Bspline surface with open end conditions,	     M
*                     representing the same geometry as Srf.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtFloat2OpenSrf, conversion                                           M
*****************************************************************************/
CagdSrfStruct *CnvrtFloat2OpenSrf(CagdSrfStruct *Srf)
{
    CagdRType UMin, UMax, VMin, VMax;
    CagdSrfStruct *TSrf1, *TSrf2;

    if (!CAGD_IS_BSPLINE_SRF(Srf)) {
	CAGD_FATAL_ERROR(CAGD_ERR_BSP_SRF_EXPECT);
	return NULL;
    }

    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);

    TSrf1 = CagdSrfRegionFromSrf(Srf, UMin, UMax, CAGD_CONST_U_DIR);
    TSrf2 = CagdSrfRegionFromSrf(TSrf1, VMin, VMax, CAGD_CONST_V_DIR);
    CagdSrfFree(TSrf1);
    return TSrf2;
}

