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 }