/*****************************************************************************
*   A multi resolution curve editor.					     *
*									     *
* Written by:  Gershon Elber				Ver 0.1, June 1993.  *
*****************************************************************************/

#include "irit_sm.h"
#include "allocate.h"
#include "attribut.h"
#include "iritprsr.h"
#include "cagd_lib.h"
#include "symb_lib.h"
#include "user_lib.h"
#include "geomat3d.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes a decomposition of surface Srf into regions that are visible    M
* with normals deviating in a cone of size as set via ConeAngle.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:         To decompose into visibility cones.                         M
*   Resolution:  Of polygonal approximation of surface Srf.                  M
*   ConeAngle:   Of coverage over the unit sphere.  ConeAngle prescribes the M
*		 opening angle of the cone. In Radians.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:  A list of trimmed surface regions of surface Srf that M
*		       are visible from a selected direction (available as   M
*		       "VIewDir" attribute on them) ad that the union of the M
*		       list covers the entire surface Srf.		     M
*                                                                            *
* SEE ALSO:                                                                  M
*   UserViewingConeSrfDomains, UserVisibilityClassify                        M
*                                                                            *
* KEYWORDS:                                                                  M
*   UserSrfVisibConeDecomp                                                   M
*****************************************************************************/
IPObjectStruct *UserSrfVisibConeDecomp(CagdSrfStruct *Srf,
				       CagdRType Resolution,
				       CagdRType ConeAngle)
{
    CagdRType
	ConeSize = sin(ConeAngle);
    CagdSrfStruct
        *NSrf = SymbSrfNormalSrf(Srf);
    IPObjectStruct *PTmp,
	*RetObjs = NULL,
        *SpherePtsDist = CGPointCoverOfUnitHemiSphere(ConeSize),
        *SrfDomains = UserViewingConeSrfDomains(Srf,
						NSrf,
						SpherePtsDist -> U.Pl,
						Resolution,
						ConeAngle,
						FALSE);

    for (PTmp = SrfDomains; PTmp != NULL; PTmp = PTmp -> Pnext) {
	TrimSrfStruct
	    *TSrfList = TrimSrfsFromContours(Srf, PTmp -> U.Pl);
	IPObjectStruct *PObj,
	    *ViewDirObj = AttrGetObjectObjAttrib(PTmp, "ViewDir");

	/* Eliminate the trimmed surfaces in PObj that are outside the     */
	/* viewing cone, in place.					   */
	TSrfList = UserVisibilityClassify(
				    AttrGetObjectObjAttrib(PTmp, "SclrSrf"),
				    TSrfList);
	if (TSrfList != NULL) {
	    PObj = GenTRIMSRFObject(TSrfList);
	    if (ViewDirObj != NULL)
	        AttrSetObjectObjAttrib(PObj, "ViewDir", ViewDirObj, TRUE);
	    LIST_PUSH(PObj, RetObjs);
	}
    }

    CagdSrfFree(NSrf);
    IPFreeObjectList(SrfDomains);
    IPFreeObject(SpherePtsDist);

    return RetObjs;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Given a decomposition of surface Srf into a set of TrimmedSrfs into      M
* regions that are inside the normal viewing cone as prescribed by SclrSrf,  M
* eliminate all regions (trimmed surfaces) that are outside the viewing      M
*  cone in the given list, TrimmedSrfs.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   SclrSrf:	   The scalar surface computed symbollically in		     M
*		   UserViewingConeSrfDomains, forming this decomosition.     M
*   TrimmedSrfs:   To verify they are within the viewing code of ViewingDir. M
*		   Trimmed surfaces inside this list that are outside the    M
*		   viewing cone are eliminated, in place.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   TrimSrfStruct *: The trimmed regions inside the viewing cone, as a       M
*		   subset of TrimmedSrfs.			             M
*                                                                            *
* SEE ALSO:                                                                  M
*   UserViewingConeSrfDomains, UserSrfVisibConeDecomp                        M
*                                                                            *
* KEYWORDS:                                                                  M
*   UserVisibilityClassify                                                   M
*****************************************************************************/
TrimSrfStruct *UserVisibilityClassify(IPObjectStruct *SclrSrf,
				      TrimSrfStruct *TrimmedSrfs)
{
    TrimSrfStruct *TSrf,
	*PrevTSrf = NULL;

    if (SclrSrf == NULL || !IP_IS_SRF_OBJ(SclrSrf)) {
	UserFatalError(USER_ERR_MISSING_ATTRIB);
	return NULL;
    }

    for (TSrf = TrimmedSrfs; TSrf != NULL; ) {
	CagdRType
	    *UVInside = TrimPointInsideTrimmedCrvs(TSrf -> TrimCrvList, TSrf),
	    *R = CagdSrfEval(SclrSrf -> U.Srfs, UVInside[0], UVInside[1]);

	if (R[1] > 0) {
	    if (PrevTSrf != NULL) {
		PrevTSrf -> Pnext = TSrf -> Pnext;
		TrimSrfFree(TSrf);
		TSrf = PrevTSrf -> Pnext;
	    }
	    else {
		TrimmedSrfs = TrimmedSrfs -> Pnext;
		TrimSrfFree(TSrf);
		TSrf = TrimmedSrfs;
	    }
	}
	else {
	    PrevTSrf = TSrf;
	    TSrf = TSrf ->Pnext;
	}
    }

    return TrimmedSrfs;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes the domain in given surface Srf that is inside the given cones  M
* of directions ConeDirs and opening angle ConeAngle.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:           To compute its visibility.                                M
*   NSrf:	   The normal surface computed symbollically for Srf.        M
*   ConeDirs:      Direction of cone's axes.                                 M
*   Resolution:    Of Cone - normal surface intersection, for polygonal      M
*		   approximation accuracy.				     M
*   ConeAngle:     The opening angle of the cone, in Radians.		     M
*   Euclidean:     Contours are in Euclidean (TRUE) or parametric (FALSE)    M
*		   space of Srf.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPPolygonStruct *:   Set of piecewise linear contours, approximating the M
*		   intersections of the normal viewing cones and the         M
*		   original surface.					     M
*                                                                            *
* SEE ALSO:                                                                  M
*   UserVisibilityClassify, UserSrfVisibConeDecomp                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   UserViewingConeSrfDomains                                                M
*****************************************************************************/
IPObjectStruct *UserViewingConeSrfDomains(CagdSrfStruct *Srf,
					  CagdSrfStruct *NSrf,
					  IPPolygonStruct *ConeDirs,
					  CagdRType Resolution,
					  CagdRType ConeAngle,
					  CagdRType Euclidean)
{
    int OldCirc;
    IPVertexStruct
        *VPts = ConeDirs -> PVertex;
    IPObjectStruct
	*CntrObjList = NULL;

    if (ConeAngle < 0.0 || ConeAngle >= M_PI / (2 + IRIT_EPS)) {
	UserFatalError(USER_ERR_WRONG_ANGLE);
	return NULL;
    }

    for ( ; VPts != NULL; VPts = VPts -> Pnext) {
	static PlaneType
	    Plane = { 1.0, 0.0, 0.0, 0.0 };
	int i;
	CagdBType HasPosVals, HasNegVals;
	CagdRType *R;
	MatrixType Mat, InvMat;
	CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *TmpSrf, *SclrSrf;

	CGGenMatrixZ2Dir(Mat, VPts -> Coord);
	MatInverseMatrix(Mat, InvMat);
	TmpSrf = CagdSrfCopy(NSrf);
	CagdSrfMatTransform(TmpSrf, InvMat);

	SymbSrfSplitScalar(TmpSrf, &SrfW, &SrfX, &SrfY, &SrfZ);
	CagdSrfFree(TmpSrf);
	if (SrfW != NULL) {
	    fprintf(stderr, "No support for rationals, at this time.\n");
	    exit(1);
	}

	TmpSrf = SymbSrfMult(SrfX, SrfX);
	CagdSrfFree(SrfX);
	SrfX = TmpSrf;

	TmpSrf = SymbSrfMult(SrfY, SrfY);
	CagdSrfFree(SrfY);
	SrfY = TmpSrf;

	TmpSrf = SymbSrfMult(SrfZ, SrfZ);
	CagdSrfFree(SrfZ);
	SrfZ = SymbSrfScalarScale(TmpSrf, SQR(tan(ConeAngle)));
	CagdSrfFree(TmpSrf);

	TmpSrf = SymbSrfAdd(SrfX, SrfY);
	CagdSrfFree(SrfX);
	CagdSrfFree(SrfY);

	SclrSrf = SymbSrfSub(TmpSrf, SrfZ);
	CagdSrfFree(TmpSrf);
	CagdSrfFree(SrfZ);

	R = SclrSrf -> Points[1];
	HasPosVals = HasNegVals = FALSE;
	for (i = SclrSrf -> ULength * SclrSrf -> VLength; i >= 0; i--) {
	    if (*R > 0.0)
	        HasPosVals = TRUE;
	    if (*R++ < 0.0)
	        HasNegVals = TRUE;
	}
	if (HasPosVals && HasNegVals) {
	    IPPolygonStruct *Cntrs, *Cntr;

	    OldCirc = IritPrsrSetPolyListCirc(TRUE);
	    Cntrs = UserCntrSrfWithPlane(SclrSrf, Plane, Resolution);
	    IritPrsrSetPolyListCirc(OldCirc);

	    if (Cntrs != NULL) {
		IPObjectStruct *CntrObj, *VecObj, *SclrSrfObj;
		IPVertexStruct *V;

		if (Euclidean) {
		    for (Cntr = Cntrs; Cntr != NULL; Cntr = Cntr -> Pnext) {
			for (V = Cntr -> PVertex; V != NULL; V = V -> Pnext) {
			    CagdRType
			        *P = CagdSrfEval(Srf, V -> Coord[1],
						      V -> Coord[2]);

			    CagdCoerceToE3(V -> Coord, &P, -1, Srf -> PType);
			}
		    }
		}
		else {
		    for (Cntr = Cntrs; Cntr != NULL; Cntr = Cntr -> Pnext) {
			for (V = Cntr -> PVertex; V != NULL; V = V -> Pnext) {
			    V -> Coord[0] = V -> Coord[1];
			    V -> Coord[1] = V -> Coord[2];
			    V -> Coord[2] = 0.0;
			}
		    }
		}

		CntrObj = GenPOLYObject(Cntrs);
		IP_SET_POLYLINE_OBJ(CntrObj);

		/* Save the viewing direction as an attribute. */
		VecObj = GenVECObject(&VPts -> Coord[0],
				      &VPts -> Coord[1],
				      &VPts -> Coord[2]);
		AttrSetObjectObjAttrib(CntrObj, "ViewDir", VecObj, FALSE);

		SclrSrfObj = GenSRFObject(SclrSrf);
		AttrSetObjectObjAttrib(CntrObj, "SclrSrf", SclrSrfObj, FALSE);

		LIST_PUSH(CntrObj, CntrObjList);
	    }
	    else {
	        CagdSrfFree(SclrSrf);
	    }
	}
	else {
	    CagdSrfFree(SclrSrf);
	}
    }

    return CntrObjList;
}

#ifdef DEBUG_VISIBILITY

void main(int argc, char **argv)
{
    IPObjectStruct *PObj;
    RealType
	Size = 0.6,
	Resolution = 20;

    if (argc > 1 && strcmp(argv[1], "-s") == 0) {
	Size = atof(argv[2]);
	argv += 2;
	argc -= 2;
    }
    if (argc > 1 && strcmp(argv[1], "-r") == 0) {
	Resolution = atof(argv[2]);
	argv += 2;
	argc -= 2;
    }

    fprintf(stderr, "Doing decomposition of size %f, resolution %f\n",
	    Size, Resolution);

    if ((PObj = IritPrsrGetDataFiles(&argv[1], 1, TRUE, TRUE)) != NULL &&
	IP_IS_SRF_OBJ(PObj)) {
	int FileIndex = 1;
	IPObjectStruct *PTmp,
	    *PTrimSrfObjs = UserSrfVisibConeDecomp(PObj -> U.Srfs,
						   Resolution, Size);

	for (PTmp = PTrimSrfObjs; PTmp != NULL; PTmp = PTmp -> Pnext) {
	    int Handler;
	    char FileName[LINE_LEN];

	    sprintf(FileName, "ViewDir%d.dat", FileIndex++);
	    if ((Handler = IritPrsrOpenDataFile(FileName,
						FALSE, FALSE)) >= 0) {
		IritPrsrPutObjectToHandler(Handler, PTmp);
		IritPrsrCloseStream(Handler, TRUE);
	    }
	}
	IPFreeObjectList(PTrimSrfObjs);
    }
    else {
	fprintf(stderr, "Expecting a surface");
    }

    exit(0);
}

#endif /* DEBUG_VISAIBILITY */