CSG2d optimizations (in-place operators, search tree)

libsrc/geom2d/csg2d.cpp

@@ -673,6 +673,7 @@ void AddIntersectionPoint(Edge edgeP, Edge edgeQ, IntersectionType i, double alp
 
 void ComputeIntersections(Solid2d & sp, Solid2d & sq)
 {
+  static Timer tall("ComputeIntersections"); RegionTimer rtall(tall);
   auto & PP = sp.polys;
   auto & QQ = sq.polys;
 
@@ -1246,6 +1247,7 @@ void CleanUpResult(Solid2d & sr)
 
 void RemoveDuplicates(Solid2d & sr)
 {
+  static Timer tall("RemoveDuplicates"); RegionTimer rtall(tall);
   for(auto & poly : sr.polys)
   {
     if(poly.first==nullptr) continue;
@@ -1308,12 +1310,142 @@ void AddIntersectionPoints ( Solid2d & s1, Solid2d & s2 )
   RemoveDuplicates(s2);
 }
 
-Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect)
+Solid2d ClipSolids ( const Solid2d & s1, const Solid2d & s2, char op)
 {
+  return ClipSolids(Solid2d{s1}, Solid2d{s2}, op);
+}
+
+Solid2d ClipSolids ( const Solid2d & s1, Solid2d && s2, char op)
+{
+  return ClipSolids(Solid2d{s1}, std::move(s2), op);
+}
+
+Solid2d ClipSolids ( Solid2d && s1, const Solid2d & s2, char op)
+{
+  return ClipSolids(std::move(s1), Solid2d{s2}, op);
+}
+
+Solid2d ClipSolids ( Solid2d && s1, Solid2d && s2, char op)
+{
+  static Timer tall("ClipSolids"); RegionTimer rtall(tall);
+  static Timer t0("copy");
+  static Timer t02("tree");
+  static Timer t03("search intersections");
+  static Timer t01("prepare");
   static Timer t1("intersection");
   static Timer t2("label");
   static Timer t3("cut");
   static Timer t4("cleanup");
+  static Timer t6("trivial union");
+
+  bool intersect = (op=='*' || op=='-');
+
+  Solid2d res;
+  res.name = s1.name;
+
+  t0.Start();
+  // Try to quickly handle parts of both solids that cannot intersect with the other one
+  int n1 = s1.polys.Size();
+  int n2 = s2.polys.Size();
+  Array<Loop> res_polys(n1+n2);
+  res_polys.SetSize(0);
+
+  t02.Start();
+  netgen::Box<2> s1_box(netgen::Box<2>::EMPTY_BOX);
+  netgen::BoxTree <2, int> tree1(s1.GetBoundingBox());
+
+  for(auto li : IntRange(n1))
+  {
+    auto box = s1.polys[li].GetBoundingBox();
+    s1_box.Add(box.PMin());
+    s1_box.Add(box.PMax());
+    tree1.Insert(box, li);
+  }
+
+  netgen::Box<2> s2_box(netgen::Box<2>::EMPTY_BOX);
+  netgen::BoxTree <2, int> tree2(s2.GetBoundingBox());
+
+  for(auto li : IntRange(n2))
+  {
+    auto box = s2.polys[li].GetBoundingBox();
+    s2_box.Add(box.PMin());
+    s2_box.Add(box.PMax());
+    tree2.Insert(box, li);
+  }
+  t02.Stop();
+
+  t03.Start();
+
+  for(auto li : IntRange(n1))
+  {
+    bool have_intersections = false;
+    auto & poly = s1.polys[li];
+    auto box = poly.GetBoundingBox();
+    tree2.GetFirstIntersecting(box.PMin(), box.PMax(), [&] (int li2)
+        {
+          return have_intersections = true;
+        });
+    if(!have_intersections)
+    {
+      if(op=='+' || op=='-')
+        res_polys.Append(std::move(poly));
+      else
+        poly.Clear();
+    }
+  }
+  t03.Stop();
+
+  for(auto li: IntRange(n1))
+    while(s1.polys.Size()>li && s1.polys[li].Size()==0)
+      s1.polys.DeleteElement(li);
+
+  t03.Start();
+  for(auto li : IntRange(n2))
+  {
+    bool have_intersections = false;
+    auto & poly = s2.polys[li];
+    auto box = poly.GetBoundingBox();
+    tree1.GetFirstIntersecting(box.PMin(), box.PMax(), [&] (int li2)
+        {
+          return have_intersections = true;
+        });
+    if(!have_intersections)
+    {
+      if(op=='+')
+        res_polys.Append(std::move(poly));
+      else
+        poly.Clear();
+    }
+  }
+  t03.Stop();
+
+  for(auto li: IntRange(n2))
+    while(s2.polys.Size()>li && s2.polys[li].Size()==0)
+      s2.polys.DeleteElement(li);
+
+  t0.Stop();
+
+  if(s1.polys.Size()==0 || s2.polys.Size()==0)
+  {
+    res.polys = std::move(res_polys);
+    return res;
+  }
+  RegionTracer rt_(0, tall, s1.polys.Size()*10000 + s2.polys.Size());
+
+  t01.Start();
+
+  if(op=='-')
+  {
+    // take complement of s2 by adding loop around everything
+    auto box = s1_box;
+    box.Add(s2_box.PMin());
+    box.Add(s2_box.PMax());
+    box.Increase(2);
+    auto pmin = box.PMin();
+    auto pmax = box.PMax();
+    s2.Append(RectanglePoly(pmin[0], pmax[0], pmin[1], pmax[1], "JUST_FOR_CLIPPING"));
+  }
+
 
   for(auto & poly : s1.polys)
     for(auto v : poly.Vertices(ALL))
@@ -1337,8 +1469,7 @@ Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect)
       v->enex = NEITHER;
     }
 
-  Solid2d res;
+  t01.Stop();
-  res.name = s1.name;
 
   t1.Start();
   ComputeIntersections(s1, s2);
@@ -1357,7 +1488,9 @@ Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect)
   RemoveDuplicates(res);
   t4.Stop();
 
-  return res;
+  res.polys.Append(std::move(res_polys));
+
+  return std::move(res);
 }
 
 bool Loop :: IsInside( Point<2> r ) const
@@ -1396,6 +1529,16 @@ bool Loop :: IsInside( Point<2> r ) const
   return ( (w % 2) != 0 );
 }
 
+netgen::Box<2> Loop :: GetBoundingBox() const
+{
+  // static Timer tall("Loop::GetBoundingBox"); RegionTimer rtall(tall);
+  netgen::Box<2> box(netgen::Box<2>::EMPTY_BOX);
+  for(auto v : Vertices(ALL))
+    box.Add(*v);
+  return box;
+}
+
+
 Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_, string bc)
     : name(name_)
 {
@@ -1422,47 +1565,38 @@ Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, strin
 
 Solid2d Solid2d :: operator+(const Solid2d & other) const
 {
-  if(polys.Size()==0)
+  static Timer t("Solid2d::operator+"); RegionTimer rt(t);
-    return other;
+  return ClipSolids(*this, other, '+');
-
-  auto res = ClipSolids(*this, other, false);
-  res.name = name;
-  return res;
 }
 
 Solid2d Solid2d :: operator*(const Solid2d & other) const
 {
-  auto res = ClipSolids(*this, other, true);
+  static Timer t("Solid2d::operator*"); RegionTimer rt(t);
-  res.name = name;
+  return ClipSolids(*this, other, '*');
-  return res;
 }
 
-Solid2d Solid2d :: operator-(const Solid2d & other_) const
+Solid2d Solid2d :: operator-(const Solid2d & other) const
 {
-  // TODO: Check dimensions of solids with bounding box
+  static Timer t("Solid2d::operator-"); RegionTimer rt(t);
-  Solid2d other = other_;
+  return ClipSolids(*this, other, '-');
-  other.Append(RectanglePoly(-1e8, 1e8, -1e8, 1e8, "JUST_FOR_CLIPPING"));
-  auto res = ClipSolids(*this, other);
-
-  res.name = name;
-  return res;
 }
 
 Solid2d & Solid2d :: operator+=(const Solid2d & other)
 {
-  *this = *this + other;
+  static Timer t("Solid2d::operator+="); RegionTimer rt(t);
+  *this = ClipSolids(std::move(*this), other, '+');
   return *this;
 }
 
 Solid2d & Solid2d :: operator*=(const Solid2d & other)
 {
-  *this = *this * other;
+  *this = ClipSolids(std::move(*this), other, '*');
   return *this;
 }
 
 Solid2d & Solid2d :: operator-=(const Solid2d & other)
 {
-  *this = *this - other;
+  *this = ClipSolids(std::move(*this), other, '-');
   return *this;
 }
 
@@ -1504,6 +1638,42 @@ bool Solid2d :: IsInside( Point<2> r ) const
   return ( (w % 2) != 0 );
 }
 
+bool Loop :: IsLeftInside( const Vertex & p0 )
+{
+  auto & p1 = *p0.next;
+  if(p0.spline)
+  {
+    auto s = *p0.spline;
+    auto v = s.GetTangent(0.5);
+    auto n = Vec<2>{-v[1], v[0]};
+    auto q = s.GetPoint(0.5) + 1e-6*n;
+    return IsInside(q);
+  }
+  auto v = p1-p0;
+  auto n = Vec<2>{-v[1], v[0]};
+  auto q = p0 + 0.5*v + 1e-6*n;
+  
+  return IsInside(q);
+}
+
+bool Loop :: IsRightInside( const Vertex & p0 )
+{
+  auto & p1 = *p0.next;
+  if(p0.spline)
+  {
+    auto s = *p0.spline;
+    auto v = s.GetTangent(0.5);
+    auto n = Vec<2>{v[1], -v[0]};
+    auto q = s.GetPoint(0.5) + 1e-6*n;
+    return IsInside(q);
+  }
+
+  auto v = p1-p0;
+  auto n = Vec<2>{v[1], -v[0]};
+  auto q = p0 + 0.5*v + 1e-6*n;
+  return IsInside(q);
+}
+
 bool Solid2d :: IsLeftInside( const Vertex & p0 )
 {
   auto & p1 = *p0.next;
@@ -1540,20 +1710,27 @@ bool Solid2d :: IsRightInside( const Vertex & p0 )
   return IsInside(q);
 }
 
-netgen::Box<2> Solid2d :: GetBoundingBox()
+netgen::Box<2> Solid2d :: GetBoundingBox() const
 {
+  static Timer tall("Solid2d::GetBoundingBox"); RegionTimer rtall(tall);
   netgen::Box<2> box(netgen::Box<2>::EMPTY_BOX);
   for(auto & poly : polys)
-      for(auto v : poly.Vertices(ALL))
+  {
-          box.Add(*v);
+    auto pbox = poly.GetBoundingBox();
+    box.Add(pbox.PMin());
+    box.Add(pbox.PMax());
+  }
 
   return box;
 }
 
 shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
 {
-  static Timer t_intersections("CSG2d - AddIntersections()");
   static Timer tall("CSG2d - GenerateSplineGeometry()");
+  static Timer t_points("add points");
+  static Timer t_segments_map("build segments map");
+  static Timer t_segments("add segments");
+  static Timer t_intersections("add intersections");
   RegionTimer rt(tall);
 
   struct Seg
@@ -1618,6 +1795,7 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
     return res;
   };
 
+  t_points.Start();
   auto insertPoint = [&](Point<2> p )
   {
     int pi = getPoint(p);
@@ -1640,9 +1818,11 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
         if(v->spline)
           insertPoint(v->spline->TangentPoint());
       }
+  t_points.Stop();
 
 
   // Generate segments from polygon edges and find left/right domain of each segment
+  t_segments_map.Start();
   int dom = 0;
   int bc = 1;
   for(auto & s : solids)
@@ -1651,6 +1831,10 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
     bool is_solid_degenerated = true; // Don't create new domain for degenerated solids
     for(auto & poly : s.polys)
     {
+      bool first = true;
+      bool is_poly_left_inside = false;
+      bool is_poly_right_inside = false;
+
       for(auto v : poly.Vertices(ALL))
       {
         auto & p0 = *v;
@@ -1675,8 +1859,15 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
           Swap(pi1,pi0);
         }
 
-        auto li = s.IsLeftInside(p0);
+        if(first)
-        auto ri = s.IsRightInside(p0);
+        {
+          is_poly_left_inside = s.IsLeftInside(p0);
+          is_poly_right_inside = s.IsRightInside(p0);
+          first = true;
+        }
+
+        auto li = is_poly_left_inside; // == poly.IsLeftInside(p0);
+        auto ri = is_poly_right_inside; // == poly.IsRightInside(p0);
 
         auto & ls = seg_map[{pi0,pi1,pi2}];
         ls.p0 = pi0;
@@ -1694,7 +1885,7 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
 
         if(li!=ri)
         {
-          if(s.IsLeftInside(p0) != flip)
+          if(li != flip)
             ls.left = dom;
           else
             ls.right = dom;
@@ -1708,10 +1899,12 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
     else
       dom--; // degenerated solid, use same domain index again
   }
+  t_segments_map.Stop();
 
   for(auto bc : Range(bcnames))
     geo->SetBCName(bc+1, bcnames[bc]);
 
+  t_segments.Start();
   for(auto const &m : seg_map)
   {
     auto ls = m.second;
@@ -1737,6 +1930,7 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
     seg->hmax = ls.maxh;
     geo->AppendSegment(seg);
   }
+  t_segments.Stop();
   return geo;
 }
 

libsrc/geom2d/csg2d.hpp

@@ -95,6 +95,12 @@ struct EdgeInfo
 struct Vertex : Point<2>
 {
   Vertex (Point<2> p) : Point<2>(p) {}
+  Vertex (const Vertex & v) : Point<2>(v)
+  {
+    spline = v.spline;
+    info = v.info;
+    is_source = true;
+  }
 
   Vertex * prev = nullptr;
   Vertex * next = nullptr;
@@ -399,15 +405,49 @@ struct Loop
       AppendVertex(*v);
   }
 
+  Loop(Loop && p) = default;
+
+  Loop & operator=(Loop && p) = default;
+
   Loop & operator=(const Loop & p)
   {
+    // static Timer t("Loop::operator="); RegionTimer rt(t);
     first = nullptr;
     if(p.first)
+    {
+      size_t n = p.Size();
+      Array<unique_ptr<Vertex>> new_verts(n);
+      {
+        size_t i = 0;
         for(const auto v : p.Vertices(ALL))
-        AppendVertex(*v);
+          new_verts[i++] = make_unique<Vertex>(*v);
+      }
+
+      for(auto i : IntRange(n-1))
+      {
+        Vertex * v  = new_verts[i].get();
+        Vertex * vn = new_verts[i+1].get();
+        v->next = vn;
+        vn->prev = v;
+      }
+      Vertex * vfirst = new_verts[0].get();
+      Vertex * vlast = new_verts[n-1].get();
+      vfirst->prev = vlast;
+      vlast->next = vfirst;
+
+      for(auto i : IntRange(1,n))
+        new_verts[n-1-i]->pnext = std::move(new_verts[n-i]);
+
+      first = std::move(new_verts[0]);
+    }
     return *this;
   }
 
+  void Clear()
+  {
+    first = nullptr;
+  }
+
   Vertex & AppendVertex(const Vertex & v)
   {
     auto & vnew = Append( static_cast<Point<2>>(v), true );
@@ -448,6 +488,8 @@ struct Loop
   }
 
   bool IsInside( Point<2> r ) const;
+  bool IsLeftInside( const Vertex & p0 );
+  bool IsRightInside( const Vertex & p0 );
 
   EdgeIterator Edges(IteratorType iterType) const
   {
@@ -551,6 +593,8 @@ struct Loop
 
     return cnt;
   }
+
+  netgen::Box<2> GetBoundingBox() const;
 };
 
 
@@ -563,12 +607,15 @@ struct Solid2d
   Solid2d() = default;
   Solid2d(string name_) : name(name_) {}
   DLL_HEADER Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_=MAT_DEFAULT, string bc_=BC_DEFAULT);
+  Solid2d(Solid2d && other) = default;
+  Solid2d(const Solid2d & other) = default;
 
   DLL_HEADER Solid2d operator+(const Solid2d & other) const;
   DLL_HEADER Solid2d operator*(const Solid2d & other) const;
   DLL_HEADER Solid2d operator-(const Solid2d & other) const;
 
-  DLL_HEADER Solid2d& operator=(const Solid2d & other) = default;
+  Solid2d& operator=(Solid2d && other) = default;
+  Solid2d& operator=(const Solid2d & other) = default;
   DLL_HEADER Solid2d& operator+=(const Solid2d & other);
   DLL_HEADER Solid2d& operator*=(const Solid2d & other);
   DLL_HEADER Solid2d& operator-=(const Solid2d & other);
@@ -632,7 +679,7 @@ struct Solid2d
     return *this;
   }
 
-  netgen::Box<2> GetBoundingBox();
+  netgen::Box<2> GetBoundingBox() const;
 };
 
 
@@ -654,7 +701,12 @@ DLL_HEADER Solid2d Circle( Point<2> center, double r, string name="", string bc=
 DLL_HEADER Solid2d Rectangle( Point<2> p0, Point<2> p1, string mat=MAT_DEFAULT, string bc=BC_DEFAULT );
 
 DLL_HEADER void AddIntersectionPoints ( Solid2d & s1, Solid2d & s2 );
-DLL_HEADER Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect=true );
+DLL_HEADER Solid2d ClipSolids ( const Solid2d & s1, const Solid2d & s2, char op);
+DLL_HEADER Solid2d ClipSolids ( const Solid2d & s1, Solid2d && s2, char op);
+DLL_HEADER Solid2d ClipSolids ( Solid2d && s1, const Solid2d & s2, char op);
+DLL_HEADER Solid2d ClipSolids ( Solid2d && s1, Solid2d && s2, char op);
+
+DLL_HEADER IntersectionType intersect(const Point<2> P1, const Point<2> P2, const Point<2> Q1, const Point<2> Q2, double& alpha, double& beta);
 
 }
 #endif // NETGEN_CSG2D_HPP_INCLUDED