/******************************************************************************
* 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.					      *
******************************************************************************/

#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++)
	CAGD_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++)
	CAGD_GEN_COPY(&NewPoints[i][Index + 1],
		      &Points[i][Index],
		      sizeof(CagdRType) * (Len - Index));

    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);
	    IritFree((VoidPtr) RefSrf -> UKnotVector);
	    IritFree((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);
	    IritFree((VoidPtr) RefSrf -> UKnotVector);
	    IritFree((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;
}