/****************************************************************************** * SCnvHul.c - computation of convex hull of a freeform surface. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Joon_Kyung Seong, Nov 2002 * ******************************************************************************/ #include "attribut.h" #include "iritgrap.h" #include "allocate.h" #include "mvar_lib.h" #include "ip_cnvrt.h" #include "geom_lib.h" #include "user_lib.h" #include "symb_loc.h" static IPObjectStruct *CHTwoSurface(CagdSrfStruct *Srf, CagdRType Fineness); static IPObjectStruct *CHThreeSurface(CagdSrfStruct *Srf, CagdRType Fineness); static double DistanceBoundary(double Pt1, double Pt2); static MvarPtStruct **ProcessUVDomainAux(CagdSrfStruct *Srf, MvarPtStruct *InPts, int *Cluster, int **Size, int TriTangents); static MvarPtStruct *CnvxSolveConstraints(CagdSrfStruct *Srf, CagdRType Fineness, int TriTangents); static MvarPtStruct *CnvxSolveConstraintsAux(CagdSrfStruct *Srf, CagdRType Fineness); static IPObjectStruct *MakeTriTangentSurface(CagdSrfStruct *Srf, MvarPtStruct **ConPts1, MvarPtStruct **ConPts2, MvarPtStruct **ConPts3, int *Size1, int *Size2, int *Size3, MvarPtStruct *TriPts); static IPObjectStruct *CombineCurves2Surfaces(CagdSrfStruct *Srf, CagdCrvStruct *CHCrvs, CagdCrvStruct *TCrv, int Parameter[3][4], PointType Pts[6], MvarPtStruct *TriPts); static MvarPtStruct *ProcessTriSolutions(CagdSrfStruct *Srf, MvarPtStruct *TriPts, double Tol); /***************************************************************************** * DESCRIPTION: M * Computes the Convex Hull of free-form surfaces. M * * * PARAMETERS: M * Srf: The two or three surfaces. M * Fineness: Accuracy of zero set computation as part of the solution. M * Value of 0.2 is a good start. M * * * RETURN VALUE: M * IPObjectStruct *: M * For two-surfaces-case: Developable surfaces which are the M * convex hull of input surfaces. For each developable M * surface, it has two trimming curves as an attribute, M * which is defined in the parameter domain. M * For three-surfaces-case: Developable surfaces which are the M * convex hull of input surfaces. For each developable M * surface, it has two trimming curves as an attribute, M * which is defined in the parameter domain. M * * * KEYWORDS: M * SymbSrfCnvxHull M *****************************************************************************/ IPObjectStruct *SymbSrfCnvxHull(CagdSrfStruct *Srf, CagdRType Fineness) { CagdSrfStruct *TSrf, *CpSrfs = NULL; IPObjectStruct *RetVal; if (Srf == NULL || Srf -> Pnext == NULL) { SYMB_FATAL_ERROR(SYMB_ERR_ILLEGAL_PARAMETERS); return NULL; } for ( ; Srf != NULL; Srf = Srf -> Pnext) { TSrf = CagdSrfCopy(Srf); if (CAGD_IS_BSPLINE_SRF(TSrf)) CagdSrfSetDomain(TSrf, 0.0, 1.0, 0.0, 1.0); LIST_PUSH(TSrf, CpSrfs); } if (CpSrfs -> Pnext -> Pnext == NULL) RetVal = CHTwoSurface(CpSrfs, Fineness); else RetVal = CHThreeSurface(CpSrfs, Fineness); CagdSrfFreeList(CpSrfs); return RetVal; } /***************************************************************************** * DESCRIPTION: M * Computes the bitangent developables between two free-form surfaces. M * * * PARAMETERS: M * Srf1: The first surface. M * Srf2: The second surface. M * Fineness: Accuracy of zero set computation as part of the solution. M * Value of 0.2 is a good start. M * * * RETURN VALUE: M * IPObjectStruct *: Developable bitangent surfaces. Actually, these bi- M * tangent surfaces are the boundary of the convex hull M * of two surfaces. M * * * KEYWORDS: M * SymbTwoSrfTangencies M *****************************************************************************/ IPObjectStruct *SymbTwoSrfTangencies(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2, CagdRType Fineness) { IPObjectStruct *CH; CagdSrfStruct *Srfs = CagdSrfCopy(Srf1); Srfs -> Pnext = CagdSrfCopy(Srf2); if (CAGD_IS_BSPLINE_SRF(Srfs)) CagdSrfSetDomain(Srfs, 0.0, 1.0, 0.0, 1.0); if (CAGD_IS_BSPLINE_SRF(Srfs -> Pnext)) CagdSrfSetDomain(Srfs -> Pnext, 0.0, 1.0, 0.0, 1.0); CH = CHTwoSurface(Srfs, Fineness); CagdSrfFreeList(Srfs); return CH; } /***************************************************************************** * DESCRIPTION: M * Computes the all the common tangencies of three surfaces. M * * * PARAMETERS: M * Srf1: The first surface. M * Srf2: The second surface. M * Srf3: The third surface. M * Fineness: Accuracy of zero set computation as part of the solution. M * Value of 0.4 is a good start. M * * * RETURN VALUE: M * IPObjectStruct *: All common tangencies of three surfaces. This object M * is a linked list of triangles which are common M * tangent to all the three surfaces. M * * * KEYWORDS: M * SymbThreeSrfTangencies M *****************************************************************************/ IPObjectStruct *SymbThreeSrfTangencies(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2, CagdSrfStruct *Srf3, CagdRType Fineness) { IPObjectStruct *Tangent, Tangents, *RetVal = IPGenLISTObject(NULL); MvarPtStruct *TriPts, *TriPt; CagdSrfStruct *Srfs = CagdSrfCopy(Srf1); int Count = 0; Srfs -> Pnext = CagdSrfCopy(Srf2); Srfs -> Pnext -> Pnext = CagdSrfCopy(Srf3); if (CAGD_IS_BSPLINE_SRF(Srfs)) CagdSrfSetDomain(Srfs, 0.0, 1.0, 0.0, 1.0); if (CAGD_IS_BSPLINE_SRF(Srfs -> Pnext)) CagdSrfSetDomain(Srfs -> Pnext, 0.0, 1.0, 0.0, 1.0); if (CAGD_IS_BSPLINE_SRF(Srfs -> Pnext -> Pnext)) CagdSrfSetDomain(Srfs -> Pnext -> Pnext, 0.0, 1.0, 0.0, 1.0); TriPts = CnvxSolveConstraintsAux(Srfs, Fineness); TriPts = ProcessTriSolutions(Srfs, TriPts, 0.001); Tangent = &Tangents; for (TriPt = TriPts; TriPt; TriPt = TriPt -> Pnext) { double *R; PointType Pt1, Pt2, Pt3; IPPolygonStruct *Pl = IPAllocPolygon(0, NULL, NULL); Pl -> PVertex = IPAllocVertex2(NULL); Pl -> PVertex -> Pnext = IPAllocVertex2(NULL); Pl -> PVertex -> Pnext -> Pnext = IPAllocVertex2(NULL); /* Evaluate Tri tangents solution points. */ R = CagdSrfEval(Srf1, TriPt -> Pt[0], TriPt -> Pt[1]); CagdCoerceToE3(Pt1, &R, -1, Srf1 -> PType); R = CagdSrfEval(Srf2, TriPt -> Pt[2], TriPt -> Pt[3]); CagdCoerceToE3(Pt2, &R, -1, Srf2 -> PType); R = CagdSrfEval(Srf3, TriPt -> Pt[4], TriPt -> Pt[5]); CagdCoerceToE3(Pt3, &R, -1, Srf3 -> PType); PT_COPY(Pl -> PVertex -> Coord, Pt1); AttrSetUVAttrib(&Pl -> PVertex -> Attr, "uvvals", TriPt -> Pt[0], TriPt -> Pt[1]); PT_COPY(Pl -> PVertex -> Pnext -> Coord, Pt2); AttrSetUVAttrib(&Pl -> PVertex -> Pnext -> Attr, "uvvals", TriPt -> Pt[2], TriPt -> Pt[3]); PT_COPY(Pl -> PVertex -> Pnext -> Pnext -> Coord, Pt3); AttrSetUVAttrib(&Pl -> PVertex -> Pnext -> Pnext -> Attr, "uvvals", TriPt -> Pt[4], TriPt -> Pt[5]); Tangent -> Pnext = IPGenPOLYObject(Pl); Tangent = Tangent -> Pnext; } Tangent -> Pnext = NULL; CagdSrfFreeList(Srfs); /* make a list of objects */ for (Tangent = Tangents.Pnext; Tangent != NULL; Tangent = Tangent -> Pnext) IPListObjectInsert(RetVal, Count++, Tangent); IPListObjectInsert(RetVal, Count, NULL); return RetVal; } /***************************************************************************** * DESCRIPTION: * * Computes the Convex Hull of two free-form surfaces. * * * * PARAMETERS: * * Srf: The two surfaces. * * Fineness: Accuracy of zero set computation as part of the solution. * * Value of 0.2 is a good start. * * * * RETURN VALUE: * * IPObjectStruct *: Developable surfaces which form the boundary of * * the convex hull of two surfaces. For each developable * * surface, it has two trimming curves as an attribute, * * which is defined in the parameter domain. * *****************************************************************************/ static IPObjectStruct *CHTwoSurface(CagdSrfStruct *Srf, CagdRType Fineness) { int *Size, i, j, Cluster = 0, Count = 0; MvarPtStruct *CHPts, *ConPt, **ConPts; CagdCrvStruct *CHCrv, *TrimCrv, CHCrvs, TrimCrvs; IPObjectStruct *CnvxHull, CnvxHulls, *RetVal = IPGenLISTObject(NULL); /* Solve the constraints. */ CHPts = CnvxSolveConstraints(Srf, Fineness, 1); /* Resolve the solution points. */ ConPts = ProcessUVDomainAux(Srf, CHPts, &Cluster, &Size, 1); /* Connect the solution points into curves. */ CHCrv = &CHCrvs; TrimCrv = &TrimCrvs; for (i = 0; i < Cluster; i++) { CagdRType **ConvexHull1, **ConvexHull2, **Trim1, **Trim2, Cp[4]; CagdCrvStruct *Crv2, *Crv1, *TCrv1, *TCrv2; CagdSrfStruct *Srf1, *Srf2; Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext); Crv1 = BspCrvNew(Size[i], 3, CAGD_PT_E3_TYPE); LIST_PUSH(Crv1, CHCrv -> Pnext); ConvexHull1 = CHCrv -> Pnext -> Points; BspKnotUniformOpen(Size[i], 3, CHCrv -> Pnext -> KnotVector); CHCrv = CHCrv -> Pnext; Crv2 = BspCrvNew(Size[i], 3, CAGD_PT_E3_TYPE); LIST_PUSH(Crv2, CHCrv -> Pnext); ConvexHull2 = CHCrv -> Pnext -> Points; BspKnotUniformOpen(Size[i], 3, CHCrv -> Pnext -> KnotVector); /* Create new curves for trimed curve. */ TCrv1 = BspCrvNew(Size[i], 2, CAGD_PT_E2_TYPE); LIST_PUSH(TCrv1, TrimCrv -> Pnext); Trim1 = TrimCrv -> Pnext -> Points; BspKnotUniformOpen(Size[i], 2, TrimCrv -> Pnext -> KnotVector); TrimCrv = TrimCrv -> Pnext; TCrv2 = BspCrvNew(Size[i], 2, CAGD_PT_E2_TYPE); LIST_PUSH(TCrv2, TrimCrv -> Pnext); Trim2 = TrimCrv -> Pnext -> Points; BspKnotUniformOpen(Size[i], 2, TrimCrv -> Pnext -> KnotVector); ConPt = ConPts[i]; for (j = 0; j < Size[i]; j++) { double *R; PointType Pt1; Cp[0] = ConPt -> Pt[0]; Cp[1] = ConPt -> Pt[1]; Cp[2] = ConPt -> Pt[2]; Cp[3] = ConPt -> Pt[3]; ConPt = ConPt -> Pnext; R = CagdSrfEval(Srf1, Cp[0], Cp[1]); CagdCoerceToE3(Pt1, &R, -1, Srf1 -> PType); ConvexHull1[1][j] = Pt1[0]; ConvexHull1[2][j] = Pt1[1]; ConvexHull1[3][j] = Pt1[2]; R = CagdSrfEval(Srf2, Cp[2], Cp[3]); CagdCoerceToE3(Pt1, &R, -1, Srf2 -> PType); ConvexHull2[1][j] = Pt1[0]; ConvexHull2[2][j] = Pt1[1]; ConvexHull2[3][j] = Pt1[2]; Trim1[1][j] = Cp[0]; Trim1[2][j] = Cp[1]; Trim2[1][j] = Cp[2]; Trim2[2][j] = Cp[3]; } CHCrv = CHCrv -> Pnext; TrimCrv = TrimCrv -> Pnext; } CHCrv -> Pnext = NULL; TrimCrv -> Pnext = NULL; /* Add a list of developable surfaces. */ CnvxHull = &CnvxHulls; for (i = 0, CHCrv = CHCrvs.Pnext, TrimCrv = TrimCrvs.Pnext; i < Cluster && CHCrv != NULL; i++) { CagdCrvStruct *Ruled1, *Ruled2, *TrimC1, *TrimC2; IPAttributeStruct *Attr; Ruled1 = CagdCrvCopy(CHCrv); Ruled1 -> Pnext = NULL; Ruled2 = CagdCrvCopy(CHCrv -> Pnext); Ruled2 -> Pnext = NULL; TrimC1 = CagdCrvCopy(TrimCrv); TrimC1 -> Pnext = NULL; TrimC2 = CagdCrvCopy(TrimCrv -> Pnext); TrimC2 -> Pnext = NULL; CnvxHull -> Pnext = IPGenSRFObject(CagdRuledSrf(Ruled1, Ruled2, 2, 2)); Attr = _AttrMallocAttribute("trim1", IP_ATTR_OBJ); Attr -> Pnext = CnvxHull -> Pnext -> Attr; CnvxHull -> Pnext -> Attr = Attr; Attr -> U.PObj = IPGenCRVObject(TrimC1); Attr -> Pnext = _AttrMallocAttribute("trim2", IP_ATTR_OBJ); Attr -> Pnext -> Pnext = NULL; Attr -> Pnext -> U.PObj = IPGenCRVObject(TrimC2); CnvxHull = CnvxHull -> Pnext; CHCrv = CHCrv -> Pnext -> Pnext; TrimCrv = TrimCrv -> Pnext -> Pnext; } CnvxHull -> Pnext = NULL; /* Make a list of objects. */ for (CnvxHull = CnvxHulls.Pnext; CnvxHull != NULL; CnvxHull = CnvxHull -> Pnext) IPListObjectInsert(RetVal, Count++, CnvxHull); IPListObjectInsert(RetVal, Count, NULL); CagdCrvFreeList(CHCrv); CagdCrvFreeList(TrimCrv); return RetVal; } /***************************************************************************** * DESCRIPTION: * * Compute distance between the point and the boundary. * * * * PARAMETERS: * * Pt1: The first uv point. * * Pt2: The second uv point. * * * * RETURN VALUE: * * double: Calculated distance between the given point and the boundary.* *****************************************************************************/ static double DistanceBoundary(double Pt1, double Pt2) { double L1, L2; L1 = (Pt1 < (1.0 - Pt1) ? Pt1 : 1.0 - Pt1); L2 = (Pt2 < (1.0 - Pt2) ? Pt2 : 1.0 - Pt2); return L1 < L2 ? L1 : L2; } /***************************************************************************** * DESCRIPTION: * * Cluster uv boundary segments and connect them. * * * * PARAMETERS: * * Srf: The surface which is defined on this uv domain. * * InPts: Solution points which are from the constraint solver. * * Cluster: The number of cluster whichi will be returned. * * Size: Size of each cluster. * * TriTangents: Index which tells the target surfaces. * * if 1 then, first and second surface will be considered * * else if 2 then, second and third surface will be considered * * else if 3 then, third and first surface will be considered * * * * RETURN VALUE: * * MvarPtStruct **: Linked list of resolved points. * *****************************************************************************/ static MvarPtStruct **ProcessUVDomainAux(CagdSrfStruct *Srf, MvarPtStruct *InPts, int *Cluster, int **Size, int TriTangents) { MvarPtStruct *CHPt, *Marker, *ConPt, *ConPt2, **RtnPts, *CHPts = MvarPtCopyList(InPts), **ConPts = (MvarPtStruct **) IritMalloc(sizeof(MvarPtStruct *) * 50); CagdSrfStruct *Srf1, *Srf2; int i, j, TmpCluster, LCluster = 0, Skip[6] = {0, 0, 0, 0, 0, 0}; (*Size) = (int *) IritMalloc(sizeof(int) * 50); if (TriTangents == 1) { Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext); } else if (TriTangents == 2) { Srf1 = CagdSrfCopy(Srf -> Pnext); Srf2 = CagdSrfCopy(Srf -> Pnext -> Pnext); } else if (TriTangents == 3) { Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext -> Pnext); } else { Srf1 = CagdSrfCopy(Srf); if (Srf -> Pnext == NULL) Srf2 = CagdSrfCopy(Srf); else Srf2 = CagdSrfCopy(Srf -> Pnext); } /* Test for inner-cross solution. */ for (CHPt = CHPts; CHPt != NULL; CHPt = CHPt -> Pnext) { CagdVecStruct *N; double N1[3], N2[3], Prod; N = CagdSrfNormal(Srf1, CHPt -> Pt[0], CHPt -> Pt[1], TRUE); N1[0] = N -> Vec[0]; N1[1] = N -> Vec[1]; N1[2] = N -> Vec[2]; N = CagdSrfNormal(Srf2, CHPt -> Pt[2], CHPt -> Pt[3], TRUE); N2[0] = N ->Vec[0]; N2[1] = N ->Vec[1]; N2[2] = N ->Vec[2]; Prod = DOT_PROD(N1, N2); if (! APX_EQ_EPS(Prod, 1.0, 0.001)) CHPt -> Dim = -1; } /* Connect the solution points into curves. */ for (CHPt = CHPts; CHPt != NULL; ) { double Small = 10.0; MvarPtStruct *TargetPt; if (CHPt -> Dim == -1) { /* We already visited this point. */ CHPt = CHPt -> Pnext; continue; } /* We should find the point which is closest to the boundary */ /* and it becomes the first solution. */ for (ConPt = CHPt; ConPt != NULL; ConPt = ConPt -> Pnext) { double NearBoundary = 0.0; if (ConPt -> Dim != -1) { NearBoundary = DistanceBoundary(ConPt -> Pt[0], ConPt -> Pt[1]); if (Small > NearBoundary) { Small = NearBoundary; TargetPt = ConPt; } } } ConPt = ConPts[LCluster] = MvarPtCopy(TargetPt); TargetPt -> Dim = -1; (*Size)[LCluster] = 1; Marker = CHPt; while (TRUE) { double STCheck; Small = 10.0; TargetPt = NULL; for (ConPt2 = CHPt; ConPt2 != NULL; ConPt2 = ConPt2 -> Pnext) { if (ConPt2 -> Dim != -1) { double Distance = (ConPt2 -> Pt[0] - ConPt -> Pt[0]) * (ConPt2 -> Pt[0] - ConPt -> Pt[0]) + (ConPt2 -> Pt[1] - ConPt -> Pt[1]) * (ConPt2 -> Pt[1] - ConPt -> Pt[1]); STCheck = (ConPt2 -> Pt[2] - ConPt -> Pt[2]) * (ConPt2 -> Pt[2] - ConPt -> Pt[2]) + (ConPt2 -> Pt[3] - ConPt -> Pt[3]) * (ConPt2 -> Pt[3] - ConPt -> Pt[3]); if (ConPt2 -> Pt[2] < 0.001) STCheck = 0.1; if (Small > Distance && STCheck < 0.2) { Small = Distance; TargetPt = ConPt2; } } } /* If distance is somewhat large, then disconnect this. */ if (Small > 0.2 || Small == 10.0) break; ConPt -> Pnext = MvarPtCopy(TargetPt); (*Size)[LCluster]++; ConPt = ConPt -> Pnext; ConPt -> Pnext = NULL; TargetPt -> Dim = -1; } LCluster++; CHPt = Marker; } /* If size of connection is smaller than 2 then pulge it out. */ TmpCluster = LCluster; RtnPts = (MvarPtStruct **)IritMalloc(sizeof(MvarPtStruct *) * TmpCluster); for (i = 0, j = 0; i < LCluster; i++) { if ((*Size)[i] >= 2 && ! Skip[i]) { int TmpSize = (*Size)[i]; RtnPts[j] = MvarPtCopyList(ConPts[i]); (*Size)[j] = TmpSize; j++; } else TmpCluster--; } LCluster = TmpCluster; *Cluster = LCluster; return RtnPts; } /***************************************************************************** * DESCRIPTION: * * Solve three constraints for two surfaces simultaneously. * * (1) F(u,v,s,t) = det[ S(u,v)-S(s,t), Su, Sv ] = 0, * * (2) Fs = 0, * * (3) Ft = 0. * * the solution points are belonged to the boundary of the convex hull. * * * * PARAMETERS: * * Srf: The two surfaces. * * Fineness: Accuracy of zero set computation as part of the solution. * * Value of 0.2 is a good start. * * TriTangents: Index which tells the target surfaces * * if 1 then, first and second surface will be considered * * else if 2 then, second and third surface will be considered * * else if 3 then, third and first surface will be considered. * * * * RETURN VALUE: * * MvarPtStruct *: Linked list of solution points. * *****************************************************************************/ static MvarPtStruct *CnvxSolveConstraints(CagdSrfStruct *Srf, CagdRType Fineness, int TriTangents) { CagdSrfStruct *Srf1, *Srf2; MvarConstraintType Constraints[3] = { MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO }; MvarMVStruct *SrfMV1, *SrfMV2, *MV1, *MV2, *MVF, **Split1, **Split2, **Split3, *A[3], *B[3], *C[3], *MVTmp1, *MVTmp2, *MVTmp3, *Finals[3]; MvarPtStruct *Pts; double SubdivTol = Fineness, NumericTol = -0.01; if (TriTangents == 1) { Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext); } else if (TriTangents == 2) { Srf1 = CagdSrfCopy(Srf -> Pnext); Srf2 = CagdSrfCopy(Srf -> Pnext -> Pnext); } else if (TriTangents == 3) { Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext -> Pnext); } else { Srf1 = CagdSrfCopy(Srf); if (Srf -> Pnext == NULL) Srf2 = CagdSrfCopy(Srf); else Srf2 = CagdSrfCopy(Srf -> Pnext); } Srf1 -> Pnext = NULL; Srf2 -> Pnext = NULL; SrfMV1 = MvarSrfToMV(Srf1); SrfMV2 = MvarSrfToMV(Srf2); MV1 = MvarPromoteMVToMV2(SrfMV1, 4, 0); /* Four variate at axes 0, 1. */ MV2 = MvarPromoteMVToMV2(SrfMV2, 4, 2); /* Four variate at axes 2, 3. */ MvarMakeMVsCompatible(&MV1, &MV2, 1, 1); /* Compute f(u,v,s,t) multivariate. */ MVTmp1 = MvarMVDerive(MV1, 0); MVTmp2 = MvarMVDerive(MV1, 1); MVTmp3 = MvarMVSub(MV1, MV2); MvarMVFree(MV1); MvarMVFree(MV2); Split1 = MvarMVSplitScalar(MVTmp1); MvarMVFree(MVTmp1); A[0] = MvarMVCopy(Split1[1]); A[1] = MvarMVCopy(Split1[2]); A[2] = MvarMVCopy(Split1[3]); Split2 = MvarMVSplitScalar(MVTmp2); MvarMVFree(MVTmp2); B[0] = MvarMVCopy(Split2[1]); B[1] = MvarMVCopy(Split2[2]); B[2] = MvarMVCopy(Split2[3]); Split3 = MvarMVSplitScalar(MVTmp3); MvarMVFree(MVTmp3); C[0] = MvarMVCopy(Split3[1]); C[1] = MvarMVCopy(Split3[2]); C[2] = MvarMVCopy(Split3[3]); MVF = MvarMVDeterminant3(A[0], A[1], A[2], B[0], B[1], B[2], C[0], C[1], C[2]); Finals[0] = MvarMVCopy(MVF); Finals[1] = MvarMVDerive(MVF, 2); Finals[2] = MvarMVDerive(MVF, 3); Pts = MvarMVsZeros(Finals, Constraints, 3, SubdivTol, NumericTol); MvarMVFree(Finals[0]); MvarMVFree(Finals[1]); MvarMVFree(Finals[2]); return Pts; } /***************************************************************************** * DESCRIPTION: * * Computes the Convex Hull of three free-form surfaces. * * * * PARAMETERS: * * Srf: The three surfaces. * * Fineness: Accuracy of zero set computation as part of the solution. * * Value of 0.2 is a good start. * * * * RETURN VALUE: * * IPObjectStruct *: Developable surfaces which form the boundary of * * the convex hull of three surfaces. For each * * developable surface, it has two trimming curves * * as an attribute, which is defined in the parameter * * domain. * *****************************************************************************/ static IPObjectStruct *CHThreeSurface(CagdSrfStruct *Srf, CagdRType Fineness) { MvarPtStruct *CHPts1, *CHPts2, *CHPts3, *TriPts, **ConPts1, **ConPts2, **ConPts3; int *Size1, *Size2, *Size3, Cluster1 = 0, Cluster2 = 0, Cluster3 = 0; /* Solve the constraints for each pair of surfaces. */ CHPts1 = CnvxSolveConstraints(Srf, Fineness * 0.5, 1); CHPts2 = CnvxSolveConstraints(Srf, Fineness * 0.5, 2); CHPts3 = CnvxSolveConstraints(Srf, Fineness * 0.5, 3); /* And solve common tangent to the three surfaces. */ TriPts = CnvxSolveConstraintsAux(Srf, Fineness); /* Resolve the solution points for each solutions. */ ConPts1 = ProcessUVDomainAux(Srf, CHPts1, &Cluster1, &Size1, 1); ConPts2 = ProcessUVDomainAux(Srf, CHPts2, &Cluster2, &Size2, 2); ConPts3 = ProcessUVDomainAux(Srf, CHPts3, &Cluster3, &Size3, 3); /* And also for the tri-tangents. */ TriPts = ProcessTriSolutions(Srf, TriPts, 0.001); /* Make the convex hull and return it. */ return MakeTriTangentSurface(Srf, ConPts1, ConPts2, ConPts3, Size1, Size2, Size3, TriPts); } /***************************************************************************** * DESCRIPTION: * * Solve six constraints for three surfaces simultaneously. * * (1) F(u,v,s,t) = det[ S(u,v)-S(s,t), Su, Sv ] = 0, * * (2) Fs = 0, * * (3) Ft = 0. * * (4) G(u,v,m,n) = det[ S(u,v)-S(m,n), Su, Sv ] = 0, * * (5) Gm = 0, * * (6) Gn = 0. * * the solution points are belonged to the boundary of the convex hull * * * * PARAMETERS: * * Srf: The three surfaces. * * Fineness: Accuracy of zero set computation as part of the solution. * * Value of 0.2 is a good start. * * * * RETURN VALUE: * * MvarPtStruct *: Linked list of solution points. * *****************************************************************************/ static MvarPtStruct *CnvxSolveConstraintsAux(CagdSrfStruct *Srf, CagdRType Fineness) { CagdSrfStruct *Srf1, *Srf2, *Srf3; MvarConstraintType Constraints[6] = { MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO }; MvarMVStruct *SrfMV1, *SrfMV2, *SrfMV3, *MV1, *MV2, *MV3, *MVF, **Split1, **Split2, **Split3, **Split4, *A[3], *B[3], *C[3], *D[3], *MVTmp1, *MVTmp2, *MVTmp3, *MVTmp4, *Finals[6]; MvarPtStruct *Pts; double SubdivTol = Fineness, NumericTol = -0.01; Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext); Srf3 = CagdSrfCopy(Srf -> Pnext -> Pnext); Srf1 -> Pnext = NULL; Srf2 -> Pnext = NULL; Srf3 -> Pnext = NULL; SrfMV1 = MvarSrfToMV(Srf1); SrfMV2 = MvarSrfToMV(Srf2); SrfMV3 = MvarSrfToMV(Srf3); MV1 = MvarPromoteMVToMV2(SrfMV1, 6, 0); /* Six variate at axes 0, 1. */ MV2 = MvarPromoteMVToMV2(SrfMV2, 6, 2); /* Six variate at axes 2, 3. */ MV3 = MvarPromoteMVToMV2(SrfMV3, 6, 4); /* Six variate at axes 4, 5. */ /* Compute f(u,v,s,t) multivariate. */ MVTmp1 = MvarMVDerive(MV1, 0); MVTmp2 = MvarMVDerive(MV1, 1); MVTmp3 = MvarMVSub(MV1, MV2); Split1 = MvarMVSplitScalar(MVTmp1); A[0] = MvarMVCopy(Split1[1]); A[1] = MvarMVCopy(Split1[2]); A[2] = MvarMVCopy(Split1[3]); Split2 = MvarMVSplitScalar(MVTmp2); B[0] = MvarMVCopy(Split2[1]); B[1] = MvarMVCopy(Split2[2]); B[2] = MvarMVCopy(Split2[3]); Split3 = MvarMVSplitScalar(MVTmp3); MvarMVFree(MVTmp3); C[0] = MvarMVCopy(Split3[1]); C[1] = MvarMVCopy(Split3[2]); C[2] = MvarMVCopy(Split3[3]); MVF = MvarMVDeterminant3(A[0], A[1], A[2], B[0], B[1], B[2], C[0], C[1], C[2]); Finals[0] = MvarMVCopy(MVF); Finals[1] = MvarMVDerive(MVF, 2); Finals[2] = MvarMVDerive(MVF, 3); MvarMVFree(MVF); /* Make G(u,v,m,n). */ MVTmp4 = MvarMVSub(MV1, MV3); Split4 = MvarMVSplitScalar(MVTmp4); MvarMVFree(MVTmp4); D[0] = MvarMVCopy(Split4[1]); D[1] = MvarMVCopy(Split4[2]); D[2] = MvarMVCopy(Split4[3]); MVF = MvarMVDeterminant3(A[0], A[1], A[2], B[0], B[1], B[2], D[0], D[1], D[2]); Finals[3] = MvarMVCopy(MVF); Finals[4] = MvarMVDerive(MVF, 4); Finals[5] = MvarMVDerive(MVF, 5); Pts = MvarMVsZeros(Finals, Constraints, 6, SubdivTol, NumericTol); MvarMVFree(Finals[0]); MvarMVFree(Finals[1]); MvarMVFree(Finals[2]); MvarMVFree(Finals[3]); MvarMVFree(Finals[4]); MvarMVFree(Finals[5]); return Pts; } /***************************************************************************** * DESCRIPTION: * * Edit curve segments to connect each pair of curves of three surfaces * * and get new trimmed curve segments. * * * * PARAMETERS: * * Srf: The three surfaces, in a linked list. * * ConPts1, ConPts2, ConPts3: Three pairs of solution points which are * * generated from two surface pair. * * Size1, Size2, Size3: Size of the corresponding solution points. * * TriPts: Solution points from the constraints solver of * * three surfaces. * * * * RETURN VALUE: * * IPObjectStruct *: Developable surfaces forming the boundary of the * * convex hull of three surfaces. * *****************************************************************************/ static IPObjectStruct *MakeTriTangentSurface(CagdSrfStruct *Srf, MvarPtStruct **ConPts1, MvarPtStruct **ConPts2, MvarPtStruct **ConPts3, int *Size1, int *Size2, int *Size3, MvarPtStruct *TriPts) { MvarPtStruct *ConPt, *TmpPt; CagdCrvStruct *CHCrv, *TrimCrv, *CHCrvs = (CagdCrvStruct *) IritMalloc(sizeof(CagdCrvStruct)), *TrimCrvs = (CagdCrvStruct *) IritMalloc(sizeof(CagdCrvStruct)); int i, j, CurSize[3], KP[3][4], Parameter[3][4]; double *RR; PointType Pt1, Pts[6]; CurSize[0] = Size1[0]; CurSize[1] = Size2[0]; CurSize[2] = Size3[0]; /* Evaluate Tri tangents solution points. */ RR = CagdSrfEval(Srf, TriPts -> Pt[0], TriPts -> Pt[1]); CagdCoerceToE3(Pt1, &RR, -1, Srf -> PType); Pts[0][0] = Pt1[0]; Pts[0][1] = Pt1[1]; Pts[0][2] = Pt1[2]; RR = CagdSrfEval(Srf -> Pnext, TriPts -> Pt[2], TriPts -> Pt[3]); CagdCoerceToE3(Pt1, &RR, -1, Srf -> Pnext -> PType); Pts[1][0] = Pt1[0]; Pts[1][1] = Pt1[1]; Pts[1][2] = Pt1[2]; RR = CagdSrfEval(Srf -> Pnext -> Pnext, TriPts -> Pt[4], TriPts -> Pt[5]); CagdCoerceToE3(Pt1, &RR, -1, Srf -> Pnext -> Pnext -> PType); Pts[2][0] = Pt1[0]; Pts[2][1] = Pt1[1]; Pts[2][2] = Pt1[2]; RR = CagdSrfEval(Srf, TriPts -> Pnext -> Pt[0], TriPts -> Pnext -> Pt[1]); CagdCoerceToE3(Pt1, &RR, -1, Srf -> PType); Pts[3][0] = Pt1[0]; Pts[3][1] = Pt1[1]; Pts[3][2] = Pt1[2]; RR = CagdSrfEval(Srf -> Pnext, TriPts -> Pnext -> Pt[2], TriPts -> Pnext -> Pt[3]); CagdCoerceToE3(Pt1, &RR, -1, Srf -> Pnext -> PType); Pts[4][0] = Pt1[0]; Pts[4][1] = Pt1[1]; Pts[4][2] = Pt1[2]; RR = CagdSrfEval(Srf -> Pnext -> Pnext, TriPts -> Pnext -> Pt[4], TriPts -> Pnext -> Pt[5]); CagdCoerceToE3(Pt1, &RR, -1, Srf -> Pnext -> Pnext -> PType); Pts[5][0] = Pt1[0]; Pts[5][1] = Pt1[1]; Pts[5][2] = Pt1[2]; /* First, make full curve segments. */ CHCrv = CHCrvs; TrimCrv = TrimCrvs; for (i = 0; i < 3; i++) { CagdRType **ConvexHull1, **ConvexHull2, **Trim1, **Trim2, Thres1 = 10.0, Thres11 = 10.0, Thres2 = 10.0, Thres22 = 10.0; CagdCrvStruct *Crv2, *Crv1, *TCrv1, *TCrv2; CagdSrfStruct *Srf1, *Srf2; PointType P1, P2, P3, P4; MvarPtStruct *PrevPt; double CP[4]; if (i == 0) { Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext); ConPt = ConPts1[0]; P1[0] = TriPts -> Pt[0]; P1[1] = TriPts -> Pt[1]; P2[0] = TriPts -> Pt[2]; P2[1] = TriPts -> Pt[3]; P3[0] = TriPts -> Pnext -> Pt[0]; P3[1] = TriPts -> Pnext -> Pt[1]; P4[0] = TriPts -> Pnext -> Pt[2]; P4[1] = TriPts -> Pnext -> Pt[3]; } else if (i == 1) { Srf1 = CagdSrfCopy(Srf -> Pnext); Srf2 = CagdSrfCopy(Srf -> Pnext -> Pnext); ConPt = ConPts2[0]; P1[0] = TriPts -> Pt[2]; P1[1] = TriPts -> Pt[3]; P2[0] = TriPts -> Pt[4]; P2[1] = TriPts -> Pt[5]; P3[0] = TriPts -> Pnext -> Pt[2]; P3[1] = TriPts -> Pnext -> Pt[3]; P4[0] = TriPts -> Pnext -> Pt[4]; P4[1] = TriPts -> Pnext -> Pt[5]; } else if (i == 2) { Srf1 = CagdSrfCopy(Srf); Srf2 = CagdSrfCopy(Srf -> Pnext -> Pnext); ConPt = ConPts3[0]; P1[0] = TriPts -> Pt[0]; P1[1] = TriPts -> Pt[1]; P2[0] = TriPts -> Pt[4]; P2[1] = TriPts -> Pt[5]; P3[0] = TriPts -> Pnext -> Pt[0]; P3[1] = TriPts -> Pnext -> Pt[1]; P4[0] = TriPts -> Pnext -> Pt[4]; P4[1] = TriPts -> Pnext -> Pt[5]; } for (TmpPt = ConPt, j = 0; TmpPt != NULL; TmpPt = TmpPt -> Pnext, j++) { double Dist1, Dist11, Dist2, Dist22; CP[0] = TmpPt -> Pt[0]; CP[1] = TmpPt -> Pt[1]; CP[2] = TmpPt -> Pt[2]; CP[3] = TmpPt -> Pt[3]; /* Find the closest point to the solution of tri-tangents. */ switch(i) { case 0: Dist1 = (CP[0] - TriPts -> Pt[0]) * (CP[0] - TriPts -> Pt[0]) + (CP[1] - TriPts -> Pt[1]) * (CP[1] - TriPts -> Pt[1]); Dist2 = (CP[2] - TriPts -> Pt[2]) * (CP[2] - TriPts -> Pt[2]) + (CP[3] - TriPts -> Pt[3]) * (CP[3] - TriPts -> Pt[3]); Dist11 = (CP[0] - TriPts -> Pnext -> Pt[0]) * (CP[0] - TriPts -> Pnext -> Pt[0]) + (CP[1] - TriPts -> Pnext -> Pt[1]) * (CP[1] - TriPts -> Pnext -> Pt[1]); Dist22 = (CP[2] - TriPts -> Pnext -> Pt[2]) * (CP[2] - TriPts -> Pnext -> Pt[2]) + (CP[3] - TriPts -> Pnext -> Pt[3]) * (CP[3] - TriPts -> Pnext -> Pt[3]); break; case 1: Dist1 = (CP[0] - TriPts -> Pt[2]) * (CP[0] - TriPts -> Pt[2]) + (CP[1] - TriPts -> Pt[3]) * (CP[1] - TriPts -> Pt[3]); Dist2 = (CP[2] - TriPts -> Pt[4]) * (CP[2] - TriPts -> Pt[4]) + (CP[3] - TriPts -> Pt[5]) * (CP[3] - TriPts -> Pt[5]); Dist11 = (CP[0] - TriPts -> Pnext -> Pt[2]) * (CP[0] - TriPts -> Pnext -> Pt[2]) + (CP[1] - TriPts -> Pnext -> Pt[3]) * (CP[1] - TriPts -> Pnext -> Pt[3]); Dist22 = (CP[2] - TriPts -> Pnext -> Pt[4]) * (CP[2] - TriPts -> Pnext -> Pt[4]) + (CP[3] - TriPts -> Pnext -> Pt[5]) * (CP[3] - TriPts -> Pnext -> Pt[5]); break; case 2: Dist1 = (CP[0] - TriPts -> Pt[0]) * (CP[0] - TriPts -> Pt[0]) + (CP[1] - TriPts -> Pt[1]) * (CP[1] - TriPts -> Pt[1]); Dist2 = (CP[2] - TriPts -> Pt[4]) * (CP[2] - TriPts -> Pt[4]) + (CP[3] - TriPts -> Pt[5]) * (CP[3] - TriPts -> Pt[5]); Dist11 = (CP[0] - TriPts -> Pnext -> Pt[0]) * (CP[0] - TriPts -> Pnext -> Pt[0]) + (CP[1] - TriPts -> Pnext -> Pt[1]) * (CP[1] - TriPts -> Pnext -> Pt[1]); Dist22 = (CP[2] - TriPts -> Pnext -> Pt[4]) * (CP[2] - TriPts -> Pnext -> Pt[4]) + (CP[3] - TriPts -> Pnext -> Pt[5]) * (CP[3] - TriPts -> Pnext -> Pt[5]); break; } if (Dist1 < Thres1) { Thres1 = Dist1; KP[i][0] = j; } if (Dist2 < Thres2) { Thres2 = Dist2; KP[i][1] = j; } if (Dist11 < Thres11) { Thres11 = Dist11; KP[i][2] = j; } if (Dist22 < Thres22) { Thres22 = Dist22; KP[i][3] = j; } } /* Create new curves. */ Crv1 = BspCrvNew(CurSize[i] + 2, 2, CAGD_PT_E3_TYPE); LIST_PUSH(Crv1, CHCrv -> Pnext); ConvexHull1 = CHCrv -> Pnext -> Points; BspKnotUniformOpen(CurSize[i] + 2, 2, CHCrv -> Pnext -> KnotVector); CHCrv = CHCrv -> Pnext; Crv2 = BspCrvNew(CurSize[i] + 2, 2, CAGD_PT_E3_TYPE); LIST_PUSH(Crv2, CHCrv -> Pnext); ConvexHull2 = CHCrv -> Pnext -> Points; BspKnotUniformOpen(CurSize[i] + 2, 2, CHCrv -> Pnext -> KnotVector); /* Create new curves for trimed curve. */ TCrv1 = BspCrvNew(CurSize[i] + 2, 2, CAGD_PT_E2_TYPE); LIST_PUSH(TCrv1, TrimCrv -> Pnext); Trim1 = TrimCrv -> Pnext -> Points; BspKnotUniformOpen(CurSize[i] + 2, 2, TrimCrv -> Pnext -> KnotVector); TrimCrv = TrimCrv -> Pnext; TCrv2 = BspCrvNew(CurSize[i] + 2, 2, CAGD_PT_E2_TYPE); LIST_PUSH(TCrv2, TrimCrv -> Pnext); Trim2 = TrimCrv -> Pnext -> Points; BspKnotUniformOpen(CurSize[i] + 2, 2, TrimCrv -> Pnext -> KnotVector); /* Increase index becuase of new inserted point. */ if (KP[i][0] < KP[i][2]) KP[i][2]++; else if (KP[i][0] > KP[i][2]) KP[i][0]++; for (j = 0, TmpPt = PrevPt = ConPt; j < CurSize[i] + 2 && TmpPt != NULL; j++, TmpPt = TmpPt -> Pnext) { double *R; PointType Pt; int After = 0; if (j != 0 && PrevPt -> Pnext != TmpPt) PrevPt = PrevPt -> Pnext; CP[0] = TmpPt -> Pt[0]; CP[1] = TmpPt -> Pt[1]; CP[2] = TmpPt -> Pt[2]; CP[3] = TmpPt -> Pt[3]; if (j == KP[i][0] || j == KP[i][2]) { double Prod1; if (j == KP[i][0]) { if (TmpPt -> Pnext == NULL) Prod1 = (P1[0] - PrevPt -> Pt[0]) * (P1[0] - CP[0]) + (P1[1] - PrevPt -> Pt[1]) * (P1[1] - CP[1]); else Prod1 = (TmpPt -> Pnext -> Pt[0] - CP[0]) * (P1[0] - CP[0]) + (TmpPt -> Pnext -> Pt[1] - CP[1]) * (P1[1] - CP[1]); if (Prod1 > 0) After = 1; else { R = CagdSrfEval(Srf1, P1[0], P1[1]); CagdCoerceToE3(Pt, &R, -1, Srf1 -> PType); ConvexHull1[1][j] = Pt[0]; ConvexHull1[2][j] = Pt[1]; ConvexHull1[3][j] = Pt[2]; R = CagdSrfEval(Srf2, P2[0], P2[1]); CagdCoerceToE3(Pt, &R, -1, Srf2 -> PType); ConvexHull2[1][j] = Pt[0]; ConvexHull2[2][j] = Pt[1]; ConvexHull2[3][j] = Pt[2]; Parameter[i][0] = j; Trim1[1][j] = P1[0]; Trim1[2][j] = P1[1]; Trim2[1][j] = P2[0]; Trim2[2][j] = P2[1]; j++; } } else if (j == KP[i][2]) { if (TmpPt -> Pnext == NULL) Prod1 = (P3[0] - PrevPt -> Pt[0]) * (P3[0] - CP[0]) + (P3[1] - PrevPt -> Pt[1]) * (P3[1] - CP[1]); else Prod1 = (TmpPt -> Pnext -> Pt[0] - CP[0]) * (P3[0] - CP[0]) + (TmpPt -> Pnext -> Pt[1] - CP[1]) * (P3[1] - CP[1]); if (Prod1 > 0) After = 1; else { R = CagdSrfEval(Srf1, P3[0], P3[1]); CagdCoerceToE3(Pt, &R, -1, Srf1 -> PType); ConvexHull1[1][j] = Pt[0]; ConvexHull1[2][j] = Pt[1]; ConvexHull1[3][j] = Pt[2]; R = CagdSrfEval(Srf2, P4[0], P4[1]); CagdCoerceToE3(Pt, &R, -1, Srf2 -> PType); ConvexHull2[1][j] = Pt[0]; ConvexHull2[2][j] = Pt[1]; ConvexHull2[3][j] = Pt[2]; Parameter[i][2] = j; Trim1[1][j] = P3[0]; Trim1[2][j] = P3[1]; Trim2[1][j] = P4[0]; Trim2[2][j] = P4[1]; j++; } } } R = CagdSrfEval(Srf1, CP[0], CP[1]); CagdCoerceToE3(Pt, &R, -1, Srf1 -> PType); ConvexHull1[1][j] = Pt[0]; ConvexHull1[2][j] = Pt[1]; ConvexHull1[3][j] = Pt[2]; R = CagdSrfEval(Srf2, CP[2], CP[3]); CagdCoerceToE3(Pt, &R, -1, Srf2 -> PType); ConvexHull2[1][j] = Pt[0]; ConvexHull2[2][j] = Pt[1]; ConvexHull2[3][j] = Pt[2]; Trim1[1][j] = CP[0]; Trim1[2][j] = CP[1]; Trim2[1][j] = CP[2]; Trim2[2][j] = CP[3]; if (After) { j++; if (j == (KP[i][0] + 1)) { R = CagdSrfEval(Srf1, P1[0], P1[1]); CagdCoerceToE3(Pt, &R, -1, Srf1 -> PType); ConvexHull1[1][j] = Pt[0]; ConvexHull1[2][j] = Pt[1]; ConvexHull1[3][j] = Pt[2]; R = CagdSrfEval(Srf2, P2[0], P2[1]); CagdCoerceToE3(Pt, &R, -1, Srf2 -> PType); ConvexHull2[1][j] = Pt[0]; ConvexHull2[2][j] = Pt[1]; ConvexHull2[3][j] = Pt[2]; Parameter[i][0] = j; Trim1[1][j] = P1[0]; Trim1[2][j] = P1[1]; Trim2[1][j] = P2[0]; Trim2[2][j] = P2[1]; } else if (j == (KP[i][2] + 1)) { R = CagdSrfEval(Srf1, P3[0], P3[1]); CagdCoerceToE3(Pt, &R, -1, Srf1 -> PType); ConvexHull1[1][j] = Pt[0]; ConvexHull1[2][j] = Pt[1]; ConvexHull1[3][j] = Pt[2]; R = CagdSrfEval(Srf2, P4[0], P4[1]); CagdCoerceToE3(Pt, &R, -1, Srf2 -> PType); ConvexHull2[1][j] = Pt[0]; ConvexHull2[2][j] = Pt[1]; ConvexHull2[3][j] = Pt[2]; Parameter[i][2] = j; Trim1[1][j] = P3[0]; Trim1[2][j] = P3[1]; Trim2[1][j] = P4[0]; Trim2[2][j] = P4[1]; } } } CHCrv = CHCrv -> Pnext; TrimCrv = TrimCrv -> Pnext; CagdSrfFree(Srf1); CagdSrfFree(Srf2); } CHCrv -> Pnext = NULL; TrimCrv -> Pnext = NULL; /* According to TriTangents solution, extract wanted subregion of each */ /* curve and return generated developable surfaces. */ return CombineCurves2Surfaces(Srf, CHCrvs -> Pnext, TrimCrvs -> Pnext, Parameter, Pts, TriPts); } /***************************************************************************** * DESCRIPTION: * * Cut six curves at TriTangents solutions and combine those curves * * into three developable surfaces. * * * * PARAMETERS: * * Srf: The three surfaces. * * CHCrvs: Six curves each generated by two surfaces pair. * * TrimCrv: Curves in UV domain. * * Parameter: Curve parameters for Pts. * * Pts: Six points which is on the common tangent plane to the * * three surfaces. * * TriPts: Solution points from constraints solver. * * * * RETURN VALUE: * * IPObjectStruct *: Developable surfaces forming the boundary of the * * convex hull of three surfaces. * *****************************************************************************/ static IPObjectStruct *CombineCurves2Surfaces(CagdSrfStruct *Srf, CagdCrvStruct *CHCrvs, CagdCrvStruct *TrimCrv, int Parameter[3][4], PointType Pts[6], MvarPtStruct *TriPts) { IPObjectStruct *CnvxHull, CnvxHulls, *RetVal = IPGenLISTObject(NULL); CagdCrvStruct *CHCrv, *RetCrv, RetCrvs, *RetTrim, RetTrims; int i, NP = 0, NTrim = 0, Count = 0; double *R, *Node; PointType Pt1, P[12]; RetCrv = & RetCrvs; RetTrim = & RetTrims; /* Cut curves at Tri tangents solutions. */ for (CHCrv = CHCrvs, i = 0; i < 3; CHCrv = CHCrv -> Pnext -> Pnext, i++) { int LParam, TParam; PointType Dir; CagdCrvStruct *CrvDeriv = CagdCrvDerive(CHCrv); double LP, TP, LPt, TPt, Inner = 0.0; switch(i) { case 0: LParam = (Parameter[0][0] < Parameter[0][2] ? Parameter[0][0] : Parameter[0][2]); TParam = (Parameter[0][0] < Parameter[0][2] ? Parameter[0][2] : Parameter[0][0]); Dir[0] = Pts[2][0] - Pts[0][0]; Dir[1] = Pts[2][1] - Pts[0][1]; Dir[2] = Pts[2][2] - Pts[0][2]; break; case 1: LParam = (Parameter[1][0] < Parameter[1][2] ? Parameter[1][0] : Parameter[1][2]); TParam = (Parameter[1][0] < Parameter[1][2] ? Parameter[1][2] : Parameter[1][0]); Dir[0] = Pts[0][0] - Pts[1][0]; Dir[1] = Pts[0][1] - Pts[1][1]; Dir[2] = Pts[0][2] - Pts[1][2]; break; case 2: LParam = (Parameter[2][0] < Parameter[2][2] ? Parameter[2][0] : Parameter[2][2]); TParam = (Parameter[2][0] < Parameter[2][2] ? Parameter[2][2] : Parameter[2][0]); Dir[0] = Pts[1][0] - Pts[0][0]; Dir[1] = Pts[1][1] - Pts[0][1]; Dir[2] = Pts[1][2] - Pts[0][2]; break; } /* Cut at the tri-tangents solution acording to the curve direction */ /* and to the tangent plane. */ Node = CagdCrvNodes(CHCrv); LP = Node[LParam]; TP = Node[TParam]; Node = CagdCrvNodes(TrimCrv); LPt = Node[LParam]; TPt = Node[TParam]; R = CagdCrvEval(CrvDeriv, LP); CagdCoerceToE3(Pt1, &R, -1, CrvDeriv -> PType); Inner = DOT_PROD(Pt1, Dir); R = CagdCrvEval(TrimCrv, LP); CagdCoerceToE2(Pt1, &R, -1, TrimCrv -> PType); if (Inner > 0) { RetCrv -> Pnext = CagdCrvRegionFromCrv(CHCrv, 0.0, LP); RetCrv = RetCrv -> Pnext; RetCrv -> Pnext = CagdCrvRegionFromCrv(CHCrv -> Pnext, 0.0, LP); RetCrv = RetCrv -> Pnext; RetCrv -> Pnext = CagdCrvRegionFromCrv(CHCrv, TP, 1.0); RetCrv = RetCrv -> Pnext; RetCrv -> Pnext = CagdCrvRegionFromCrv(CHCrv -> Pnext, TP, 1.0); /* For trimming curves. */ RetTrim -> Pnext = CagdCrvRegionFromCrv(TrimCrv, 0.0, LPt); RetTrim = RetTrim -> Pnext; RetTrim -> Pnext = CagdCrvRegionFromCrv(TrimCrv -> Pnext, 0.0, LPt); RetTrim = RetTrim -> Pnext; RetTrim -> Pnext = CagdCrvRegionFromCrv(TrimCrv, TPt, 1.0); RetTrim = RetTrim -> Pnext; RetTrim -> Pnext = CagdCrvRegionFromCrv(TrimCrv -> Pnext, TPt, 1.0); } else { RetCrv -> Pnext = CagdCrvRegionFromCrv(CHCrv, LP, TP); RetCrv = RetCrv -> Pnext; RetCrv -> Pnext = CagdCrvRegionFromCrv(CHCrv -> Pnext, LP, TP); RetTrim -> Pnext = CagdCrvRegionFromCrv(TrimCrv, LPt, TPt); RetTrim = RetTrim -> Pnext; RetTrim -> Pnext = CagdCrvRegionFromCrv(TrimCrv -> Pnext, LPt, TPt); } R = CagdCrvEval(CHCrv, LP); CagdCoerceToE3(P[NP++], &R, -1, CHCrv -> PType); R = CagdCrvEval(CHCrv -> Pnext, LP); CagdCoerceToE3(P[NP++], &R, -1, CHCrv -> Pnext -> PType); R = CagdCrvEval(CHCrv, TP); CagdCoerceToE3(P[NP++], &R, -1, CHCrv -> PType); R = CagdCrvEval(CHCrv -> Pnext, TP); CagdCoerceToE3(P[NP++], &R, -1, CHCrv -> Pnext -> PType); TrimCrv = TrimCrv -> Pnext -> Pnext; RetCrv = RetCrv -> Pnext; RetTrim = RetTrim -> Pnext; } /* Now, add a list of developable surfaces. */ CnvxHull = &CnvxHulls; for (RetCrv = RetCrvs.Pnext, RetTrim = RetTrims.Pnext; RetCrv != NULL; ) { CagdCrvStruct *Ruled1, *Ruled2, *TrimC1, *TrimC2; IPAttributeStruct *Attr; Ruled1 = CagdCrvCopy(RetCrv); Ruled1 -> Pnext = NULL; Ruled2 = CagdCrvCopy(RetCrv -> Pnext); Ruled2 -> Pnext = NULL; TrimC1 = CagdCrvCopy(RetTrim); TrimC1 -> Pnext = NULL; TrimC2 = CagdCrvCopy(RetTrim -> Pnext); TrimC2 -> Pnext = NULL; CnvxHull -> Pnext = IPGenSRFObject(CagdRuledSrf(Ruled1, Ruled2, 2, 2)); Attr = _AttrMallocAttribute("trim1", IP_ATTR_OBJ); Attr -> Pnext = CnvxHull -> Pnext -> Attr; CnvxHull -> Pnext -> Attr = Attr; Attr -> U.PObj = IPGenCRVObject(TrimC1); Attr -> Pnext = _AttrMallocAttribute("trim2", IP_ATTR_OBJ); Attr -> Pnext -> Pnext = NULL; Attr -> Pnext -> U.PObj = IPGenCRVObject(TrimC2); CnvxHull = CnvxHull -> Pnext; RetCrv = RetCrv -> Pnext -> Pnext; RetTrim = RetTrim -> Pnext -> Pnext; } CnvxHull -> Pnext = NULL; /* Make a list of objects. */ for (CnvxHull = CnvxHulls.Pnext; CnvxHull != NULL; CnvxHull = CnvxHull -> Pnext) IPListObjectInsert(RetVal, Count++, CnvxHull); IPListObjectInsert(RetVal, Count, NULL); CagdCrvFreeList(CHCrv); CagdCrvFreeList(TrimCrv); CagdCrvFreeList(RetCrv); CagdCrvFreeList(RetTrim); return RetVal; } /***************************************************************************** * DESCRIPTION: * * Handle duplicate solution points and unappropriate solution points. * * * * PARAMETERS: * * Srf: The three surfaces. * * TriPts: Solution points from the constraints solver. * * Tol: Tolerence to decide whether two points are same or not. * * * * RETURN VALUE: * * MvarPtStruct *: Resolved solution points. * *****************************************************************************/ static MvarPtStruct *ProcessTriSolutions(CagdSrfStruct *Srf, MvarPtStruct *TriPts, double Tol) { MvarPtStruct *Pt, *RtnPt, *RtnPts = (MvarPtStruct *)IritMalloc(sizeof(MvarPtStruct) * 1); RtnPt = RtnPts; for (Pt = TriPts; Pt != NULL; Pt = Pt -> Pnext) { CagdVecStruct *N; double N1[3], N2[3], N3[3]; if (Pt -> Dim == -1) continue; /* Surface normal should be the same direction. */ N = CagdSrfNormal(Srf, Pt -> Pt[0], Pt -> Pt[1], TRUE); N1[0] = N -> Vec[0]; N1[1] = N -> Vec[1]; N1[2] = N -> Vec[2]; N = CagdSrfNormal(Srf -> Pnext, Pt -> Pt[2], Pt -> Pt[3], TRUE); N2[0] = N -> Vec[0]; N2[1] = N -> Vec[1]; N2[2] = N -> Vec[2]; N = CagdSrfNormal(Srf -> Pnext -> Pnext, Pt -> Pt[4], Pt -> Pt[5], TRUE); N3[0] = N -> Vec[0]; N3[1] = N -> Vec[1]; N3[2] = N -> Vec[2]; if (DOT_PROD(N1, N2) > 0 && DOT_PROD(N1, N3) > 0) { MvarPtStruct *TmpPt; RtnPt -> Pnext = MvarPtCopy(Pt); RtnPt = RtnPt -> Pnext; for (TmpPt = Pt -> Pnext; TmpPt != NULL; TmpPt = TmpPt -> Pnext) { if (APX_EQ_EPS(Pt -> Pt[0], TmpPt -> Pt[0], Tol) && APX_EQ_EPS(Pt -> Pt[1], TmpPt -> Pt[1], Tol)) TmpPt -> Dim = -1; } } } return RtnPts -> Pnext; }