0028830: HalfSpace command chooses the wrong side of the given shell
authormsv <msv@opencascade.com>
Fri, 9 Jun 2017 07:46:34 +0000 (10:46 +0300)
committerbugmaster <bugmaster@opencascade.com>
Thu, 15 Jun 2017 08:31:39 +0000 (11:31 +0300)
Improve the algorithm BRepPrimAPI_MakeHalfSpace. Earlier it made projection of the point only on faces. If the nearest point does not conform to normal projection criterion the result is wrong. The fix includes search of projection on edges and vertices. This makes the algorithm robust for half spaces with boundaries.

src/BRepPrimAPI/BRepPrimAPI_MakeHalfSpace.cxx
tests/bugs/modalg_7/bug28830 [new file with mode: 0644]

index c397825..60df649 100644 (file)
@@ -18,7 +18,7 @@
 #include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
 #include <BRepBuilderAPI_MakeVertex.hxx>
-#include <BRepExtrema_ExtPF.hxx>
+#include <BRepExtrema_DistShapeShape.hxx>
 #include <BRepLProp_SLProps.hxx>
 #include <BRepPrimAPI_MakeHalfSpace.hxx>
 #include <gp.hxx>
 #include <TopoDS_Vertex.hxx>
 
 //=======================================================================
-//function : FindExtrema
-//purpose  : This finction is called to find the nearest normal projection
-//           of a point <aPnt> on a face <aFace>.
-//           1) return true if extrema are found.
-//           2) Set in:
-//             - Dist : The lower distance found.
-//             - anOppositePnt : The corresponding point lying on the face
-//             - U,V : The parameters of <anOppositePnt> on the face <aFace>
-//=======================================================================
-static Standard_Real FindExtrema(const gp_Pnt&        aPnt,
-                                const TopoDS_Face&   aFace,
-                                      Standard_Real& Dist,
-                                      gp_Pnt&        anOppositePnt,
-                                      Standard_Real& U,
-                                      Standard_Real& V)
+//function : getNormalOnFace
+//purpose  : 
+//=======================================================================
+
+static gp_Dir getNormalOnFace(const TopoDS_Face& theFace,
+                              const Standard_Real theU,
+                              const Standard_Real theV)
 {
-  Standard_Integer intvalue=0;
-  Dist = RealLast();
-  TopoDS_Vertex RefVertex = BRepBuilderAPI_MakeVertex(aPnt);
-  
-  BRepExtrema_ExtPF ext(RefVertex,aFace);
-  
-  if (ext.IsDone() && ext.NbExt() > 0) {
-    // the point projection exist
-    Standard_Integer nbext = ext.NbExt();
-    for (Standard_Integer iext = 1; iext <= nbext; iext++) {
-      if (ext.SquareDistance(iext) < Dist) {
-       Dist     = ext.SquareDistance(iext);
-       intvalue = iext;
+  Standard_Real aPrec = gp::Resolution();
+  BRepLProp_SLProps aProps(BRepAdaptor_Surface(theFace), theU, theV, 2, aPrec);
+  gp_Dir aNormal = aProps.Normal();
+  if (theFace.Orientation() == TopAbs_REVERSED)
+    aNormal.Reverse();
+  return aNormal;
+}
+
+//=======================================================================
+//function : getNormalFromEdge
+//purpose  : Get average normal at the point with the given parameter on the edge
+//=======================================================================
+
+static Standard_Boolean  getNormalFromEdge(const TopoDS_Shape& theShape,
+                                           const TopoDS_Edge& theEdge,
+                                           const Standard_Real thePar,
+                                           gp_Dir& theNormal)
+{
+  gp_XYZ aSum;
+  TopExp_Explorer ex(theShape, TopAbs_FACE);
+  for (; ex.More(); ex.Next()) {
+    const TopoDS_Face& aF = TopoDS::Face(ex.Current());
+    TopExp_Explorer ex1(aF, TopAbs_EDGE);
+    for (; ex1.More(); ex1.Next()) {
+      if (ex1.Current().IsSame(theEdge)) {
+        Standard_Real f, l;
+        Handle(Geom2d_Curve) aC2d = BRep_Tool::CurveOnSurface(theEdge, aF, f, l);
+        gp_Pnt2d aP2d = aC2d->Value(thePar);
+        gp_Dir aNorm = getNormalOnFace(aF, aP2d.X(), aP2d.Y());
+        aSum += aNorm.XYZ();
       }
     }
-    Dist = sqrt(Dist);
-  } else {
-    // compute the min distance with the face vertices
-    TopExp_Explorer explo(aFace, TopAbs_VERTEX);
-    for(; explo.More(); explo.Next()) {
-      const TopoDS_Vertex& vertex = TopoDS::Vertex(explo.Current());
-      gp_Pnt2d             param  = BRep_Tool::Parameters(vertex, aFace);
-      gp_Pnt               Pnt    = BRep_Tool::Pnt(vertex);
-      
-      Standard_Real        d2      = Pnt.SquareDistance(aPnt);
-      if(d2 < Dist) {
-       Dist          = d2;
-       anOppositePnt = Pnt;
-       param.Coord(U, V);
+  }
+  if (aSum.SquareModulus() > gp::Resolution()) {
+    theNormal = aSum;
+    return Standard_True;
+  }
+  return Standard_False;
+}
+
+//=======================================================================
+//function : getNormalFromVertex
+//purpose  : Get average normal at the point of the vertex
+//=======================================================================
+
+static Standard_Boolean  getNormalFromVertex(const TopoDS_Shape& theShape,
+                                             const TopoDS_Vertex& theVer,
+                                             gp_Dir& theNormal)
+{
+  gp_XYZ aSum;
+  TopExp_Explorer ex(theShape, TopAbs_FACE);
+  for (; ex.More(); ex.Next()) {
+    const TopoDS_Face& aF = TopoDS::Face(ex.Current());
+    TopExp_Explorer ex1(aF, TopAbs_VERTEX);
+    for (; ex1.More(); ex1.Next()) {
+      if (ex1.Current().IsSame(theVer)) {
+        gp_Pnt2d aP2d = BRep_Tool::Parameters(theVer, aF);
+        gp_Dir aNorm = getNormalOnFace(aF, aP2d.X(), aP2d.Y());
+        aSum += aNorm.XYZ();
       }
     }
-    Dist = sqrt(Dist);
+  }
+  if (aSum.SquareModulus() > gp::Resolution()) {
+    theNormal = aSum;
     return Standard_True;
   }
+  return Standard_False;
+}
 
-  // Final result.
-  anOppositePnt = ext.Point(intvalue);
-  ext.Parameter(intvalue,U,V);
+//=======================================================================
+//function : FindExtrema
+//purpose  : This finction is called to find the nearest normal projection
+//           of a point <aPnt> on a shape <aShape>.
+//           1) return true if extrema is found.
+//           2) Set in:
+//             - theMinPnt : The solution point
+//             - theNormal : The normal direction to the shape at projection point
+//=======================================================================
+static Standard_Boolean FindExtrema(const gp_Pnt&        thePnt,
+                                    const TopoDS_Shape&  theShape,
+                                    gp_Pnt&              theMinPnt,
+                                    gp_Dir&              theNormal)
+{
+  TopoDS_Vertex aRefVertex = BRepBuilderAPI_MakeVertex(thePnt);
+  
+  BRepExtrema_DistShapeShape ext(aRefVertex, theShape);
+  
+  if (!ext.IsDone() || ext.NbSolution() == 0)
+    return Standard_False;
 
-  return Standard_True;
+  // the point projection exist
+  Standard_Integer nbext = ext.NbSolution();
+  // try to find a projection on face
+  for (Standard_Integer iext = 1; iext <= nbext; iext++) {
+    if (ext.SupportTypeShape2(iext) == BRepExtrema_IsInFace) {
+      TopoDS_Face aF = TopoDS::Face(ext.SupportOnShape2(iext));
+      theMinPnt = ext.PointOnShape2(iext);
+      Standard_Real aU, aV;
+      ext.ParOnFaceS2(iext, aU, aV);
+      theNormal = getNormalOnFace(aF, aU, aV);
+      return Standard_True;
+    }
+  }
+
+  // if not found then take any edge or vertex solution
+  for (Standard_Integer iext = 1; iext <= nbext; iext++) {
+    if (ext.SupportTypeShape2(iext) == BRepExtrema_IsOnEdge) {
+      theMinPnt = ext.PointOnShape2(iext);
+      Standard_Real aPar;
+      ext.ParOnEdgeS2(iext, aPar);
+      TopoDS_Edge aE = TopoDS::Edge(ext.SupportOnShape2(iext));
+      if (getNormalFromEdge(theShape, aE, aPar, theNormal))
+        return Standard_True;
+    }
+    else if (ext.SupportTypeShape2(iext) == BRepExtrema_IsVertex) {
+      theMinPnt = ext.PointOnShape2(iext);
+      TopoDS_Vertex aV = TopoDS::Vertex(ext.SupportOnShape2(iext));
+      if (getNormalFromVertex(theShape, aV, theNormal))
+        return Standard_True;
+    }
+  }
+  return Standard_False;
 }
 
+//=======================================================================
+//function : isOutside
+//purpose  : 
+//=======================================================================
 
+static Standard_Boolean isOutside(const gp_Pnt&      thePnt,
+                                  const gp_Pnt&      thePonF,
+                                  const gp_Dir&      theNormal)
+{
+  gp_Dir anOppRef(thePnt.XYZ() - thePonF.XYZ());
+  Standard_Real aSca = theNormal * anOppRef;
+  // outside if same directions
+  return aSca > 0.;
+}
 
 //=======================================================================
 //function : BRepPrimAPI_MakeHalfSpace
 //purpose  : 
 //=======================================================================
 
-BRepPrimAPI_MakeHalfSpace::BRepPrimAPI_MakeHalfSpace(const TopoDS_Face& Face,
-                                            const gp_Pnt&      RefPnt)
+BRepPrimAPI_MakeHalfSpace::BRepPrimAPI_MakeHalfSpace(const TopoDS_Face& theFace,
+                                                     const gp_Pnt&      theRefPnt)
 {
   // Set the flag is <IsDone> to False.
   NotDone();
 
-  TopoDS_Shell Shell;
-
-//  gp_Pnt CurPnt, MinPnt;
-  gp_Pnt MinPnt;
-  Standard_Real U, V;
-
-  Standard_Real MinDist;
-  if (FindExtrema(RefPnt,Face,MinDist,MinPnt,U,V)) {
-    BRep_Builder myBuilder;
-    myBuilder.MakeShell(Shell);
-    myBuilder.Add(Shell,Face);
-    
-    // Normal, scalair product and direction.
-    Standard_Real Prec = gp::Resolution();
-//    BRepLProp_SLProps Props(BRepAdaptor_Surface(Face),U,V,2,Prec);
-    BRepLProp_SLProps Props = BRepLProp_SLProps(BRepAdaptor_Surface(Face),U,V,2,Prec);
-    gp_Dir Normale    = Props.Normal();
-    gp_Dir OppRef(RefPnt.XYZ()-MinPnt.XYZ());
-    Standard_Real Sca = Normale*OppRef;
-    
+  TopoDS_Shell aShell;
+
+  gp_Pnt aMinPnt;
+  gp_Dir aNormal;
+  if (FindExtrema(theRefPnt, theFace, aMinPnt, aNormal)) {
+    Standard_Boolean toReverse = isOutside(theRefPnt, aMinPnt, aNormal);
+
     // Construction of the open solid.
-    myBuilder.MakeSolid(mySolid);
-    if (Sca > 0.) {
-      // Same directions: inverted case.
-      Shell.Reverse();
+    BRep_Builder().MakeShell(aShell);
+    BRep_Builder().Add(aShell, theFace);
+    BRep_Builder().MakeSolid(mySolid);
+    if (toReverse) {
+      aShell.Reverse();
     }
-    
-    myBuilder.Add(mySolid,Shell);
+    BRep_Builder().Add(mySolid, aShell);
     Done();
   }
 }
@@ -142,54 +216,25 @@ BRepPrimAPI_MakeHalfSpace::BRepPrimAPI_MakeHalfSpace(const TopoDS_Face& Face,
 //purpose  : 
 //=======================================================================
 
-BRepPrimAPI_MakeHalfSpace::BRepPrimAPI_MakeHalfSpace(const TopoDS_Shell& Shell,
-                                            const gp_Pnt&       RefPnt)
+BRepPrimAPI_MakeHalfSpace::BRepPrimAPI_MakeHalfSpace(const TopoDS_Shell& theShell,
+                                                     const gp_Pnt&       theRefPnt)
 {
   // Set the flag is <IsDone> to False.
   NotDone();
 
-  gp_Pnt CurPnt, MinPnt;
-  TopoDS_Face CurFace, MinFace;
-  Standard_Real MinDist = RealLast();
-  Standard_Real CurDist, U, V, MinU=0, MinV=0;
-
   // Find the point of the skin closest to the reference point.
-  Standard_Boolean YaExt = Standard_False;
-
-  TopoDS_Shell aShell = Shell;
-  TopExp_Explorer Exp(aShell,TopAbs_FACE);
-  while (Exp.More()) {
-    CurFace = TopoDS::Face(Exp.Current());
-    if ( FindExtrema(RefPnt,CurFace,CurDist,CurPnt,U,V)) {
-      YaExt = Standard_True;
-      if (CurDist < MinDist) {
-       MinDist = CurDist;
-       MinPnt  = CurPnt;
-       MinU    = U;
-       MinV    = V;
-       MinFace = CurFace;
-      }
-    }
-    Exp.Next();
-  }
+  gp_Pnt aMinPnt;
+  gp_Dir aNormal;
+  if (FindExtrema(theRefPnt, theShell, aMinPnt, aNormal)) {
+    Standard_Boolean toReverse = isOutside(theRefPnt, aMinPnt, aNormal);
 
-  if ( YaExt) {
-    // Normal, scalar product and direction.
-    BRep_Builder myBuilder;
-    Standard_Real Prec = gp::Resolution();
-//    BRepLProp_SLProps Props(BRepAdaptor_Surface(MinFace),MinU,MinV,2,Prec);
-    BRepLProp_SLProps Props = BRepLProp_SLProps(BRepAdaptor_Surface(MinFace),MinU,MinV,2,Prec);
-    gp_Dir Normale    = Props.Normal();
-    gp_Dir OppRef(RefPnt.XYZ()-MinPnt.XYZ());
-    Standard_Real Sca = Normale*OppRef;
-    
     // Construction of the open solid.
-    myBuilder.MakeSolid(mySolid);
-    if (Sca > 0.) {
-      // Same directions: inverted case.
+    TopoDS_Shell aShell = theShell;
+    BRep_Builder().MakeSolid(mySolid);
+    if (toReverse) {
       aShell.Reverse();
     }
-    myBuilder.Add(mySolid,aShell);
+    BRep_Builder().Add(mySolid, aShell);
     Done();
   }
 }
@@ -217,5 +262,3 @@ BRepPrimAPI_MakeHalfSpace::operator TopoDS_Solid() const
 {
   return Solid();
 }
-
-
diff --git a/tests/bugs/modalg_7/bug28830 b/tests/bugs/modalg_7/bug28830
new file mode 100644 (file)
index 0000000..d980445
--- /dev/null
@@ -0,0 +1,19 @@
+puts "========"
+puts "OCC28830"
+puts "========"
+puts ""
+########################################
+# HalfSpace command chooses the wrong side of the given shell
+########################################
+
+# Restore the initial shape
+restore [locate_data_file bug28830_halfspace.brep] sh
+
+point p 9.30222203002736 0.87421058209264 1.54257060749683
+halfspace h sh 9.30222203002736 0.87421058209264 1.54257060749683
+
+if {![regexp "IN" [bclassify h p]]} {
+  puts "Error: halfspace is wrong"
+} else {
+  puts "OK: halfspace is good"
+}