/******************************************************************************
* Curvature.c - curvature computation of curves and surfaces.		      *
*******************************************************************************
* Written by Gershon Elber, March 93.					      *
******************************************************************************/

#include "symb_loc.h"

#define MAX_POS_REF_ITERATION 20

/******************************************************************************
* DESCRIPTION:                                                               M
* Computes a scalar curve representing the curvature of a planar curve.      M
*   The given curve is assumed to be planar and only its x and y coordinates M
* are considered.							     M
*   Then the curvature k is equal to					     M
*      .  ..    .  ..							     V
*      X  Y  -  Y  X							     V
* k =  -------------							     V
*         .2   .2  3/2 							     V
*      (  X  + Y  )							     V
*   Since we cannot represent k because of the square root, we compute and   M
* represent k^2.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       To compute the square of teh curvature field for.             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   The square of the curvature field of Crv.             M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv2DCurvatureSqr, curvature                                         M
*****************************************************************************/
CagdCrvStruct *SymbCrv2DCurvatureSqr(CagdCrvStruct *Crv)
{
    CagdBType
        IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
		  *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
		  *CTmp1, *CTmp2, *CTmp3, *Numer, *Denom, *CurvatureSqr,
	*Crv1Deriv = CagdCrvDerive(Crv),
	*Crv2Deriv = CagdCrvDerive(Crv1Deriv);

    SymbCrvSplitScalar(Crv1Deriv, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
    SymbCrvSplitScalar(Crv2Deriv, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);
    CagdCrvFree(Crv1Deriv);
    CagdCrvFree(Crv2Deriv);

    CTmp1 = SymbCrvMult(Crv1X, Crv2Y);
    CTmp2 = SymbCrvMult(Crv2X, Crv1Y);
    CTmp3 = SymbCrvSub(CTmp1, CTmp2);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp2);
    Numer = SymbCrvMult(CTmp3, CTmp3);
    CagdCrvFree(CTmp3);

    CTmp1 = SymbCrvMult(Crv1X, Crv1X);
    CTmp2 = SymbCrvMult(Crv1Y, Crv1Y);
    CTmp3 = SymbCrvAdd(CTmp1, CTmp2);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp2);
    CTmp1 = SymbCrvMult(CTmp3, CTmp3);
    Denom = SymbCrvMult(CTmp1, CTmp3);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp3);
    if (IsRational) {
	CTmp1 = SymbCrvMult(Crv1W, Crv1W);
	CTmp2 = SymbCrvMult(CTmp1, CTmp1);
	CTmp3 = SymbCrvMult(CTmp2, Numer);
	CagdCrvFree(CTmp1);
	CagdCrvFree(CTmp2);
	CagdCrvFree(Numer);
	Numer = CTmp3;

	CTmp1 = SymbCrvMult(Crv2W, Crv2W);
	CTmp2 = SymbCrvMult(CTmp1, Denom);
	CagdCrvFree(CTmp1);
	CagdCrvFree(Denom);
	Denom = CTmp2;
    }
    if (CAGD_IS_BSPLINE_CRV(Denom)) {
	CTmp1 = SymbMakePosCrvCtlPolyPos(Denom);
	CagdCrvFree(Denom);
	Denom = CTmp1;
    }

    CagdMakeCrvsCompatible(&Denom, &Numer, TRUE, TRUE);
    CurvatureSqr = SymbCrvMergeScalar(Denom, Numer, NULL, NULL);
    CagdCrvFree(Denom);
    CagdCrvFree(Numer);

    CagdCrvFree(Crv1X);
    CagdCrvFree(Crv1Y);
    CagdCrvFree(Crv2X);
    CagdCrvFree(Crv2Y);
    if (Crv1Z)
	CagdCrvFree(Crv1Z);
    if (Crv2Z)
	CagdCrvFree(Crv2Z);
    if (Crv1W)
	CagdCrvFree(Crv1W);
    if (Crv2W)
	CagdCrvFree(Crv2W);

    return CurvatureSqr;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Computes a vector field curve representing the curvature of a curve, in    M
* the binormal direction, that is kB, and square it:			     M
*									     M
*       .   .. 			        .   ..	2			     V
*       C x C 			      ( C x C )				     V
* kB =  -----  		and	k^2 = --------				     V
*         .  3 			          .  6				     V
*       | C |  			        | C |				     V
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To compute vector field curvatrue for,                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Computed vector field curvature of Crv                M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv3DCurvatureSqr, curvature                                         M
*****************************************************************************/
CagdCrvStruct *SymbCrv3DCurvatureSqr(CagdCrvStruct *Crv)
{
    CagdCrvStruct
	*CTmp, *CTmp2, *Numer, *Denom, *CurvatureSqr,
	*Crv1Deriv = CagdCrvDerive(Crv),
	*Crv2Deriv = CagdCrvDerive(Crv1Deriv);

    CTmp = SymbCrvCrossProd(Crv1Deriv, Crv2Deriv);
    CagdCrvFree(Crv2Deriv);
    Numer = SymbCrvDotProd(CTmp, CTmp);
    CagdCrvFree(CTmp);

    CTmp = SymbCrvDotProd(Crv1Deriv, Crv1Deriv);
    CagdCrvFree(Crv1Deriv);
    CTmp2 = SymbCrvMult(CTmp, CTmp);
    Denom = SymbCrvMult(CTmp, CTmp2);
    CagdCrvFree(CTmp);
    CagdCrvFree(CTmp2);

    if (CAGD_IS_RATIONAL_CRV(Denom) || CAGD_IS_RATIONAL_CRV(Numer)) {
	CTmp = SymbCrvInvert(Denom);
	CurvatureSqr = SymbCrvMult(CTmp, Numer);
	CagdCrvFree(CTmp);
    }
    else {
	CagdCrvStruct *PCrvW, *PCrvX, *PCrvY, *PCrvZ;

	SymbCrvSplitScalar(Numer, &PCrvW, &PCrvX, &PCrvY, &PCrvZ);
	CagdMakeCrvsCompatible(&Denom, &PCrvX, TRUE, TRUE);
	CagdMakeCrvsCompatible(&Denom, &PCrvY, TRUE, TRUE);
	CagdMakeCrvsCompatible(&Denom, &PCrvZ, TRUE, TRUE);
	CurvatureSqr = SymbCrvMergeScalar(Denom, PCrvX, PCrvY, PCrvZ);
	CagdCrvFree(PCrvX);
	CagdCrvFree(PCrvY);
	CagdCrvFree(PCrvZ);
    }

    CagdCrvFree(Denom);
    CagdCrvFree(Numer);

    return CurvatureSqr;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Computes a vector field curve representing the radius (1/curvature) of a   M
* curve, in the normal direction, that is N / k:			     M
*									     M
*                   .   ..    .	     .  6        .   ..    .     .  2	     V
*          k N     (C x C ) x C	   | C |      ( (C x C ) x C ) | C |	     V
* N / k =  ----- = ------------ . --------- = -----------------------	     V
*            2          .  4       .   .. 2           .   .. 2		     V
*           k         | C |       (C x C )           (C x C )		     V
*                                                                            M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To compute the normal field with radius as magnitude.        M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Computed normal field with 1 / k as magnitude.        M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv3DRadiusNormal, curvature                                         M
*****************************************************************************/
CagdCrvStruct *SymbCrv3DRadiusNormal(CagdCrvStruct *Crv)
{
    CagdCrvStruct
	*PCrvW, *PCrvX, *PCrvY, *PCrvZ,
	*CTmp, *CTmp2, *Numer, *Denom, *RadiusNormal,
	*Crv1Deriv = CagdCrvDerive(Crv),
	*Crv2Deriv = CagdCrvDerive(Crv1Deriv);

    CTmp = SymbCrvCrossProd(Crv1Deriv, Crv2Deriv);
    CagdCrvFree(Crv2Deriv);
    Denom = SymbCrvDotProd(CTmp, CTmp);
    CTmp2 = SymbCrvCrossProd(CTmp, Crv1Deriv);
    CagdCrvFree(CTmp);
    CTmp = SymbCrvDotProd(Crv1Deriv, Crv1Deriv);
    Numer = SymbCrvMultScalar(CTmp2, CTmp);
    CagdCrvFree(CTmp2);
    CagdCrvFree(CTmp);

    if (CAGD_IS_RATIONAL_CRV(Denom) || CAGD_IS_RATIONAL_CRV(Numer)) {
	CTmp = SymbCrvInvert(Denom);
	RadiusNormal = SymbCrvMult(CTmp, Numer);
	CagdCrvFree(CTmp);
    }
    else {
	SymbCrvSplitScalar(Numer, &PCrvW, &PCrvX, &PCrvY, &PCrvZ);
	CagdMakeCrvsCompatible(&Denom, &PCrvX, TRUE, TRUE);
	CagdMakeCrvsCompatible(&Denom, &PCrvY, TRUE, TRUE);
	CagdMakeCrvsCompatible(&Denom, &PCrvZ, TRUE, TRUE);
	RadiusNormal = SymbCrvMergeScalar(Denom, PCrvX, PCrvY, PCrvZ);
	CagdCrvFree(PCrvX);
	CagdCrvFree(PCrvY);
	CagdCrvFree(PCrvZ);
    }

    CagdCrvFree(Denom);
    CagdCrvFree(Numer);

    return RadiusNormal;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes a vector field curve representing the curvature of a curve, in    M
* the normal direction, that is kN.					     M
*                .   ..      .       .   ..     .			     V
*                C x C       C     ( C x C  ) x C     			     V
* kN = kB x T =  -----  x  ----- = --------------			     V
*                  .  3      .           .  4				     V
*                | C |     | C |       | C |				     V
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To compute the normal curvature field.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Computed normal curvature field.                      M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv3DCurvatureNormal, curvature                                      M
*****************************************************************************/
CagdCrvStruct *SymbCrv3DCurvatureNormal(CagdCrvStruct *Crv)
{
    CagdBType
        IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ,
	*CTmp, *Numer, *Denom, *CurvatureNormal,
	*Crv1Deriv = CagdCrvDerive(Crv),
	*Crv2Deriv = CagdCrvDerive(Crv1Deriv);

    CTmp = SymbCrvCrossProd(Crv1Deriv, Crv2Deriv);
    CagdCrvFree(Crv2Deriv);
    Numer = SymbCrvCrossProd(CTmp, Crv1Deriv);
    CagdCrvFree(CTmp);
    SymbCrvSplitScalar(Numer, &CrvW, &CrvX, &CrvY, &CrvZ);
    CagdCrvFree(Numer);

    CTmp = SymbCrvDotProd(Crv1Deriv, Crv1Deriv);
    CagdCrvFree(Crv1Deriv);
    Denom = SymbCrvMult(CTmp, CTmp);

    if (IsRational) {
	CagdCrvStruct *DenomCrvW, *DenomCrvX, *DenomCrvY, *DenomCrvZ;

	SymbCrvSplitScalar(Denom,
			   &DenomCrvW, &DenomCrvX, &DenomCrvY, &DenomCrvZ);
	CagdCrvFree(Denom);

	CTmp = SymbCrvMult(CrvW, DenomCrvX);
	CagdCrvFree(CrvW);
	CrvW = CTmp;

	CTmp = SymbCrvMult(CrvX, DenomCrvW);
	CagdCrvFree(CrvX);
	CrvX = CTmp;

	CTmp = SymbCrvMult(CrvY, DenomCrvW);
	CagdCrvFree(CrvY);
	CrvY = CTmp;

	CTmp = SymbCrvMult(CrvZ, DenomCrvW);
	CagdCrvFree(CrvZ);
	CrvZ = CTmp;

	CagdCrvFree(DenomCrvW);
	CagdCrvFree(DenomCrvX);

	CurvatureNormal = SymbCrvMult(CTmp, Numer);
	CagdCrvFree(CTmp);
    }
    else {
	CagdMakeCrvsCompatible(&Denom, &CrvX, TRUE, TRUE);
	CagdMakeCrvsCompatible(&Denom, &CrvY, TRUE, TRUE);
	CagdMakeCrvsCompatible(&Denom, &CrvZ, TRUE, TRUE);
	CrvW = Denom;
    }
    CurvatureNormal = SymbCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ);
    CagdCrvFree(CrvX);
    CagdCrvFree(CrvY);
    CagdCrvFree(CrvZ);
    CagdCrvFree(CrvW);

    return CurvatureNormal;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Computes a scalar curve representing the curvature sign of a planar curve. M
*   The given curve is assumed to be planar and only its x and y coordinates M
* are considered.							     M
*   Then the curvature sign is equal to					     M
*      .  ..    .  ..							     V
* s =  X  Y  -  Y  X							     V
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:                                                                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:                                                         M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv2DCurvatureSign, curvature                                        M
*****************************************************************************/
CagdCrvStruct *SymbCrv2DCurvatureSign(CagdCrvStruct *Crv)
{
    CagdBType
        IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z,
		  *Crv2W, *Crv2X, *Crv2Y, *Crv2Z,
		  *CTmp1, *CTmp2, *Numer, *Denom, *CurvatureSign,
	*Crv1Deriv = CagdCrvDerive(Crv),
	*Crv2Deriv = CagdCrvDerive(Crv1Deriv);

    SymbCrvSplitScalar(Crv1Deriv, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z);
    SymbCrvSplitScalar(Crv2Deriv, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z);
    CagdCrvFree(Crv1Deriv);
    CagdCrvFree(Crv2Deriv);

    CTmp1 = SymbCrvMult(Crv1X, Crv2Y);
    CTmp2 = SymbCrvMult(Crv2X, Crv1Y);
    Numer = SymbCrvSub(CTmp1, CTmp2);
    CagdCrvFree(CTmp1);
    CagdCrvFree(CTmp2);

    if (IsRational) {
	Denom = SymbCrvMult(Crv1W, Crv2W);

	CagdMakeCrvsCompatible(&Denom, &Numer, TRUE, TRUE);
	CurvatureSign = SymbCrvMergeScalar(Denom, Numer, NULL, NULL);
	CagdCrvFree(Denom);
	CagdCrvFree(Numer);
    }
    else {
	CurvatureSign = Numer;
    }

    CagdCrvFree(Crv1X);
    CagdCrvFree(Crv1Y);
    CagdCrvFree(Crv2X);
    CagdCrvFree(Crv2Y);
    if (Crv1Z)
	CagdCrvFree(Crv1Z);
    if (Crv2Z)
	CagdCrvFree(Crv2Z);
    if (Crv1W)
	CagdCrvFree(Crv1W);
    if (Crv2W)
	CagdCrvFree(Crv2W);

    return CurvatureSign;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Given a scalar curve that is positive, refine it until all its control     M
* points has positive coefficients. Always returns a Bspline curve.          M
*                                                                            *
* PARAMETERS:                                                                M
*   OrigCrv:    To refine until all its control points are non negative.     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:    Refined positive curve with positive control points. M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbMakePosCrvCtlPolyPos, refinement                                     M
*****************************************************************************/
CagdCrvStruct *SymbMakePosCrvCtlPolyPos(CagdCrvStruct *OrigCrv)
{
    int l;
    CagdCrvStruct
	*RefCrv = NULL;

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

    for (l = 0; l < MAX_POS_REF_ITERATION; l++) {
	int i, j,
	    Len = RefCrv -> Length,
	    Order = RefCrv -> Order;
	CagdRType
	    *KV = RefCrv -> KnotVector,
	    *Nodes = BspKnotNodes(KV, Len + Order, Order),
	    *Pts = RefCrv -> Points[1];

	for (i = j = 0; i < Len; i++) {
	    if (ABS(Pts[i]) < SQR(EPSILON))
		Pts[i] = 0.0;
	    if (Pts[i] < 0)	    /* To refine at negative control points. */
	        Nodes[j++] = Nodes[i];
	}
	if (j == 0) {
	    IritFree(Nodes);
	    break;
	}
	else {
	    CagdCrvStruct
	        *NewCrv = CagdCrvRefineAtParams(RefCrv, FALSE, Nodes, j);

	    CagdCrvFree(RefCrv);
	    RefCrv = NewCrv;
	    IritFree(Nodes);
	}
    }

    return RefCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a planar curve, finds all its inflection points by finding the zero  M
* set of the sign of the curvature function of the curve.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To find all its inflection points.                           M
*   Epsilon:    Accuracy control.                                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPtStruct *:    A list of parameter values on Crv that are inflection M
*                      points.                                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv2DInflectionPts, curvature, inflection points                     M
*****************************************************************************/
CagdPtStruct *SymbCrv2DInflectionPts(CagdCrvStruct *Crv, CagdRType Epsilon)
{
    CagdCrvStruct
	*CrvtrSign = SymbCrv2DCurvatureSign(Crv);
    CagdPtStruct
	*InflectionPts = SymbCrvZeroSet(CrvtrSign, 1, Epsilon);

    CagdCrvFree(CrvtrSign);

    return InflectionPts;    
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a planar curve, finds all its extreme curvatrue points by finding    M
* the set of extreme locations on  the curvature function of Crv.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       To find all int extrem curvature locations.                   M
*   Epsilon:   Accuracy control.                                             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPtStruct *:   A list of parameter values on Crv that have extrem     M
*                     curvature values.                                      M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrv2DExtremCrvtrPts, curvature                                       M
*****************************************************************************/
CagdPtStruct *SymbCrv2DExtremCrvtrPts(CagdCrvStruct *Crv, CagdRType Epsilon)
{
    if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 2) {
	CagdCrvStruct
	    *CrvtrSqrCrv = SymbCrv2DCurvatureSqr(Crv);
	CagdPtStruct
	    *ExtremCrvtrPts = SymbCrvExtremSet(CrvtrSqrCrv, 1, Epsilon);

	CagdCrvFree(CrvtrSqrCrv);

	return ExtremCrvtrPts;    
    }
    else if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 3) {
	CagdCrvStruct
	    *CrvtrNormalCrv = SymbCrv3DCurvatureNormal(Crv),
	    *CrvtrMag = SymbCrvDotProd(CrvtrNormalCrv, CrvtrNormalCrv);
	CagdPtStruct
	    *ExtremCrvtrPts = SymbCrvExtremSet(CrvtrMag, 1, Epsilon);

	CagdCrvFree(CrvtrNormalCrv);
	CagdCrvFree(CrvtrMag);

	return ExtremCrvtrPts;    
    }
    else {
	SYMB_FATAL_ERROR(SYMB_ERR_ONLY_2D_OR_3D);
	return NULL;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes coefficients of the first fundamental form of given surface Srf.  M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:         Do compute the coefficients of the FFF for.                 M
*   DuSrf:       First derivative of Srf with respect to U goes to here.     M
*   DvSrf:       First derivative of Srf with respect to V goes to here.     M
*   FffG11:      FFF G11 scalar field.                                       M
*   FffG12:      FFF G12 scalar field.                                       M
*   FffG22:      FFF G22 scalar field.                                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfFff, first fundamental form                                       M
*****************************************************************************/
void SymbSrfFff(CagdSrfStruct *Srf,
		CagdSrfStruct **DuSrf,
		CagdSrfStruct **DvSrf,
		CagdSrfStruct **FffG11,
		CagdSrfStruct **FffG12,
		CagdSrfStruct **FffG22)
{
    *DuSrf = CagdSrfDerive(Srf, CAGD_CONST_U_DIR);
    *DvSrf = CagdSrfDerive(Srf, CAGD_CONST_V_DIR);

    *FffG11 = SymbSrfDotProd(*DuSrf, *DuSrf);
    *FffG12 = SymbSrfDotProd(*DuSrf, *DvSrf);
    *FffG22 = SymbSrfDotProd(*DvSrf, *DvSrf);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes coefficients of the first fundamental form of given surface Srf.  M
*   These coefficients are using non normalized normal that is also returned.M
*                                                                            *
* PARAMETERS:                                                                M
*   DuSrf:       First derivative of Srf with respect to U.		     M
*   DvSrf:       First derivative of Srf with respect to V.		     M
*   SffL11:      SFF L11 scalar field returned herein.                       M
*   SffL12:      SFF L12 scalar field returned herein.                       M
*   SffL22:      SFF L22 scalar field returned herein.                       M
*   SNormal:     Unnormalized normal vector field returned herein.           M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfFff, second fundamental form                                      M
*****************************************************************************/
void SymbSrfSff(CagdSrfStruct *DuSrf,
		CagdSrfStruct *DvSrf,
		CagdSrfStruct **SffL11,
		CagdSrfStruct **SffL12,
		CagdSrfStruct **SffL22,
		CagdSrfStruct **SNormal)
{
    CagdSrfStruct
	*DuuSrf = CagdSrfDerive(DuSrf, CAGD_CONST_U_DIR),
	*DuvSrf = CagdSrfDerive(DuSrf, CAGD_CONST_V_DIR),
	*DvvSrf = CagdSrfDerive(DvSrf, CAGD_CONST_V_DIR);

    *SNormal = SymbSrfCrossProd(DvSrf, DuSrf);
    *SffL11 = SymbSrfDotProd(DuuSrf, *SNormal);
    *SffL12 = SymbSrfDotProd(DuvSrf, *SNormal);
    *SffL22 = SymbSrfDotProd(DvvSrf, *SNormal);

    CagdSrfFree(DuuSrf);
    CagdSrfFree(DuvSrf);
    CagdSrfFree(DvvSrf);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes the expression of Srf11 * Srf22 - Srf12 * Srf21, which is a	     M
* determinant of a 2 by 2 matrix.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf11, Srf12, Srf21, Srf22:  The four factors of the determinant.        M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *: A scalar field representing the determinant computation.M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrf2DDeterminant, determinant                                        M
*****************************************************************************/
CagdSrfStruct *SymbSrf2DDeterminant(CagdSrfStruct *Srf11,
				    CagdSrfStruct *Srf12,
				    CagdSrfStruct *Srf21,
				    CagdSrfStruct *Srf22)
{
    CagdSrfStruct
	*Prod1 = SymbSrfMult(Srf11, Srf22),
	*Prod2 = SymbSrfMult(Srf21, Srf12),
	*Add12 = SymbSrfSub(Prod1, Prod2);

    CagdSrfFree(Prod1);
    CagdSrfFree(Prod2);
    return Add12;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes curvature upper bound as Xi = k1^2 + k2^2, where k1 and k2 are    M
* the principal curvatures. 						     M
* Gij are the coefficients of the first fundamental form and Lij are of the  M
* second, using non unit normal n,					     M
*									     M
*      ( G11 L22 + G22 L11 - 2 G12 L12 )^2 - 2 |G| |L|			     V
* Xi = -----------------------------------------------			     V
*			     |G|^2 ||n||^2				     V
*									     M
* See: "Second Order Surface Analysis Using Hybrid of Symbolic and Numeric   M
* Operators", By Gershon Elber and Elaine Cohen, Transaction on graphics,    M
* Vol. 12, No. 2, pp 160-178, April 1993.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       Surface to compute curvature bound for.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A scalar field representing the curvature bound.      M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfCurvatureUpperBound, curvature                                    M
*****************************************************************************/
CagdSrfStruct *SymbSrfCurvatureUpperBound(CagdSrfStruct *Srf)
{
    CagdSrfStruct *DuSrf, *DvSrf, *SNormal, *STmp1, *STmp2, *STmp3, *STmp4,
	*Numer, *FffDeterminant, *SffDeterminant, *CurvatureBound,
	*Denom, *FffG11, *FffG12, *FffG22, *SffL11, *SffL12, *SffL22;

    SymbSrfFff(Srf, &DuSrf, &DvSrf, &FffG11, &FffG12, &FffG22);
    SymbSrfSff(DuSrf, DvSrf, &SffL11, &SffL12, &SffL22, &SNormal);
    CagdSrfFree(DuSrf);
    CagdSrfFree(DvSrf);

    STmp1 = SymbSrfMult(FffG11, SffL22);
    STmp2 = SymbSrfMult(FffG22, SffL11);
    STmp3 = SymbSrfMult(FffG12, SffL12);
    STmp4 = SymbSrfScalarScale(STmp3, 2.0);
    CagdSrfFree(STmp3);
    STmp3 = SymbSrfAdd(STmp1, STmp2);
    CagdSrfFree(STmp1);
    CagdSrfFree(STmp2);
    STmp1 = SymbSrfSub(STmp3, STmp4);
    CagdSrfFree(STmp3);
    CagdSrfFree(STmp4);
    STmp2 = SymbSrfMult(STmp1, STmp1);
    CagdSrfFree(STmp1);

    FffDeterminant = SymbSrf2DDeterminant(FffG11, FffG12, FffG12, FffG22);
    SffDeterminant = SymbSrf2DDeterminant(SffL11, SffL12, SffL12, SffL22);
    CagdSrfFree(FffG11);
    CagdSrfFree(FffG12);
    CagdSrfFree(FffG22);
    CagdSrfFree(SffL11);
    CagdSrfFree(SffL12);
    CagdSrfFree(SffL22);

    STmp1 = SymbSrfMult(FffDeterminant, SffDeterminant);
    STmp3 = SymbSrfScalarScale(STmp1, 2.0);
    CagdSrfFree(STmp1);
    Numer = SymbSrfSub(STmp2, STmp3);
    CagdSrfFree(STmp2);
    CagdSrfFree(STmp3);

    STmp1 = SymbSrfDotProd(SNormal, SNormal);
    CagdSrfFree(SNormal);

    STmp2 = SymbSrfMult(FffDeterminant, FffDeterminant);
    Denom = SymbSrfMult(STmp1, STmp2);
    CagdSrfFree(STmp1);
    CagdSrfFree(STmp2);

    CagdMakeSrfsCompatible(&Denom, &Numer, TRUE, TRUE, TRUE, TRUE);
    CurvatureBound = SymbSrfMergeScalar(Denom, Numer, NULL, NULL);
    CagdSrfFree(Denom);
    CagdSrfFree(Numer);

    return CurvatureBound;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes normal curvature bound in given isoparametric direction.          M
*   This turns out to be (L11 . n) / G11 for u and (L22 . n) / G22 for v.    M
*   Herein the square of these equations is computed symbolically and        M
* returned.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      To compute normal curvature in an isoparametric direction Dir. M
*   Dir:      Direction to compute normal curvature. Either U or V.          M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A scalar field representing the normal curvature       M
*                     square of Srf in dirction Dir.                         M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbSrfIsoDirNormalCurvatureBound, curvature                             M
*****************************************************************************/
CagdSrfStruct *SymbSrfIsoDirNormalCurvatureBound(CagdSrfStruct *Srf,
						 CagdSrfDirType Dir)
{
    CagdSrfStruct *DuSrf, *DvSrf, *SNormal, *STmp1, *STmp2, *STmp3, *STmp4,
	*SNormalSize, *FffG11, *FffG12, *FffG22, *SffL11, *SffL12, *SffL22,
	*CurvatureBound;

    SymbSrfFff(Srf, &DuSrf, &DvSrf, &FffG11, &FffG12, &FffG22);
    SymbSrfSff(DuSrf, DvSrf, &SffL11, &SffL12, &SffL22, &SNormal);
    CagdSrfFree(DuSrf);
    CagdSrfFree(DvSrf);

    SNormalSize = SymbSrfDotProd(SNormal, SNormal);

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    STmp2 = SymbSrfMult(SffL11, SffL11);
	    STmp3 = SymbSrfMult(FffG11, FffG11);
	    STmp4 = SymbSrfMult(SNormalSize, STmp3);
	    CagdSrfFree(STmp3);
	    STmp1 = SymbSrfInvert(STmp4);
	    CagdSrfFree(STmp4);
	    break;
	case CAGD_CONST_V_DIR:
	    STmp2 = SymbSrfMult(SffL22, SffL22);
	    STmp3 = SymbSrfMult(FffG22, FffG22);
	    STmp4 = SymbSrfMult(SNormalSize, STmp3);
	    CagdSrfFree(STmp3);
	    STmp1 = SymbSrfInvert(STmp4);
	    CagdSrfFree(STmp4);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
	    STmp1 = STmp2 = NULL;
    }

    CurvatureBound = SymbSrfMult(STmp1, STmp2);
    CagdSrfFree(STmp1);
    CagdSrfFree(STmp2);
    CagdSrfFree(SNormalSize);
    CagdSrfFree(FffG11);
    CagdSrfFree(FffG12);
    CagdSrfFree(FffG22);
    CagdSrfFree(SffL11);
    CagdSrfFree(SffL12);
    CagdSrfFree(SffL22);

    return CurvatureBound;
}
