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

#include "cagd_loc.h"

/******************************************************************************
* Given two bezier curves - multiply them coordinatewise.		      *
* The two curves are promoted to same point type before the multiplication    *
* can take place.							      *
* Return a bezier curve representing their product.			      *
******************************************************************************/
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))
	FATAL_ERROR(CAGD_ERR_BZR_CRV_EXPECT);

    Crv1 = CagdCrvCopy(Crv1);
    Crv2 = CagdCrvCopy(Crv2);
    if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE))
	FATAL_ERROR(CAGD_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;
}

/******************************************************************************
* Given two bezier curve lists - multiply them one at a time.		      *
* Return a bezier curve lists representing their products.		      *
******************************************************************************/
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;
}

/******************************************************************************
* Given two bezier surfaces - multiply them coordinatewise.		      *
* The two surfaces are promoted to same point type before the multiplication  *
* can take place.							      *
* Return a bezier surface representing their product.			      *
******************************************************************************/
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))
	FATAL_ERROR(CAGD_ERR_BZR_SRF_EXPECT);

    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);
    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
	FATAL_ERROR(CAGD_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;
}

/******************************************************************************
* Given a rational bezier curve - computes its derivative curve using the     *
* quotient rule for differentiation.					      *
******************************************************************************/
CagdCrvStruct *BzrCrvDeriveRational(CagdCrvStruct *Crv)
{
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *DCrvW, *DCrvX, *DCrvY, *DCrvZ,
	*TCrv1, *TCrv2, *DeriveCrv;

    CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    if (CrvW)
	DCrvW = BzrCrvDerive(CrvW);
    else {
    	DCrvW = NULL;
	FATAL_ERROR(CAGD_ERR_RATIONAL_EXPECTED);
    }
    
    if (CrvX) {
	DCrvX = BzrCrvDerive(CrvX);

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

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

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

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

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

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

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

	CagdCrvFree(CrvZ);
	CagdCrvFree(DCrvZ);
	CrvZ = CagdCrvSub(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))
	FATAL_ERROR(CAGD_ERR_CRV_FAIL_CMPT);

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

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

    return DeriveCrv;
}

/******************************************************************************
* Given a rational bezier surface - computes its derivative curve using the   *
* quotient rule for differentiation.					      *
******************************************************************************/
CagdSrfStruct *BzrSrfDeriveRational(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ,
	*TSrf1, *TSrf2, *DeriveSrf;

    CagdSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ);
    if (SrfW)
	DSrfW = BzrSrfDerive(SrfW, Dir);
    else {
    	DSrfW = NULL;
	FATAL_ERROR(CAGD_ERR_RATIONAL_EXPECTED);
    }

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

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

	CagdSrfFree(SrfX);
	CagdSrfFree(DSrfX);
	SrfX = CagdSrfSub(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 = CagdSrfSub(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 = CagdSrfSub(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))
	FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);

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

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

    return DeriveSrf;
}

/******************************************************************************
* Given a bezier curve - convert it to (possibly) piecewise cubic.	      *
* If the curve is							      *
* 1. A cubic - a copy if it is returned.				      *
* 2. Lower than cubic - a degree raised (to a cubic) curve is returned.       *
* 3. Higher than cubic - a C^1 piecewise cubic approximation is computed.     *
*    Note in this case a list of polynomial cubic curves is returned. Tol is  *
*    then used for the distance tolerance error measure for the approximation *
*   If, however, NoRational is set, rational curves of any order will also be *
* approximated using cubic polynomials. Furthermore if the total length of    *
* control polygon is more than MaxLen, the curve is subdivided until this is  *
* not the case.								      *
******************************************************************************/
CagdCrvStruct *BzrApproxBzrCrvAsCubics(CagdCrvStruct *Crv, CagdRType Tol,
				       CagdRType MaxLen, CagdBType NoRational)
{
    CagdCrvStruct *TCrv, *CubicCrvs,
	*CubicCrvsMaxLen = NULL;

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

    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)
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    CubicCrvs = BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol);
    }

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

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

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

    CagdCrvFreeList(CubicCrvs);
    return CubicCrvsMaxLen;
}

/******************************************************************************
* Given a bezier curve with order larger than cubic, approximate it using     *
* piecewise cubic curves. A C^1 continuous approximation is computed so that  *
* the approximation is at most sqrt(Tol2) away from the given curve.	      *
*   Input curve can be rational, although output is always a polynomial.      *
******************************************************************************/
CagdCrvStruct *BzrApproxBzrCrvAsCubicPoly(CagdCrvStruct *Crv, CagdRType Tol2)
{
    CagdBType
	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;

    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 = CagdCrvSub(Crv, CubicBzr);
    DistSqrCrv = CagdCrvDotProd(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)
	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;
	return Crv1Approx;	
    }
}