/******************************************************************************
* NrmlCone.c - normal cone bounds for normal fields.			      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, October 95.					      *
******************************************************************************/

#include "symb_loc.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes a tangent cone for a given curve, by examine the control polygon  M
* of the curve and deriving its angular span.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:    To compute a tangent cone for.                                   M
*   Planar: If TRUE, only the X and Y coefficients are considered.           M
*                                                                            *
* RETURN VALUE:                                                              M
*   SymbNormalConeStruct *:   The computed tangent cone.                     M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbNormalConeOverlap, SymbNormalConeForSrf                              M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbTangentConeForCrv, tangents, tangent bound                           M
*****************************************************************************/
SymbNormalConeStruct *SymbTangentConeForCrv(CagdCrvStruct *Crv, int Planar)
{
    int i, j, Length;
    CagdRType **Points, t, ConeAngle;
    CagdVType ConeAxis, Tan;
    SymbNormalConeStruct
	*TangentCone = (SymbNormalConeStruct *)
			 	    IritMalloc(sizeof(SymbNormalConeStruct));

    Crv = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE);

    Points = Crv -> Points;
    Length = Crv -> Length;
    
    PT_RESET(ConeAxis);

    /* Make sure coefficients of nrmlCrv are all unit length normals.       */
    /* Also compute the average vector at the same time.		    */
    for (i = 1; i < Length; i++) {
	if (Planar) {
	    for (j = 0; j < 2; j++)
	        Tan[j] = Points[j + 1][i] - Points[j + 1][i - 1];
	    Tan[2] = 0.0;
	}
	else {
	    for (j = 0; j < 3; j++)
	        Tan[j] = Points[j + 1][i] - Points[j + 1][i - 1];
	}
	if ((t = PT_LENGTH(Tan)) > IRIT_EPS) {
	    t = 1.0 / t;
	    PT_SCALE(Tan, t);
	}

	PT_ADD(ConeAxis, ConeAxis, Tan);
    }
    if ((t = PT_LENGTH(ConeAxis)) < IRIT_UEPS) {
        ConeAxis[0] = ConeAxis[1] = 1.0;
	ConeAxis[2] = 0.0;
    }
    else {
        t = 1.0 / t;
        PT_SCALE(ConeAxis, t);
    }

    /* Find the maximal angle between ConeAxis and the vector in mesh. */
    ConeAngle = 1.0;
    for (i = 1; i < Length; i++) {
        CagdRType InnerProd;

	if (Planar) {
	    for (j = 0; j < 2; j++)
	        Tan[j] = Points[j + 1][i] - Points[j + 1][i - 1];
	    Tan[2] = 0.0;
	}
	else {
	    for (j = 0; j < 3; j++)
	        Tan[j] = Points[j + 1][i] - Points[j + 1][i - 1];
	}
	if ((t = PT_LENGTH(Tan)) > IRIT_UEPS) {
	    t = 1.0 / t;
	    PT_SCALE(Tan, t);

	    InnerProd = DOT_PROD(ConeAxis, Tan);

	    if (ConeAngle > InnerProd)
	        ConeAngle = InnerProd;
	}
    }
    	
    PT_COPY(TangentCone -> ConeAxis, ConeAxis);
    TangentCone -> ConeAngle = acos(ConeAngle);
    return TangentCone;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes a normal cone for a given surface, by computing the normal field  M
* of the surface and deriving the angular span of this normal field by       M
* testing the angular span of all control vector in the normal field.        M
*   A normal field is searched for as "_NormalSrf" attribute in Srf or       M
* computed locally of no such attribute is found.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:    To compute a normal cone for.                                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   SymbNormalConeStruct *:   The computed normal cone.                      M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbNormalConeOverlap, SymbTangentConeForCrv                             M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbNormalConeForSrf, normals, normal bound                              M
*****************************************************************************/
SymbNormalConeStruct *SymbNormalConeForSrf(CagdSrfStruct *Srf)
{
    CagdBType LocalNrmlSrf;
    int i, MeshSize;
    CagdSrfStruct
	*NrmlSrf = (CagdSrfStruct *) AttrGetPtrAttrib(Srf -> Attr,
						      "_NormalSrf");
    CagdRType **Points, *XPts, *YPts, *ZPts, ConeAngle;
    CagdVType ConeAxis;
    SymbNormalConeStruct
	*NormalCone = (SymbNormalConeStruct *)
			 	    IritMalloc(sizeof(SymbNormalConeStruct));

    if (NrmlSrf == NULL) {
	NrmlSrf = SymbSrfNormalSrf(Srf);
	LocalNrmlSrf = TRUE;
    }
    else
    	LocalNrmlSrf = FALSE;
    if (NrmlSrf -> PType != CAGD_PT_E3_TYPE) {
	CagdSrfStruct
	    *TSrf = CagdCoerceSrfTo(NrmlSrf, CAGD_PT_E3_TYPE);

	if (LocalNrmlSrf)
	    CagdSrfFree(NrmlSrf);
	NrmlSrf = TSrf;
	LocalNrmlSrf = TRUE;
    }

    Points = NrmlSrf -> Points;
    XPts = Points[1];
    YPts = Points[2];
    ZPts = Points[3];

    MeshSize = NrmlSrf -> ULength * NrmlSrf -> VLength;
    
    PT_RESET(ConeAxis);

    /* Make sure coefficients of nrmlSrf are all unit length normals.       */
    /* Also compute the average vector at the same time.		    */
    for (i = 0; i < MeshSize; i++) {
	CagdRType
	    Len = sqrt(SQR(XPts[i]) + SQR(YPts[i]) + SQR(ZPts[i]));

	if (Len != 0.0) {
	    XPts[i] /= Len;
	    YPts[i] /= Len;
	    ZPts[i] /= Len;
	}

	ConeAxis[0] += XPts[i];
	ConeAxis[1] += YPts[i];
	ConeAxis[2] += ZPts[i];
    }
    PT_SCALE(ConeAxis, (1.0 / MeshSize));

    /* Find the maximal angle between ConeAxis and the vector in mesh. */
    ConeAngle = 1.0;
    for (i = 0; i < MeshSize; i++) {
    	CagdRType
    	    InnerProd = ConeAxis[0] * XPts[i] +
			ConeAxis[1] * YPts[i] +
			ConeAxis[2] * ZPts[i];

	if (ConeAngle > InnerProd)
	    ConeAngle = InnerProd;
    }

    if (LocalNrmlSrf)
	CagdSrfFree(NrmlSrf);
    	
    PT_COPY(NormalCone -> ConeAxis, ConeAxis);
    NormalCone -> ConeAngle = acos(ConeAngle);
    return NormalCone;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Tests if the given two normal cones overlap or not.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   NormalCone1, NormalCone2:  The two normal cones to test for angular      M
		overlap.						     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdBType:   TRUE if overlap, FALSE otherwise.                           M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbNormalConeOverlap                                                    M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbNormalConeOverlap, normals, normal bound                             M
*****************************************************************************/
CagdBType SymbNormalConeOverlap(SymbNormalConeStruct *NormalCone1,
				SymbNormalConeStruct *NormalCone2)
{
    CagdRType
    	Angle = acos(DOT_PROD(NormalCone1 -> ConeAxis,
			      NormalCone2 -> ConeAxis));

    return Angle < NormalCone1 -> ConeAngle + NormalCone2 -> ConeAngle;
}