/******************************************************************************
* 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 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                                                    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);

	CagdSrfFree(NrmlSrf);
	NrmlSrf = TSrf;
    }

    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;
}