/******************************************************************************
* Offset.c - computes offset approximation to curves and surfaces.            *
*******************************************************************************
* Written by Gershon Elber, March. 93.					      *
******************************************************************************/

#include "symb_loc.h"
#include "extra_fn.h"

#define MAX_OFFSET_IMPROVE_ITERS	20
#define NORMAL_PERTURB			1e-4
#define OFFSET_TRIM_EPS			1.25

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve and an offset amount OffsetDist, returns an approximation to M
* the offset curve by offseting the control polygon in the normal direction. M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:          To approximate its offset curve with distance OffsetDist.  M
*   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
*   BezInterp:    If TRUE, control points are interpolated when the curve is M
*		  reduced to a Bezier form. Otherwise, control points are    M
*		  translated OffsetDist amount only, under estimating the    M
*		  Offset. 						     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   An approximation to the offset curve.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvOffset, offset                                                    M
*****************************************************************************/
CagdCrvStruct *SymbCrvOffset(CagdCrvStruct *Crv,
			     CagdRType OffsetDist,
			     CagdBType BezInterp)
{
    CagdBType
	IsBezierCurve = FALSE,
	IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    int i, j,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType),
	Order = Crv -> Order,
	Length = Crv -> Length;
    CagdBType
	HasNewKV = FALSE;
    CagdVecStruct *N;
    CagdCrvStruct
	*OffsetCrv = CagdCrvCopy(Crv);
    CagdRType *Nodes, *NodePtr,
	**Points = OffsetCrv -> Points,
	*KV = NULL;

    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	    HasNewKV = TRUE;
	    KV = BspKnotUniformOpen(Length, Order, NULL);
	    IsBezierCurve = TRUE;
	    break;
	case CAGD_CBSPLINE_TYPE:
	    KV = Crv -> KnotVector;
	    IsBezierCurve = Crv -> Length == Crv -> Order;
	    break;
	case CAGD_CPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    return NULL;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    return NULL;
    }

    Nodes = BspKnotNodes(KV, Length + Order, Order);
    NodePtr = Nodes;

    /* Interpolate the compute control points instead of under estimating */
    /* This offset.						          */
    if (BezInterp && IsBezierCurve) {
	CagdCrvStruct *TCrv;
	if (IsRational) {
	    TCrv = CagdCoerceCrvTo(OffsetCrv, CAGD_MAKE_PT_TYPE(FALSE,
				   CAGD_NUM_OF_PT_COORD(OffsetCrv -> PType)));

	    CagdCrvFree(OffsetCrv);
	    OffsetCrv = TCrv;
	    Points = OffsetCrv -> Points;
	}

	for (j = 0; j < Length; j++, NodePtr++) {
	    CagdRType
		*R = CagdCrvEval(Crv, *NodePtr);

	    if ((N = CagdCrvNormal(Crv, *NodePtr)) == NULL &&
		(N = CagdCrvNormal(Crv, NodePtr == Nodes ?
				         *NodePtr + NORMAL_PERTURB :
				         *NodePtr - NORMAL_PERTURB)) == NULL) {
		SYMB_FATAL_ERROR(SYMB_ERR_CANNOT_COMP_NORMAL);
		return NULL;
	    }

	    for (i = 1; i <= MaxCoord; i++)
		Points[i][j] = R[i] / (IsRational ? R[0] : 1.0) +
			       N -> Vec[i - 1] * OffsetDist;
	}

	TCrv = CagdCrvCopy(OffsetCrv);
	for (i = 1; i <= MaxCoord; i++)
	    BzrCrvInterp(TCrv -> Points[i], OffsetCrv -> Points[i], Length);

	CagdCrvFree(OffsetCrv);
	OffsetCrv = TCrv;
    }
    else {
	for (j = 0; j < Length; j++, NodePtr++) {
	    if ((N = CagdCrvNormal(Crv, *NodePtr)) == NULL &&
		(N = CagdCrvNormal(Crv, NodePtr == Nodes ?
				         *NodePtr + NORMAL_PERTURB :
				         *NodePtr - NORMAL_PERTURB)) == NULL) {
		SYMB_FATAL_ERROR(SYMB_ERR_CANNOT_COMP_NORMAL);
		return NULL;
	    }

	    for (i = 1; i <= MaxCoord; i++)
		Points[i][j] +=
		    N -> Vec[i - 1] * OffsetDist *
			(IsRational ? Points[W][j] : 1.0);
	}
    }

    if (HasNewKV)
	IritFree((VoidPtr) KV);
    IritFree((VoidPtr) Nodes);

    return OffsetCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve and an offset amount OffsetDist, returns an approximation to M
* the offset curve by offseting the control polygon in the normal direction. M
*   If resulting offset is not satisfying the required tolerance the curve   M
* is subdivided and the algorithm recurses on both parts.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:          To approximate its offset curve with distance OffsetDist.  M
*   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
*   Tolerance:    Accuracy control.                                          M
*   BezInterp:    If TRUE, control points are interpolated when the curve is M
*		  reduced to a Bezier form. Otherwise, control points are    M
*		  translated OffsetDist amount only, under estimating the    M
*		  Offset. 						     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   An approximation to the offset curve, to within       M
*                      Tolerance.					     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvSubdivOffset, offset                                              M
*****************************************************************************/
CagdCrvStruct *SymbCrvSubdivOffset(CagdCrvStruct *Crv,
				   CagdRType OffsetDist,
				   CagdRType Tolerance,
				   CagdBType BezInterp)
{
    CagdCrvStruct
	*OffCrv = SymbCrvOffset(Crv, OffsetDist, BezInterp),
	*Dist = SymbCrvSub(Crv, OffCrv),
	*DistSqr = SymbCrvDotProd(Dist, Dist);
    CagdRType *R, MinVal, MaxVal, TMin, TMax;

    CagdCrvFree(Dist);

    R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> Length, TRUE);
    MinVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);
    R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> Length, FALSE);
    MaxVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);

    CagdCrvFree(DistSqr);

    CagdCrvDomain(Crv, &TMin, &TMax);
    if ((ABS(MinVal - ABS(OffsetDist)) > Tolerance ||
	 ABS(MaxVal - ABS(OffsetDist)) > Tolerance) &&
	TMax - TMin > NORMAL_PERTURB * 10.0) {
	CagdBType NewCrv;
	CagdCrvStruct *Crv1, *Crv2, *OffCrv1, *OffCrv2;

	if (CAGD_IS_BEZIER_CRV(Crv)) {
	    /* Promote to a Bspline so we can keep track of parametrization. */
	    Crv = CnvrtBezier2BsplineCrv(Crv);
	    NewCrv = TRUE;
	}
	else
	    NewCrv = FALSE;

	if (TMax - TMin > NORMAL_PERTURB * 10.0) {
	    CagdCrvFree(OffCrv);

	    Crv1 = CagdCrvSubdivAtParam(Crv, (TMin + TMax) / 2.0);
	    Crv2 = Crv1 -> Pnext;
	    Crv1 -> Pnext = NULL;
	    OffCrv1 = SymbCrvSubdivOffset(Crv1, OffsetDist, Tolerance, BezInterp);
	    OffCrv2 = SymbCrvSubdivOffset(Crv2, OffsetDist, Tolerance, BezInterp);
	    CagdCrvFree(Crv1);
	    CagdCrvFree(Crv2);

	    OffCrv = CagdMergeCrvCrv(OffCrv1, OffCrv2, TRUE);
	    CagdCrvFree(OffCrv1);
	    CagdCrvFree(OffCrv2);
	}

	if (NewCrv)
	    CagdCrvFree(Crv);
    }

    return OffCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a surface and an offset amount OffsetDist, returns an approximation  M
* to the offset surface by offseting the control mesh in the normal	     M
* direction.                                                                 M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:         To approximate its offset surface with distance OffsetDist. M
*   OffsetDist:  Amount of offset. Negative denotes other offset direction.  M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   An approximation to the offset surface.               M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfOffset, offset                                                    M
*****************************************************************************/
CagdSrfStruct *SymbSrfOffset(CagdSrfStruct *Srf, CagdRType OffsetDist)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int	i, Row, Col,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType),
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength;
    CagdBType
	HasNewKV = FALSE;
    CagdVecStruct *N;
    CagdSrfStruct
	*OffsetSrf = CagdSrfCopy(Srf);
    CagdRType *UNodes, *VNodes, *UNodePtr, *VNodePtr,
	**Points = OffsetSrf -> Points,
	*UKV = NULL,
	*VKV = NULL;

    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	    HasNewKV = TRUE;
	    UKV = BspKnotUniformOpen(ULength, UOrder, NULL);
	    VKV = BspKnotUniformOpen(VLength, VOrder, NULL);
	    break;
	case CAGD_SBSPLINE_TYPE:
	    UKV = Srf -> UKnotVector;
	    VKV = Srf -> VKnotVector;
	    break;
	case CAGD_SPOWER_TYPE:
	    SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT);
	    return NULL;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    return NULL;
    }

    UNodes = BspKnotNodes(UKV, ULength + UOrder, UOrder);
    VNodes = BspKnotNodes(VKV, VLength + VOrder, VOrder);

    if (IsNotRational)
	for (Row = 0, VNodePtr = VNodes; Row < VLength; Row++, VNodePtr++)
	    for (Col = 0, UNodePtr = UNodes; Col < ULength; Col++, UNodePtr++) {
	    	N = CagdSrfNormal(Srf, *UNodePtr, *VNodePtr);
	    	for (i = 1; i <= MaxCoord; i++)
		    Points[i][CAGD_MESH_UV(OffsetSrf, Col, Row)] +=
			N -> Vec[i - 1] * OffsetDist;
	    }
    else
	for (Row = 0, VNodePtr = VNodes; Row < VLength; Row++, VNodePtr++)
	    for (Col = 0, UNodePtr = UNodes; Col < ULength; Col++, UNodePtr++) {
	    	N = CagdSrfNormal(Srf, *UNodePtr, *VNodePtr);
	    	for (i = 1; i <= MaxCoord; i++)
	    	    Points[i][CAGD_MESH_UV(OffsetSrf, Col, Row)] +=
			N -> Vec[i - 1] * OffsetDist *
			    Points[W][CAGD_MESH_UV(OffsetSrf, Col, Row)];
	    }

    if (HasNewKV) {
	IritFree((VoidPtr) UKV);
	IritFree((VoidPtr) VKV);
    }

    IritFree((VoidPtr) UNodes);
    IritFree((VoidPtr) VNodes);

    return OffsetSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a surface and an offset amount OffsetDist, returns an approximation  M
* to the offset surface by offseting the control mesh in the normal	     M
* direction.                                                                 M
*   If resulting offset is not satisfying the required tolerance the surface M
* is subdivided and the algorithm recurses on both parts.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:         To approximate its offset surface with distance OffsetDist. M
*   OffsetDist:  Amount of offset. Negative denotes other offset direction.  M
*   Tolerance:    Accuracy control.                                          M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   An approximation to the offset surface, to within     M
*                      Tolerance.				             M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfSubdivOffset, offset                                              M
*****************************************************************************/
CagdSrfStruct *SymbSrfSubdivOffset(CagdSrfStruct *Srf,
				   CagdRType OffsetDist,
				   CagdRType Tolerance)
{
    CagdSrfStruct
	*OffSrf = SymbSrfOffset(Srf, OffsetDist),
	*Dist = SymbSrfSub(Srf, OffSrf),
	*DistSqr = SymbSrfDotProd(Dist, Dist);
    CagdRType *R, MinVal, MaxVal;
    CagdSrfDirType Dir;

    CagdSrfFree(Dist);

    R = SymbExtremumCntPtVals(DistSqr -> Points,
			      DistSqr -> ULength * DistSqr -> VLength, TRUE);
    MinVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);
    R = SymbExtremumCntPtVals(DistSqr -> Points,
			      DistSqr -> ULength * DistSqr -> VLength, FALSE);
    MaxVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]);

    CagdSrfFree(DistSqr);

    if (ABS(MinVal - ABS(OffsetDist)) > Tolerance ||
	ABS(MaxVal - ABS(OffsetDist)) > Tolerance) {
	CagdSrfStruct *Srf1, *Srf2, *OffSrf1, *OffSrf2;
	CagdRType UMin, UMax, VMin, VMax, t;
	CagdBType NewSrf;

	/* Make sure it is a Bspline surface. */
	if (CAGD_IS_BEZIER_SRF(Srf)) {
	    /* Promote to a Bspline so we can keep track of parametrization. */
	    Srf = CnvrtBezier2BsplineSrf(Srf);
	    NewSrf = TRUE;
	}
	else
	    NewSrf = FALSE;

	CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
	if (MAX(UMax - UMin, VMax - VMin) > NORMAL_PERTURB * 10.0) {
	    CagdSrfFree(OffSrf);

	    if (UMax - UMin > VMax - VMin) {
		Dir = CAGD_CONST_U_DIR;
		t = (UMin + UMax) / 2.0;
	    }
	    else {
		Dir = CAGD_CONST_V_DIR;
		t = (VMin + VMax) / 2.0;
	    }
	    Srf1 = CagdSrfSubdivAtParam(Srf, t, Dir);
	    Srf2 = Srf1 -> Pnext;
	    Srf1 -> Pnext = NULL;
	    OffSrf1 = SymbSrfSubdivOffset(Srf1, OffsetDist, Tolerance);
	    OffSrf2 = SymbSrfSubdivOffset(Srf2, OffsetDist, Tolerance);
	    CagdSrfFree(Srf1);
	    CagdSrfFree(Srf2);

	    OffSrf = CagdMergeSrfSrf(OffSrf1, OffSrf2, Dir, TRUE, TRUE);
	    CagdSrfFree(OffSrf1);
	    CagdSrfFree(OffSrf2);
	}

	if (NewSrf)
	    CagdSrfFree(Srf);
    }

    return OffSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve and an offset amount OffsetDist, returns an approximation to M
* the offset curve by offseting the control polygon in the normal direction. M
*   This function computes an approximation to the offset using		     M
* OffsetAprxFunc, measure the error and use it to refine and decrease the    M
* error adaptively.							     M
*   Bezier curves are promoted to Bsplines curves.			     M
*   See also: Gershon Elber and Elaine Cohen, "Error Bounded Variable	     M
* Distance Offset Operator for Free Form Curves and Surfaces". International M
* Journal of Computational Geometry & Applications, Vol. 1, Num. 1, March    M
* 1991, pp 67-78.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   OrigCrv:      To approximate its offset curve with distance OffsetDist.  M
*   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
*   OffsetError:  Tolerance control.                                         M
*   OffsetAprxFunc:  A function that can be used to approximate an offset    M
*                    of a curve. If NULL SymbCrvOffset function is selected. M
*   BezInterp:    If TRUE, control points are interpolated when the curve is M
*		  reduced to a Bezier form. Otherwise, control points are    M
*		  translated OffsetDist amount only, under estimating the    M
*		  Offset. 						     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   An approximation to the offset curve, to within       M
*                      OffsetError.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvAdapOffset, offset                                                M
*****************************************************************************/
CagdCrvStruct *SymbCrvAdapOffset(CagdCrvStruct *OrigCrv,
				 CagdRType OffsetDist,
				 CagdRType OffsetError,
				 SymbOffCrvFuncType OffsetAprxFunc,
				 CagdBType BezInterp)
{
    int i;
    CagdBType
	IsRational = CAGD_IS_RATIONAL_CRV(OrigCrv);
    CagdRType Min, Max, TMin, TMax,
	OffsetDistSqr = SQR(OffsetDist);
    CagdCrvStruct *OffsetCrv, *Crv;

    switch (OrigCrv -> GType) {
	case CAGD_CBEZIER_TYPE:
	    Crv = CnvrtBezier2BsplineCrv(OrigCrv);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    Crv = CagdCrvCopy(OrigCrv);
	    break;
	case CAGD_CPOWER_TYPE:
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    Crv = NULL;
	    break;
    }

    if (OffsetAprxFunc == NULL)
	OffsetAprxFunc = SymbCrvOffset;

    CagdCrvDomain(Crv, &TMin, &TMax);

    for (i = 0; i < MAX_OFFSET_IMPROVE_ITERS; i++) {
	CagdCrvStruct *DiffCrv, *DistSqrCrv;

	OffsetCrv = OffsetAprxFunc(Crv, OffsetDist, BezInterp);
	DiffCrv = SymbCrvSub(OffsetCrv, Crv);
	DistSqrCrv = SymbCrvDotProd(DiffCrv, DiffCrv);
	CagdCrvFree(DiffCrv);

	CagdCrvMinMax(DistSqrCrv, 1, &Min, &Max);

	if (OffsetDistSqr - Min < OffsetError &&
	    Max - OffsetDistSqr < OffsetError) {
	    /* Error is within bounds - returns this offset approximation. */
	    CagdCrvFree(DistSqrCrv);
	    break;
	}
	else {
	    /* Refine in regions where the error is too large. */
	    int j, k,
	        Length = DistSqrCrv -> Length,
	        Order = DistSqrCrv -> Order,
	        KVLen = Length + Order;
	    CagdRType
		*KV = DistSqrCrv -> KnotVector,
		*Nodes = BspKnotNodes(KV, KVLen, Order),
		*RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * 2 * Length);

	    for (j = k = 0; j < Length; j++) {
		CagdRType
		    *Pt = CagdCrvEval(DistSqrCrv, Nodes[j]),
		    V = OffsetDistSqr - (IsRational ? Pt[1] / Pt[0] : Pt[1]);

		if (ABS(V) > OffsetError) {
		    int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]);

		    if (APX_EQ(KV[Index], Nodes[j])) {
			if (j > 0)
			    RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0;
			if (j < Length - 1)
			    RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0;
		    }
		    else
			RefKV[k++] = Nodes[j];
		}
	    }
	    CagdCrvFree(DistSqrCrv);
	    IritFree((VoidPtr) Nodes);

	    if (k == 0) {
		/* No refinement was found needed - return current curve. */
		IritFree((VoidPtr) RefKV);
		break;
	    }
	    else {
		CagdCrvStruct
		    *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k);

		IritFree((VoidPtr) RefKV);
		CagdCrvFree(Crv);
		Crv = CTmp;
	    }
	}
    }

    CagdCrvFree(Crv);
    return OffsetCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Same function as CagdCrvAdapOffset, but trims the self intersection loops. M
*   See also: Gershon Elber and Elaine Cohen, "Error Bounded Variable	     M
* Distance Offset Operator for Free Form Curves and Surfaces". International M
* Journal of Computational Geometry & Applications, Vol. 1, Num. 1, March    M
* 1991, pp 67-78.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   OrigCrv:      To approximate its offset curve with distance OffsetDist.  M
*   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
*   OffsetError:  Tolerance control.                                         M
*   OffsetAprxFunc:  A function that can be used to approximate an offset    M
*                    of a curve. If NULL SymbCrvOffset function is selected. M
*		     Third parameter of SymbOffCrvFuncType is optional.	     M
*   BezInterp:    If TRUE, control points are interpolated when the curve is M
*		  reduced to a Bezier form. Otherwise, control points are    M
*		  translated OffsetDist amount only, under estimating the    M
*		  Offset. 						     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   An approximation to the offset curve, to within       M
*                      OffsetError.                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvAdapOffsetTrim, offset                                            M
*****************************************************************************/
CagdCrvStruct *SymbCrvAdapOffsetTrim(CagdCrvStruct *OrigCrv,
				     CagdRType OffsetDist,
				     CagdRType OffsetError,
				     SymbOffCrvFuncType OffsetAprxFunc,
				     CagdBType BezInterp)
{
    CagdBType
        IsRational = CAGD_IS_RATIONAL_CRV(OrigCrv);
    CagdCrvStruct *CrvtrSqr, *CrvtrSign, *Crv, *Crvs, *LastCrv,
	*CrvOffsetList = NULL;
    CagdPtStruct *Pts, *PtsHead, *LastPts;
    CagdRType
        OffsetDistSqr1 = 1.0 / SQR(OffsetDist);

    switch (OrigCrv -> GType) {
	case CAGD_CBEZIER_TYPE:
	    Crv = CnvrtBezier2BsplineCrv(OrigCrv);
	    break;
	case CAGD_CBSPLINE_TYPE:
	    Crv = CagdCrvCopy(OrigCrv);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV);
	    Crv = NULL;
	    break;
    }

    if (OffsetDist == 0.0)
	return Crv;

    CrvtrSqr = SymbCrv2DCurvatureSqr(Crv);
    CrvtrSign = SymbCrv2DCurvatureSign(Crv);

    PtsHead = SymbCrvConstSet(CrvtrSqr, 1, SQR(OffsetError),
			      OffsetDistSqr1 / SQR(OFFSET_TRIM_EPS));

    for (Pts = PtsHead, LastPts = NULL; Pts != NULL;) {
	CagdRType *CrvtrSignPt, CrvtrSignVal;

	CrvtrSignPt = CagdCrvEval(CrvtrSign, Pts -> Pt[0]);
	CrvtrSignVal = IsRational ? CrvtrSignPt[1] / CrvtrSignPt[0]
				  : CrvtrSignPt[1];

	if (CrvtrSignVal * OffsetDist > 0.0) {
	    /* Remove point from list. */
	    if (LastPts != NULL) {
		LastPts -> Pnext = Pts -> Pnext;
		CagdPtFree(Pts);
		Pts = LastPts -> Pnext;
	    }
	    else {
		PtsHead = PtsHead -> Pnext;
		CagdPtFree(Pts);
		Pts = PtsHead;
	    }
	}
	else {
	    LastPts = Pts;
	    Pts = Pts -> Pnext;
	}
    }

    for (Pts = PtsHead;
	 Pts != NULL || Crv != NULL;
	 Pts = (Pts != NULL ? Pts -> Pnext : NULL)) {
	CagdCrvStruct
	    *Crv1 = Pts ? CagdCrvSubdivAtParam(Crv, Pts -> Pt[0])
			: CagdCrvCopy(Crv),
	    *Crv2 = Pts ? Crv1 -> Pnext : NULL;
	CagdRType Min, Max, *CrvtrSqrPt, CrvtrSqrVal,
						*CrvtrSignPt, CrvtrSignVal;

	CagdCrvDomain(Crv1, &Min, &Max);

	CrvtrSqrPt = CagdCrvEval(CrvtrSqr, (Min + Max) / 2);
	CrvtrSqrVal = CrvtrSqrPt[1] / CrvtrSqrPt[0];
	CrvtrSignPt = CagdCrvEval(CrvtrSign, (Min + Max) / 2);
	CrvtrSignVal = IsRational ? CrvtrSignPt[1] / CrvtrSignPt[0]
				  : CrvtrSignPt[1];

	if (CrvtrSqrVal < OffsetDistSqr1 / SQR(OFFSET_TRIM_EPS) ||
	    CrvtrSignVal * OffsetDist > 0.0) {
	    CagdCrvStruct
		*Crv1Off = SymbCrvAdapOffset(Crv1, OffsetDist, OffsetError,
					     OffsetAprxFunc, BezInterp);

	    LIST_PUSH(Crv1Off, CrvOffsetList);
	}

	CagdCrvFree(Crv1);
	CagdCrvFree(Crv);
	Crv = Crv2;
    }

    Crvs = CagdListReverse(CrvOffsetList);
    CrvOffsetList = NULL;
    LastCrv = Crvs;
    Crvs = Crvs -> Pnext;
    for (Crv = Crvs; Crv != NULL; Crv = Crv -> Pnext) {
	CagdCrvStruct *CTmp, *CTmp2;
	CagdSrfStruct
	    *DistSrf = SymbSrfDistCrvCrv(LastCrv, Crv);
	CagdPtStruct
	    *IPts = SymbSrfDistFindPoints(DistSrf, OffsetError, FALSE);

	CagdSrfFree(DistSrf);

	if (IPts != NULL) {
	    if (IPts -> Pnext) {
		SYMB_FATAL_ERROR(SYMB_ERR_TOO_COMPLEX);
		return NULL;
	    }
	    else {
		CagdRType TMin, TMax;

		CagdCrvDomain(LastCrv, &TMin, &TMax);
		CTmp = CagdCrvRegionFromCrv(LastCrv, TMin, IPts -> Pt[0]);
		CagdCrvFree(LastCrv);
		LastCrv = CTmp;

		CagdCrvDomain(Crv, &TMin, &TMax);
		CTmp = CagdCrvRegionFromCrv(Crv, IPts -> Pt[1], TMax);

		CTmp2 = CagdMergeCrvCrv(LastCrv, CTmp, FALSE);
		CagdCrvFree(CTmp);
		CagdCrvFree(LastCrv);
		LastCrv = CTmp2;
	    }
	}
	else {
	    /* Simply chain the pieces together. */
	    CTmp = CagdMergeCrvCrv(LastCrv, Crv, FALSE);
	    CagdCrvFree(LastCrv);
	    LastCrv = CTmp;
	}
    }

    CagdPtFreeList(PtsHead);
    CagdCrvFree(CrvtrSqr);
    CagdCrvFree(CrvtrSign);

    return LastCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a curve and an offset amount OffsetDist, returns an approximation to M
* the offset curve by least square fitting a curve to samples taken on the   M
* offset curve.								     M
*   Resulting curve of order Order (degree of Crv if Order == 0) will have   M
* NumOfDOF control points that least sqaure fit NumOfSamples samples on the  M
* offset curve.								     M
*   Tolerance will be updated to hold an error distance measure.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:          To approximate its offset curve with distance OffsetDist.  M
*   OffsetDist:   Amount of offset. Negative denotes other offset direction. M
*   NumOfSamples: Number of samples to sample the offset curve at.           M
*   NumOfDOF:     Number of degrees of freedom on the newly computed offset  M
*		  approximation. This is thesame as the number of control    M
*		  points the new curve will have.			     M
*   Order:        Of the newly constructed offset approximation. If equal to M
*		  zero, the order of Crv will be used.			     M
*   Tolerance:    To return an error estimate in the L-infinity norm.        M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   An approximation to the offset curve.		     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvLeastSquarOffset, offset                                          M
*****************************************************************************/
CagdCrvStruct *SymbCrvLeastSquarOffset(CagdCrvStruct *Crv,
				       CagdRType OffsetDist,
				       int NumOfSamples,
				       int NumOfDOF,
				       int Order,
				       CagdRType *Tolerance)
{
    int i;
    CagdRType TMin, TMax, t, dt;
    CagdVType Tangent;
    CagdCrvStruct *OffApprox,
	*TangentCrv = CagdCrvDerive(Crv);
    CagdPtStruct
	*Pt = NULL,
	*PtList = NULL;

    /* Create NumOfSamples samples on the offset curve. */
    CagdCrvDomain(Crv, &TMin, &TMax);
    dt = (TMax - TMin) / (NumOfSamples - 1);
    for (i = 0, t = TMin; i < NumOfSamples; i++, t += dt) {
	CagdRType *R;

	if (PtList == NULL)
	    PtList = Pt = CagdPtNew();
	else {
	    Pt -> Pnext = CagdPtNew();
	    Pt = Pt -> Pnext;
	}

	R = CagdCrvEval(Crv, t);
	CagdCoerceToE3(Pt -> Pt, &R, -1, Crv -> PType);

	R = CagdCrvEval(TangentCrv, t);
	CagdCoerceToE2(Tangent, &R, -1, TangentCrv -> PType);
	Tangent[2] = 0.0;
	PT_NORMALIZE(Tangent);
	
	Pt -> Pt[0] += Tangent[1] * OffsetDist;
	Pt -> Pt[1] -= Tangent[0] * OffsetDist;
    }

    OffApprox = BspCrvInterpPts(PtList,
				Order == 0 ? Crv -> Order : Order,
				MIN(NumOfDOF, NumOfSamples),
				CAGD_UNIFORM_PARAM);

    *Tolerance = BspCrvInterpPtsError(OffApprox, PtList, CAGD_UNIFORM_PARAM);

    CagdPtFreeList(PtList);
    CagdCrvFree(TangentCrv);

    return OffApprox;
}
