/******************************************************************************
* CagdZero.c - computes the zeros and extremes of a given object.	      *
*******************************************************************************
* Written by Gershon Elber, March 93.					      *
******************************************************************************/

#include "cagd_loc.h"

#define SET_DEFAULT_EPSILON	1e-3
#define SET_DEFAULT_EPSILON	1e-3
#define ZERO_APX_EQ(x, y)		(ABS(x - y) < GlblSetEpsilon * 10)
#define MID_RATIO			0.5123456789

static CagdPtStruct
    *GlblPtList = NULL;
static CagdRType CrvTMin, CrvTMax,
    GlblSetEpsilon = SET_DEFAULT_EPSILON;

static void CagdScalarCrvConstSet(CagdCrvStruct *Crv, CagdRType ConstVal);
static void CagdScalarCrvExtremSet(CagdCrvStruct *Crv);
static void CagdInsertNewParam(CagdRType t);

/******************************************************************************
* Computes the zero set of a given curve, in the given axis (1-3 for X-Z).    *
*   Returned is a list of the zero set points holding the parameter values at *
* Pt[0] of each point.							      *
******************************************************************************/
CagdPtStruct *CagdCrvZeroSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon)
{
    return CagdCrvConstSet(Crv, Axis, Epsilon, 0.0);
}

/******************************************************************************
* Computes the extreme set of a given curve, in given axis (1-3 for X-Z).     *
*   Returned is a list of the extreme set points holding the parameter        *
* values at Pt[0] of each point.					      *
*   One could compute the derivative of the curve and find its zero set.      *
* However, for rational curves, this will double the degree and slow down the *
* computation considerably.						      *
******************************************************************************/
CagdPtStruct *CagdCrvExtremSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon)
{
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *ScalarCrv, *TCrv;

    GlblSetEpsilon = Epsilon;

    CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);

    switch (Axis) {
	case 1:
	    if (CrvX)
		ScalarCrv = CrvX;
	    else
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    break;
	case 2:
	    if (CrvY)
		ScalarCrv = CrvY;
	    else
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    break;
	case 3:
	    if (CrvZ)
		ScalarCrv = CrvZ;
	    else
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
    }

    Crv = CagdCrvMergeScalar(CrvW, ScalarCrv, NULL, NULL);
    if (CrvW)
	CagdCrvFree(CrvW);
    if (CrvX)
	CagdCrvFree(CrvX);
    if (CrvY)
	CagdCrvFree(CrvY);
    if (CrvZ)
	CagdCrvFree(CrvZ);
    
    GlblPtList = NULL;
    if (CAGD_IS_BEZIER_CRV(Crv)) {
	/* We need to save domains. */
	TCrv = CnvrtBezier2BsplineCrv(Crv);
	CagdCrvFree(Crv);
	Crv = TCrv;
    }

    CagdCrvDomain(Crv, &CrvTMin, &CrvTMax);
    CagdScalarCrvExtremSet(Crv);
    CagdCrvFree(Crv);

    return GlblPtList;
}

/******************************************************************************
* Computes the constant set of a given curve, in the given axis (1-3 for X-Z).*
* Returned is a list of the constant set points holding the parameter values  *
* at Pt[0] of each point.						      *
******************************************************************************/
CagdPtStruct *CagdCrvConstSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon,
							    CagdRType ConstVal)
{
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *ScalarCrv, *TCrv;

    GlblSetEpsilon = Epsilon;

    CagdCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);

    switch (Axis) {
	case 1:
	    if (CrvX)
		ScalarCrv = CrvX;
	    else
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    break;
	case 2:
	    if (CrvY)
		ScalarCrv = CrvY;
	    else
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    break;
	case 3:
	    if (CrvZ)
		ScalarCrv = CrvZ;
	    else
		FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE);
    }

    Crv = CagdCrvMergeScalar(CrvW, ScalarCrv, NULL, NULL);
    if (CrvW)
	CagdCrvFree(CrvW);
    if (CrvX)
	CagdCrvFree(CrvX);
    if (CrvY)
	CagdCrvFree(CrvY);
    if (CrvZ)
	CagdCrvFree(CrvZ);
    
    GlblPtList = NULL;
    if (CAGD_IS_BEZIER_CRV(Crv)) {
	/* We need to save domains. */
	TCrv = CnvrtBezier2BsplineCrv(Crv);
	CagdCrvFree(Crv);
	Crv = TCrv;
    }

    CagdCrvDomain(Crv, &CrvTMin, &CrvTMax);
    CagdScalarCrvConstSet(Crv, ConstVal);
    CagdCrvFree(Crv);

    return GlblPtList;
}

/******************************************************************************
* Computes the zero set of a scalar curve. Curve might be rational.	      *
* Assumes open end condition on the curve parametric domain.		      *
* Zero set points are append to the GlblPtList point list.		      *
******************************************************************************/
static void CagdScalarCrvConstSet(CagdCrvStruct *Crv, CagdRType ConstVal)
{
    int i,
	Len = Crv -> Length,
        Sign = 0;
    CagdRType
	*ScalarPts = Crv -> Points[1],
	*ScalarWPts = Crv -> Points[0];

    if (CagdCrvPosNegWeights(Crv)) {
	i = 0; 				 /* Force subdivision of the curve. */
    }
    else {
	for (i = 0; i < Len; i++) {
	    CagdRType
		V = (ScalarWPts ? ScalarPts[i] / ScalarWPts[i] :
                                  ScalarPts[i]) - ConstVal;
	    int NewSign = SIGN(V);

	    if (Sign * NewSign < 0)
		break;
	    else if (NewSign)
		Sign = NewSign;
	}
    }

    if (i < Len) {
	CagdRType TMin, TMax, TMid;

	/* Control polygon is both positive and negative. */
	CagdCrvDomain(Crv, &TMin, &TMax);
	TMid = TMin * MID_RATIO + TMax * (1.0 - MID_RATIO);

	if (TMax - TMin < GlblSetEpsilon / 5.0) {    /* Small enough - stop. */
	    CagdInsertNewParam((TMax + TMin) / 2.0);
	}
	else { 						 /* Needs to subdiv. */
	    CagdCrvStruct
	        *Crv1 = CagdCrvSubdivAtParam(Crv, TMid),
		*Crv2 = Crv1 -> Pnext;

	    CagdScalarCrvConstSet(Crv1, ConstVal);
	    CagdScalarCrvConstSet(Crv2, ConstVal);

	    CagdCrvFree(Crv1);
	    CagdCrvFree(Crv2);
	}
    }
}

/******************************************************************************
* Computes the extrem set of a scalar curve. Curve might be rational.	      *
* Assumes open end condition on the curve parametric domain.		      *
* Extreme set points are append to the GlblPtList point list.		      *
******************************************************************************/
static void CagdScalarCrvExtremSet(CagdCrvStruct *Crv)
{
    int i,
	Len = Crv -> Length,
        Sign = 0;
    CagdRType
	*ScalarPts = Crv -> Points[1],
	*ScalarWPts = Crv -> Points[0];

    if (CagdCrvPosNegWeights(Crv)) {
	i = 0; 				 /* Force subdivision of the curve. */
    }
    else {
	CagdRType
	    PrevV = ScalarWPts ? ScalarPts[0] / ScalarWPts[0] : ScalarPts[0];

	for (i = 0; i < Len; i++) {
	    CagdRType
		V = (ScalarWPts ? ScalarPts[i] / ScalarWPts[i] : ScalarPts[i]),
		DV = V - PrevV;

	    int NewSign = SIGN(DV);

	    if (Sign * NewSign < 0)
		break;
	    else if (NewSign)
		Sign = NewSign;

	    PrevV = V;
	}
    }

    if (i < Len) {
	CagdRType TMin, TMax, TMid;

	/* Control polygon is both increasing and decreasing. */
	CagdCrvDomain(Crv, &TMin, &TMax);
	TMid = TMin * MID_RATIO + TMax * (1.0 - MID_RATIO);

	if (TMax - TMin < GlblSetEpsilon / 5.0) {    /* Small enough - stop. */
	    CagdInsertNewParam((TMax + TMin) / 2.0);
	}
	else { 						 /* Needs to subdiv. */
	    CagdCrvStruct
	        *Crv1 = CagdCrvSubdivAtParam(Crv, TMid),
		*Crv2 = Crv1 -> Pnext;

	    CagdScalarCrvExtremSet(Crv1);
	    CagdScalarCrvExtremSet(Crv2);

	    CagdCrvFree(Crv1);
	    CagdCrvFree(Crv2);
	}
    }
    else {
	CagdRType V1, V2, TMin, TMax;

	CagdCrvDomain(Crv, &TMin, &TMax);

	/* This segment is monotone. Test if end points are extremes. */
	V1 = ScalarWPts ? ScalarPts[0] / ScalarWPts[0] : ScalarPts[0];
	V2 = ScalarWPts ? ScalarPts[1] / ScalarWPts[1] : ScalarPts[1];
	if (APX_EQ(V1, V2))
	    CagdInsertNewParam(TMin);

	V1 = ScalarWPts ? ScalarPts[Len - 2] / ScalarWPts[Len - 2]
			: ScalarPts[Len - 2];
	V2 = ScalarWPts ? ScalarPts[Len - 1] / ScalarWPts[Len - 1]
			: ScalarPts[Len - 1];
	if (APX_EQ(V1, V2))
	    CagdInsertNewParam(TMax);
    }
}

/******************************************************************************
* Returns TRUE iff the Crv is not rational or rational with weights that are  *
* entirely positive or entirely negative.				      *
******************************************************************************/
CagdBType CagdCrvPosNegWeights(CagdCrvStruct *Crv)
{
    int i;
    CagdBType HasNeg, HasPos;
    CagdRType
	*Weights = Crv -> Points[0];

    if (Weights == NULL)
	return FALSE;				   /* Curve is not rational. */

    for (HasNeg = HasPos = FALSE, i = Crv -> Length - 1; i >= 0; i--) {
	HasNeg |= *Weights < 0.0;
	HasPos |= *Weights++ > 0.0;
    }

    return HasNeg && HasPos;
}

/******************************************************************************
* Insert a single t value into existing GlblPtList, provided no equal t value *
* exists already in the list. List is order incrementally.		      *
******************************************************************************/
static void CagdInsertNewParam(CagdRType t)
{
    CagdPtStruct *PtTmp, *PtLast, *Pt;

    if (APX_EQ(t, CrvTMin) || APX_EQ(t, CrvTMax))
	return;

    Pt = IritMalloc(sizeof(CagdPtStruct));
    Pt -> Pt[0] = t;
    Pt -> Pnext = NULL;

    if (GlblPtList) {
	for (PtTmp = GlblPtList, PtLast = NULL;
	     PtTmp != NULL;
	     PtLast = PtTmp, PtTmp = PtTmp -> Pnext) {
	    if (ZERO_APX_EQ(PtTmp -> Pt[0], t)) {
	        IritFree(Pt);
		return;
	    }
	    if (PtTmp -> Pt[0] > t)
	        break;
	}
	if (PtTmp) {
	    /* Insert the new point in the middle of the list. */
	    Pt -> Pnext = PtTmp;
	    if (PtLast)
		PtLast -> Pnext = Pt;
	    else
		GlblPtList = Pt;
	}
	else {
	    /* Insert the new point as the last point in the list. */
	    PtLast -> Pnext = Pt;
	}
    }
    else
        GlblPtList = Pt;
}