/****************************************************************************** * Prisa.c - piecewise ruled srf approximation and layout (prisa). * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * 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, SYMB_CRV_APPROX_UNIFORM, TRUE), *Poly2 = SymbCrv2Polyline(Crv2, SamplesPerCurve, SYMB_CRV_APPROX_UNIFORM, 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; }