/******************************************************************************
* Adap_iso.c - adaptive isoline surface extraction algorithm.		      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, Aug. 92.					      *
******************************************************************************/

#include "symb_loc.h"

#define SIMILAR_PARAM_VAL	1e-4

#define CONVEX_BLEND(v1, v2, b1, b2, t) { \
		CagdRType \
		    Blend = -b1 / ( b2 - b1 ); \
		t = v2 * Blend + v1 * (1.0 - Blend); \
	    }

static int MinSubdivLevel = 1;

static CagdCrvStruct *SymbAdapIsoExtractAux(int Level,
					    CagdSrfStruct *Srf,
					    CagdSrfStruct *NSrf,
					    SymbAdapIsoDistSqrFuncType
					                      AdapIsoDistFunc,
					    CagdCrvStruct *Crv1,
					    CagdCrvStruct *NCrv1,
					    CagdCrvStruct *Crv2,
					    CagdCrvStruct *NCrv2,
					    CagdRType Crv1Param,
					    CagdRType Crv2Param,
					    CagdSrfDirType Dir,
					    CagdRType Eps2,
					    CagdBType FullIso,
					    CagdBType SinglePath);
static CagdCrvStruct *CopyRegionFromCrv(CagdCrvStruct *Crv,
					CagdRType TMin,
					CagdRType TMax);
static CagdRType ComputeMidParam(CagdRType Param1,
				 CagdRType Param2,
				 CagdSrfDirType Dir,
				 CagdSrfStruct *Srf);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Sets minimum level of subdivision forced in the adaptive iso extraction.   M
*                                                                            *
* PARAMETERS:                                                                M
*   MinLevel:    At least that many subdivision will occur.                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSetAdapIsoExtractMinLevel, adaptive isocurves                        M
*****************************************************************************/
void SymbSetAdapIsoExtractMinLevel(int MinLevel)
{
    MinSubdivLevel = MinLevel;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Extracts a valid coverage set of isolines from the given surface in the    M
* given direction and epsilon.						     M
*   If FullIso is TRUE, all extracted isocurves are spanning the entire      M
* parametric domain.							     M
*   If SinglePath is TRUE, the entire coverage is going to be a single       M
* curve.								     M
*   If NSrf != NULL, every second curve will be a vector field curve         M
* representing the unnormalized normal for the previous Euclidean curve.     M
* This mode disable the SinglePath mode.				     M
*   See also function SymbSetAdapIsoExtractMinLevel.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To compute adaptive isocurve coverage form                   M
*   NSrf:       Normal vector field defining the normals of Srf.             M
*   AdapIsoDistFunc: Optional function to invoke with the two adjacent       M
*		isoparametric curves of the coverage to evaluate the         M
*		distance between them.					     M
*   Dir:        Direction of adaptive isocurve extraction. Either U or V.    M
*   Eps:        Tolerance of adaptive isocurve cuverage. For every point P   M
*               on Srf there will be a point Q in one of the extracted       M
*               isocurves such the |P - Q| < Eps.			     M
*   FullIso:    Do we want all isocurves to span the entire domain?          M
*   SinglePath: Do we want a single curve through them all?                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A list of curves representing the computed adaptive    M
*                     isocurve coverage for surface Srf. If normal field,    M
*                     NSrf, is prescribed, normal curves are concatenated    M
*                     alternatingly in this list.                            M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbAdapIsoExtract, adaptive isocurves                                   M
*****************************************************************************/
CagdCrvStruct *SymbAdapIsoExtract(CagdSrfStruct *Srf,
				  CagdSrfStruct *NSrf,
				  SymbAdapIsoDistSqrFuncType AdapIsoDistFunc,
				  CagdSrfDirType Dir,
				  CagdRType Eps,
				  CagdBType FullIso,
				  CagdBType SinglePath)
{
    CagdBType
	SrfBezier = FALSE;
    CagdRType Crv1Param, Crv2Param;
    CagdCrvStruct *Crv1, *Crv2, *NCrv1, *NCrv2, *AllAdapIso, *TCrv;

    if (NSrf != NULL)
	SinglePath = FALSE;

    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	    Srf = CnvrtBezier2BsplineSrf(Srf);
	    SrfBezier = TRUE;
	    break;
	case CAGD_SBSPLINE_TYPE:
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_WRONG_SRF);
	    break;
    }

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    Crv1Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> UKnotVector[0] + IRIT_EPS : IRIT_EPS;
	    Crv2Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> UKnotVector[Srf -> ULength +
				   Srf -> UOrder - 1] - IRIT_EPS
		: 1.0 - IRIT_EPS;
	    break;
	case CAGD_CONST_V_DIR:
	    Crv1Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> VKnotVector[0] + IRIT_EPS : IRIT_EPS;
	    Crv2Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> VKnotVector[Srf -> VLength +
				   Srf -> VOrder - 1] - IRIT_EPS
		: 1.0 - IRIT_EPS;
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
	    Crv1Param = 0.0;
	    Crv2Param = 1.0;
	    break;
    }

    Crv1 = CagdCrvFromSrf(Srf, Crv1Param, Dir);
    Crv2 = CagdCrvFromSrf(Srf, Crv2Param, Dir);
    if (NSrf != NULL) {
	NCrv1 = CagdCrvFromSrf(NSrf, Crv1Param, Dir);
	NCrv2 = CagdCrvFromSrf(NSrf, Crv2Param, Dir);
    }
    else
	NCrv1 = NCrv2 = NULL;

    /* Compute the adaptive iso curves. */
    AllAdapIso = SymbAdapIsoExtractAux(0, Srf, NSrf, AdapIsoDistFunc,
				       Crv1, NCrv1, Crv2, NCrv2,
				       Crv1Param, Crv2Param,
				       Dir, Eps * Eps, FullIso, SinglePath);

    /* Chain first and last iso curves that always span the entire domain. */
    if (AllAdapIso != NULL) {
	Crv1 -> Pnext = AllAdapIso;
	TCrv = CagdListLast(AllAdapIso);
	TCrv -> Pnext = Crv2;
    }
    else
	Crv1 -> Pnext = Crv2;

    if (NSrf != NULL) {
	NCrv1 -> Pnext = Crv1 -> Pnext;
	Crv1 -> Pnext = NCrv1;
	NCrv2 -> Pnext = Crv2 -> Pnext;
	Crv2 -> Pnext = NCrv2;
    }

    if (SrfBezier)
	CagdSrfFree(Srf);

    return Crv1;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* An auxiliary function of SymbAdapIsoExtract. Computes the distance square  *
* between the given two curves, extract the regions that are far than that   *
* and recursively invoke this function with the sub-curves.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   Level:         Of recursion.                                             *
*   Srf:           This adaptive isocurve coverage is computed for.          *
*   NSrf:          Normal vectorfield of Srf.                                *
*   AdapIsoDistFunc: Optional function to invoke with the two adjacent       *
*		   isoparametric curves of the coverage to evaluate the      *
*		   distance between them.				     *
*   Crv1:          First curve to handle.                                    *
*   NCrv1:         Normal vector field of first curve Crv1.                  *
*   Crv2:          Second curve to handle.                                   *
*   NCrv2:         Normal vector field of second curve Crv2.                 *
*   Crv1Param:     Parameter along the other direction of Srf, of Crv1.      *
*   Crv2Param:     Parameter along the other direction of Srf, of Crv2.      *
*   Dir:           Direction of adaptive isocurve extraction.                *
*   Eps2:          Tolerance of adaptive isocurve extraction.                *
*   FullIso:       Do we want all isocurves to span all parametric domain?   *
*   SinglePath:    Do we want everything in a single path?                   *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:  A list of isocurves covering Srf to within Eps2.       *
*****************************************************************************/
static CagdCrvStruct *SymbAdapIsoExtractAux(int Level,
					    CagdSrfStruct *Srf,
					    CagdSrfStruct *NSrf,
					    SymbAdapIsoDistSqrFuncType
					                      AdapIsoDistFunc,
					    CagdCrvStruct *Crv1,
					    CagdCrvStruct *NCrv1,
					    CagdCrvStruct *Crv2,
					    CagdCrvStruct *NCrv2,
					    CagdRType Crv1Param,
					    CagdRType Crv2Param,
					    CagdSrfDirType Dir,
					    CagdRType Eps2,
					    CagdBType FullIso,
					    CagdBType SinglePath)
{
    int i, KVLen;
    CagdPointType PType;
    CagdRType Dist, *KV, LastT, DistCrv2Min, DistCrv2Max,
	**Points, *PtsW, *PtsX, LastDist,
	*Nodes = NULL;
    CagdBType CloseEnough, LastCloseEnough;
    CagdCrvStruct *TCrv, *DistCrv2,
	*AllAdapIsoTail = NULL,
	*AllAdapIso = NULL;

    if (AdapIsoDistFunc != NULL) { /* Call function to compute distance sqr. */
	DistCrv2 = AdapIsoDistFunc(Level, Crv1, NCrv1, Crv2, NCrv2);
    }
    else {	       /* Symbolically compute the distance square function. */
	CagdCrvStruct
	    *DiffCrv = SymbCrvSub(Crv1, Crv2);

	DistCrv2 = SymbCrvDotProd(DiffCrv, DiffCrv);
	CagdCrvFree(DiffCrv);
    }
    PType = DistCrv2 -> PType;

    /* Simple heuristic how much to refine DistCrv2: */
    KVLen = DistCrv2 -> Length * 2;
    DistCrv2Min = DistCrv2 -> KnotVector[0];
    DistCrv2Max = DistCrv2 -> KnotVector[DistCrv2 -> Order +
					 DistCrv2 -> Length - 1];
    KV = (CagdRType *) IritMalloc(sizeof(CagdRType) * KVLen);
    for (i = 0; i < KVLen; i++)
	KV[i] = DistCrv2Min + (i + 1) * (DistCrv2Max - DistCrv2Min) /
								(KVLen + 1);

    TCrv = BspCrvKnotInsertNDiff(DistCrv2, FALSE, KV, KVLen);
    IritFree((VoidPtr) KV);
    CagdCrvFree(DistCrv2);
    DistCrv2 = TCrv;

    Nodes = CagdCrvNodes(DistCrv2);
    LastT = Nodes[0];
    Points = DistCrv2 -> Points,
    PtsW = Points[W],
    PtsX = Points[X],
    LastDist = (PType == CAGD_PT_E1_TYPE ? PtsX[0]
					 : (PtsX[0] / PtsW[0])) - Eps2;
    LastCloseEnough = LastDist < 0.0 && MinSubdivLevel <= Level;

    for (i = 1; i < DistCrv2 -> Length; i++) {
	Dist = (PType == CAGD_PT_E1_TYPE ? PtsX[i]
					 : (PtsX[i] / PtsW[i])) - Eps2;
	CloseEnough = Dist < 0.0 && MinSubdivLevel <= Level;

	if (CloseEnough != LastCloseEnough ||
	    (i == DistCrv2 -> Length - 1 && !LastCloseEnough)) {
	    CagdRType t;
	    int KnotIndex,
	        KVLength = Crv1 -> Length + Crv1 -> Order;

	    if (CloseEnough == LastCloseEnough)
		t = Nodes[DistCrv2 -> Length - 1];
	    else
		CONVEX_BLEND(Nodes[i - 1], Nodes[i], LastDist, Dist, t);

	    /* If parameter t is close "enough" to an existing knot. use it. */
	    KnotIndex = BspKnotLastIndexLE(Crv1 -> KnotVector, KVLength, t);
	    if (KnotIndex >= 0 &&
		t - Crv1 -> KnotVector[KnotIndex] < SIMILAR_PARAM_VAL)
		t = Crv1 -> KnotVector[KnotIndex];
	    KnotIndex = BspKnotFirstIndexG(Crv1 -> KnotVector, KVLength, t);
	    if (KnotIndex < KVLength &&
		Crv1 -> KnotVector[KnotIndex] - t < SIMILAR_PARAM_VAL)
		t = Crv1 -> KnotVector[KnotIndex];

	    if (t - LastT > SIMILAR_PARAM_VAL &&
		(CloseEnough ||
		 (i == DistCrv2 -> Length - 1 && !LastCloseEnough))) {
		/* We are at the end of a region that is not close enough. */
		if (SinglePath) {
		    CagdRType
		        MidParam1 = (Crv1Param * 2.0 + Crv2Param) / 3.0,
			MidParam2 = (Crv1Param + Crv2Param * 2.0) / 3.0;
		    CagdCrvStruct
		        *AdapIso1, *AdapIso2, *AdapIso3, *AdapIso,
			*MidCrv1 = CagdCrvFromSrf(Srf, MidParam1, Dir),
			*MidCrv2 = CagdCrvFromSrf(Srf, MidParam2, Dir);

		    if (FullIso) {
			AdapIso1 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 Crv1, NULL,
							 MidCrv1, NULL,
							 Crv1Param, MidParam1,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 MidCrv1, NULL,
							 MidCrv2, NULL,
							 MidParam1, MidParam2,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso3 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 MidCrv2, NULL,
							 Crv2, NULL,
							 MidParam2, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidCrv1;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidCrv1; 
			if (AdapIso2 != NULL)
			    MidCrv1 -> Pnext = AdapIso2;

			if (AdapIso2 != NULL) {
			    for (TCrv = AdapIso2;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidCrv2;
			}
			else
			    MidCrv1 -> Pnext = MidCrv2;
			if (AdapIso2 != NULL)
			    MidCrv2 -> Pnext = AdapIso3;

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidCrv2;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);

			break; /* Only one (entire domain) region! */
		    }
		    else {
			CagdCrvStruct *Region1, *Region2,
			    *MidRegion1, *MidRegion2;

			Region1 = CopyRegionFromCrv(Crv1, LastT, t),
			Region2 = CopyRegionFromCrv(Crv2, LastT, t);
			MidRegion1 = CopyRegionFromCrv(MidCrv1, LastT, t);
			MidRegion2 = CopyRegionFromCrv(MidCrv2, LastT, t);
			CagdCrvFree(MidCrv1);
			CagdCrvFree(MidCrv2);

			AdapIso1 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 Region1, NULL,
							 MidRegion1, NULL,
							 Crv1Param, MidParam1,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 MidRegion1, NULL,
							 MidRegion2, NULL,
							 MidParam1, MidParam2,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso3 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 MidRegion2, NULL,
							 Region2, NULL,
							 MidParam2, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidRegion1;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidRegion1; 
			if (AdapIso2 != NULL)
			    MidRegion1 -> Pnext = AdapIso2;

			if (AdapIso2 != NULL) {
			    for (TCrv = AdapIso2;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidRegion2;
			}
			else
			    MidRegion1 -> Pnext = MidRegion2;
			if (AdapIso2 != NULL)
			    MidRegion2 -> Pnext = AdapIso3;

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidRegion2;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);
		    }
		}
		else { /* Return all the isocurves as a list of curves. */
		    CagdRType
		        MidParam = ComputeMidParam(Crv1Param, Crv2Param,
						   Dir, Srf);
		    CagdCrvStruct *AdapIso1, *AdapIso2, *AdapIso,
			*MidCrv = CagdCrvFromSrf(Srf, MidParam, Dir),
			*MidNCrv = NSrf ? CagdCrvFromSrf(NSrf, MidParam, Dir)
					: NULL;

		    if (FullIso) {
			AdapIso1 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 Crv1, NCrv1,
							 MidCrv, MidNCrv,
							 Crv1Param, MidParam,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 MidCrv, MidNCrv,
							 Crv2, NCrv2,
							 MidParam, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidCrv;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidCrv;

			if (NSrf != NULL) {
			    MidCrv -> Pnext = MidNCrv;
			    MidCrv = MidNCrv;
			}

			if (AdapIso2 != NULL)
			    MidCrv -> Pnext = AdapIso2;

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidCrv;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);

			break; /* Only one (entire domain) region! */
		    }
		    else {
			CagdCrvStruct *Region1, *Region2, *MidRegion,
			    *NRegion1 = NULL,
			    *NRegion2 = NULL,
			    *MidNRegion = NULL;

			Region1 = CopyRegionFromCrv(Crv1, LastT, t);
			Region2 = CopyRegionFromCrv(Crv2, LastT, t);
			MidRegion = CopyRegionFromCrv(MidCrv, LastT, t);
			if (NSrf) {
			    NRegion1 = CopyRegionFromCrv(NCrv1, LastT, t);
			    NRegion2 = CopyRegionFromCrv(NCrv2, LastT, t);
			    MidNRegion = CopyRegionFromCrv(MidNCrv, LastT,
							      t);
			    CagdCrvFree(MidNCrv);
			}
			CagdCrvFree(MidCrv);

			AdapIso1 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 Region1, NRegion1,
							 MidRegion, MidNRegion,
							 Crv1Param, MidParam,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = SymbAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 AdapIsoDistFunc,
							 MidRegion, MidNRegion,
							 Region2, NRegion2,
							 MidParam, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidRegion;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidRegion;

			if (NSrf != NULL) {
			    MidRegion -> Pnext = MidNRegion;
			    MidRegion = MidNRegion;
			}

			if (AdapIso2 != NULL)
			    MidRegion -> Pnext = AdapIso2;

			CagdCrvFree(Region1);
			CagdCrvFree(Region2);
			if (NSrf != NULL) {
			    CagdCrvFree(NRegion1);
			    CagdCrvFree(NRegion2);
			}

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidRegion;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);
		    }
		}
	    }

	    /* We are at the beginning of a region not close enough. */
	    LastT = t;
	}

	LastDist = Dist;
	LastCloseEnough = CloseEnough;
    }

    CagdCrvFree(DistCrv2);
    IritFree((VoidPtr) Nodes);

    if (SinglePath) { /* Add connecting isocurves into the existing curves. */
	/* Not implemented yet. */
    }

    return AllAdapIso;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Extracts a sub region of the given curve, but also copies its attributes.*
*                                                                            *
* PARAMETERS:                                                                *
*   Crv:        To extracts a region from.                                   *
*   TMin, TMax: Domain of region.                                            *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:   Extracted region.                                     *
*****************************************************************************/
static CagdCrvStruct *CopyRegionFromCrv(CagdCrvStruct *Crv,
					CagdRType TMin,
					CagdRType TMax)
{
    CagdCrvStruct
	*Region = CagdCrvRegionFromCrv(Crv, TMin, TMax);

    Region -> Attr = AttrCopyAttributes(Crv -> Attr);

    return Region;
}


/*****************************************************************************
* DESCRIPTION:                                                               *
*   Estimates a parameter between Param1 and Param2 inSrf in Dir. Selects    *
* interior knots before selecting an average between parameters.             *
*                                                                            *
* PARAMETERS:                                                                *
*   Param1, Param2:  Two parameters to find in between.                      *
*   Dir:             Direction in Srf of Param1/2.                           *
*   Srf:             Surface to explore.                                     *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdRType:       A parameter between Param1 and Param2.                  *
*****************************************************************************/
static CagdRType ComputeMidParam(CagdRType Param1,
				 CagdRType Param2,
				 CagdSrfDirType Dir,
				 CagdSrfStruct *Srf)
{
    CagdRType t, *KV;
    int KVLen, Index1, Index2;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    KV = Srf -> UKnotVector;
	    KVLen = Srf -> UOrder + Srf -> ULength;
	    break;
	case CAGD_CONST_V_DIR:
	    KV = Srf -> VKnotVector;
	    KVLen = Srf -> VOrder + Srf -> VLength;
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
	    KV = NULL;
	    KVLen = 0;
    }

    Index1 = BspKnotLastIndexLE(KV, KVLen, Param1);
    Index2 = BspKnotLastIndexLE(KV, KVLen, Param2);

    if (Index1 < 0 || Index2 < 0 || Index1 >= KVLen || Index2 >= KVLen) {
	SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE);
	return -IRIT_INFNTY;
    }

    t = KV[(Index1 + Index2) / 2];
    if (t < Param1 || t > Param2 || APX_EQ(t, Param1) || APX_EQ(t, Param2))
	t = (Param1 + Param2) / 2.0;

    return t;
}