/******************************************************************************
* SBsp_Sym.c - Bspline surface symbolic computation.			      *
*******************************************************************************
* Written by Gershon Elber, Sep. 92.					      *
******************************************************************************/

#include "cagd_loc.h"

static CagdCrvStruct *BspCrvMultAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2);
static CagdSrfStruct *BspSrfMultAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2);

/******************************************************************************
* Given two bspline curves - multiply them coordinatewise.		      *
* The two curves are promoted to same point type before the multiplication    *
* can take place. The bspline curves are subdivided into piecewise bezier     *
* curves, each is multiplied individually and then all curves are merged.     *
* Return a bspline curves representing their product.			      *
******************************************************************************/
CagdCrvStruct *BspCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *ProdCrv, *TCrv;

    Crv1 = CagdCrvCopy(Crv1);
    Crv2 = CagdCrvCopy(Crv2);
    if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE))
	FATAL_ERROR(CAGD_ERR_CRV_FAIL_CMPT);

    TCrv = BspCrvMultAux(Crv1, Crv2);

    if (CAGD_IS_BEZIER_CRV(TCrv)) {
	ProdCrv = CnvrtBezier2BsplineCrv(TCrv);
	CagdCrvFree(TCrv);
    }
    else
	ProdCrv = TCrv;

    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    return ProdCrv;
}

/******************************************************************************
* Auxiliary routine. Subdivides the curves into Bezier curves, multiply       *
* the Bezier curves and merge them back. All is done simultaneously.          *
******************************************************************************/
static CagdCrvStruct *BspCrvMultAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *Crv1a, *Crv1b, *Crv2a, *Crv2b, *CrvA, *CrvB, *ProdCrv;

    if (Crv1 -> Order != Crv1 -> Length ||
	Crv2 -> Order != Crv2 -> Length) {
	CagdRType SubdivVal = Crv1 -> Order != Crv1 -> Length ?
				Crv1 -> KnotVector[(Crv1 -> Length +
						    Crv1 -> Order) / 2] :
				Crv2 -> KnotVector[(Crv2 -> Length +
						    Crv2 -> Order) / 2];

	/* Subdivide. */
	Crv1a = BspCrvSubdivAtParam(Crv1, SubdivVal);
	Crv1b = Crv1a -> Pnext;
	Crv1a -> Pnext = NULL;

	Crv2a = BspCrvSubdivAtParam(Crv2, SubdivVal);
	Crv2b = Crv2a -> Pnext;
	Crv2a -> Pnext = NULL;

	CrvA = BspCrvMultAux(Crv1a, Crv2a);
	CrvB = BspCrvMultAux(Crv1b, Crv2b);
	CagdCrvFree(Crv1a);
	CagdCrvFree(Crv1b);
	CagdCrvFree(Crv2a);
	CagdCrvFree(Crv2b);

	ProdCrv = CagdMergeCrvCrv(CrvA, CrvB);
	CagdCrvFree(CrvA);
	CagdCrvFree(CrvB);
    }
    else {
	int i;
	CagdRType TMin, TMax;
	CagdCrvStruct
	    *Crv1Bzr = CnvrtBspline2BezierCrv(Crv1),
	    *Crv2Bzr = CnvrtBspline2BezierCrv(Crv2),
	    *ProdCrvAux = BzrCrvMult(Crv1Bzr, Crv2Bzr);

	CagdCrvDomain(Crv1, &TMin, &TMax);
	ProdCrv = CnvrtBezier2BsplineCrv(ProdCrvAux);
	for (i = 0; i < ProdCrv -> Order; i++) {
	    ProdCrv -> KnotVector[i] = TMin;
	    ProdCrv -> KnotVector[i + ProdCrv -> Order] = TMax;
	}

	CagdCrvFree(Crv1Bzr);
	CagdCrvFree(Crv2Bzr);
	CagdCrvFree(ProdCrvAux);
    }

    return ProdCrv;
}

/******************************************************************************
* Given two bspline surfaces - multiply them coordinatewise.		      *
* The two surfaces are promoted to same point type before the multiplication  *
* can take place. The bspline surfaces are subdivided into piecewise bezier   *
* patches, each is multiplied individually and then all patches are merged.   *
* Return a bspline surface representing their product.			      *
******************************************************************************/
CagdSrfStruct *BspSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *ProdSrf, *TSrf;

    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);
    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE))
	FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);

    TSrf = BspSrfMultAux(Srf1, Srf2);

    if (CAGD_IS_BEZIER_SRF(TSrf)) {
	ProdSrf = CnvrtBezier2BsplineSrf(TSrf);
	CagdSrfFree(TSrf);
    }
    else
	ProdSrf = TSrf;

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return ProdSrf;
}

/******************************************************************************
* Auxiliary routine. Subdivides the surfaces into Bezier patches, multiply    *
* the Bezier patches and merge them back. All is done simultaneously.         *
******************************************************************************/
static CagdSrfStruct *BspSrfMultAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *Srf1a, *Srf1b, *Srf2a, *Srf2b, *SrfA, *SrfB, *ProdSrf;

    if (Srf1 -> UOrder != Srf1 -> ULength ||
	Srf2 -> UOrder != Srf2 -> ULength) {
	CagdRType SubdivVal = Srf1 -> UOrder != Srf1 -> ULength ?
				Srf1 -> UKnotVector[(Srf1 -> ULength +
						     Srf1 -> UOrder) / 2] :
				Srf2 -> UKnotVector[(Srf2 -> ULength +
						     Srf2 -> UOrder) / 2];

	/* Subdivide along U. */
	Srf1a = BspSrfSubdivAtParam(Srf1, SubdivVal, CAGD_CONST_U_DIR);
	Srf1b = Srf1a -> Pnext;
	Srf1a -> Pnext = NULL;

	Srf2a = BspSrfSubdivAtParam(Srf2, SubdivVal, CAGD_CONST_U_DIR);
	Srf2b = Srf2a -> Pnext;
	Srf2a -> Pnext = NULL;

	SrfA = BspSrfMultAux(Srf1a, Srf2a);
	SrfB = BspSrfMultAux(Srf1b, Srf2b);
	CagdSrfFree(Srf1a);
	CagdSrfFree(Srf1b);
	CagdSrfFree(Srf2a);
	CagdSrfFree(Srf2b);

	ProdSrf = CagdMergeSrfSrf(SrfA, SrfB, CAGD_CONST_U_DIR, TRUE);
	CagdSrfFree(SrfA);
	CagdSrfFree(SrfB);
    }
    else if (Srf1 -> VOrder != Srf1 -> VLength) {
	CagdRType SubdivVal = Srf1 -> VOrder != Srf1 -> VLength ?
				Srf1 -> VKnotVector[(Srf1 -> VLength +
						     Srf1 -> VOrder) / 2] :
				Srf2 -> VKnotVector[(Srf2 -> VLength +
						     Srf1 -> VOrder) / 2];

	/* Subdivide along V. */
	/* Subdivide along U. */
	Srf1a = BspSrfSubdivAtParam(Srf1, SubdivVal, CAGD_CONST_V_DIR);
	Srf1b = Srf1a -> Pnext;
	Srf1a -> Pnext = NULL;

	Srf2a = BspSrfSubdivAtParam(Srf2, SubdivVal, CAGD_CONST_V_DIR);
	Srf2b = Srf2a -> Pnext;
	Srf2a -> Pnext = NULL;

	SrfA = BspSrfMultAux(Srf1a, Srf2a);
	SrfB = BspSrfMultAux(Srf1b, Srf2b);
	CagdSrfFree(Srf1a);
	CagdSrfFree(Srf1b);
	CagdSrfFree(Srf2a);
	CagdSrfFree(Srf2b);

	ProdSrf = CagdMergeSrfSrf(SrfA, SrfB, CAGD_CONST_V_DIR, TRUE);
	CagdSrfFree(SrfA);
	CagdSrfFree(SrfB);
    }
    else {
	int i;
	CagdRType UMin, UMax, VMin, VMax;
	CagdSrfStruct
	    *Srf1Bzr = CnvrtBspline2BezierSrf(Srf1),
	    *Srf2Bzr = CnvrtBspline2BezierSrf(Srf2),
	    *ProdSrfAux = BzrSrfMult(Srf1Bzr, Srf2Bzr);

	CagdSrfDomain(Srf1, &UMin, &UMax, &VMin, &VMax);
	ProdSrf = CnvrtBezier2BsplineSrf(ProdSrfAux);
	for (i = 0; i < ProdSrf -> UOrder; i++) {
	    ProdSrf -> UKnotVector[i] = UMin;
	    ProdSrf -> UKnotVector[i + ProdSrf -> UOrder] = UMax;
	}
	for (i = 0; i < ProdSrf -> VOrder; i++) {
	    ProdSrf -> VKnotVector[i] = VMin;
	    ProdSrf -> VKnotVector[i + ProdSrf -> VOrder] = VMax;
	}

	CagdSrfFree(Srf1Bzr);
	CagdSrfFree(Srf2Bzr);
	CagdSrfFree(ProdSrfAux);
    }

    return ProdSrf;
}

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

    CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    if (CrvW)
	DCrvW = BspCrvDerive(CrvW);
    else
	FATAL_ERROR(CAGD_ERR_RATIONAL_EXPECTED);

    if (CrvX) {
	DCrvX = BspCrvDerive(CrvX);

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

	CagdCrvFree(CrvX);
	CrvX = CagdCrvSub(TCrv1, TCrv2);
	CagdCrvFree(TCrv1);
	CagdCrvFree(TCrv2);
    }
    if (CrvY) {
	DCrvY = BspCrvDerive(CrvY);

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

	CagdCrvFree(CrvY);
	CrvY = CagdCrvSub(TCrv1, TCrv2);
	CagdCrvFree(TCrv1);
	CagdCrvFree(TCrv2);
    }
    if (CrvZ) {
	DCrvZ = BspCrvDerive(CrvZ);

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

	CagdCrvFree(CrvZ);
	CrvZ = CagdCrvSub(TCrv1, TCrv2);
	CagdCrvFree(TCrv1);
	CagdCrvFree(TCrv2);
    }
    
    TCrv1 = BspCrvMult(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);
	CagdCrvFree(DCrvX);
    }
    if (CrvY) {
	CagdCrvFree(CrvY);
	CagdCrvFree(DCrvY);
    }
    if (CrvZ) {
	CagdCrvFree(CrvZ);
	CagdCrvFree(DCrvZ);
    }
    if (CrvW) {
	CagdCrvFree(CrvW);
	CagdCrvFree(DCrvW);
    }

    return DeriveCrv;
}

/******************************************************************************
* Given a rational bezier surface - computes its derivative curve using the   *
* quotient rule for differentiation.					      *
******************************************************************************/
CagdSrfStruct *BspSrfDeriveRational(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 = BspSrfDerive(SrfW, Dir);
    else
	FATAL_ERROR(CAGD_ERR_RATIONAL_EXPECTED);

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

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

	CagdSrfFree(SrfX);
	SrfX = CagdSrfSub(TSrf1, TSrf2);
	CagdSrfFree(TSrf1);
	CagdSrfFree(TSrf2);
    }
    if (SrfY) {
	DSrfY = BspSrfDerive(SrfY, Dir);

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

	CagdSrfFree(SrfY);
	SrfY = CagdSrfSub(TSrf1, TSrf2);
	CagdSrfFree(TSrf1);
	CagdSrfFree(TSrf2);
    }
    if (SrfZ) {
	DSrfZ = BspSrfDerive(SrfZ, Dir);

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

	CagdSrfFree(SrfZ);
	SrfZ = CagdSrfSub(TSrf1, TSrf2);
	CagdSrfFree(TSrf1);
	CagdSrfFree(TSrf2);
    }
    
    TSrf1 = BspSrfMult(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);
	CagdSrfFree(DSrfX);
    }
    if (SrfY) {
	CagdSrfFree(SrfY);
	CagdSrfFree(DSrfY);
    }
    if (SrfZ) {
	CagdSrfFree(SrfZ);
	CagdSrfFree(DSrfZ);
    }
    if (SrfW) {
	CagdSrfFree(SrfW);
	CagdSrfFree(DSrfW);
    }

    return DeriveSrf;
}