0026338: STL export (especially binary) needs a lot of time if selected export path is not local

Method StlAPI_Writer::Write() is reimplemented to write triangulation directly, without conversion to StlMesh_Mesh.

New DRAW command "tessellate" is added to generate rapidly triangulation of prescribed size (on surface).

Command "tricheck" is protected to deal correctly with triangulation without UV data.

New tests added: perf de bug26338_1 and _2; bugs stlvrml bug26338

Correction of testing environment

src/MeshTest/MeshTest.cxx

@@ -232,6 +232,155 @@ options:\n\
   return 0;
 }
 
+//=======================================================================
+//function : tessellate
+//purpose  : 
+//=======================================================================
+static Standard_Integer tessellate (Draw_Interpretor& /*di*/, Standard_Integer nbarg, const char** argv)
+{
+  if (nbarg != 5)
+  {
+    std::cerr << "Builds regular triangulation with specified number of triangles\n"
+                 "    Usage: tessellate result {surface|face} nbu nbv\n"
+                 "    Triangulation is put into the face with natural bounds (result);\n"
+                 "    it will have 2*nbu*nbv triangles and (nbu+1)*(nbv+1) nodes";
+    return 1;
+  }
+
+  const char *aResName = argv[1];
+  const char *aSrcName = argv[2];
+  int aNbU = Draw::Atoi (argv[3]);
+  int aNbV = Draw::Atoi (argv[4]);
+
+  if (aNbU <= 0 || aNbV <= 0)
+  {
+    std::cerr << "Error: Arguments nbu and nbv must be both greater than 0\n";
+    return 1;
+  }
+
+  Handle(Geom_Surface) aSurf = DrawTrSurf::GetSurface(aSrcName);
+  double aUMin, aUMax, aVMin, aVMax;
+  if (! aSurf.IsNull())
+  {
+    aSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
+  }
+  else
+  {
+    TopoDS_Shape aShape = DBRep::Get(aSrcName);
+    if (aShape.IsNull() || aShape.ShapeType() != TopAbs_FACE)
+    {
+      std::cerr << "Error: " << aSrcName << " is not a face\n";
+      return 1;
+    }
+    TopoDS_Face aFace = TopoDS::Face (aShape);
+    aSurf = BRep_Tool::Surface (aFace);
+    if (aSurf.IsNull())
+    {
+      std::cerr << "Error: Face " << aSrcName << " has no surface\n";
+      return 1;
+    }
+
+    BRepTools::UVBounds (aFace, aUMin, aUMax, aVMin, aVMax);
+  }
+  if (Precision::IsInfinite (aUMin) || Precision::IsInfinite (aUMax) || 
+      Precision::IsInfinite (aVMin) || Precision::IsInfinite (aVMax))
+  {
+    std::cerr << "Error: surface has infinite parametric range, aborting\n";
+    return 1;
+  }
+
+  BRepBuilderAPI_MakeFace aFaceMaker (aSurf, aUMin, aUMax, aVMin, aVMax, Precision::Confusion());
+  if (! aFaceMaker.IsDone())
+  {
+    std::cerr << "Error: cannot build face with natural bounds, aborting\n";
+    return 1;
+  }
+  TopoDS_Face aFace = aFaceMaker;
+
+  // create triangulation
+  int aNbNodes = (aNbU + 1) * (aNbV + 1);
+  int aNbTriangles = 2 * aNbU * aNbV;
+  Handle(Poly_Triangulation) aTriangulation =
+    new Poly_Triangulation (aNbNodes, aNbTriangles, Standard_False);
+
+  // fill nodes
+  TColgp_Array1OfPnt &aNodes = aTriangulation->ChangeNodes();
+  GeomAdaptor_Surface anAdSurf (aSurf);
+  double aDU = (aUMax - aUMin) / aNbU;
+  double aDV = (aVMax - aVMin) / aNbV;
+  for (int iU = 0, iShift = 1; iU <= aNbU; iU++, iShift += aNbV + 1)
+  {
+    double aU = aUMin + iU * aDU;
+    for (int iV = 0; iV <= aNbV; iV++)
+    {
+      double aV = aVMin + iV * aDV;
+      gp_Pnt aP = anAdSurf.Value (aU, aV);
+      aNodes.SetValue (iShift + iV, aP);
+    }
+  }
+
+  // fill triangles
+  Poly_Array1OfTriangle &aTriangles = aTriangulation->ChangeTriangles();
+  for (int iU = 0, iShift = 1, iTri = 0; iU < aNbU; iU++, iShift += aNbV + 1)
+  {
+    for (int iV = 0; iV < aNbV; iV++)
+    {
+      int iBase = iShift + iV;
+      Poly_Triangle aTri1 (iBase, iBase + aNbV + 2, iBase + 1);
+      Poly_Triangle aTri2 (iBase, iBase + aNbV + 1, iBase + aNbV + 2);
+      aTriangles.SetValue (++iTri, aTri1);
+      aTriangles.SetValue (++iTri, aTri2);
+    }
+  }
+
+  // put triangulation to face
+  BRep_Builder B;
+  B.UpdateFace (aFace, aTriangulation);
+
+  // fill edge polygons
+  TColStd_Array1OfInteger aUMinIso (1, aNbV + 1), aUMaxIso (1, aNbV + 1);
+  for (int iV = 0; iV <= aNbV; iV++)
+  {
+    aUMinIso.SetValue (1 + iV, 1 + iV);
+    aUMaxIso.SetValue (1 + iV, 1 + iV + aNbU * (1 + aNbV));
+  }
+  TColStd_Array1OfInteger aVMinIso (1, aNbU + 1), aVMaxIso (1, aNbU + 1);
+  for (int iU = 0; iU <= aNbU; iU++)
+  {
+    aVMinIso.SetValue (1 + iU,  1 + iU  * (1 + aNbV));
+    aVMaxIso.SetValue (1 + iU, (1 + iU) * (1 + aNbV));
+  }
+  Handle(Poly_PolygonOnTriangulation) aUMinPoly = new Poly_PolygonOnTriangulation (aUMinIso);
+  Handle(Poly_PolygonOnTriangulation) aUMaxPoly = new Poly_PolygonOnTriangulation (aUMaxIso);
+  Handle(Poly_PolygonOnTriangulation) aVMinPoly = new Poly_PolygonOnTriangulation (aVMinIso);
+  Handle(Poly_PolygonOnTriangulation) aVMaxPoly = new Poly_PolygonOnTriangulation (aVMaxIso);
+  for (TopExp_Explorer exp (aFace, TopAbs_EDGE); exp.More(); exp.Next())
+  {
+    TopoDS_Edge anEdge = TopoDS::Edge (exp.Current());
+    Standard_Real aFirst, aLast;
+    Handle(Geom2d_Curve) aC = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
+    gp_Pnt2d aPFirst = aC->Value (aFirst);
+    gp_Pnt2d aPLast  = aC->Value (aLast);
+    if (Abs (aPFirst.X() - aPLast.X()) < 0.1 * (aUMax - aUMin)) // U=const
+    {
+      if (BRep_Tool::IsClosed (anEdge, aFace))
+        B.UpdateEdge (anEdge, aUMinPoly, aUMaxPoly, aTriangulation);
+      else
+        B.UpdateEdge (anEdge, (aPFirst.X() < 0.5 * (aUMin + aUMax) ? aUMinPoly : aUMaxPoly), aTriangulation);
+    }
+    else // V=const
+    {
+      if (BRep_Tool::IsClosed (anEdge, aFace))
+        B.UpdateEdge (anEdge, aVMinPoly, aVMaxPoly, aTriangulation);
+      else
+        B.UpdateEdge (anEdge, (aPFirst.Y() < 0.5 * (aVMin + aVMax) ? aVMinPoly : aVMaxPoly), aTriangulation);
+    }
+  }
+
+  DBRep::Set (aResName, aFace);
+  return 0;
+}
+
 //=======================================================================
 //function : MemLeakTest
 //purpose  : 
@@ -1620,6 +1769,7 @@ void  MeshTest::Commands(Draw_Interpretor& theCommands)
   g = "Mesh Commands";
 
   theCommands.Add("incmesh","Builds triangular mesh for the shape, run w/o args for help",__FILE__, incrementalmesh, g);
+  theCommands.Add("tessellate","Builds triangular mesh for the surface, run w/o args for help",__FILE__, tessellate, g);
   theCommands.Add("MemLeakTest","MemLeakTest",__FILE__, MemLeakTest, g);
   theCommands.Add("fastdiscret","fastdiscret shape deflection",__FILE__, fastdiscret, g);
   theCommands.Add("mesh","mesh result Shape deflection",__FILE__, triangule, g);

src/MeshTest/MeshTest_PluginCommands.cxx

@@ -366,7 +366,6 @@ static Standard_Integer tricheck (Draw_Interpretor& di, int n, const char ** a)
       TopLoc_Location aLoc;
       Handle(Poly_Triangulation) aT = BRep_Tool::Triangulation(aFace, aLoc);
       const TColgp_Array1OfPnt& aPoints = aT->Nodes();
-      const TColgp_Array1OfPnt2d& aPoints2d = aT->UVNodes();
       const gp_Trsf& trsf = aLoc.Transformation();
 
       TColgp_Array1OfPnt pnts(1,2);
@@ -381,12 +380,16 @@ static Standard_Integer tricheck (Draw_Interpretor& di, int n, const char ** a)
         DrawTrSurf::Set (name, poly);
         DrawTrSurf::Set (name, pnts(1));
         DrawTrSurf::Set (name, pnts(2));
-        pnts2d(1) = aPoints2d(n1);
-        pnts2d(2) = aPoints2d(n2);
-        Handle(Poly_Polygon2D) poly2d = new Poly_Polygon2D (pnts2d);
-        DrawTrSurf::Set (name, poly2d);
-        DrawTrSurf::Set (name, pnts2d(1));
-        DrawTrSurf::Set (name, pnts2d(2));
+        if (aT->HasUVNodes())
+        {
+          const TColgp_Array1OfPnt2d& aPoints2d = aT->UVNodes();
+          pnts2d(1) = aPoints2d(n1);
+          pnts2d(2) = aPoints2d(n2);
+          Handle(Poly_Polygon2D) poly2d = new Poly_Polygon2D (pnts2d);
+          DrawTrSurf::Set (name, poly2d);
+          DrawTrSurf::Set (name, pnts2d(1));
+          DrawTrSurf::Set (name, pnts2d(2));
+        }
       }
       di << "\n";
     }
@@ -428,10 +431,12 @@ static Standard_Integer tricheck (Draw_Interpretor& di, int n, const char ** a)
       TopLoc_Location aLoc;
       Handle(Poly_Triangulation) aT = BRep_Tool::Triangulation(aFace, aLoc);
       const TColgp_Array1OfPnt& aPoints = aT->Nodes();
-      const TColgp_Array1OfPnt2d& aPoints2d = aT->UVNodes();
       const gp_Trsf& trsf = aLoc.Transformation();
       DrawTrSurf::Set (name, aPoints(inode).Transformed(trsf));
-      DrawTrSurf::Set (name, aPoints2d(inode));
+      if (aT->HasUVNodes())
+      {
+        DrawTrSurf::Set (name, aT->UVNodes()(inode));
+      }
 
       di << "{" << iface << " " << inode << "} ";
     }
@@ -520,7 +525,6 @@ static Standard_Integer tricheck (Draw_Interpretor& di, int n, const char ** a)
       di << "Not connected mesh inside face " << aFaceId << "\n";
 
       const TColgp_Array1OfPnt& aPoints = aT->Nodes();
-      const TColgp_Array1OfPnt2d& aPoints2d = aT->UVNodes();
       const gp_Trsf& trsf = aLoc.Transformation();
 
       TColgp_Array1OfPnt pnts(1,2);
@@ -536,12 +540,16 @@ static Standard_Integer tricheck (Draw_Interpretor& di, int n, const char ** a)
         DrawTrSurf::Set (name, poly);
         DrawTrSurf::Set (name, pnts(1));
         DrawTrSurf::Set (name, pnts(2));
-        pnts2d(1) = aPoints2d(aLink.FirstNode());
-        pnts2d(2) = aPoints2d(aLink.LastNode());
-        Handle(Poly_Polygon2D) poly2d = new Poly_Polygon2D (pnts2d);
-        DrawTrSurf::Set (name, poly2d);
-        DrawTrSurf::Set (name, pnts2d(1));
-        DrawTrSurf::Set (name, pnts2d(2));
+        if (aT->HasUVNodes())
+        {
+          const TColgp_Array1OfPnt2d& aPoints2d = aT->UVNodes();
+          pnts2d(1) = aPoints2d(aLink.FirstNode());
+          pnts2d(2) = aPoints2d(aLink.LastNode());
+          Handle(Poly_Polygon2D) poly2d = new Poly_Polygon2D (pnts2d);
+          DrawTrSurf::Set (name, poly2d);
+          DrawTrSurf::Set (name, pnts2d(1));
+          DrawTrSurf::Set (name, pnts2d(2));
+        }
       }
       di << "\n";
     }

src/StlAPI/StlAPI_ErrorStatus.hxx

@@ -23,7 +23,8 @@ enum StlAPI_ErrorStatus
 {
 StlAPI_StatusOK,
 StlAPI_MeshIsEmpty,
-StlAPI_CannotOpenFile
+StlAPI_CannotOpenFile,
+StlAPI_WriteError
 };
 
 #endif // _StlAPI_ErrorStatus_HeaderFile

src/StlAPI/StlAPI_Writer.cxx

@@ -15,11 +15,16 @@
 #include <Bnd_Box.hxx>
 #include <BRepBndLib.hxx>
 #include <OSD_Path.hxx>
+#include <OSD_OpenFile.hxx>
 #include <RWStl.hxx>
 #include <StlAPI_Writer.hxx>
 #include <StlMesh_Mesh.hxx>
 #include <StlTransfer.hxx>
-#include <TopoDS_Shape.hxx>
+#include <BRep_Tool.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopExp_Explorer.hxx>
+#include <Poly_Triangulation.hxx>
 
 StlAPI_Writer::StlAPI_Writer()
 {
@@ -32,26 +37,172 @@ Standard_Boolean& StlAPI_Writer::ASCIIMode()
   return theASCIIMode;
 }
 
-StlAPI_ErrorStatus StlAPI_Writer::Write(const TopoDS_Shape& theShape, const Standard_CString theFileName)
+// Auxiliary tools
+namespace
 {
-  OSD_Path aFile(theFileName);
-  StlTransfer::RetrieveMesh(theShape, theStlMesh);
-
-  if (theStlMesh.IsNull() || theStlMesh->IsEmpty())
-    return StlAPI_MeshIsEmpty;
-
-  // Write the built mesh
-  Standard_Boolean wasFileOpened = Standard_False;
-  if (theASCIIMode) {
-    wasFileOpened = RWStl::WriteAscii(theStlMesh, aFile);
-    }  
-  else {
-    wasFileOpened = RWStl::WriteBinary(theStlMesh, aFile);
+  // Tool to get triangles from triangulation taking into account face
+  // orientation and location
+  class TriangleAccessor
+  {
+  public:
+    TriangleAccessor (const TopoDS_Face& aFace)
+    {
+      TopLoc_Location aLoc;
+      myPoly = BRep_Tool::Triangulation (aFace, aLoc);
+      myTrsf = aLoc.Transformation();
+      myNbTriangles = (myPoly.IsNull() ? 0 : myPoly->Triangles().Length());
+      myInvert = (aFace.Orientation() == TopAbs_REVERSED);
+      if (myTrsf.IsNegative())
+        myInvert = ! myInvert;
     }
 
-  if (!wasFileOpened)
-    return StlAPI_CannotOpenFile;
+    int NbTriangles () const { return myNbTriangles; } 
 
-  return StlAPI_StatusOK;
+    // get i-th triangle and outward normal
+    void GetTriangle (int iTri, gp_Vec &theNormal, gp_Pnt &thePnt1, gp_Pnt &thePnt2, gp_Pnt &thePnt3)
+    {
+      // get positions of nodes
+      int iNode1, iNode2, iNode3;
+      myPoly->Triangles()(iTri).Get (iNode1, iNode2, iNode3);
+      thePnt1 = myPoly->Nodes()(iNode1);
+      thePnt2 = myPoly->Nodes()(myInvert ? iNode3 : iNode2);
+      thePnt3 = myPoly->Nodes()(myInvert ? iNode2 : iNode3);
+
+      // apply transormation if not identity
+      if (myTrsf.Form() != gp_Identity)
+      {
+        thePnt1.Transform (myTrsf);
+        thePnt2.Transform (myTrsf);
+        thePnt3.Transform (myTrsf);
+      }
+
+      // calculate normal
+      theNormal = (thePnt2.XYZ() - thePnt1.XYZ()) ^ (thePnt3.XYZ() - thePnt1.XYZ());
+      Standard_Real aNorm = theNormal.Magnitude();
+      if (aNorm > gp::Resolution())
+        theNormal /= aNorm;
+    }
+
+  private:
+    Handle(Poly_Triangulation) myPoly;
+    gp_Trsf myTrsf;
+    int myNbTriangles;
+    bool myInvert;
+  };
+
+  // convert to float and, on big-endian platform, to little-endian representation
+  inline float convertFloat (Standard_Real aValue)
+  {
+#ifdef OCCT_BINARY_FILE_DO_INVERSE
+    return OSD_BinaryFile::InverseShortReal ((float)aValue);
+#else
+    return (float)aValue;
+#endif
+  }
+}
+
+StlAPI_ErrorStatus StlAPI_Writer::Write(const TopoDS_Shape& theShape, const Standard_CString theFileName)
+{
+  // open file
+  FILE* aFile = OSD_OpenFile (theFileName, "wb");
+  if (!aFile)
+    return StlAPI_CannotOpenFile;
+
+  // write
+  if (theASCIIMode)
+  {
+    // header 
+    Fprintf (aFile, "solid shape, STL ascii file, created with Open CASCADE Technology\n");
+
+    // facets
+    for (TopExp_Explorer exp (theShape, TopAbs_FACE); exp.More(); exp.Next())
+    {
+      TriangleAccessor aTool (TopoDS::Face (exp.Current()));
+      for (int iTri = 1; iTri <= aTool.NbTriangles(); iTri++)
+      {
+        gp_Vec aNorm;
+        gp_Pnt aPnt1, aPnt2, aPnt3;
+        aTool.GetTriangle (iTri, aNorm, aPnt1, aPnt2, aPnt3);
+
+        Fprintf (aFile,
+          " facet normal %12e %12e %12e\n"
+          "   outer loop\n"
+          "     vertex %12e %12e %12e\n"
+          "     vertex %12e %12e %12e\n"
+          "     vertex %12e %12e %12e\n"
+          "   endloop\n"
+          " endfacet\n",
+          aNorm.X(), aNorm.Y(), aNorm.Z(),
+          aPnt1.X(), aPnt1.Y(), aPnt1.Z(),
+          aPnt2.X(), aPnt2.Y(), aPnt2.Z(),
+          aPnt3.X(), aPnt3.Y(), aPnt3.Z());
+      }
+    }
+
+    // footer
+    Fprintf (aFile, "endsolid shape\n");
+  }
+  else
+  {
+    // header block (meaningless 80 bytes)
+    Fprintf (aFile, "%-80.80s", "STL binary file, created with Open CASCADE Technology");
+
+    // number of facets
+    int32_t aNbTri = 0;
+    for (TopExp_Explorer exp (theShape, TopAbs_FACE); exp.More(); exp.Next())
+    {
+      TopLoc_Location aLoc;
+      Handle(Poly_Triangulation) aPoly =
+        BRep_Tool::Triangulation (TopoDS::Face (exp.Current()), aLoc);
+      if (! aPoly.IsNull())
+        aNbTri += aPoly->NbTriangles();
+    }
+    // suppose that number of triangles must be little endian...
+#ifdef OCCT_BINARY_FILE_DO_INVERSE
+    aNbTri = OSD_BinaryFile::InverseInteger (aNbTri);
+#endif
+    fwrite (&aNbTri, sizeof(int32_t), 1, aFile);
+
+    // facets
+    struct Facet {
+      float nx, ny, nz;
+      float x1, y1, z1;
+      float x2, y2, z2;
+      float x3, y3, z3;
+      uint16_t dummy;
+    } f;
+    f.dummy = 0;
+    for (TopExp_Explorer exp (theShape, TopAbs_FACE); exp.More(); exp.Next())
+    {
+      TriangleAccessor aTool (TopoDS::Face (exp.Current()));
+      for (int iTri = 1; iTri <= aTool.NbTriangles(); iTri++)
+      {
+        gp_Vec aNorm;
+        gp_Pnt aPnt1, aPnt2, aPnt3;
+        aTool.GetTriangle (iTri, aNorm, aPnt1, aPnt2, aPnt3);
+
+        f.nx = convertFloat (aNorm.X());
+        f.ny = convertFloat (aNorm.Y());
+        f.nz = convertFloat (aNorm.Z());
+
+        f.x1 = convertFloat (aPnt1.X());
+        f.y1 = convertFloat (aPnt1.Y());
+        f.z1 = convertFloat (aPnt1.Z());
+
+        f.x2 = convertFloat (aPnt2.X());
+        f.y2 = convertFloat (aPnt2.Y());
+        f.z2 = convertFloat (aPnt2.Z());
+
+        f.x3 = convertFloat (aPnt3.X());
+        f.y3 = convertFloat (aPnt3.Y());
+        f.z3 = convertFloat (aPnt3.Z());
+
+        fwrite (&f, 50 /* 50 bytes per facet */, 1, aFile);
+      }
+    }
+  }
+
+  fclose (aFile);
+  return ferror(aFile) ? StlAPI_WriteError : StlAPI_StatusOK;
 }
 

tests/bugs/stlvrml/bug26338

@@ -0,0 +1,28 @@
+puts "========"
+puts "0026338: STL export (especially binary) needs a lot of time if selected export path is not local"
+puts "========"
+puts ""
+
+pload MODELING XSDRAW
+
+# check that export to STL correctly takes into account shape locations
+box b1 5 5 5 
+box b2 5 5 5 
+ttranslate b2 10 10 10 
+compound b1 b2 comp
+incmesh comp 1.
+
+# write to binary STL
+writestl comp $imagedir/${casename}.stl 1
+
+# load STL
+readstl result $imagedir/${casename}.stl
+
+# check that bounding box is 
+set bnd [boundingstr result]
+checkreal "XMin" [lindex $bnd 0]  0. 1e-5 0.
+checkreal "YMin" [lindex $bnd 1]  0. 1e-5 0.
+checkreal "ZMin" [lindex $bnd 2]  0. 1e-5 0.
+checkreal "XMax" [lindex $bnd 3] 15. 1e-5 0.
+checkreal "YMax" [lindex $bnd 4] 15. 1e-5 0.
+checkreal "ZMax" [lindex $bnd 5] 15. 1e-5 0.

tests/perf/de/bug26338_1

@@ -0,0 +1,16 @@
+puts "========"
+puts "0026338: STL export (especially binary) needs a lot of time if selected export path is not local"
+puts "========"
+puts ""
+
+pload MODELING XSDRAW
+
+# make sphere triangulated with 2M triangles
+sphere s 10
+tessellate result s 1000 1000
+trinfo result
+
+# write to binary STL
+chrono s reset; chrono s start
+writestl result $imagedir/${casename}-binary.stl 1
+chrono s stop; chrono s show

tests/perf/de/bug26338_2

@@ -0,0 +1,16 @@
+puts "========"
+puts "0026338: STL export (especially binary) needs a lot of time if selected export path is not local"
+puts "========"
+puts ""
+
+pload MODELING XSDRAW
+
+# make sphere triangulated with 2M triangles
+sphere s 10
+tessellate result s 1000 1000
+trinfo result
+
+# write to ascii STL
+chrono s reset; chrono s start
+writestl result $imagedir/${casename}-ascii.stl 0
+chrono s stop; chrono s show

tests/perf/grids.list

@@ -2,3 +2,4 @@
 002 ncollection
 003 bspline
 004 fclasses
+005 de