/******************************************************************************
* BspBoehm.c - Implements Boehm single knot insertion routine.		      *
* Based on:								      *
* "Recursive proof of Boehm's knot insertion technique", by Phillip J Barry   *
* Ronald N Goldman, CAD, Volume 20 number 4 1988, pp 181-182.		      *
* The original paper by Boehm W. is also sighted there.			      *
*******************************************************************************
* Written by Gershon Elber, Aug. 90.					      *
******************************************************************************/

#ifdef __MSDOS__
#include <stdlib.h>
#endif /* __MSDOS__ */

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

/******************************************************************************
* Returns a new curve refined at t (t is inserted as a new knot in Crv).      *
* If however the multiplicity of t in the current knot vector is equal	      *
* (or greater!?) to the degree or t is not in the curve parametric domain, no *
* new knot is insert and NULL is returned instead.			      *
* Control mesh is updated as follows (P is old ctl polygon, Q is new):	      *
* Let Index be the last knot in old knot vector less than t and		      *
* let j be j = Index - order + 1. Also let k be the curve order.	      *
*									      *
*  Case 1: Q(i) = P(i), i <= j						      *
*									      *
*		  t     -    t(i)        t(i+k-1)  -   t		      *
*  case 2: Q(i) = --------------- P(i) + --------------- P(i-1), j<i<=Index   *
*		  t(i+k-1) - t(i)        t(i+k-1) - t(i)		      *
*									      *
*  case 3: Q(i) = P(i-1), Index < i					      *
*									      *
*  Note: Altough works, this is not the optimal way to insert many knot!      *
******************************************************************************/
CagdCrvStruct *BspCrvKnotInsert(CagdCrvStruct *Crv, CagdRType t)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    CagdRType
	*KnotVector = Crv -> KnotVector;
    int i, j,
	k = Crv -> Order,
	Len = Crv -> Length,
	KVLen = k + Len,
	Index = BspKnotLastIndexL(KnotVector, KVLen, t),
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct
	*RefinedCrv = CagdCrvNew(Crv -> GType, Crv -> PType, Len + 1);
    CagdRType
	**Points = Crv -> Points,
	**NewPoints = RefinedCrv -> Points;

    if (!BspKnotParamInDomain(KnotVector, Len, k, t))
	FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

    /* Update the rest of the RefinedCrv data structure (easy part). */
    RefinedCrv -> Order = k;
    RefinedCrv -> KnotVector = BspKnotInsertOne(Crv -> KnotVector, k, Len, t);

    /* Case 1: Copy all points upto (Index - Order + 1) as is. */
    for (i = IsNotRational; i <= MaxCoord; i++)
	GEN_COPY(NewPoints[i], Points[i], sizeof(CagdRType) * (Index - k + 2));

    /* Case 2: Convex blend of exactly 2 points. */
    for (j = Index - k + 2; j <= Index; j++)
	for (i = IsNotRational; i <= MaxCoord; i++)
	    NewPoints[i][j] =
	       ((t - KnotVector[j]) * Points[i][j] +
		(KnotVector[j + k - 1] - t) * Points[i][j - 1]) /
				      (KnotVector[j + k - 1] - KnotVector[j]);

    /* Case 3: Copy all points upto the end. */
    for (i = IsNotRational; i <= MaxCoord; i++)
	GEN_COPY(&NewPoints[i][Index + 1],
		 &Points[i][Index],
		 sizeof(CagdRType) * (Len - Index));

    return RefinedCrv;
}

/******************************************************************************
*  Insert n knot all with the value t. In no case will the multiplicity of    *
* knot be greater or equal to the curve order.				      *
*  Note: Altough works, this is not the optimal way to insert many knot!      *
******************************************************************************/
CagdCrvStruct *BspCrvKnotInsertNSame(CagdCrvStruct *Crv, CagdRType t, int n)
{
    int Mult = BspKnotFindMult(Crv -> KnotVector, Crv -> Order,
							    Crv -> Length, t);
    CagdCrvStruct *RefinedCrv,
        *TempCrv = Crv;

    n = MIN(n, Crv -> Order - Mult - 1);

    if (n > 0)
	while (n--) {
	    RefinedCrv = BspCrvKnotInsert(TempCrv, t);
	    if (TempCrv != Crv)
		CagdCrvFree(TempCrv);
	    TempCrv = RefinedCrv;
	}
    else
	RefinedCrv = CagdCrvCopy(Crv);

    return RefinedCrv;
}

/******************************************************************************
*  Insert n knot with different values as defined by t. If however Replace is *
* TRUE, the knot are simply replacing the current ones.			      *
******************************************************************************/
CagdCrvStruct *BspCrvKnotInsertNDiff(CagdCrvStruct *Crv, CagdBType Replace,
							   CagdRType *t, int n)
{
    int i;
    CagdCrvStruct *RefinedCrv,
        *TempCrv = Crv;

    if (Replace) {
	if (Crv -> Order + Crv -> Length != n)
	    FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
	for (i = 1; i < n; i++) if (t[i] < t[i - 1])
	    FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);

    	RefinedCrv = CagdCrvCopy(Crv);
	for (i = 0; i < n; i++) RefinedCrv -> KnotVector[i] = *t++;
    }
    else {
    	for (i = 0; i < n; i++) {
    	    RefinedCrv = BspCrvKnotInsert(TempCrv, t[i]);
    	    if (TempCrv != Crv) CagdCrvFree(TempCrv);
    	    TempCrv = RefinedCrv;
    	}
    }

    return RefinedCrv;
}

/******************************************************************************
* Returns a new srf refined at t (t is inserted as a new knot in Crv) in Dir. *
******************************************************************************/
CagdSrfStruct *BspSrfKnotInsert(CagdSrfStruct *BspSrf, CagdSrfDirType Dir,
								   CagdRType t)
{
    int Row, Col,
	ULength = BspSrf -> ULength,
	VLength = BspSrf -> VLength,
	UOrder = BspSrf -> UOrder,
	VOrder = BspSrf -> VOrder;
    CagdCrvStruct *Crv, *RefCrv;
    CagdSrfStruct
	*RefSrf = NULL;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    RefSrf = BspSrfNew(ULength + 1, VLength, UOrder, VOrder,
							      BspSrf -> PType);
	    CagdFree((VoidPtr) RefSrf -> UKnotVector);
	    CagdFree((VoidPtr) RefSrf -> VKnotVector);
	    RefSrf -> VKnotVector = BspKnotCopy(BspSrf -> VKnotVector,
					 BspSrf -> VLength + BspSrf -> VOrder);

	    for (Row = 0; Row < VLength; Row++) {
		Crv = BspSrfCrvFromMesh(BspSrf, Row, CAGD_CONST_V_DIR);
		RefCrv = BspCrvKnotInsert(Crv, t);

		if (Row == 0) {		  /* Figure out refined knot vector. */
		    RefSrf -> UKnotVector = BspKnotCopy(RefCrv -> KnotVector,
					   RefCrv -> Length + RefCrv -> Order);
		}

		CagdCrvToMesh(RefCrv, Row, CAGD_CONST_V_DIR, RefSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(RefCrv);
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    RefSrf = BspSrfNew(ULength, VLength + 1, UOrder, VOrder,
							     BspSrf -> PType);
	    CagdFree((VoidPtr) RefSrf -> UKnotVector);
	    CagdFree((VoidPtr) RefSrf -> VKnotVector);
	    RefSrf -> UKnotVector = BspKnotCopy(BspSrf -> UKnotVector,
					BspSrf -> ULength + BspSrf -> UOrder);

	    for (Col = 0; Col < ULength; Col++) {
		Crv = BspSrfCrvFromMesh(BspSrf, Col, CAGD_CONST_U_DIR);
		RefCrv = BspCrvKnotInsert(Crv, t);

		if (Col == 0) {		  /* Figure out refined knot vector. */
		    RefSrf -> VKnotVector = BspKnotCopy(RefCrv -> KnotVector,
					   RefCrv -> Length + RefCrv -> Order);
		}

		CagdCrvToMesh(RefCrv, Col, CAGD_CONST_U_DIR, RefSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(RefCrv);
	    }
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return RefSrf;
}

/******************************************************************************
*  Insert n knot all with the value t in direction Dir. In no case will the   *
* multiplicity of knot be greater or equal to the curve order.		      *
*  Note: Altough works, this is not the optimal way to insert many knot!      *
******************************************************************************/
CagdSrfStruct *BspSrfKnotInsertNSame(CagdSrfStruct *BspSrf, CagdSrfDirType Dir,
							    CagdRType t, int n)
{
    int Row, Col,
	ULength = BspSrf -> ULength,
	VLength = BspSrf -> VLength,
	UOrder = BspSrf -> UOrder,
	VOrder = BspSrf -> VOrder;
    CagdCrvStruct *Crv, *RefCrv;
    CagdSrfStruct
	*RefSrf = NULL;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    RefSrf = BspSrfNew(ULength + n, VLength, UOrder, VOrder,
							      BspSrf -> PType);
	    CagdFree((VoidPtr) RefSrf -> UKnotVector);
	    CagdFree((VoidPtr) RefSrf -> VKnotVector);
	    RefSrf -> VKnotVector = BspKnotCopy(BspSrf -> VKnotVector,
					 BspSrf -> VLength + BspSrf -> VOrder);

	    for (Row = 0; Row < VLength; Row++) {
		Crv = BspSrfCrvFromMesh(BspSrf, Row, CAGD_CONST_V_DIR);
		RefCrv = BspCrvKnotInsertNSame(Crv, t, n);

		if (Row == 0) {		  /* Figure out refined knot vector. */
		    RefSrf -> UKnotVector = BspKnotCopy(RefCrv -> KnotVector,
					   RefCrv -> Length + RefCrv -> Order);
		}

		CagdCrvToMesh(RefCrv, Row, CAGD_CONST_V_DIR, RefSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(RefCrv);
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    RefSrf = BspSrfNew(ULength, VLength + n, UOrder, VOrder,
							     BspSrf -> PType);
	    CagdFree((VoidPtr) RefSrf -> UKnotVector);
	    CagdFree((VoidPtr) RefSrf -> VKnotVector);
	    RefSrf -> UKnotVector = BspKnotCopy(BspSrf -> UKnotVector,
					BspSrf -> ULength + BspSrf -> UOrder);

	    for (Col = 0; Col < ULength; Col++) {
		Crv = BspSrfCrvFromMesh(BspSrf, Col, CAGD_CONST_U_DIR);
		RefCrv = BspCrvKnotInsertNSame(Crv, t, n);

		if (Col == 0) {		  /* Figure out refined knot vector. */
		    RefSrf -> VKnotVector = BspKnotCopy(RefCrv -> KnotVector,
					   RefCrv -> Length + RefCrv -> Order);
		}

		CagdCrvToMesh(RefCrv, Col, CAGD_CONST_U_DIR, RefSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(RefCrv);
	    }
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_SRF);
	    break;
    }

    return RefSrf;
}

/******************************************************************************
*  Insert n knot with different values as defined by t. If however Replace is *
* TRUE, the knot are simply replacing the current ones.			      *
******************************************************************************/
CagdSrfStruct *BspSrfKnotInsertNDiff(CagdSrfStruct *BspSrf, CagdSrfDirType Dir,
					      int Replace, CagdRType *t, int n)
{
    int i, Row, Col,
	ULength = BspSrf -> ULength,
	VLength = BspSrf -> VLength,
	UOrder = BspSrf -> UOrder,
	VOrder = BspSrf -> VOrder;
    CagdCrvStruct *Crv, *RefCrv;
    CagdSrfStruct
	*RefSrf = NULL;

    if (Replace) {
	for (i = 1; i < n; i++) if (t[i] < t[i - 1])
	    FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);

    	switch (Dir) {
	    case CAGD_CONST_U_DIR:
		if (BspSrf -> UOrder + BspSrf -> ULength != n)
		    FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);

		RefSrf = CagdSrfCopy(BspSrf);
		for (i = 0; i < n; i++) RefSrf -> UKnotVector[i] = *t++;
		break;
	    case CAGD_CONST_V_DIR:
		if (BspSrf -> VOrder + BspSrf -> VLength != n)
		    FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);

		RefSrf = CagdSrfCopy(BspSrf);
		for (i = 0; i < n; i++) RefSrf -> VKnotVector[i] = *t++;
		break;
	    default:
		FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
		break;

	}

	return RefSrf;
    }

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    RefSrf = BspSrfNew(ULength + n, VLength, UOrder, VOrder,
							      BspSrf -> PType);
	    CagdFree((VoidPtr) RefSrf -> UKnotVector);
	    CagdFree((VoidPtr) RefSrf -> VKnotVector);
	    RefSrf -> VKnotVector = BspKnotCopy(BspSrf -> VKnotVector,
					 BspSrf -> VLength + BspSrf -> VOrder);

	    for (Row = 0; Row < VLength; Row++) {
		Crv = BspSrfCrvFromMesh(BspSrf, Row, CAGD_CONST_V_DIR);
		RefCrv = BspCrvKnotInsertNDiff(Crv, FALSE, t, n);

		if (Row == 0) {		  /* Figure out refined knot vector. */
		    RefSrf -> UKnotVector = BspKnotCopy(RefCrv -> KnotVector,
					   RefCrv -> Length + RefCrv -> Order);
		}

		CagdCrvToMesh(RefCrv, Row, CAGD_CONST_V_DIR, RefSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(RefCrv);
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    RefSrf = BspSrfNew(ULength, VLength + n, UOrder, VOrder,
							     BspSrf -> PType);
	    CagdFree((VoidPtr) RefSrf -> UKnotVector);
	    CagdFree((VoidPtr) RefSrf -> VKnotVector);
	    RefSrf -> UKnotVector = BspKnotCopy(BspSrf -> UKnotVector,
					BspSrf -> ULength + BspSrf -> UOrder);

	    for (Col = 0; Col < ULength; Col++) {
		Crv = BspSrfCrvFromMesh(BspSrf, Col, CAGD_CONST_U_DIR);
		RefCrv = BspCrvKnotInsertNDiff(Crv, FALSE, t, n);

		if (Col == 0) {		  /* Figure out refined knot vector. */
		    RefSrf -> VKnotVector = BspKnotCopy(RefCrv -> KnotVector,
					   RefCrv -> Length + RefCrv -> Order);
		}

		CagdCrvToMesh(RefCrv, Col, CAGD_CONST_U_DIR, RefSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(RefCrv);
	    }
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return RefSrf;
}