/******************************************************************************
* CagdAIso.c - adaptive isoline surface extraction algorithm.		      *
*******************************************************************************
* Written by Gershon Elber, Aug. 92.					      *
******************************************************************************/

#include "cagd_loc.h"

#define CONVEX_BLEND(v1, v2, b1, b2, t) { \
		CagdRType \
		    Blend = -b1 / ( b2 - b1 ); \
		t = v2 * Blend + v1 * (1.0 - Blend); \
	    }

static int MinSubdivLevel = 1;

static CagdCrvStruct *CagdAdapIsoExtractAux(int Level,
				     CagdSrfStruct *Srf, CagdSrfStruct *NSrf,
				     CagdCrvStruct *Crv1, CagdCrvStruct *NCrv1,
				     CagdCrvStruct *Crv2, CagdCrvStruct *NCrv2,
				     CagdRType Crv1Param, CagdRType Crv2Param,
				     CagdSrfDirType Dir, CagdRType Eps2,
				     CagdBType FullIso, CagdBType SinglePath);

/******************************************************************************
* Sets minimum level of subdivision forced in the adaptive iso extraction.    *
******************************************************************************/
void CagdSetAdapIsoExtractMinLevel(int MinLevel)
{
    MinSubdivLevel = MinLevel;
}

/******************************************************************************
* Extract a valid coverage set of isolines from the given surface in the      *
* given direction and epsilon.						      *
*   If FullIso is TRUE, all extracted isocurves are spanning the entire       *
* parametric domain.							      *
*   If SinglePath is TRUE, the entire coverage is going to be a single curve  *
*   If NSrf != NULL, every second curve will be a vector field curve rep. the *
* unnormalized normal for the previous Euclidean curve.	This mode disable the *
* SinglePath mode.							      *
******************************************************************************/
CagdCrvStruct *CagdAdapIsoExtract(CagdSrfStruct *Srf, CagdSrfStruct *NSrf,
				  CagdSrfDirType Dir,
				  CagdRType Eps, CagdBType FullIso,
				  CagdBType SinglePath)
{
    CagdBType
	SrfBezier = FALSE;
    CagdRType Crv1Param, Crv2Param;
    CagdCrvStruct *Crv1, *Crv2, *NCrv1, *NCrv2, *AllAdapIso, *TCrv;

    if (NSrf != NULL)
	SinglePath = FALSE;

    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	    Srf = CnvrtBezier2BsplineSrf(Srf);
	    SrfBezier = TRUE;
	    break;
	case CAGD_SBSPLINE_TYPE:
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_WRONG_SRF);
	    break;
    }

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    Crv1Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> UKnotVector[0] + EPSILON : EPSILON;
	    Crv2Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> UKnotVector[Srf -> ULength +
				   Srf -> UOrder - 1] - EPSILON
		: 1.0 - EPSILON;
	    break;
	case CAGD_CONST_V_DIR:
	    Crv1Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> VKnotVector[0] + EPSILON : EPSILON;
	    Crv2Param = Srf -> GType == CAGD_SBSPLINE_TYPE ? 
		Srf -> VKnotVector[Srf -> VLength +
				   Srf -> VOrder - 1] - EPSILON
		: 1.0 - EPSILON;
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    Crv1 = CagdCrvFromSrf(Srf, Crv1Param, Dir);
    Crv2 = CagdCrvFromSrf(Srf, Crv2Param, Dir);
    if (NSrf != NULL) {
	NCrv1 = CagdCrvFromSrf(NSrf, Crv1Param, Dir);
	NCrv2 = CagdCrvFromSrf(NSrf, Crv2Param, Dir);
    }
    else
	NCrv1 = NCrv2 = NULL;

    /* Compute the adaptive iso curves. */
    AllAdapIso = CagdAdapIsoExtractAux(0, Srf, NSrf,
				       Crv1, NCrv1, Crv2, NCrv2,
				       Crv1Param, Crv2Param,
				       Dir, Eps * Eps, FullIso, SinglePath);

    /* Chain first and last iso curves that always span the entire domain. */
    if (AllAdapIso != NULL) {
	Crv1 -> Pnext = AllAdapIso;
	for (TCrv = AllAdapIso; TCrv -> Pnext != NULL; TCrv = TCrv -> Pnext);
	TCrv -> Pnext = Crv2;
    }
    else
	Crv1 -> Pnext = Crv2;

    if (NSrf != NULL) {
	NCrv1 -> Pnext = Crv1 -> Pnext;
	Crv1 -> Pnext = NCrv1;
	NCrv2 -> Pnext = Crv2 -> Pnext;
	Crv2 -> Pnext = NCrv2;
    }

    if (SrfBezier)
	CagdSrfFree(Srf);

    return Crv1;
}

/******************************************************************************
* An auxiliary function of CagdAdapIsoExtract. Computes the distance square   *
* between the given two curves, extract the regions that are far than that    *
* and recursively invoke this function with them.			      *
******************************************************************************/
static CagdCrvStruct *CagdAdapIsoExtractAux(int Level,
				     CagdSrfStruct *Srf, CagdSrfStruct *NSrf,
				     CagdCrvStruct *Crv1, CagdCrvStruct *NCrv1,
				     CagdCrvStruct *Crv2, CagdCrvStruct *NCrv2,
				     CagdRType Crv1Param, CagdRType Crv2Param,
				     CagdSrfDirType Dir, CagdRType Eps2,
				     CagdBType FullIso, CagdBType SinglePath)
{
    int i, KVLen;
    CagdCrvStruct
	*DiffCrv = CagdCrvSub(Crv1, Crv2),
	*DistCrv2 = CagdCrvDotProd(DiffCrv, DiffCrv);  /* Dist.^2 function. */
    CagdPointType
	PType = DistCrv2 -> PType;
    CagdRType Dist, *KV, LastT, DistCrv2Min, DistCrv2Max,
	**Points, *PtsW, *PtsX, LastDist,
	*Nodes = NULL;
    CagdBType CloseEnough, LastCloseEnough;
    CagdCrvStruct *TCrv,
	*AllAdapIsoTail = NULL,
	*AllAdapIso = NULL;

    CagdCrvFree(DiffCrv);

    /* Simple heuristic how much to refine DistCrv2: */
    KVLen = DistCrv2 -> Length * 2;
    DistCrv2Min = DistCrv2 -> KnotVector[0];
    DistCrv2Max = DistCrv2 -> KnotVector[DistCrv2 -> Order +
					 DistCrv2 -> Length - 1];
    KV = (CagdRType *) IritMalloc(sizeof(CagdRType) * KVLen);
    for (i = 0; i < KVLen; i++)
	KV[i] = DistCrv2Min + (i + 1) * (DistCrv2Max - DistCrv2Min) /
								(KVLen + 1);

    TCrv = BspCrvKnotInsertNDiff(DistCrv2, FALSE, KV, KVLen);
    IritFree((VoidPtr) KV);
    CagdCrvFree(DistCrv2);
    DistCrv2 = TCrv;

    Nodes = BspKnotNodes(DistCrv2 -> KnotVector,
			 DistCrv2 -> Length + DistCrv2 -> Order,
			 DistCrv2 -> Order);
    LastT = Nodes[0];
    Points = DistCrv2 -> Points,
    PtsW = Points[W],
    PtsX = Points[X],
    LastDist = (PType == CAGD_PT_E1_TYPE ? PtsX[0]
					 : (PtsX[0] / PtsW[0])) - Eps2;
    LastCloseEnough = LastDist < 0.0 && MinSubdivLevel <= Level;

    for (i = 1; i < DistCrv2 -> Length; i++) {
	Dist = (PType == CAGD_PT_E1_TYPE ? PtsX[i]
					 : (PtsX[i] / PtsW[i])) - Eps2;
	CloseEnough = Dist < 0.0 && MinSubdivLevel <= Level;

	if (CloseEnough != LastCloseEnough ||
	    (i == DistCrv2 -> Length - 1 && !LastCloseEnough)) {
	    CagdRType t;

	    if (CloseEnough == LastCloseEnough)
		t = Nodes[DistCrv2 -> Length - 1];
	    else
		CONVEX_BLEND(Nodes[i - 1], Nodes[i], LastDist, Dist, t);

	    if (CloseEnough ||
		(i == DistCrv2 -> Length - 1 && !LastCloseEnough)) {
		/* We are at the end of a region that is not close enough. */
		if (SinglePath) {
		    CagdRType
		        MidParam1 = (Crv1Param * 2.0 + Crv2Param) / 3.0,
			MidParam2 = (Crv1Param + Crv2Param * 2.0) / 3.0;
		    CagdCrvStruct
		        *AdapIso1, *AdapIso2, *AdapIso3, *AdapIso,
			*MidCrv1 = CagdCrvFromSrf(Srf, MidParam1, Dir),
			*MidCrv2 = CagdCrvFromSrf(Srf, MidParam2, Dir);

		    if (FullIso) {
			AdapIso1 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 Crv1, NULL,
							 MidCrv1, NULL,
							 Crv1Param, MidParam1,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 MidCrv1, NULL,
							 MidCrv2, NULL,
							 MidParam1, MidParam2,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso3 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 MidCrv2, NULL,
							 Crv2, NULL,
							 MidParam2, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidCrv1;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidCrv1; 
			if (AdapIso2 != NULL)
			    MidCrv1 -> Pnext = AdapIso2;

			if (AdapIso2 != NULL) {
			    for (TCrv = AdapIso2;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidCrv2;
			}
			else
			    MidCrv1 -> Pnext = MidCrv2;
			if (AdapIso2 != NULL)
			    MidCrv2 -> Pnext = AdapIso3;

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidCrv2;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);

			break; /* Only one (entire domain) region! */
		    }
		    else {
			CagdCrvStruct *Region1, *Region2,
			    *MidRegion1, *MidRegion2;

			Region1 = CagdCrvRegionFromCrv(Crv1, LastT, t),
			Region2 = CagdCrvRegionFromCrv(Crv2, LastT, t);
			MidRegion1 = CagdCrvRegionFromCrv(MidCrv1, LastT, t);
			MidRegion2 = CagdCrvRegionFromCrv(MidCrv2, LastT, t);
			CagdCrvFree(MidCrv1);
			CagdCrvFree(MidCrv2);

			AdapIso1 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 Region1, NULL,
							 MidRegion1, NULL,
							 Crv1Param, MidParam1,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 MidRegion1, NULL,
							 MidRegion2, NULL,
							 MidParam1, MidParam2,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso3 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 MidRegion2, NULL,
							 Region2, NULL,
							 MidParam2, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidRegion1;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidRegion1; 
			if (AdapIso2 != NULL)
			    MidRegion1 -> Pnext = AdapIso2;

			if (AdapIso2 != NULL) {
			    for (TCrv = AdapIso2;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidRegion2;
			}
			else
			    MidRegion1 -> Pnext = MidRegion2;
			if (AdapIso2 != NULL)
			    MidRegion2 -> Pnext = AdapIso3;

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidRegion2;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);
		    }
		}
		else { /* Return all the isocurves as a list of curves. */
		    CagdRType
		        MidParam = (Crv1Param + Crv2Param) / 2.0;
		    CagdCrvStruct *AdapIso1, *AdapIso2, *AdapIso,
			*MidCrv = CagdCrvFromSrf(Srf, MidParam, Dir),
			*MidNCrv = NSrf ? CagdCrvFromSrf(NSrf, MidParam, Dir)
					: NULL;

		    if (FullIso) {
			AdapIso1 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 Crv1, NCrv1,
							 MidCrv, MidNCrv,
							 Crv1Param, MidParam,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 MidCrv, MidNCrv,
							 Crv2, NCrv2,
							 MidParam, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidCrv;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidCrv;

			if (NSrf != NULL) {
			    MidCrv -> Pnext = MidNCrv;
			    MidCrv = MidNCrv;
			}

			if (AdapIso2 != NULL)
			    MidCrv -> Pnext = AdapIso2;

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidCrv;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);

			break; /* Only one (entire domain) region! */
		    }
		    else {
			CagdCrvStruct *Region1, *Region2, *NRegion1, *NRegion2,
			    *MidRegion, *MidNRegion;

			Region1 = CagdCrvRegionFromCrv(Crv1, LastT, t);
			Region2 = CagdCrvRegionFromCrv(Crv2, LastT, t);
			MidRegion = CagdCrvRegionFromCrv(MidCrv, LastT, t);
			if (NSrf) {
			    NRegion1 = CagdCrvRegionFromCrv(NCrv1, LastT, t);
			    NRegion2 = CagdCrvRegionFromCrv(NCrv2, LastT, t);
			    MidNRegion = CagdCrvRegionFromCrv(MidNCrv, LastT,
							      t);
			    CagdCrvFree(MidNCrv);
			}
			CagdCrvFree(MidCrv);

			AdapIso1 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 Region1, NRegion1,
							 MidRegion, MidNRegion,
							 Crv1Param, MidParam,
							 Dir, Eps2,
							 FullIso, SinglePath);
			AdapIso2 = CagdAdapIsoExtractAux(Level + 1, Srf, NSrf,
							 MidRegion, MidNRegion,
							 Region2, NRegion2,
							 MidParam, Crv2Param,
							 Dir, Eps2,
							 FullIso, SinglePath);

			if (AdapIso1 != NULL) {
			    for (TCrv = AdapIso1;
				 TCrv -> Pnext != NULL;
				 TCrv = TCrv -> Pnext);
			    TCrv -> Pnext = MidRegion;
			    AdapIso = AdapIso1;
			}
			else
			    AdapIso = MidRegion;

			if (NSrf != NULL) {
			    MidRegion -> Pnext = MidNRegion;
			    MidRegion = MidNRegion;
			}

			if (AdapIso2 != NULL)
			    MidRegion -> Pnext = AdapIso2;

			CagdCrvFree(Region1);
			CagdCrvFree(Region2);
			if (NSrf != NULL) {
			    CagdCrvFree(NRegion1);
			    CagdCrvFree(NRegion2);
			}

			/* Chain these isolines to the entire set.*/
			if (AllAdapIso)
			    AllAdapIsoTail -> Pnext = AdapIso;
			else
			    AllAdapIso = AdapIso;

			for (AllAdapIsoTail = MidRegion;
			     AllAdapIsoTail -> Pnext != NULL;
			     AllAdapIsoTail = AllAdapIsoTail -> Pnext);
		    }
		}
	    }
	    else {
		/* We are at the beginning of a region not close enough. */
		LastT = t;
	    }
	}

	LastDist = Dist;
	LastCloseEnough = CloseEnough;
    }

    CagdCrvFree(DistCrv2);
    IritFree((VoidPtr) Nodes);

    if (SinglePath) { /* Add connecting isocurves into the existing curves. */
	/* Not implemented yet. */
    }

    return AllAdapIso;
}