/******************************************************************************
* Prisa.c - piecewise ruled srf approximation and layout (prisa).	      *
*******************************************************************************
* Written by Gershon Elber, Apr. 93.					      *
******************************************************************************/

#include "symb_loc.h"

static CagdSrfStruct *SymbPiecewiseRuledSrfAux(CagdSrfStruct *Srf,
					       CagdBType ConsistentDir,
					       CagdRType Epsilon,
					       CagdSrfDirType Dir);
static void SymbInterCircCirc(CagdRType Center1[3],
			      CagdRType Radius1,
			      CagdRType Center2[3],
			      CagdRType Radius2,
			      CagdRType Inter1[3],
			      CagdRType Inter2[3]);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Computes a piecewise ruled surface approximation to a given set of         M
* surfaces with given Epsilon, and lay them out "nicely" onto the XY plane,  M
* by approximating each ruled surface as a developable surface with          M
* SamplesPerCurve samples.						     M
*   Dir controls the direction of ruled approximation, SpaceScale and        M
* Offset controls the placement of the different planar pieces.		     M
*   Prisa is the hebrew word for the process of flattening out a three       M
* dimensional surface. I have still to find an english word for it.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srfs:            To approximate and faltten out.                         M
*   SamplesPerCurve: During the approximation of a ruled surface as a        M
*		     developable surface.				     M
*   Epsilon:         Accuracy control for the piecewise ruled surface        M
*		     approximation.					     M
*   Dir:             Direction of ruled/developable surface approximation.   M
*		     Either U or V.					     M
*   Space:           A vector in the XY plane to denote the amount of        M
*		     translation from one flattened out surface to the next. M
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A list of planar surfaces denoting the layout (prisa) M
*		       of the given Srfs to the accuracy requested.	     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbAllPrisaSrfs, layout, prisa                                          M
*****************************************************************************/
CagdSrfStruct *SymbAllPrisaSrfs(CagdSrfStruct *Srfs,
				int SamplesPerCurve,
			        CagdRType Epsilon,
				CagdSrfDirType Dir,
			        CagdVType Space)
{
    int SrfIndex = 1;
    CagdSrfStruct *Srf,
	*PrisaSrfsAll = NULL;
    CagdVType Offset;

    for (Srf = Srfs; Srf != NULL; Srf = Srf -> Pnext, SrfIndex++) {
	CagdSrfStruct *TSrf, *RSrf, *RuledSrfs;

	if (Epsilon > 0) {
	    RuledSrfs = SymbPiecewiseRuledSrfApprox(Srf, FALSE, Epsilon, Dir);

	    Offset[0] = SrfIndex * Space[0];
	    Offset[1] = 0.0;
	    Offset[2] = 0.0;

	    for (RSrf = RuledSrfs; RSrf != NULL; RSrf = RSrf -> Pnext) {
		TSrf = SymbPrisaRuledSrf(RSrf, SamplesPerCurve,
					 Space[1], Offset);
		TSrf -> Pnext = PrisaSrfsAll;
		PrisaSrfsAll = TSrf;
	    }

	    CagdSrfFreeList(RuledSrfs);
	}
	else {
	    /* Return the 3D ruled surface approximation. */
	    RuledSrfs = SymbPiecewiseRuledSrfApprox(Srf, FALSE, -Epsilon, Dir);

	    for (TSrf = RuledSrfs;
		 TSrf -> Pnext != NULL;
		 TSrf = TSrf -> Pnext);
	    TSrf -> Pnext = PrisaSrfsAll;
	    PrisaSrfsAll = RuledSrfs;
	}
    }

    return PrisaSrfsAll;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Constructs a piecewise ruled surface approximation to the given surface,   M
* Srf, in the given direction, Dir, that is close to the surface to within   M
* Epsilon.								     M
*   If ConsitentDir then ruled surface parametrization is set to be the      M
* same as original surface Srf. Otherwise, ruling dir is always		     M
* CAGD_CONST_V_DIR.							     M
*   Surface is assumed to have point types E3 or P3 only.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:           To approximate using piecewise ruled surfaces.            M
*   ConsistentDir: Do we want parametrization to be the same as Srf?         M
*   Epsilon:       Accuracy of piecewise ruled surface approximation.        M
*   Dir:           Direction of piecewise ruled surface approximation.       M
*		   Either U or V.					     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   A list of ruled surfaces approximating Srf to within  M
*                      Epsilon in direction Dir.                             M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbPiecewiseRuledSrfApprox, layout, prisa, ruled surface approximation  M
*****************************************************************************/
CagdSrfStruct *SymbPiecewiseRuledSrfApprox(CagdSrfStruct *Srf,
					   CagdBType ConsistentDir,
					   CagdRType Epsilon,
					   CagdSrfDirType Dir)
{
    CagdSrfStruct *RuledSrfs;

    switch (Srf -> PType) {
	case CAGD_PT_E3_TYPE:
	case CAGD_PT_P3_TYPE:
	    Srf = CagdSrfCopy(Srf);
	    break;
	case CAGD_PT_P1_TYPE:
	case CAGD_PT_P2_TYPE:
	case CAGD_PT_P4_TYPE:
	case CAGD_PT_P5_TYPE:
	    Srf = CagdCoerceSrfTo(Srf, CAGD_PT_P3_TYPE);
	    break;
	case CAGD_PT_E1_TYPE:
	case CAGD_PT_E2_TYPE:
	case CAGD_PT_E4_TYPE:
	case CAGD_PT_E5_TYPE:
	    Srf = CagdCoerceSrfTo(Srf, CAGD_PT_E3_TYPE);
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT);
	    Srf = NULL;
	    break;
    }

    RuledSrfs = SymbPiecewiseRuledSrfAux(Srf, ConsistentDir, Epsilon, Dir);
    CagdSrfFree(Srf);
    return RuledSrfs;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Auxiliary function to function SymbPiecewiseRuledSrfApprox		     *
*****************************************************************************/
static CagdSrfStruct *SymbPiecewiseRuledSrfAux(CagdSrfStruct *Srf,
					       CagdBType ConsistentDir,
					       CagdRType Epsilon,
					       CagdSrfDirType Dir)
{
    CagdSrfStruct *DiffSrf, *DistSqrSrf,
	*RuledSrf = CagdSrfCopy(Srf);
    CagdRType *XPts, *WPts, UMin, UMax, VMin, VMax,
	**Points = RuledSrf -> Points,
	MaxError = 0.0,
	t = 0.0;
    int i, j, k, Index,
	ULength = RuledSrf -> ULength,
	VLength = RuledSrf -> VLength;
    PointType E3PtStart, E3PtEnd, E3Pt;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    for (j = 0; j < VLength; j++) {
		/* First order approximation to the ratios of the   */
		/* distance of interior point to the end points.    */
		CagdCoerceToE3(E3PtStart, Points,
			       CAGD_MESH_UV(RuledSrf, ULength / 2, 0),
			       Srf -> PType);
		CagdCoerceToE3(E3PtEnd, Points,
			       CAGD_MESH_UV(RuledSrf, ULength / 2,
					              VLength - 1),
			       Srf -> PType);
		CagdCoerceToE3(E3Pt, Points,
			       CAGD_MESH_UV(RuledSrf, ULength / 2, j),
			       Srf -> PType);
		PT_SUB(E3PtStart, E3PtStart, E3PtEnd);
		PT_SUB(E3Pt, E3Pt, E3PtEnd);
		t = PT_LENGTH(E3PtStart);
		t = APX_EQ(t, 0.0) ? 0.5 : PT_LENGTH(E3Pt) / t;

		for (i = 0; i < ULength; i++) {
		    CagdCoerceToE3(E3PtStart, Points,
				   CAGD_MESH_UV(RuledSrf, i, 0),
				   Srf -> PType);
		    CagdCoerceToE3(E3PtEnd, Points,
				   CAGD_MESH_UV(RuledSrf, i, VLength - 1),
				   Srf -> PType);
		    PT_BLEND(E3Pt, E3PtStart, E3PtEnd, t);

		    Index = CAGD_MESH_UV(RuledSrf, i, j);
		    if (CAGD_IS_RATIONAL_PT(RuledSrf -> PType)) {
			for (k = 0; k < 3; k++)
			    Points[k + 1][Index] = 
				E3Pt[k] * Points[0][Index];
		    }
		    else {
			for (k = 0; k < 3; k++)
			    Points[k + 1][Index] = E3Pt[k];
		    }
		}
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    for (i = 0; i < ULength; i++) {
		/* First order approximation to the ratios of the   */
		/* distance of interior point to the end points.    */
		CagdCoerceToE3(E3PtStart, Points,
			       CAGD_MESH_UV(RuledSrf, 0, VLength / 2),
			       Srf -> PType);
		CagdCoerceToE3(E3PtEnd, Points,
			       CAGD_MESH_UV(RuledSrf, ULength - 1,
						      VLength / 2),
			       Srf -> PType);
		CagdCoerceToE3(E3Pt, Points,
			       CAGD_MESH_UV(RuledSrf, i, VLength / 2),
			       Srf -> PType);
		PT_SUB(E3PtStart, E3PtStart, E3PtEnd);
		PT_SUB(E3Pt, E3Pt, E3PtEnd);
		t = PT_LENGTH(E3PtStart);
		t = APX_EQ(t, 0.0) ? 0.5 : PT_LENGTH(E3Pt) / t;

		for (j = 0; j < VLength; j++) {
		    CagdCoerceToE3(E3PtStart, Points,
				   CAGD_MESH_UV(RuledSrf, 0, j),
				   Srf -> PType);
		    CagdCoerceToE3(E3PtEnd, Points,
				   CAGD_MESH_UV(RuledSrf, ULength - 1, j),
				   Srf -> PType);
		    PT_BLEND(E3Pt, E3PtStart, E3PtEnd, t);

		    Index = CAGD_MESH_UV(RuledSrf, i, j);
		    if (CAGD_IS_RATIONAL_PT(RuledSrf -> PType)) {
			for (k = 0; k < 3; k++)
			    Points[k + 1][Index] = 
				E3Pt[k] * Points[0][Index];
		    }
		    else {
			for (k = 0; k < 3; k++)
			    Points[k + 1][Index] = E3Pt[k];
		    }
		}
	    }
	    break;
	default:
	    SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    DiffSrf = SymbSrfSub(Srf, RuledSrf);
    CagdSrfFree(RuledSrf);
    DistSqrSrf = SymbSrfDotProd(DiffSrf, DiffSrf);
    CagdSrfFree(DiffSrf);
    XPts = DistSqrSrf -> Points[1];
    WPts = CAGD_IS_RATIONAL_PT(DistSqrSrf -> PType) ?
					DistSqrSrf -> Points[0] : NULL;

    for (i = DistSqrSrf -> ULength * DistSqrSrf -> VLength; i > 0; i--) {
	CagdRType
	    V = WPts != NULL ? *XPts++ / *WPts++ : *XPts++;

	if (MaxError < V)
	    MaxError = V;
    }
    CagdSrfFree(DistSqrSrf);

    if (MaxError > SQR(Epsilon)) {
	/* Subdivide and try again. */
	CagdSrfStruct *Srf1, *Srf2, *RuledSrf1, *RuledSrf2;

	CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
	t = Dir == CAGD_CONST_V_DIR ? (UMax + UMin) / 2 : (VMax + VMin) / 2;
	Srf1 = CagdSrfSubdivAtParam(Srf, t, CAGD_OTHER_DIR(Dir));
	Srf2 = Srf1 -> Pnext;
	Srf1 -> Pnext = NULL;

	RuledSrf1 = SymbPiecewiseRuledSrfAux(Srf1, ConsistentDir, Epsilon, Dir);
	RuledSrf2 = SymbPiecewiseRuledSrfAux(Srf2, ConsistentDir, Epsilon, Dir);
	CagdSrfFree(Srf1);
	CagdSrfFree(Srf2);

	for (Srf1 = RuledSrf1; Srf1 -> Pnext != NULL; Srf1 = Srf1 -> Pnext);
	Srf1 -> Pnext = RuledSrf2;
	return RuledSrf1;
    }
    else {
	/* Returns the ruled surface as a linear in the ruled direction. */
	CagdCrvStruct
	    *Crv1 = CagdCrvFromMesh(Srf, 0, CAGD_OTHER_DIR(Dir)),
	    *Crv2 = CagdCrvFromMesh(Srf,
				    Dir == CAGD_CONST_V_DIR ? ULength - 1
							    : VLength - 1,
				    CAGD_OTHER_DIR(Dir));

	RuledSrf = CagdRuledSrf(Crv1, Crv2, 2, 2);
	if (ConsistentDir && Dir == CAGD_CONST_V_DIR) {
	    /* Needs to reverse the ruled surface so it matches Srf. */
	    CagdSrfStruct
		*TSrf = CagdSrfReverse2(RuledSrf);

	    CagdSrfFree(RuledSrf);
	    RuledSrf = TSrf;
	}

	if (CAGD_IS_BSPLINE_SRF(Srf)) {
	    /* Updates the knot vector domain. */
	    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
	    if (Dir == CAGD_CONST_V_DIR)
	        BspKnotAffineTrans(RuledSrf -> UKnotVector,
				   RuledSrf -> ULength + RuledSrf -> UOrder,
				   UMin, UMax - UMin);
	    else
	        BspKnotAffineTrans(RuledSrf -> VKnotVector,
				   RuledSrf -> VLength + RuledSrf -> VOrder,
				   VMin, VMax - VMin);
	}
	CagdCrvFree(Crv1);
	CagdCrvFree(Crv2);
	return RuledSrf;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Layout a single ruled surface, by approximating it as a set of polygons.   M
*   The given ruled surface might be non-developable, in which case	     M
* approximation will be of a surface with no twist.			     M
*   The ruled surface is assumed to be constructed using CagdRuledSrf and    M
* that the ruled direction is consisnt and is always CAGD_CONST_V_DIR.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:             A ruled surface to layout flat on the XY plane.         M
*   SamplesPerCurve: During the approximation of a ruled surface as a        M
*		     developable surface.				     M
*   Space:	     Increment on Y on the offset vector, after this         M
*		     surface was placed in the XY plane.		     M
*   Offset:	     A vector in the XY plane to denote the amount of        M
*		     translation for the flatten surface in the XY plane.    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A planar surface in the XY plane approximating the     M
*                     falttening process of Srf.			     M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbPrisaRuledSrf, layout, prisa                                         M
*****************************************************************************/
CagdSrfStruct *SymbPrisaRuledSrf(CagdSrfStruct *Srf,
				 int SamplesPerCurve,
				 CagdRType Space,
				 CagdVType Offset)
{
    int i, j,
	VLength = Srf -> VLength;
    CagdRType Angle, PtTmp1[3], PtTmp2[3], PtTmp3[3];
    CagdCrvStruct
	*Crv1 = CagdCrvFromMesh(Srf, 0, CAGD_CONST_V_DIR),
	*Crv2 = CagdCrvFromMesh(Srf, VLength - 1, CAGD_CONST_V_DIR);
    CagdPolylineStruct
	*Poly1 = SymbCrv2Polyline(Crv1, SamplesPerCurve, FALSE, TRUE),
	*Poly2 = SymbCrv2Polyline(Crv2, SamplesPerCurve, FALSE, TRUE),
	*Poly1Prisa = CagdPolylineNew(Poly1 -> Length),
	*Poly2Prisa = CagdPolylineNew(Poly2 -> Length);
    CagdPolylnStruct
        *Pt1 = Poly1 -> Polyline,
        *Pt2 = Poly2 -> Polyline,
        *Pt1Prisa = Poly1Prisa -> Polyline,
        *Pt2Prisa = Poly2Prisa -> Polyline;
    CagdMType Mat1, Mat;
    CagdBBoxStruct BBox;

    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    /* Anchor the location of the first line. */
    for (j = 0; j < 3; j++) {
	Pt1Prisa -> Pt[j] = 0.0;
	Pt2Prisa -> Pt[j] = 0.0;
    }
    PT_SUB(PtTmp1, Pt1 -> Pt, Pt2 -> Pt);
    Pt2Prisa -> Pt[1] = PT_LENGTH(PtTmp1);

    /* Move alternately along Poly1 and Poly2 and anchor the next point of */
    /* the next triangle using the distances to current Pt1 and Pt2.	   */
    for (i = 2; i < Poly1 -> Length + Poly2 -> Length; i++) {
	CagdRType Dist1, Dist2, Inter1[3], Inter2[3],
	    *NextPt = i & 0x01 ? Pt1[1].Pt : Pt2[1].Pt;

	/* Compute distance from previous two locations to new location. */ 
	PT_SUB(PtTmp1, Pt1 -> Pt, NextPt);
	Dist1 = PT_LENGTH(PtTmp1);
	PT_SUB(PtTmp1, Pt2 -> Pt, NextPt);
	Dist2 = PT_LENGTH(PtTmp1);

	/* Find the (two) intersection points of circles with radii Dist?. */
	SymbInterCircCirc(Pt1Prisa -> Pt, Dist1, Pt2Prisa -> Pt, Dist2,
			  Inter1, Inter2);

	/* Find which of the two intersection points is the "good" one. */
	PT_SUB(PtTmp1, Inter1, Pt1Prisa -> Pt);
	PT_SUB(PtTmp2, Inter1, Pt2Prisa -> Pt);
	CROSS_PROD(PtTmp3, PtTmp1, PtTmp2);
	if (i & 0x01) {
	    Pt1Prisa++;
	    for (j = 0; j < 3; j++)
		Pt1Prisa -> Pt[j] = (PtTmp3[2] > 0 ? Inter2[j] : Inter1[j]);
	    Pt1++;
	}
	else {
	    Pt2Prisa++;
	    for (j = 0; j < 3; j++)
		Pt2Prisa -> Pt[j] = (PtTmp3[2] > 0 ? Inter2[j] : Inter1[j]);
	    Pt2++;
	}
    }

    /* Save centering location so we can orient the resulting data nicely. */
    PT_COPY(PtTmp1, Poly1Prisa -> Polyline[Poly1Prisa -> Length / 2].Pt);
    PT_COPY(PtTmp2, Poly2Prisa -> Polyline[Poly2Prisa -> Length / 2].Pt);
    PT_SUB(PtTmp3, PtTmp2, PtTmp1);

    Crv1 = CnvrtPolyline2LinBsplineCrv(Poly1Prisa);
    Crv2 = CnvrtPolyline2LinBsplineCrv(Poly2Prisa);
    CagdPolylineFree(Poly1);
    CagdPolylineFree(Poly2);
    CagdPolylineFree(Poly1Prisa);
    CagdPolylineFree(Poly2Prisa);

    Srf = CagdRuledSrf(Crv1, Crv2, 2, 2);
    CagdCrvFree(Crv1);
    CagdCrvFree(Crv2);

    /* Translate PtTmp1 to the origin. */
    MatGenMatTrans(-PtTmp1[0], -PtTmp1[1], -PtTmp1[2], Mat);

    /* Rotate PtTmp3 direction to the +Y direction. */
    Angle = atan2(PtTmp3[1], PtTmp3[0]);
    MatGenMatRotZ1(M_PI / 2 - Angle, Mat1);

    MatMultTwo4by4(Mat, Mat, Mat1);

    CagdSrfMatTransform(Srf, Mat);

    /* Translate by the Offset. */
    CagdSrfBBox(Srf, &BBox);
    MatGenMatTrans(Offset[0], Offset[1] - BBox.Min[1], Offset[2], Mat);
    Offset[1] += (BBox.Max[1] - BBox.Min[1]) + Space;

    CagdSrfMatTransform(Srf, Mat);

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Finds the two intersection points of the given two planar circles. It is   *
* assumed intersection exists.						     *
*                                                                            *
* PARAMETERS:                                                                *
*   Center1, Radius1:  Geometry of first circle.                             *
*   Center2, Radius2:  Geometry of second circle.                            *
*   Inter1, Inter2:    Where the two intersection locations will be placed.  *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SymbInterCircCirc(CagdRType Center1[3],
			      CagdRType Radius1,
			      CagdRType Center2[3],
			      CagdRType Radius2,
			      CagdRType Inter1[3],
			      CagdRType Inter2[3])
{
    CagdRType A, B, C, Delta,
	a = Center2[0] - Center1[0],
	b = Center2[1] - Center1[1],
	c = (SQR(Radius1) - SQR(Radius2) +
	     SQR(Center2[0]) - SQR(Center1[0]) +
	     SQR(Center2[1]) - SQR(Center1[1])) / 2;

    if (PT_APX_EQ(Center1, Center2)) {
	Inter1[0] = Inter2[0] = Radius1;
	Inter1[1] = Inter2[1] = 0.0;
    }
    else if (ABS(a) > ABS(b)) {
	/* Solve for Y first. */
	A = SQR(b) / SQR(a) + 1;
	B = 2 * (b * Center1[0] / a - b * c / SQR(a) - Center1[1]);
	C = SQR(c / a) + SQR(Center1[0]) + SQR(Center1[1])
				-2 * c * Center1[0] / a - SQR(Radius1);
	Delta = SQR(B) - 4 * A * C;
	if (Delta < 0)  /* If no solution, do something almost reasonable. */
	    Delta = 0;
	Inter1[1] = (-B + sqrt(Delta)) / (2 * A);
	Inter2[1] = (-B - sqrt(Delta)) / (2 * A);

	Inter1[0] = (c - b * Inter1[1]) / a;
	Inter2[0] = (c - b * Inter2[1]) / a;
    }
    else {
	/* Solve for X first. */
	A = SQR(a) / SQR(b) + 1;
	B = 2 * (a * Center1[1] / b - a * c / SQR(b) - Center1[0]);
	C = SQR(c / b) + SQR(Center1[0]) + SQR(Center1[1])
				-2 * c * Center1[1] / b - SQR(Radius1);
	Delta = SQR(B) - 4 * A * C;
	if (Delta < 0)  /* If no solution, do something almost reasonable. */
	    Delta = 0;
	Inter1[0] = (-B + sqrt(Delta)) / (2 * A);
	Inter2[0] = (-B - sqrt(Delta)) / (2 * A);

	Inter1[1] = (c - a * Inter1[0]) / b;
	Inter2[1] = (c - a * Inter2[0]) / b;
    }

    Inter1[2] = Inter2[2] = 0.0;
}
