/******************************************************************************
* CagdCmps.c - Composition computation (symbolic).			      *
*******************************************************************************
* Written by Gershon Elber, Apr. 93.					      *
******************************************************************************/

#include <string.h>
#include "cagd_loc.h"

#define ZERO_SET_EPS 0.0001

static CagdCrvStruct *CagdComposeCrvCrvAux(CagdCrvStruct *Crv1,
					   CagdCrvStruct *Crv2);
static CagdCrvStruct *CagdComposeCrvCrvAux2(CagdCrvStruct *Crv1,
					    CagdCrvStruct *Crv2);
static CagdCrvStruct *CagdComposeSrfCrvAux(CagdSrfStruct *Srf,
					   CagdCrvStruct *Crv);
static CagdCrvStruct **CagdComputeCurvePowers(CagdCrvStruct *Crv, int Order);

/******************************************************************************
* Given two curves, Crv1 and Crv2, computes the composition Crv1(Crv2(t)).    *
******************************************************************************/
CagdCrvStruct *CagdComposeCrvCrv(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *CmpsCrv;
    CagdRType TMax, TMin, CTMax, CTMin;

    switch (Crv1 -> GType) {
	case CAGD_CBEZIER_TYPE:
	case CAGD_CBSPLINE_TYPE:
	    break;
	case CAGD_CPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }

    switch (Crv2 -> GType) {
	case CAGD_CBEZIER_TYPE:
	case CAGD_CBSPLINE_TYPE:
	    break;
	case CAGD_CPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }

    CmpsCrv = CagdComposeCrvCrvAux(Crv1, Crv2);

    /* Now map the domain of the curve to be crv2 domain. */
    CagdCrvDomain(Crv2, &TMin, &TMax);
    CagdCrvDomain(CmpsCrv, &CTMin, &CTMax);
    if (CmpsCrv -> GType == CAGD_CBEZIER_TYPE) {
	CagdCrvStruct
	    *CTmp = CnvrtBezier2BsplineCrv(CmpsCrv);

	CagdCrvFree(CmpsCrv);
	CmpsCrv = CTmp;
    }
    BspKnotAffineTrans(CmpsCrv -> KnotVector,
		       CmpsCrv -> Length + CmpsCrv -> Order,
		       TMin - CTMin,
		       (TMax - TMin) / (CTMax - CTMin));

    return CmpsCrv;
}

/******************************************************************************
* Aux. function. Subdivides Crv2 until it is a Bezier curve, invokes the      *
* Bezier composition code on them, and merges them back to a complete curve.  *
* At this point, curves can be either Bezier or Bspline only.		      *
* Curves are both assumed to have open end condition.			      *
******************************************************************************/
static CagdCrvStruct *CagdComposeCrvCrvAux(CagdCrvStruct *Crv1,
					   CagdCrvStruct *Crv2)
{
    CagdCrvStruct *CmpsCrv;

    if (Crv2 -> Length > Crv2 -> Order) {
	CagdRType t;

	/* Crv2 is not a Bezier segment. Subdivide, compute for each segment */
	/* and merge back.						     */
	CagdCrvStruct *Crv2a, *Crv2b, *CmpsCrvA, *CmpsCrvB;

	t = Crv2 -> KnotVector[(Crv2 -> Order + Crv2 -> Length) / 2];

	Crv2a = CagdCrvSubdivAtParam(Crv2, t);
	Crv2b = Crv2a -> Pnext;

	CmpsCrvA = CagdComposeCrvCrvAux(Crv1, Crv2a);
	CmpsCrvB = CagdComposeCrvCrvAux(Crv1, Crv2b);
	CagdCrvFree(Crv2a);
	CagdCrvFree(Crv2b);

	CmpsCrv = CagdMergeCrvCrv(CmpsCrvA, CmpsCrvB);
	CagdCrvFree(CmpsCrvA);
	CagdCrvFree(CmpsCrvB);
    }
    else {
    	/* Crv2 is a Bezier segment - compute its composition. */
    	CmpsCrv = CagdComposeCrvCrvAux2(Crv1, Crv2);
    }

    return CmpsCrv;    
}

/******************************************************************************
* Aux. function. Subdivides Crv1 until it is a Bezier curve, invokes the      *
* Bezier composition code on them, and merges them back to a complete curve.  *
* At this point, Crv1 is a Bezier curve, Crv2 can be either Bezier or Bspline.*
* Curves are both assumed to have open end condition.			      *
******************************************************************************/
static CagdCrvStruct *CagdComposeCrvCrvAux2(CagdCrvStruct *Crv1,
					    CagdCrvStruct *Crv2)
{
    CagdRType t, TMin, TMax, CTMax, CTMin;
    CagdCrvStruct *CmpsCrv;

    if (Crv1 -> Length > Crv1 -> Order) {
	/* Needs to compose in pieces after subdividing Crv2 at the interior */
	/* knots of Crv1, by finding where Crv2(r) = t, the interior knot.   */
	CagdCrvStruct *Crv1a, *Crv1b, *Crv2a, *Crv2b, *CmpsCrvA, *CmpsCrvB;
	CagdPtStruct *Pts;

	t = Crv1 -> KnotVector[(Crv1 -> Order + Crv1 -> Length) / 2];
	Crv1a = CagdCrvSubdivAtParam(Crv1, t);
	Crv1b = Crv1a -> Pnext;

	Pts = CagdCrvConstSet(Crv2, 1, ZERO_SET_EPS, t);
	if (Pts) {
	    if (Pts -> Pnext != NULL)
		FATAL_ERROR(CAGD_ERR_REPARAM_NOT_MONOTONE);
	    t = Pts -> Pt[0];
	    CagdPtFreeList(Pts);
	    Crv2a = CagdCrvSubdivAtParam(Crv2, t);
	    Crv2b = Crv2a -> Pnext;
	}
	else
	    Crv2a = Crv2b = NULL;

	CmpsCrvA = CagdComposeCrvCrvAux2(Crv1a, Crv2a ? Crv2a : Crv2);
	CmpsCrvB = CagdComposeCrvCrvAux2(Crv1b, Crv2b ? Crv2b : Crv2);
	CagdCrvFree(Crv1a);
	CagdCrvFree(Crv1b);
	if (Crv2a)
	    CagdCrvFree(Crv2a);
	if (Crv2b)
	    CagdCrvFree(Crv2b);

	CmpsCrv = CagdMergeCrvCrv(CmpsCrvA, CmpsCrvB);
	CagdCrvFree(CmpsCrvA);
	CagdCrvFree(CmpsCrvB);
    }
    else {
	/* We can compose the curves, but first make sure the domain of */
	/* Crv1 is [0, 1] which is also the range of Crv2.		*/
	CagdCrvDomain(Crv1, &TMin, &TMax);
	if (!APX_EQ(TMin, 0.0) || !APX_EQ(TMax, 1.0)) {
	    CagdRType
		Translate = -TMin;

	    Crv2 = CagdCrvCopy(Crv2);
	    CagdCrvTransform(Crv2, &Translate, 1.0 / (TMax - TMin));
	    CmpsCrv = BzrComposeCrvCrv(Crv1, Crv2);
	    CagdCrvFree(Crv2);
	}
	else
	    CmpsCrv = BzrComposeCrvCrv(Crv1, Crv2);

	/* Now map the domain of the curve to be crv2 domain. */
	CagdCrvDomain(Crv2, &TMin, &TMax);
	CagdCrvDomain(CmpsCrv, &CTMin, &CTMax);
	if (CmpsCrv -> GType == CAGD_CBEZIER_TYPE) {
	    CagdCrvStruct
		*CTmp = CnvrtBezier2BsplineCrv(CmpsCrv);

	    CagdCrvFree(CmpsCrv);
	    CmpsCrv = CTmp;
	}
	BspKnotAffineTrans(CmpsCrv -> KnotVector,
			   CmpsCrv -> Length + CmpsCrv -> Order,
			   TMin - CTMin,
			   (TMax - TMin) / (CTMax - CTMin));
    }

    return CmpsCrv;
}

/******************************************************************************
* Given two Bezier curves, Crv1 and Crv2, computes composition Crv1(Crv2(t)). *
* See: "Freeform surfcae analysis using a hybrid of symbolic and numeric      *
* computation" by Gershon Elber, PhD thesis, University of Utah, 1992.	      *
******************************************************************************/
CagdCrvStruct *BzrComposeCrvCrv(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_CRV(Crv1);
    int i, j, k, CmpsOrder,
	Order = Crv1 -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv1 -> PType);
    CagdRType
	Translate = 0.0;
    CagdCrvStruct **Crv2Factors,
    	*CmpsCrv = NULL;

    if (CAGD_NUM_OF_PT_COORD(Crv2 -> PType) != 1)
	FATAL_ERROR(CAGD_ERR_WRONG_PT_TYPE);

    Crv2Factors = CagdComputeCurvePowers(Crv2, Order);
    CmpsCrv = BzrCrvNew(Crv2Factors[0] -> Length, Crv1 -> PType);
    CmpsOrder = CmpsCrv -> Order;

    for (j = !IsRational; j <= MaxCoord; j++) {
    	CagdRType
	    *CmpsPoints = CmpsCrv -> Points[j],
	    *Crv1Points = Crv1 -> Points[j];

	for (i = 0; i < Order; i++) {
	    CagdCrvStruct
		*CTmp = CagdCrvCopy(Crv2Factors[i]);
	    CagdRType
	        *CTmpPoints = CTmp -> Points[1];

	    CagdCrvTransform(CTmp, &Translate, *Crv1Points++);
	    if (i == 0) {
       		CAGD_GEN_COPY(CmpsPoints, CTmpPoints,
			      CmpsOrder * sizeof(CagdRType));
	    }
	    else {
	    	for (k = 0; k < CmpsOrder; k++)
	    	    CmpsPoints[k] += CTmpPoints[k];
	    }
	}
    }

    for (i = 0; i < Order; i++)
	CagdCrvFree(Crv2Factors[i]);

    if (CAGD_IS_RATIONAL_CRV(Crv2)) {
	CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *CTmp;

	CagdCrvSplitScalar(CmpsCrv, &CrvW, &CrvX, &CrvY, &CrvZ);
	CTmp = CagdCrvMergeScalar(Crv2Factors[Order], CrvX, CrvY, CrvZ);
	CagdCrvFree(CmpsCrv);
	CmpsCrv = CTmp;

	if (CrvX)
	    CagdCrvFree(CrvX);
	if (CrvY)
	    CagdCrvFree(CrvY);
	if (CrvZ)
	    CagdCrvFree(CrvZ);
	
	CagdCrvFree(Crv2Factors[Order]);
    }

    IritFree((VoidPtr) Crv2Factors);

    return CmpsCrv;
}

/******************************************************************************
* Given a curve Crv and surface Srf, computes the composition Srf(Crv(t)).    *
******************************************************************************/
CagdCrvStruct *CagdComposeSrfCrv(CagdSrfStruct *Srf, CagdCrvStruct *Crv)
{
    CagdCrvStruct *CmpsCrv;

    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	    break;
	case CAGD_SBSPLINE_TYPE:
	    FATAL_ERROR(CAGD_ERR_BSPLINE_NO_SUPPORT);
	    break;
	case CAGD_SPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }

    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	case CAGD_CBSPLINE_TYPE:
	    break;
	case CAGD_CPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }

    CmpsCrv = CagdComposeSrfCrvAux(Srf, Crv);

    return CmpsCrv;
}

/******************************************************************************
* Aux. function. Subdivides Crv until it is a Bezier curve, invokes the       *
* Bezier composition code on them, and merges them back to a complete curve.  *
* At this point, the curve can be either Bezier or Bspline only.              *
* Curve is assumed to have open end condition.				      *
******************************************************************************/
static CagdCrvStruct *CagdComposeSrfCrvAux(CagdSrfStruct *Srf,
					   CagdCrvStruct *Crv)
{
    CagdRType t, TMin, TMax;
    CagdCrvStruct *CmpsCrv;

    if (Crv -> Length > Crv -> Order) {
	/* Crv is not a Bezier segment. Subdivide, compute for each segment  */
	/* and merge back.						     */
	CagdCrvStruct *CrvA, *CrvB, *CmpsCrvA, *CmpsCrvB;

	t = Crv -> KnotVector[(Crv -> Order + Crv -> Length) / 2];

	CrvA = CagdCrvSubdivAtParam(Crv, t);
	CrvB = CrvA -> Pnext;

	CmpsCrvA = CagdComposeSrfCrvAux(Srf, CrvA);
	CmpsCrvB = CagdComposeSrfCrvAux(Srf, CrvB);
	CagdCrvFree(CrvA);
	CagdCrvFree(CrvB);

	CmpsCrv = CagdMergeCrvCrv(CmpsCrvA, CmpsCrvB);
	CagdCrvFree(CmpsCrvA);
	CagdCrvFree(CmpsCrvB);
    }
    else {
    	/* Crv is a Bezier segment - compute its composition. */
    	CmpsCrv = BzrComposeSrfCrv(Srf, Crv);

	/* Now map the domain of the composed curve to be crv's domain. */
	CagdCrvDomain(Crv, &TMin, &TMax);
	if (CmpsCrv -> GType == CAGD_CBEZIER_TYPE) {
	    CagdCrvStruct
		*CTmp = CnvrtBezier2BsplineCrv(CmpsCrv);

	    CagdCrvFree(CmpsCrv);
	    CmpsCrv = CTmp;
	}
	BspKnotAffineTrans(CmpsCrv -> KnotVector,
			   CmpsCrv -> Length + CmpsCrv -> Order,
			   TMin - CmpsCrv -> KnotVector[0],
			   TMax - TMin);
    }

    return CmpsCrv;    
}

/******************************************************************************
* Given a Bezier curve Crv and a Bezier surface Srf, computes the composition *
* Srf(Crv(t)).								      *
* See: "Freeform surfcae analysis using a hybrid of symbolic and numeric      *
* computation" by Gershon Elber, PhD thesis, University of Utah, 1992.	      *
******************************************************************************/
CagdCrvStruct *BzrComposeSrfCrv(CagdSrfStruct *Srf, CagdCrvStruct *Crv)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, k, l, CmpsOrder,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdRType
	Translate = 0.0;
    CagdCrvStruct **CrvUFactors, **CrvVFactors, *CrvW, *CrvX, *CrvY, *CrvZ,
	*CrvU, *CrvV,
    	*CmpsCrv = NULL;

    if (CAGD_NUM_OF_PT_COORD(Crv -> PType) != 2)
	FATAL_ERROR(CAGD_ERR_WRONG_PT_TYPE);
    CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    CrvU = CrvW == NULL ? CagdCrvCopy(CrvX)
		        : CagdCrvMergeScalar(CrvW, CrvX, NULL, NULL);
    CrvV = CrvW == NULL ? CagdCrvCopy(CrvY)
			: CagdCrvMergeScalar(CrvW, CrvY, NULL, NULL);
    if (CrvW)
	CagdCrvFree(CrvW);
    CagdCrvFree(CrvX);
    CagdCrvFree(CrvY);

    CrvUFactors = CagdComputeCurvePowers(CrvU, UOrder);
    CrvVFactors = CagdComputeCurvePowers(CrvV, VOrder);

    CmpsCrv = BzrCrvNew(CrvUFactors[0] -> Length + 
			CrvVFactors[0] -> Length - 1, Srf -> PType);
    CmpsOrder = CmpsCrv -> Order;

    for (k = !IsRational; k <= MaxCoord; k++) {
    	CagdRType
	    *CmpsPoints = CmpsCrv -> Points[k],
	    *SPoints = Srf -> Points[k];

	for (j = 0; j < VOrder; j++) {
	    for (i = 0; i < UOrder; i++) {
		CagdCrvStruct
		    *CTmp = CagdCrvMult(CrvUFactors[i], CrvVFactors[j]);
		CagdRType
		    *CTmpPoints = CTmp -> Points[1];

		CagdCrvTransform(CTmp, &Translate, *SPoints++);

		if (i == 0 && j == 0) {
		    CAGD_GEN_COPY(CmpsPoints, CTmpPoints,
				  CmpsOrder * sizeof(CagdRType));
		}
		else {
		    for (l = 0; l < CmpsOrder; l++)
			CmpsPoints[l] += CTmpPoints[l];
		}
	    }
	}
    }

    for (i = 0; i < UOrder; i++)
    	CagdCrvFree(CrvUFactors[i]);
    for (i = 0; i < VOrder; i++)
    	CagdCrvFree(CrvVFactors[i]);

    if (CAGD_IS_RATIONAL_CRV(Crv)) {
	CagdCrvStruct *CTmp,
	    *NewCrvW = CagdCrvMult(CrvUFactors[UOrder], CrvVFactors[VOrder]);

	CagdCrvSplitScalar(CmpsCrv, &CrvW, &CrvX, &CrvY, &CrvZ);
	CTmp = CagdCrvMergeScalar(NewCrvW, CrvX, CrvY, CrvZ);
	CagdCrvFree(NewCrvW);
	CagdCrvFree(CmpsCrv);
	CmpsCrv = CTmp;

	if (CrvX)
	    CagdCrvFree(CrvX);
	if (CrvY)
	    CagdCrvFree(CrvY);
	if (CrvZ)
	    CagdCrvFree(CrvZ);

	CagdCrvFree(CrvUFactors[UOrder]);
	CagdCrvFree(CrvVFactors[VOrder]);
    }

    IritFree((VoidPtr) CrvUFactors);
    IritFree((VoidPtr) CrvVFactors);

    CagdCrvFree(CrvU);
    CagdCrvFree(CrvV);

    return CmpsCrv;
}

/******************************************************************************
* Computes the factors of the Bernstein polynomials where t is a scalar curve.*
*									      *
*   n            n              n-i         i				      *
*  B (crv(t)) = ( ) (1 - crv(t))    (crv(t))				      *
*   i            i							      *
*									      *
* The curve crv(t) is a scalar, possibly rational curve.		      *
* The constant 1 is equal to wcrv(t) if crv(t) is rational. If rational, the  *
* returned vector, index Order will contain wcrv(t)^(Order-1).		      *
* See: "Freeform surface analysis using a hybrid of symbolic and numeric      *
* computation" by Gershon Elber, PhD thesis, University of Utah, 1992.	      *
******************************************************************************/
static CagdCrvStruct **CagdComputeCurvePowers(CagdCrvStruct *Crv, int Order)
{
    CagdBType
	IsRational = CAGD_IS_RATIONAL_CRV(Crv);
    int i;
    CagdCrvStruct *Crv_1,
	**CrvFactors1_Crv = (CagdCrvStruct **)
			    IritMalloc((Order + 1) * sizeof(CagdCrvStruct *)),
	**CrvFactorsCrv = (CagdCrvStruct **)
			    IritMalloc((Order + 1) * sizeof(CagdCrvStruct *)),
	**CrvFactors = (CagdCrvStruct **)
			    IritMalloc((Order + 1) * sizeof(CagdCrvStruct *));
    CagdRType *Points,
	Translate = 0.0;

    if (IsRational) {
	CagdCrvStruct *CrvX, *CrvY, *CrvZ, *CrvW, *CTmp;

	CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
	Crv_1 = CagdCrvSub(CrvW, CrvX);

	/* Set CrvFactors[Order] to CrvW(t)^(Order-1) if curve is rational. */
	CrvFactors[Order] = CagdCrvCopy(CrvW);
	for (i = 1; i < Order - 1; i++) {
	    CTmp = CagdCrvMult(CrvFactors[Order], CrvW);
	    CagdCrvFree(CrvFactors[Order]);
	    CrvFactors[Order] = CTmp;
	}

	CagdCrvFree(CrvW);
	CagdCrvFree(CrvX);
    }
    else {
    	Crv_1 = CagdCrvCopy(Crv);
	Points = Crv_1 -> Points[1];
    	for (i = 0; i < Crv -> Order; i++, Points++)
    	    *Points = 1.0 - *Points;

	CrvFactors[Order] = NULL;
    }

    for (i = 0; i < Order; i++) {
    	if (i == 0) {
    	    CrvFactors1_Crv[0] = NULL;
    	    CrvFactorsCrv[0] = NULL;
	}
    	else if (i == 1) {
    	    CrvFactors1_Crv[1] = Crv_1;
    	    CrvFactorsCrv[1] = CagdCrvCopy(Crv);
    	}
    	else {
    	    CrvFactors1_Crv[i] = CagdCrvMult(CrvFactors1_Crv[i - 1], Crv_1);
    	    CrvFactorsCrv[i] = CagdCrvMult(CrvFactorsCrv[i - 1], Crv);
    	}
    }

    for (i = 0; i < Order; i++) {
	if (i == 0) {
	    CrvFactors[i]= CagdCrvCopy(CrvFactors1_Crv[Order - 1]);
	}
	else if (i == Order - 1) {
	    CrvFactors[i] = CagdCrvCopy(CrvFactorsCrv[Order - 1]);
	}
	else {
	    CrvFactors[i] = CagdCrvMult(CrvFactors1_Crv[Order - 1 - i],
					CrvFactorsCrv[i]);
	}
	CagdCrvTransform(CrvFactors[i],
			 &Translate, CagdIChooseK(i, Order - 1));
    }

    for (i = 1; i < Order; i++) {
    	CagdCrvFree(CrvFactorsCrv[i]);
    	CagdCrvFree(CrvFactors1_Crv[i]);
    }
    IritFree((VoidPtr) CrvFactorsCrv);
    IritFree((VoidPtr) CrvFactors1_Crv);

    return CrvFactors;
}