/******************************************************************************
 * Blossom.c - Polar forms and Blossoming related functions.		      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, Feb. 99.					      *
******************************************************************************/

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

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes the Blossom over the given points, Pts, with knot sequence      M
* Knots, and Blossoming factors BlsmVals.  Evaluation is conducted via the   M
* Cox - De Boor algorithm with  a possibly different parameter at each       M
* iteration as prescribed via the Blossoming factors.  Note that the Bezier  M
* case is supported via the case for which the Knots are NULL.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Pts:        Coefficeints or scalar control points to blossom.            M
*   PtsStep:    Step size between coefficients, typically one.		     M
*   Order:      Order of the freeform geometry,				     M
*   Knots:      Knots of the freeform geometry.  If NULL assumed Bezier.     M
*   KnotsLen:   Length of Knots knot vectors.				     M
*   BlsmVals:   Blossoming values to consider.                               M
*   BlsmLen:    Length of BlsmVals vector.                                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType:    Evaluated Blossom                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdBlossomEval                                                          M
*****************************************************************************/
CagdRType CagdBlossomEval(CagdRType *Pts,
			  int PtsStep,
			  int Order,
			  CagdRType *Knots,
			  int KnotsLen,
			  CagdRType *BlsmVals,
			  int BlsmLen)
{
    STATIC_DATA int
	LastPtsStep = 0,
	LastBlsmLen = 0,
	LclVecLen = 0;
    STATIC_DATA CagdRType
	*LastPts = NULL,
	*LastKnots = NULL,
	LastResult = 0.0,
	*LastBlsmVals = NULL,
	*LclVec = NULL;
    int k, i, J,
        LastJ = -1;
    CagdRType
	TMin = Knots == NULL ? 0.0 : Knots[1],
	TMax = Knots == NULL ? 1.0 : Knots[KnotsLen - 2];

    /* Cache and compare last and present blossom values in case we could    */
    /* use old results.  Also verify Pts, PtsStep and Knots.		     */
    if (LastBlsmLen < BlsmLen) {
	if (LastBlsmVals != NULL)
	    IritFree(LastBlsmVals);
	LastBlsmVals = (CagdRType *) IritMalloc(sizeof(CagdRType) * BlsmLen);
    }
    else if (LastBlsmLen == BlsmLen &&
	     LastKnots == Knots && 
	     LastPts == Pts &&
	     LastPtsStep == PtsStep &&
	     GEN_CMP(LastBlsmVals, BlsmVals, sizeof(CagdRType) * BlsmLen) == 0)
	return LastResult;

    GEN_COPY(LastBlsmVals, BlsmVals, sizeof(CagdRType) * BlsmLen);
    LastBlsmLen = BlsmLen;
    LastKnots = Knots;
    LastPts = Pts;
    LastPtsStep = PtsStep;

    /* We are about to blend control points, in place - make a copy: */
    if (LclVecLen < KnotsLen + 1) {
	if (LclVec != NULL)
	    IritFree(LclVec);

	/* In order to prevent from doing this very often. */
	LclVecLen = (KnotsLen + 1) * 2;
	LclVec = (CagdRType *) IritMalloc(sizeof(CagdRType) * LclVecLen);
    }
    if (PtsStep == 1)
        GEN_COPY(LclVec, Pts, sizeof(CagdRType) * (KnotsLen - Order));
    else {
        for (i = 0; i < KnotsLen - Order; i++) {
	    LclVec[i] = *Pts;
	    Pts += PtsStep;
	}
    }

    if (Order <= BlsmLen)
	CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);

    for (k = 0; k < BlsmLen; k++) {
	CagdRType u, *R, *K,
	    BlsmVal = BlsmVals[k];
	int OrderK = Order - 1 - k;

	if (BlsmVal < TMin || BlsmVal > TMax)
	    CAGD_FATAL_ERROR(CAGD_ERR_WRONG_DOMAIN);

	if (Knots == NULL) {					  /* Bezier. */
	    u = BlsmVal;
	    J = Order - 1;

	    for (i = J, R = &LclVec[J]; i > J - OrderK && i >= 0; i--, R--)
		*R = i == 0 ? u * *R : u * *R + (1.0 - u) * R[-1];
	}
	else {							 /* Bspline. */
	    CagdRType
		ClippedBlsmVal = BlsmVal;

	    CAGD_VALIDATE_MIN_MAX_DOMAIN(ClippedBlsmVal, TMin, TMax);
	    J = BspKnotLastIndexLE(Knots, KnotsLen, ClippedBlsmVal);

	    if (J - LastJ > 1) {
		if (k == 0)
		    LastJ = J;
		else
		    J = LastJ + 1;
	    }
	    LastJ = J;

	    for (i = J, R = &LclVec[J], K = &Knots[i];
		 i > J - OrderK && i >= 0;
		 i--, R--, K--) {
		u = K[OrderK] - *K > 0 ? (BlsmVal - *K) / (K[OrderK] - *K)
		                       : 0.0;

		*R = i == 0 ? u * *R : u * *R + (1.0 - u) * R[-1];
	    }
	}
    }

    return LastResult = LclVec[J];
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes the Blossom over the given curve, Crv, and Blossoming factors   M
* BlsmVals.          					                     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        Curve to blossom.				             M
*   BlsmVals:   Blossoming values to consider.                               M
*   BlsmLen:    Length of BlsmVals vector; assumed less than curve order!    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCtlPtStruct *:    Evaluated Blossom                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdCrvBlossomEval                                                       M
*****************************************************************************/
CagdCtlPtStruct *CagdCrvBlossomEval(CagdCrvStruct *Crv,
				    CagdRType *BlsmVals,
				    int BlsmLen)
{
    int i;
    CagdCtlPtStruct
	*CtlPt = CagdCtlPtNew(Crv -> PType);

    for (i = !CAGD_IS_RATIONAL_PT(Crv -> PType);
	 i <= CAGD_NUM_OF_PT_COORD(Crv -> PType);
	 i++) {
        CtlPt -> Coords[i] =
	    CagdBlossomEval(Crv -> Points[i], 1, Crv -> Order,
			    CAGD_IS_BEZIER_CRV(Crv) ? NULL
						    : Crv -> KnotVector,
			    Crv -> Order + Crv -> Length,
			    BlsmVals, BlsmLen);
    }

    return CtlPt;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes a new curve with its degree raised to NewOrder, given curve     M
* Crv, using blossoming.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        Curve to degree raise.				             M
*   NewOrder:   New desired order of curve Crv.		    		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:    Degree raised curve, or NULL if error.               M
*                                                                            *
* SEE ALSO:                                                                  M
*   CagdCrvBlossomDegreeRaise, BspCrvDegreeRaise, BzrCrvDegreeRaise	     M
*   CagdCrvDegreeRaise, CagdSrfBlossomDegreeRaiseN		             M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdCrvBlossomDegreeRaiseN                                               M
*****************************************************************************/
CagdCrvStruct *CagdCrvBlossomDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder)
{
    int i,
	Order = Crv -> Order,
	DiffOrder = NewOrder - Order;

    if (DiffOrder <= 0) {
	CAGD_FATAL_ERROR(CAGD_ERR_WRONG_ORDER);
	return NULL;
    }

    for (i = Order; i < NewOrder; i++) {
        CagdCrvStruct
	    *RCrv = CagdCrvBlossomDegreeRaise(Crv);

	if (Crv -> Order > Order)
	    CagdCrvFree(Crv);
	Crv = RCrv;
    }

    return Crv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes a new curve with its degree raised once, given curve Crv, using M
* blossoming.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        Curve to degree raise.				             M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:    Degree raised curve, or NULL if error.               M
*                                                                            *
* SEE ALSO:                                                                  M
*   CagdCrvBlossomDegreeRaiseN, BspCrvDegreeRaise, BzrCrvDegreeRaise	     M
*   CagdCrvDegreeRaise, CagdSrfBlossomDegreeRaise                            M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdCrvBlossomDegreeRaise                                                M
*****************************************************************************/
CagdCrvStruct *CagdCrvBlossomDegreeRaise(CagdCrvStruct *Crv)
{
    CagdBType
	IsBezier = FALSE,
	NewCrv = FALSE;
    CagdPointType
	PType = Crv -> PType;
    int i, j, k, l, m, NewLen,
	Len = Crv -> Length,
	Order = Crv -> Order;
    CagdRType *NewKV, *KV, *BlossomValues;
    CagdCrvStruct *RCrv, *TCrv;

    if (CAGD_IS_BEZIER_CRV(Crv)) {            /* Convert to a Bspline curve. */
	NewCrv = TRUE;
	IsBezier = TRUE;
	Crv = CnvrtBezier2BsplineCrv(Crv);
    }
    else if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }
    if (!BspCrvHasOpenEC(Crv)) {
	TCrv = CnvrtFloat2OpenCrv(Crv);
	if (NewCrv)
	    CagdCrvFree(Crv);
	Crv = TCrv;
	NewCrv = TRUE;
    }
    KV = Crv -> KnotVector;

    /* We now have an open end Bspline curve to deal with. */

    /* Allocate a vector for the new knot vector that will be sufficient and */
    /* populate it:							     */
    NewKV = BspKnotDegreeRaisedKV(KV, Len, Order, Order + 1, &NewLen);
    NewLen -= Order + 1; /* Get the number of control points we should have. */

    RCrv = BspCrvNew(NewLen, Order + 1, PType);
    GEN_COPY(RCrv -> KnotVector, NewKV,
	     sizeof(CagdRType) * (NewLen + RCrv -> Order));

    BlossomValues = (CagdRType *) IritMalloc(sizeof(CagdRType) * Order);

    for (i = !CAGD_IS_RATIONAL_PT(PType);
	 i <= CAGD_NUM_OF_PT_COORD(PType);
	 i++) {
	CagdRType
	    *Points = Crv -> Points[i],
	    *NewPoints = RCrv -> Points[i];

	/* Evaluate the degree raised control points via the symmetric sum   */
	/* of Order blossoms of original curve's order.			     */
        for (l = 0; l < NewLen; l++) {
	    *NewPoints = 0.0;

	    for (j = 0; j < Order; j++) {
	        for (k = m = 0; k < Order; k++) {
		    if (k != j)
		        BlossomValues[m++] = NewKV[l + 1 + k];
		}

		*NewPoints += CagdBlossomEval(Points, 1, Order,
					      KV, Len + Order,
					      BlossomValues, Order - 1);
	    }
	    *NewPoints++ /= Order;
	}
    }

    if (IsBezier) {
	TCrv = CnvrtBspline2BezierCrv(RCrv);
	CagdCrvFree(RCrv);
	RCrv = TCrv;
    }

    if (NewCrv)
	CagdCrvFree(Crv);

    IritFree(BlossomValues);
    IritFree(NewKV);

    return RCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes a new surface with its degree raised to NewOrder, given surface M
* Srf, using blossoming.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Surface to degree raise.				     M
*   NewUOrder:  New U order of Srf.					     M
*   NewVOrder:  New V order of Srf.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:    Degree raised surface, or NULL if error.             M
*                                                                            *
* SEE ALSO:                                                                  M
*   CagdSrfBlossomDegreeRaise, BspSrfDegreeRaise, BzrSrfDegreeRaise	     M
*   CagdSrfDegreeRaise, CagdCrvBlossomDegreeRaiseN	                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdSrfBlossomDegreeRaiseN                                               M
*****************************************************************************/
CagdSrfStruct *CagdSrfBlossomDegreeRaiseN(CagdSrfStruct *Srf,
					  int NewUOrder,
					  int NewVOrder)
{
    int i,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	DiffUOrder = NewUOrder - UOrder,
	DiffVOrder = NewVOrder - VOrder;

    if (DiffUOrder <= 0 && DiffVOrder <= 0) {
	CAGD_FATAL_ERROR(CAGD_ERR_WRONG_ORDER);
	return NULL;
    }

    Srf = CagdSrfCopy(Srf);

    for (i = UOrder; i < NewUOrder; i++) {
        CagdSrfStruct
	    *RSrf = CagdSrfBlossomDegreeRaise(Srf, CAGD_CONST_U_DIR);

	CagdSrfFree(Srf);
	Srf = RSrf;
    }
    for (i = VOrder; i < NewVOrder; i++) {
        CagdSrfStruct
	    *RSrf = CagdSrfBlossomDegreeRaise(Srf, CAGD_CONST_V_DIR);

	CagdSrfFree(Srf);
	Srf = RSrf;
    }

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes a new surface with its degree raised once, given surface Srf,   M
* using blossoming.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Surface to degree raise.			             M
*   Dir:       Direction of degree raising. Either U or V.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:    Degree raised surface, or NULL if error.             M
*                                                                            *
* SEE ALSO:                                                                  M
*   CagdSrfBlossomDegreeRaiseN, BspSrfDegreeRaise, BzrSrfDegreeRaise	     M
*   CagdSrfDegreeRaise, CagdCrvBlossomDegreeRaise                            M
*                                                                            *
* KEYWORDS:                                                                  M
*   CagdSrfBlossomDegreeRaise                                                M
*****************************************************************************/
CagdSrfStruct *CagdSrfBlossomDegreeRaise(CagdSrfStruct *Srf,
					 CagdSrfDirType Dir)
{
    CagdBType
	IsBezier = FALSE,
	NewSrf = FALSE;
    CagdPointType
	PType = Srf -> PType;
    int i, j, k, l, m, NewLen, Len, Order, StepSize;
    CagdRType *NewKV, *KV, *BlossomValues;
    CagdSrfStruct *RSrf, *TSrf;

    if (CAGD_IS_BEZIER_SRF(Srf)) {          /* Convert to a Bspline surface. */
	NewSrf = TRUE;
	IsBezier = TRUE;
	Srf = CnvrtBezier2BsplineSrf(Srf);
    }
    else if (CAGD_IS_PERIODIC_SRF(Srf)) {
	NewSrf = TRUE;
	Srf = CnvrtPeriodic2FloatSrf(Srf);
    }
    if (!BspSrfHasOpenEC(Srf)) {
	TSrf = CnvrtFloat2OpenSrf(Srf);
	if (NewSrf)
	    CagdSrfFree(Srf);
	Srf = TSrf;
	NewSrf = TRUE;
    }

    /* We now have an open end Bspline surface to deal with. */

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    KV = Srf -> UKnotVector;
	    Len = Srf -> ULength;
	    Order = Srf -> UOrder;

	    /* Allocate a vector for the new knot vector that will be        */
	    /* sufficient and populate it:				     */
	    NewKV = BspKnotDegreeRaisedKV(KV, Len, Order, Order + 1, &NewLen);
	    NewLen -= Order + 1; /* Get num. of control pnts we should have. */

	    RSrf = BspSrfNew(NewLen, Srf -> VLength,
			     Order + 1, Srf -> VOrder, PType);
	    GEN_COPY(RSrf -> UKnotVector, NewKV,
		     sizeof(CagdRType) * (NewLen + RSrf -> UOrder));
	    GEN_COPY(RSrf -> VKnotVector, Srf -> VKnotVector,
		     sizeof(CagdRType) * (Srf -> VLength + Srf -> VOrder));
	    StepSize = CAGD_NEXT_U(Srf);
	    break;
        case CAGD_CONST_V_DIR:
	    KV = Srf -> VKnotVector;
	    Len = Srf -> VLength;
	    Order = Srf -> VOrder;

	    /* Allocate a vector for the new knot vector that will be        */
	    /* sufficient and populate it:				     */
	    NewKV = BspKnotDegreeRaisedKV(KV, Len, Order, Order + 1, &NewLen);
	    NewLen -= Order + 1; /* Get num. of control pnts we should have. */

	    RSrf = BspSrfNew(Srf -> ULength, NewLen,
			     Srf -> UOrder, Order + 1, PType);
	    GEN_COPY(RSrf -> UKnotVector, Srf -> UKnotVector,
		     sizeof(CagdRType) * (Srf -> ULength + Srf -> UOrder));
	    GEN_COPY(RSrf -> VKnotVector, NewKV,
		     sizeof(CagdRType) * (NewLen + RSrf -> VOrder));
	    StepSize = CAGD_NEXT_V(Srf);
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    BlossomValues = (CagdRType *) IritMalloc(sizeof(CagdRType) * Order);

    for (i = !CAGD_IS_RATIONAL_PT(PType);
	 i <= CAGD_NUM_OF_PT_COORD(PType);
	 i++) {
        int n;
	CagdRType *NewPoints,
	    *Points = Srf -> Points[i];

	/* Evaluate the degree raised control points via the symmetric sum   */
	/* of Order blossoms of original surface's order.		     */
	switch (Dir) {
	    case CAGD_CONST_U_DIR:
	        for (n = 0; n < Srf -> VLength; n++) {
		    NewPoints = &RSrf -> Points[i][CAGD_MESH_UV(RSrf, 0, n)];
		    for (l = 0; l < NewLen; l++) {
		        *NewPoints = 0.0;

			if (Order > 1) {
			    for (j = 0; j < Order; j++) {
			        for (k = m = 0; k < Order; k++) {
				    if (k != j)
				        BlossomValues[m++] = NewKV[l + 1 + k];
				}

				*NewPoints += CagdBlossomEval(Points, StepSize,
							      Order, KV,
							      Len + Order,
							      BlossomValues,
							      Order - 1);
			    }
			    *NewPoints /= Order;
			}
			else
			    *NewPoints = *Points;

			NewPoints += CAGD_NEXT_U(RSrf);
		    }
		    Points += CAGD_NEXT_V(Srf);
		}
		break;
	    case CAGD_CONST_V_DIR:
		for (n = 0; n < Srf -> ULength; n++) {
		    NewPoints = &RSrf -> Points[i][CAGD_MESH_UV(RSrf, n, 0)];
		    for (l = 0; l < NewLen; l++) {
		        *NewPoints = 0.0;

			if (Order > 1) {
			    for (j = 0; j < Order; j++) {
			        for (k = m = 0; k < Order; k++) {
				    if (k != j)
				        BlossomValues[m++] = NewKV[l + 1 + k];
				}

				*NewPoints += CagdBlossomEval(Points, StepSize,
							      Order, KV,
							      Len + Order,
							      BlossomValues,
							      Order - 1);
			    }
			    *NewPoints /= Order;
			}
			else
			    *NewPoints = *Points;

			NewPoints += CAGD_NEXT_V(RSrf);
		    }
		    Points += CAGD_NEXT_U(Srf);
		}
		break;
	}
    }

    if (IsBezier) {
	TSrf = CnvrtBspline2BezierSrf(RSrf);
	CagdSrfFree(RSrf);
	RSrf = TSrf;
    }

    if (NewSrf)
	CagdSrfFree(Srf);

    IritFree(BlossomValues);
    IritFree(NewKV);

    return RSrf;
}