/******************************************************************************
* Bzr_Sym.c - Bezier symbolic computation.				      *
*******************************************************************************
* Written by Gershon Elber, Sep. 92.					      *
******************************************************************************/

#include "symb_loc.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two Bezier curves - multiply them coordinatewise.		     M
*   The two curves are promoted to same point type before the multiplication M
* can take place.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:   The two curves to multiply.                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  The product Crv1 * Crv2 coordinatewise.                M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvMult, product                                                      M
*****************************************************************************/
CagdCrvStruct *BzrCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdBType IsNotRational;
    int i, j, k, MaxCoord, Order,
	Order1 = Crv1 -> Order,
	Order2 = Crv2 -> Order,
	Degree1 = Order1 - 1,
	Degree2 = Order2 - 1;
    CagdCrvStruct *ProdCrv;
    CagdRType **Points, **Points1, **Points2;

    if (!CAGD_IS_BEZIER_CRV(Crv1) || !CAGD_IS_BEZIER_CRV(Crv2))
	SYMB_FATAL_ERROR(SYMB_ERR_BZR_CRV_EXPECT);

    Crv1 = CagdCrvCopy(Crv1);
    Crv2 = CagdCrvCopy(Crv2);
    if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE))
	SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);

    ProdCrv = BzrCrvNew(Order = Order1 + Order2 - 1, Crv1 -> PType);
    IsNotRational = !CAGD_IS_RATIONAL_CRV(ProdCrv);
    MaxCoord = CAGD_NUM_OF_PT_COORD(ProdCrv -> PType);

    Points = ProdCrv -> Points;
    Points1 = Crv1 -> Points;
    Points2 = Crv2 -> Points;

    for (i = 0; i < Order; i++)
	for (k = IsNotRational; k <= MaxCoord; k++)
	    Points[k][i] = 0.0;

    for (i = 0; i < Order1; i++) {
	for (j = 0; j < Order2; j++) {
	    CagdRType
		Coef = CagdIChooseK(i, Degree1) *
		       CagdIChooseK(j, Degree2) /
		       CagdIChooseK(i + j, Degree1 + Degree2);

	    for (k = IsNotRational; k <= MaxCoord; k++)
		Points[k][i+j] += Coef * Points1[k][i] * Points2[k][j];
	}
    }

    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    return ProdCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two Bezier curve lists - multiply them one at a time.		     M
*   Return a Bezier curve lists representing their products.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1Lst:    First list of Bezier curves to multiply.                     M
*   Crv2Lst:    Second list of Bezier curves to multiply.                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A list of product curves                              M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvMultList, product                                                  M
*****************************************************************************/
CagdCrvStruct *BzrCrvMultList(CagdCrvStruct *Crv1Lst, CagdCrvStruct *Crv2Lst)
{
    CagdCrvStruct *Crv1, *Crv2,
	*ProdCrvTail = NULL,
	*ProdCrvList = NULL;

    for (Crv1 = Crv1Lst, Crv2 = Crv2Lst;
	 Crv1 != NULL && Crv2 != NULL;
	 Crv1 = Crv1 -> Pnext, Crv2 = Crv2 -> Pnext) {
	CagdCrvStruct
	    *ProdCrv = BzrCrvMult(Crv1, Crv2);

	if (ProdCrvList == NULL)
	    ProdCrvList = ProdCrvTail = ProdCrv;
	else {
	    ProdCrvTail -> Pnext = ProdCrv;
	    ProdCrvTail = ProdCrv;
	}
    }

    return ProdCrvList;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two Bezier surfaces - multiply them coordinatewise.		     M
*   The two surfaces are promoted to same point type before multiplication   M
* can take place.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:   The two surfaces to multiply.		                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  The product Srf1 * Srf2 coordinatewise.                M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfMult, product                                                      M
*****************************************************************************/
CagdSrfStruct *BzrSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdBType IsNotRational;
    int i, j, k, l, m, MaxCoord, UOrder, VOrder,
	UOrder1 = Srf1 -> UOrder,
	VOrder1 = Srf1 -> VOrder,
	UOrder2 = Srf2 -> UOrder,
	VOrder2 = Srf2 -> VOrder,
	UDegree1 = UOrder1 - 1,
	VDegree1 = VOrder1 - 1,
	UDegree2 = UOrder2 - 1,
	VDegree2 = VOrder2 - 1;
    CagdSrfStruct *ProdSrf;
    CagdRType **Points, **Points1, **Points2;

    if (!CAGD_IS_BEZIER_SRF(Srf1) || !CAGD_IS_BEZIER_SRF(Srf2))
	SYMB_FATAL_ERROR(SYMB_ERR_BZR_SRF_EXPECT);

    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);
    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    ProdSrf = BzrSrfNew(UOrder = UOrder1 + UOrder2 - 1,
			VOrder = VOrder1 + VOrder2 - 1, Srf1 -> PType);
    IsNotRational = !CAGD_IS_RATIONAL_SRF(ProdSrf);
    MaxCoord = CAGD_NUM_OF_PT_COORD(ProdSrf -> PType);

    Points = ProdSrf -> Points;
    Points1 = Srf1 -> Points;
    Points2 = Srf2 -> Points;

    for (k = IsNotRational; k <= MaxCoord; k++) {
	CagdRType
	    *Pts = Points[k];

	for (i = 0; i < UOrder * VOrder; i++)
	    *Pts++ = 0.0;
    }

    for (i = 0; i < UOrder1; i++) {
	for (j = 0; j < VOrder1; j++) {
	    for (l = 0; l < UOrder2; l++) {
		for (m = 0; m < VOrder2; m++) {
		    int Index = CAGD_MESH_UV(ProdSrf, i + l, j + m),
			Index1 = CAGD_MESH_UV(Srf1, i, j),
			Index2 = CAGD_MESH_UV(Srf2, l, m);
		    CagdRType
			Coef1 = CagdIChooseK(i, UDegree1) *
			        CagdIChooseK(l, UDegree2) /
				CagdIChooseK(i + l, UDegree1 + UDegree2),
			Coef2 = CagdIChooseK(j, VDegree1) *
			        CagdIChooseK(m, VDegree2) /
				CagdIChooseK(j + m, VDegree1 + VDegree2);

		    for (k = IsNotRational; k <= MaxCoord; k++)
			Points[k][Index] += Coef1 * Coef2 * 
			    Points1[k][Index1] * Points2[k][Index2];
		}
	    }
	}
    }

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return ProdSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a rational Bezier curve - computes its derivative curve (Hodograph) M
* using the quotient rule for differentiation.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:         Bezier curve to differentiate.                             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:    Differentiated rational Bezier curve.               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvDeriveRational, derivatives                                        M
*****************************************************************************/
CagdCrvStruct *BzrCrvDeriveRational(CagdCrvStruct *Crv)
{
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *DCrvW, *DCrvX, *DCrvY, *DCrvZ,
	*TCrv1, *TCrv2, *DeriveCrv;

    SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    if (CrvW)
	DCrvW = BzrCrvDerive(CrvW);
    else {
    	DCrvW = NULL;
	SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED);
    }
    
    if (CrvX) {
	DCrvX = BzrCrvDerive(CrvX);

	TCrv1 = BzrCrvMult(DCrvX, CrvW);
	TCrv2 = BzrCrvMult(CrvX, DCrvW);

	CagdCrvFree(CrvX);
	CagdCrvFree(DCrvX);
	CrvX = SymbCrvSub(TCrv1, TCrv2);
	CagdCrvFree(TCrv1);
	CagdCrvFree(TCrv2);
    }

    if (CrvY) {
	DCrvY = BzrCrvDerive(CrvY);

	TCrv1 = BzrCrvMult(DCrvY, CrvW);
	TCrv2 = BzrCrvMult(CrvY, DCrvW);

	CagdCrvFree(CrvY);
	CagdCrvFree(DCrvY);
	CrvY = SymbCrvSub(TCrv1, TCrv2);
	CagdCrvFree(TCrv1);
	CagdCrvFree(TCrv2);
    }

    if (CrvZ) {
	DCrvZ = BzrCrvDerive(CrvZ);

	TCrv1 = BzrCrvMult(DCrvZ, CrvW);
	TCrv2 = BzrCrvMult(CrvZ, DCrvW);

	CagdCrvFree(CrvZ);
	CagdCrvFree(DCrvZ);
	CrvZ = SymbCrvSub(TCrv1, TCrv2);
	CagdCrvFree(TCrv1);
	CagdCrvFree(TCrv2);
    }

    TCrv1 = BzrCrvMult(CrvW, CrvW);
    CagdCrvFree(CrvW);
    CrvW = TCrv1;

    if (!CagdMakeCrvsCompatible(&CrvW, &CrvX, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&CrvW, &CrvY, TRUE, TRUE) ||
	!CagdMakeCrvsCompatible(&CrvW, &CrvZ, TRUE, TRUE))
	SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);

    DeriveCrv = SymbCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ);

    if (CrvX)
	CagdCrvFree(CrvX);
    if (CrvY)
	CagdCrvFree(CrvY);
    if (CrvZ)
	CagdCrvFree(CrvZ);
    if (CrvW) {
	CagdCrvFree(CrvW);
	CagdCrvFree(DCrvW);
    }

    return DeriveCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a rational Bezier surface - computes its derivative surface in      M
* direction Dir, using the quotient rule for differentiation.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:         Bezier surface to differentiate.                           M
*   Dir:         Direction of Differentiation. Either U or V.                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:    Differentiated rational Bezier surface.             M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfDeriveRational, derivatives                                        M
*****************************************************************************/
CagdSrfStruct *BzrSrfDeriveRational(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ,
	*TSrf1, *TSrf2, *DeriveSrf;

    SymbSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ);
    if (SrfW)
	DSrfW = BzrSrfDerive(SrfW, Dir);
    else {
    	DSrfW = NULL;
	SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED);
    }

    if (SrfX) {
	DSrfX = BzrSrfDerive(SrfX, Dir);

	TSrf1 = BzrSrfMult(DSrfX, SrfW);
	TSrf2 = BzrSrfMult(SrfX, DSrfW);

	CagdSrfFree(SrfX);
	CagdSrfFree(DSrfX);
	SrfX = SymbSrfSub(TSrf1, TSrf2);
	CagdSrfFree(TSrf1);
	CagdSrfFree(TSrf2);
    }
    if (SrfY) {
	DSrfY = BzrSrfDerive(SrfY, Dir);

	TSrf1 = BzrSrfMult(DSrfY, SrfW);
	TSrf2 = BzrSrfMult(SrfY, DSrfW);

	CagdSrfFree(SrfY);
	CagdSrfFree(DSrfY);
	SrfY = SymbSrfSub(TSrf1, TSrf2);
	CagdSrfFree(TSrf1);
	CagdSrfFree(TSrf2);
    }
    if (SrfZ) {
	DSrfZ = BzrSrfDerive(SrfZ, Dir);

	TSrf1 = BzrSrfMult(DSrfZ, SrfW);
	TSrf2 = BzrSrfMult(SrfZ, DSrfW);

	CagdSrfFree(SrfZ);
	CagdSrfFree(DSrfZ);
	SrfZ = SymbSrfSub(TSrf1, TSrf2);
	CagdSrfFree(TSrf1);
	CagdSrfFree(TSrf2);
    }
    
    TSrf1 = BzrSrfMult(SrfW, SrfW);
    CagdSrfFree(SrfW);
    SrfW = TSrf1;

    if (!CagdMakeSrfsCompatible(&SrfW, &SrfX, TRUE, TRUE, TRUE, TRUE) ||
	!CagdMakeSrfsCompatible(&SrfW, &SrfY, TRUE, TRUE, TRUE, TRUE) ||
	!CagdMakeSrfsCompatible(&SrfW, &SrfZ, TRUE, TRUE, TRUE, TRUE))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    DeriveSrf = SymbSrfMergeScalar(SrfW, SrfX, SrfY, SrfZ);

    if (SrfX)
	CagdSrfFree(SrfX);
    if (SrfY)
	CagdSrfFree(SrfY);
    if (SrfZ)
	CagdSrfFree(SrfZ);
    if (SrfW) {
	CagdSrfFree(SrfW);
	CagdSrfFree(DSrfW);
    }

    return DeriveSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bezier curve - convert it to (possibly) piecewise cubic.	     M
* If the curve is							     M
* 1. A cubic - a copy if it is returned.				     M
* 2. Lower than cubic - a degree raised (to a cubic) curve is returned.      M
* 3. Higher than cubic - a C^1 continuous piecewise cubic approximation      M
*    is computed for Crv.						     M
*									     M
* In case 3 a list of polynomial cubic curves is returned. Tol is then used  M
* for the distance tolerance error measure for the approximation.	     M
*   If, however, NoRational is set, rational curves of any order will also   M
* be approximated using cubic polynomials.				     M
*   Furthermore if the total length of control polygon is more than MaxLen,  M
* the curve is subdivided until this is not the case.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:          To approximate using cubic Bezier polynomials.             M
*   Tol:          Accuracy control.                                          M
*   MaxLen:       Maximum arc length of curve.                               M
*   NoRational:   Do we want to approximate rational curves as well?         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A list of cubic Bezier polynomials approximating Crv.  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrApproxBzrCrvAsCubics, conversion, approximation                       M
*****************************************************************************/
CagdCrvStruct *BzrApproxBzrCrvAsCubics(CagdCrvStruct *Crv,
				       CagdRType Tol,
				       CagdRType MaxLen,
				       CagdBType NoRational)
{
    CagdBType
	NewCrv = FALSE;
    CagdCrvStruct *TCrv, *CubicCrvs,
	*CubicCrvsMaxLen = NULL;

    if (CAGD_IS_RATIONAL_CRV(Crv) && NoRational)
	return BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
        NewCrv = TRUE;
        Crv = CnvrtPeriodic2FloatCrv(Crv);
    }
    if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) {
	TCrv = BspCrvOpenEnd(Crv);

	if (NewCrv)
	    CagdCrvFree(Crv);
	Crv = TCrv;
	NewCrv = TRUE;
    }

    switch (Crv -> Order) {
	case 2:
	    CubicCrvs = BzrCrvDegreeRaiseN(Crv, 4);
	    break;
	case 3:
	    CubicCrvs = BzrCrvDegreeRaise(Crv);
	    break;
	case 4:
	    CubicCrvs = CagdCrvCopy(Crv);
	default:
	    if (Crv -> Order < 4)
		SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE);
	    CubicCrvs = BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);
    }

    for (TCrv = CubicCrvs; TCrv != NULL; TCrv = TCrv -> Pnext) {
	CagdCrvStruct *TCrv2,
	    *MaxLenCrv = SymbLimitCrvArcLen(TCrv, MaxLen);

	for (TCrv2 = MaxLenCrv;
	     TCrv2 -> Pnext != NULL;
	     TCrv2 = TCrv2 -> Pnext);

	TCrv2 -> Pnext = CubicCrvsMaxLen;
	CubicCrvsMaxLen = MaxLenCrv;
    }

    if (NewCrv)
	CagdCrvFree(Crv);

    CagdCrvFreeList(CubicCrvs);
    return CubicCrvsMaxLen;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bezier curve with order larger than cubic, approximate it using    M
* piecewise cubic curves.						     M
*   A C^1 continuous approximation is computed so that the approximation is  M
* at most sqrt(Tol2) away from the given curve.				     M
*   Input curve can be rational, although output is always polynomial.       M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       To approximate using cubic Bezier curves.                     M
*   Tol2:      Accuracy control.                                             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   List of cubic Bezier curves approximating Crv.        M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrApproxBzrCrvAsCubicPoly, approximation, conversion                    M
*****************************************************************************/
CagdCrvStruct *BzrApproxBzrCrvAsCubicPoly(CagdCrvStruct *Crv, CagdRType Tol2)
{
    CagdBType
	NewCrv = FALSE,
	ApproxOK = TRUE,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i,
	Order = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdPointType
	PType = Crv -> PType,
	CubicPType = CAGD_MAKE_PT_TYPE(FALSE, MaxCoord);
    CagdCrvStruct *DistSqrCrv, *DiffCrv,
	*CubicBzr = BzrCrvNew(4, CubicPType);
    CagdRType E3Pt1[3], E3Pt2[3], E3Pt3[3], E3Pt4[3], Tan1[3], Tan2[3],
	**CubicPoints = CubicBzr -> Points,
	**Points = Crv -> Points;

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
        NewCrv = TRUE;
        Crv = CnvrtPeriodic2FloatCrv(Crv);
    }
    if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) {
	CagdCrvStruct
	    *TCrv = BspCrvOpenEnd(Crv);

	if (NewCrv)
	    CagdCrvFree(Crv);
	Crv = TCrv;
	NewCrv = TRUE;
    }

    CagdCoerceToE3(E3Pt1, Points, 0, PType);
    CagdCoerceToE3(E3Pt2, Points, 1, PType);
    CagdCoerceToE3(E3Pt3, Points, Order - 2, PType);
    CagdCoerceToE3(E3Pt4, Points, Order - 1, PType);
	
    /* End of the two points must match. */
    for (i = 0; i < MaxCoord; i++) {
	CubicPoints[i + 1][0] = E3Pt1[i];
	CubicPoints[i + 1][3] = E3Pt4[i];
    }

    /* Tangents at the end of the two points must match. */
    PT_SUB(Tan1, E3Pt2, E3Pt1);
    PT_SUB(Tan2, E3Pt4, E3Pt3);
    PT_SCALE(Tan1, (Order - 1.0) / 3.0);
    PT_SCALE(Tan2, (Order - 1.0) / 3.0);

    for (i = 0; i < MaxCoord; i++) {
	CubicPoints[i + 1][1] = E3Pt1[i] + Tan1[i];
	CubicPoints[i + 1][2] = E3Pt4[i] - Tan2[i];
    }

    /* Compare the two curves by computed the distance square between */
    /* corresponding parameter values.				      */
    DiffCrv = SymbCrvSub(Crv, CubicBzr);
    DistSqrCrv = SymbCrvDotProd(DiffCrv, DiffCrv);
    CagdCrvFree(DiffCrv);
    Points = DistSqrCrv -> Points;
    if (IsNotRational) {
	CagdRType
	    *Dist = Points[1];

	for (i = DistSqrCrv -> Order - 1; i >= 0; i--) {
	    if (*Dist++ > Tol2) {
		ApproxOK = FALSE;
		break;
	    }
	}		
    }
    else {
	CagdRType
	    *Dist0 = Points[0],
	    *Dist1 = Points[1];

	for (i = DistSqrCrv -> Order - 1; i >= 0; i--) {
	    if (*Dist1++ / *Dist0++ > Tol2) {
		ApproxOK = FALSE;
		break;
	    }
	}		
    }
    CagdCrvFree(DistSqrCrv);

    if (ApproxOK) {
	if (NewCrv)
	    CagdCrvFree(Crv);
	return CubicBzr;
    }
    else {
	CagdCrvStruct
	    *Crv1 = BzrCrvSubdivAtParam(Crv, 0.5),
	    *Crv2 = Crv1 -> Pnext,
	    *Crv1Approx = BzrApproxBzrCrvAsCubicPoly(Crv1, Tol2),
	    *Crv2Approx = BzrApproxBzrCrvAsCubicPoly(Crv2, Tol2);

	CagdCrvFree(Crv1);
	CagdCrvFree(Crv2);

	for (Crv1 = Crv1Approx; Crv1 -> Pnext != NULL; Crv1 = Crv1 ->Pnext);
	Crv1 -> Pnext = Crv2Approx;
	if (NewCrv)
	    CagdCrvFree(Crv);
	return Crv1Approx;	
    }
}

