/******************************************************************************
* SBsp_Sym.c - Bspline surface symbolic computation.			      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, Sep. 92.					      *
******************************************************************************/

#include "symb_loc.h"

#define NODE_EQUAL_SHIFT 0.8

static int
    BspMultUsingInterpolation = TRUE;

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

/*****************************************************************************
* DESCRIPTION:                                                               M
* Sets method of Bspline product computation.                                M
*                                                                            *
* PARAMETERS:                                                                M
*   BspMultUsingInter:  If TRUE, Bspline product is computed by setting an   M
*                       interpolation problem. Otherwise, by decomposing the M
*                       Bspline geometry to Bezier geometry.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:         Previous setting.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspMultInterpFlag, product 		                                     M
*****************************************************************************/
int BspMultInterpFlag(int BspMultUsingInter)
{
    int i = BspMultUsingInterpolation;

    BspMultUsingInterpolation = BspMultUsingInter;
    return i;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two Bspline curves - multiply them coordinatewise.		     M
*   The two curves are promoted to same point type before the multiplication M
* can take place.							     M
*   See also BspMultInterpFlag.                                              M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv1, Crv2:   The two curves to multiply.                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  The product Crv1 * Crv2 coordinatewise.                M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvMult, product                                                      M
*****************************************************************************/
CagdCrvStruct *BspCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
{
    CagdCrvStruct *ProdCrv, *TCrv;

    Crv1 = CagdCrvCopy(Crv1);
    Crv2 = CagdCrvCopy(Crv2);

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

    if (BspMultUsingInterpolation) {
	CagdPointType
	    PType = Crv1 -> PType;
	CagdBType
	    IsRational = CAGD_IS_RATIONAL_PT(PType);
	int i, j, KVLen, ResLength,
	    NumCoords = CAGD_NUM_OF_PT_COORD(PType),
	    ResOrder = Crv1 -> Order + Crv2 -> Order - 1;
	CagdRType *R,
	    *KV = BspKnotContinuityMergeTwo(Crv1 -> KnotVector,
					    Crv1 -> Length + Crv1 -> Order,
					    Crv1 -> Order,
					    Crv2 -> KnotVector,
					    Crv2 -> Length + Crv2 -> Order,
					    Crv2 -> Order,
					    ResOrder, &KVLen),
	    *KVNodes = BspKnotNodes(KV, KVLen, ResOrder);
	CagdCtlPtStruct
	    *CtlPt = NULL,
	    *CtlPtList = NULL;

	ResLength = KVLen - ResOrder;

	/* Verify that all nodes are distinct. */
	for (i = 0, R = KVNodes; i < ResLength - 1; i++, R++) {
	    if (APX_EQ(R[0], R[1])) {
		if (i > 0)
		    R[0] = R[0] * NODE_EQUAL_SHIFT +
			   R[-1] * (1 - NODE_EQUAL_SHIFT);
	    }
	}

	/* Evaluate the multiplication at the node values. */
	for (i = 0, R = KVNodes; i < ResLength; i++, R++) {
	    CagdRType *Evl;

	    if (CtlPt == NULL)
		CtlPt = CtlPtList = CagdCtlPtNew(PType);
	    else {
		CtlPt -> Pnext = CagdCtlPtNew(PType);
		CtlPt = CtlPt -> Pnext;
	    }

	    Evl = CagdCrvEval(Crv1, *R);
	    SYMB_GEN_COPY(CtlPt -> Coords, Evl,
			  CAGD_MAX_PT_SIZE * sizeof(CagdRType));
	    Evl = CagdCrvEval(Crv2, *R);
	    for (j = !IsRational; j <= NumCoords; j++)
	        CtlPt -> Coords[j] *= Evl[j];
	}

	ProdCrv = BspCrvInterpolate(CtlPtList, ResLength, KVNodes, KV,
				    ResLength, ResOrder, FALSE);
	IritFree((VoidPtr) KVNodes);
	CagdCtlPtFreeList(CtlPtList);
    }
    else {
	TCrv = BspCrvMultAux(Crv1, Crv2);

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

    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    return ProdCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Auxiliary routine. Subdivides the curves into Bezier curves, multiply      *
* the Bezier curves and merge them back. All is done simultaneously.         *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv1, Crv2:   The two curves to multiply.                                *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:  The product Crv1 * Crv2 coordinatewise.                *
*****************************************************************************/
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, FALSE);
	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;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given two Bspline surfaces - multiply them coordinatewise.		     M
*   The two surfaces are promoted to same point type before multiplication   M
* can take place.							     M
*   See also BspMultInterpFlag.                                              M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf1, Srf2:   The two surfaces to multiply.		                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  The product Srf1 * Srf2 coordinatewise.                M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfMult, product                                                      M
*****************************************************************************/
CagdSrfStruct *BspSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2)
{
    CagdSrfStruct *ProdSrf, *TSrf;

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

/*****************************************************************************
* DESCRIPTION:                                                               *
* Auxiliary routine. Subdivides the surfaces into Bezier surfaces, multiply  *
* the Bezier surfaces and merge them back. All is done simultaneously.       *
*                                                                            *
* PARAMETERS:                                                                *
*   Srf1, Srf2:   The two surfaces to multiply.                              *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdSrfStruct *:  The product Srf1 * Srf2 coordinatewise.                *
*****************************************************************************/
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, FALSE, FALSE);
	CagdSrfFree(SrfA);
	CagdSrfFree(SrfB);
    }
    else if (Srf1 -> VOrder != Srf1 -> VLength ||
	     Srf2 -> VOrder != Srf2 -> VLength) {
	CagdRType
	    SubdivVal = Srf1 -> VOrder != Srf1 -> VLength ?
				Srf1 -> VKnotVector[(Srf1 -> VLength +
						     Srf1 -> VOrder) / 2] :
				Srf2 -> VKnotVector[(Srf2 -> VLength +
						     Srf2 -> VOrder) / 2];

	/* Subdivide along V. */
	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, FALSE, FALSE);
	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;
}

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

    SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
    if (CrvW)
	DCrvW = BspCrvDerive(CrvW);
    else
	SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED);

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

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

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

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

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

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

	CagdCrvFree(CrvZ);
	CrvZ = SymbCrvSub(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))
	SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT);

    DeriveCrv = SymbCrvMergeScalar(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;
}

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

    SymbSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ);
    if (SrfW)
	DSrfW = BspSrfDerive(SrfW, Dir);
    else
	SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED);

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

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

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

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

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

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

	CagdSrfFree(SrfZ);
	SrfZ = SymbSrfSub(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))
	SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT);

    DeriveSrf = SymbSrfMergeScalar(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;
}