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

#include "symb_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 SymbScalarCrvConstSet(CagdCrvStruct *Crv, CagdRType ConstVal);
static void SymbScalarCrvExtremSet(CagdCrvStruct *Crv);
static void SymbInsertNewParam(CagdRType t);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes the zero set of a given curve, in the given axis (1-3 for X-Z).   M
*   Returned is a list of the zero set points holding the parameter values   M
* att Pt[0] of each point.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To compute its zero set.                                     M
*   Axis:       The axis of Crv to compute zero set for.                     M
*   Epsilon:    Tolerance control.                                           M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPtStruct *:   List of parameter values form which Crv is zero in     M
*                     axis Axis.                                             M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvZeroSet, zero set, symbolic computation                           M
*****************************************************************************/
CagdPtStruct *SymbCrvZeroSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon)
{
    return SymbCrvConstSet(Crv, Axis, Epsilon, 0.0);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes the extremum set of a given curve, in given axis (1-3 for X-Z).   M
*   Returned is a list of the extreme set points holding the parameter       M
* values at Pt[0] of each point.					     M
*   One could compute the derivative of the curve and find its zero set.     M
*   However, for rational curves, this will double the degree and slow down  M
* the computation considerably.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To compute its extremum set.                                 M
*   Axis:       The axis of Crv to compute extremum set for.                 M
*   Epsilon:    Tolerance control.                                           M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPtStruct *:   List of parameter values form which Crv has an         M
*                     extremum value in axis Axis.                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvExtremSet, extremum set, symbolic computation                     M
*****************************************************************************/
CagdPtStruct *SymbCrvExtremSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon)
{
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *TCrv,
	*ScalarCrv = NULL;

    GlblSetEpsilon = Epsilon;

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

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

    Crv = SymbCrvMergeScalar(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);
    SymbScalarCrvExtremSet(Crv);
    CagdCrvFree(Crv);

    return GlblPtList;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes the constant set of a given curve, in the given axis 	     M
* (1-3 for X-Z).							     M
*   Returned is a list of the constant set points holding the parameter	     M
* values at Pt[0] of each point.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To compute its constant set.                                 M
*   Axis:       The axis of Crv to compute constant set for.                 M
*   Epsilon:    Tolerance control.                                           M
*   ConstVal:   The value at which to compute the constant set.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPtStruct *:   List of parameter values form which Crv has an         M
*                     value of ConstVal in axis Axis.                        M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvConstSet, constant set, zero set, symbolic computation            M
*****************************************************************************/
CagdPtStruct *SymbCrvConstSet(CagdCrvStruct *Crv,
			      int Axis,
			      CagdRType Epsilon,
			      CagdRType ConstVal)
{
    CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *TCrv,
	*ScalarCrv = NULL;

    GlblSetEpsilon = Epsilon;

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

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

    Crv = SymbCrvMergeScalar(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);
    SymbScalarCrvConstSet(Crv, ConstVal);
    CagdCrvFree(Crv);

    return GlblPtList;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* 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.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv:        To compute its constant set.                                 *
*   ConstVal:   Thep value at which to compute the constant set.	     *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SymbScalarCrvConstSet(CagdCrvStruct *Crv, CagdRType ConstVal)
{
    int i,
	Len = Crv -> Length,
        Sign = 0;
    CagdRType
	*ScalarPts = Crv -> Points[1],
	*ScalarWPts = Crv -> Points[0];

    if (SymbCrvPosNegWeights(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. */
	    SymbInsertNewParam((TMax + TMin) / 2.0);
	}
	else { 						 /* Needs to subdiv. */
	    CagdCrvStruct
	        *Crv1 = CagdCrvSubdivAtParam(Crv, TMid),
		*Crv2 = Crv1 -> Pnext;

	    SymbScalarCrvConstSet(Crv1, ConstVal);
	    SymbScalarCrvConstSet(Crv2, ConstVal);

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

/*****************************************************************************
* DESCRIPTION:                                                               *
* 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.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv:        To compute its extremum set.                                 *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SymbScalarCrvExtremSet(CagdCrvStruct *Crv)
{
    int i,
	Len = Crv -> Length,
        Sign = 0;
    CagdRType
	*ScalarPts = Crv -> Points[1],
	*ScalarWPts = Crv -> Points[0];

    if (SymbCrvPosNegWeights(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. */
	    SymbInsertNewParam((TMax + TMin) / 2.0);
	}
	else { 						 /* Needs to subdiv. */
	    CagdCrvStruct
	        *Crv1 = CagdCrvSubdivAtParam(Crv, TMid),
		*Crv2 = Crv1 -> Pnext;

	    SymbScalarCrvExtremSet(Crv1);
	    SymbScalarCrvExtremSet(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))
	    SymbInsertNewParam(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))
	    SymbInsertNewParam(TMax);
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns TRUE iff the Crv is not rational or rational with weights that are M
* entirely positive or entirely negative.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       To examine for same sign weights, if any.                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:   TRUE if no weights or all of same sign.                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvPosNegWeights, symbolic computation                               M
*****************************************************************************/
CagdBType SymbCrvPosNegWeights(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;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Insert a single t value into existing GlblPtList, provided no equal t      *
* value exists already in the list. List is ordered incrementally.	     *
*                                                                            *
* PARAMETERS:                                                                *
*   t:         New value to insert to global GlblPtList list.		     *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SymbInsertNewParam(CagdRType t)
{
    CagdPtStruct *PtTmp, *PtLast, *Pt;

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

    Pt = CagdPtNew();
    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;
}
