0027234: Code duplication: Convert_CompBezierCurvesToBSplineCurve* in ShapeConstruct
[occt.git] / src / ShapeConstruct / ShapeConstruct_MakeTriangulation.cxx
1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14
15 #include <BRep_Builder.hxx>
16 #include <BRep_Tool.hxx>
17 #include <BRepBuilderAPI_MakeEdge.hxx>
18 #include <BRepBuilderAPI_MakePolygon.hxx>
19 #include <Geom_Plane.hxx>
20 #include <gp_Pln.hxx>
21 #include <gp_Vec.hxx>
22 #include <Precision.hxx>
23 #include <ShapeAnalysis_Curve.hxx>
24 #include <ShapeAnalysis_Edge.hxx>
25 #include <ShapeAnalysis_Wire.hxx>
26 #include <ShapeConstruct_MakeTriangulation.hxx>
27 #include <TColgp_HArray1OfPnt.hxx>
28 #include <TColgp_SequenceOfPnt.hxx>
29 #include <TColStd_Array1OfInteger.hxx>
30 #include <TColStd_SequenceOfInteger.hxx>
31 #include <TopoDS.hxx>
32 #include <TopoDS_Edge.hxx>
33 #include <TopoDS_Face.hxx>
34 #include <TopoDS_Iterator.hxx>
35 #include <TopoDS_Shell.hxx>
36 #include <TopoDS_Wire.hxx>
37 #include <TopTools_HSequenceOfShape.hxx>
38
39 //=======================================================================
40 //function : IsRightContour (static)
41 //purpose  : 
42 //=======================================================================
43 Standard_Boolean IsRightContour (const TColgp_SequenceOfPnt& pts, const Standard_Real prec)
44 {
45   Standard_Integer len = pts.Length();
46   if (len < 4) return Standard_True;
47
48   TColgp_Array1OfPnt thePts(1,len);
49   Standard_Integer i; // svv Jan11 2000 : porting on DEC
50   for (i = 1; i <= len; i++) thePts(i) = pts(i);
51   
52   gp_XYZ Norm(0,0,0);
53   if (ShapeAnalysis_Curve::IsPlanar(thePts,Norm,prec)) {
54
55     // Make polygonal wire from points
56     BRepBuilderAPI_MakePolygon mkPoly;
57     for (i = 1; i <= len; i++) mkPoly.Add(thePts(i));
58     mkPoly.Close();
59     mkPoly.Build();
60     if (mkPoly.IsDone()) {
61
62       // Create mean plane
63       gp_XYZ center(0,0,0);
64       for (i = 1; i <= len; i++) center += thePts(i).XYZ();
65       center /= len;
66       gp_Pnt pc(center); gp_Dir pd(Norm); gp_Pln pl(pc,pd);
67       Handle(Geom_Plane) thePlane = new Geom_Plane(pl);
68
69       BRep_Builder B;
70       TopoDS_Face theFace;
71       B.MakeFace(theFace,thePlane,Precision::Confusion());
72       TopoDS_Wire theWire = mkPoly.Wire();
73       B.Add(theFace,theWire);
74       Handle(ShapeAnalysis_Wire) saw = new ShapeAnalysis_Wire(theWire,theFace,prec);
75       return !saw->CheckSelfIntersection();
76     }
77   }
78
79   return Standard_False;
80 }
81
82 //=======================================================================
83 //function : MeanNormal (static)
84 //purpose  : 
85 //=======================================================================
86
87  gp_Vec MeanNormal (const TColgp_Array1OfPnt& pts)
88 {
89   Standard_Integer len = pts.Length();
90   if (len < 3) return gp_Vec(0,0,0);
91
92   Standard_Integer i;
93   gp_XYZ theCenter(0,0,0);
94   for (i = 1; i <= len; i++) theCenter += pts(i).XYZ();
95   theCenter /= len;
96
97   gp_XYZ theNorm(0,0,0);
98   for (i = 1; i <= len; i++) {
99     gp_Vec v1(pts(i).XYZ() - theCenter);
100     gp_Vec v2(pts((i==len)? 1 : i+1).XYZ() - theCenter);
101     theNorm += (v1^v2).XYZ();
102   }
103
104   return gp_Vec(theNorm).Normalized();
105 }
106
107 //=======================================================================
108 //function : ShapeConstruct_MakeTriangulation
109 //purpose  : 
110 //=======================================================================
111
112 ShapeConstruct_MakeTriangulation::ShapeConstruct_MakeTriangulation (const TColgp_Array1OfPnt& pnts,
113                                                                     const Standard_Real prec)
114 {
115   myPrecision = (prec > 0.0)? prec : Precision::Confusion();
116   // Make polygonal wire from points
117   BRepBuilderAPI_MakePolygon mkPoly;
118   for (Standard_Integer i = pnts.Lower(); i <= pnts.Upper(); i++) mkPoly.Add(pnts(i));
119   mkPoly.Close();
120   mkPoly.Build();
121   if (mkPoly.IsDone()) {
122     myWire = mkPoly.Wire();
123     Build();
124   }
125 }
126
127 //=======================================================================
128 //function : ShapeConstruct_MakeTriangulation
129 //purpose  : 
130 //=======================================================================
131
132 ShapeConstruct_MakeTriangulation::ShapeConstruct_MakeTriangulation (const TopoDS_Wire& wire,
133                                                                     const Standard_Real prec)
134 {
135   myPrecision = (prec > 0.0)? prec : Precision::Confusion();
136   myWire = wire;
137   Build();
138 }
139
140 //=======================================================================
141 //function : Build
142 //purpose  : 
143 //=======================================================================
144
145  void ShapeConstruct_MakeTriangulation::Build()
146 {
147   if (myShape.IsNull()) {
148     // Triangulate polygonal wire
149     if (!myWire.IsNull()) Triangulate(myWire);
150   }
151 }
152
153 //=======================================================================
154 //function : IsDone
155 //purpose  : 
156 //=======================================================================
157
158  Standard_Boolean ShapeConstruct_MakeTriangulation::IsDone() const
159 {
160   return (!myShape.IsNull());
161 }
162
163 //=======================================================================
164 //function : Triangulate
165 //purpose  : 
166 //=======================================================================
167
168  void ShapeConstruct_MakeTriangulation::Triangulate (const TopoDS_Wire& wire)
169 {
170   // Fill sequence of edges
171   Handle(TopTools_HSequenceOfShape) theEdges = new TopTools_HSequenceOfShape;
172   for (TopoDS_Iterator ite(wire); ite.More(); ite.Next()) theEdges->Append(ite.Value());
173   Standard_Integer NbEdges = theEdges->Length();
174
175   if (NbEdges < 4) AddFacet(wire);
176   else {
177
178     // Fill array of end points
179     ShapeAnalysis_Edge sae;
180     Handle(TColgp_HArray1OfPnt) theAPnt = new TColgp_HArray1OfPnt(1,NbEdges);
181     for (Standard_Integer i = 1; i <= NbEdges; i++)
182       theAPnt->SetValue(i,BRep_Tool::Pnt(sae.FirstVertex(TopoDS::Edge(theEdges->Value(i)))));
183
184     TopoDS_Wire theNewWire;
185
186     // Check planarity on array of points
187     gp_XYZ Norm(0,0,0);
188     if (ShapeAnalysis_Curve::IsPlanar(theAPnt->Array1(),Norm,myPrecision)) AddFacet(wire);
189     else {
190
191       // Current contour is not planar - triangulate
192       TColStd_SequenceOfInteger theSegments;
193       Standard_Integer cindex = 1, seqFlag = 1;
194       TColStd_Array1OfInteger llimits(1,2), rlimits(1,2);
195       llimits.Init(NbEdges+1); rlimits.Init(0);
196
197       // Find all "good" planar segments
198       while (seqFlag) {
199
200         // Set up limits for current sequence
201         Standard_Integer llimit = llimits(seqFlag);
202         Standard_Integer rlimit = rlimits(seqFlag);
203         if (rlimit <= llimit) rlimit += NbEdges;
204
205         // Check stop condition
206         if ((rlimit - llimit) > (NbEdges - 2)) break;
207
208         TColgp_SequenceOfPnt theSPnt;
209
210         // Add 3 points of minimal facet
211         Standard_Integer lindex = (cindex==1? NbEdges : cindex-1);
212         Standard_Integer rindex = (cindex==NbEdges? 1 : cindex+1);
213         theSPnt.Append(theAPnt->Value(lindex));
214         theSPnt.Append(theAPnt->Value(cindex));
215         theSPnt.Append(theAPnt->Value(rindex));
216         
217         // Try prepending points
218         Standard_Boolean isGood = Standard_True;
219         while (isGood) {
220           cindex = (lindex==1? NbEdges : lindex-1);
221           if (cindex == rindex) { AddFacet(wire); return; }
222           Standard_Integer cid = cindex;
223           if (rlimit > NbEdges && cid <= llimit) cid += NbEdges;
224           if (cid > llimit && cid < rlimit) isGood = Standard_False;
225           if (isGood) {
226             theSPnt.Prepend(theAPnt->Value(cindex));
227             isGood = IsRightContour(theSPnt,myPrecision);
228             if (isGood) lindex = cindex;
229             else theSPnt.Remove(1);
230           }
231         }
232
233         // Try appending points
234         isGood = Standard_True;
235         while (isGood) {
236           cindex = (rindex==NbEdges? 1 : rindex+1);
237           if (cindex == lindex) { AddFacet(wire); return; }
238           Standard_Integer cid = cindex;
239           if (rlimit > NbEdges && cid <= llimit) cid += NbEdges;
240           if (cid > llimit && cid < rlimit) isGood = Standard_False;
241           if (isGood) {
242             theSPnt.Append(theAPnt->Value(cindex));
243             isGood = IsRightContour(theSPnt,myPrecision);
244             if (isGood) rindex = cindex;
245             else theSPnt.Remove(theSPnt.Length());
246           }
247         }
248
249         // Record new planar segment
250         theSegments.Append(lindex);
251         theSegments.Append(rindex);
252
253         // Set up new limits
254         cindex = rindex;
255         rlimits(seqFlag) = rindex;
256         if (llimit > NbEdges) llimits(seqFlag) = lindex;
257
258         // Set up sequence index
259         seqFlag = (seqFlag == 1)? 2 : 1;
260       }
261
262       cindex = theSegments.Length();
263       if (cindex%4 == 0) {
264         gp_Vec theBaseNorm = MeanNormal(theAPnt->Array1());
265         // Identify sequence of matching facets
266         Standard_Integer len = cindex/2, lindex, rindex, seqlen, j;
267         Standard_Real theC, theC1 = 0.0, theC2 = 0.0;
268         Standard_Integer ii; // svv Jan11 2000 : porting on DEC
269         // skl : change "i" to "ii"
270         for (ii = 0; ii < len; ii++) { 
271           // Compute normal collinearity for facet
272           lindex = theSegments(ii*2+1);
273           rindex = theSegments(ii*2+2);
274           seqlen = ((rindex > lindex)? rindex - lindex + 1 : NbEdges - lindex + rindex + 1);
275           TColgp_Array1OfPnt theArr(1,seqlen);
276           for (j = 1; j <= seqlen; j++) {
277             theArr(j) = theAPnt->Value(lindex);
278             lindex++;
279             if (lindex > NbEdges) lindex = 1;
280           }
281           theC = theBaseNorm*MeanNormal(theArr);
282           if (ii%2) theC2 += theC;      else theC1 += theC;
283         }
284         Standard_Integer best1 = ((theC1 > theC2)? 0 : 2);
285
286         // Add "best" planar facets
287         BRep_Builder B;
288         TopoDS_Wire theFacetWire;
289         TopoDS_Edge theLEdge, theSLEdge;
290         len = cindex/4;
291         // Check for special case - 1 shared line
292         Standard_Boolean special = (len == 2 &&
293                                     theSegments(best1+1) == theSegments(4+best1+2) &&
294                                     theSegments(best1+2) == theSegments(4+best1+1));
295         Standard_Integer pindex = theSegments((len-1)*4+best1+2);
296         for (ii = 0; ii < len; ii++) {
297           lindex = theSegments(ii*4+best1+1);
298           rindex = theSegments(ii*4+best1+2);
299           if (special && !theSLEdge.IsNull()) {
300             TopoDS_Shape aLocalShape = theSLEdge.Reversed();
301             theLEdge = TopoDS::Edge(aLocalShape);       
302             //      theLEdge = TopoDS::Edge(theSLEdge.Reversed());
303           }
304           else {
305             TopoDS_Vertex v1 = sae.FirstVertex(TopoDS::Edge(theEdges->Value(lindex)));
306             TopoDS_Vertex v2 = sae.FirstVertex(TopoDS::Edge(theEdges->Value(rindex)));
307             v1.Orientation(TopAbs_FORWARD);
308             v2.Orientation(TopAbs_REVERSED);
309             theLEdge = BRepBuilderAPI_MakeEdge(v1,v2);
310             theSLEdge = theLEdge;
311           }
312           // Create new facet wire
313           B.MakeWire(theFacetWire);
314           B.Add(theFacetWire,theLEdge.Reversed());
315           Standard_Integer cur_edge = lindex;
316           while (cur_edge != rindex) {
317             TopoDS_Edge theNEdge = TopoDS::Edge(theEdges->Value(cur_edge));
318             B.Add(theFacetWire,theNEdge);
319             if (cur_edge == NbEdges) cur_edge = 1;
320             else cur_edge++;
321           }
322           AddFacet(theFacetWire);
323           if (!special) {
324             // Add edges to the new wire
325             if (theNewWire.IsNull()) B.MakeWire(theNewWire);
326             cur_edge = pindex;
327             while (cur_edge != lindex) {
328               TopoDS_Edge theNEdge = TopoDS::Edge(theEdges->Value(cur_edge));
329               B.Add(theNewWire,theNEdge);
330               if (cur_edge == NbEdges) cur_edge = 1;
331               else cur_edge++;
332             }
333             B.Add(theNewWire,theLEdge);
334             pindex = rindex;
335           }
336         }
337       }
338
339       // Clear storage variables
340       theEdges.Nullify();
341       theAPnt.Nullify();
342
343       // Triangulate the rest of edges
344       if (!theNewWire.IsNull()) Triangulate(theNewWire);
345     }
346   }
347 }
348
349 //=======================================================================
350 //function : AddFacet
351 //purpose  : 
352 //=======================================================================
353
354  void ShapeConstruct_MakeTriangulation::AddFacet (const TopoDS_Wire& wire)
355 {
356   if (wire.IsNull()) return;
357
358   // Fill sequence of points
359   ShapeAnalysis_Edge sae;
360   TColgp_SequenceOfPnt pts;
361   for (TopoDS_Iterator ite(wire); ite.More(); ite.Next())
362     pts.Append(BRep_Tool::Pnt(sae.FirstVertex(TopoDS::Edge(ite.Value()))));
363   Standard_Integer i, nbp = pts.Length();
364   if (nbp < 3) return;
365
366   // Create mean plane
367   Standard_Real cMod, maxMod = 0.0;
368   gp_XYZ maxVec, Normal(0,0,0);
369   for (i = 1; i <= nbp; i++) {
370     gp_XYZ vb(pts(i).XYZ());
371     gp_XYZ v1(pts(i==nbp? 1 : i+1).XYZ()-vb);
372     cMod = v1.SquareModulus();
373     if (cMod == 0.0) continue;
374     else if (cMod > maxMod) { maxMod = cMod; maxVec = v1; }
375     gp_XYZ v2(pts(i==1? nbp : i-1).XYZ()-vb);
376     cMod = v2.SquareModulus();
377     if (cMod == 0.0) continue;
378     else if (cMod > maxMod) { maxMod = cMod; maxVec = v2; }
379     Normal += gp_XYZ(v1^v2);
380   }
381   if (Normal.SquareModulus() == 0.0) {
382     if (maxMod == 0.0)
383       Normal = gp_XYZ(0,0,1);
384     else if (maxVec.X() != 0.0)
385       Normal = gp_XYZ(-maxVec.Y()/maxVec.X(),1,0);
386     else if (maxVec.Y() != 0.0)
387       Normal = gp_XYZ(0,-maxVec.Z()/maxVec.Y(),1);
388     else
389       Normal = gp_XYZ(1,0,0);
390   }
391
392   gp_Pln Pln(pts(1),gp_Dir(Normal));
393   Handle(Geom_Plane) thePlane = new Geom_Plane(Pln);
394
395   // Mean plane created - build face
396   BRep_Builder B;
397   TopoDS_Face theFace;
398   B.MakeFace(theFace,thePlane,Precision::Confusion());
399   B.Add(theFace,wire);
400
401   // Add new face to the shell
402   if (myShape.IsNull()) myShape = theFace;
403   else {
404     if (myShape.ShapeType() == TopAbs_FACE) {
405       TopoDS_Face firstFace = TopoDS::Face(myShape);
406       TopoDS_Shell theShell;
407       B.MakeShell(theShell);
408       myShape = theShell;
409       B.Add(myShape,firstFace);
410     }
411     B.Add(myShape,theFace);
412   }
413 }