0033433: Shape Healing - Implement a new mode to keep initial types of curves

The implemented method modifies curves in wire so it becomes connected in vertices well with minor shifting of the curves

src/BRepTest/BRepTest_CurveCommands.cxx

@@ -42,9 +42,12 @@
 #include <Geom2d_Circle.hxx>
 #include <Geom2dAPI_Interpolate.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
+#include <Message.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Shape.hxx>
 #include <TopoDS_Wire.hxx>
+#include <ShapeFix_Wire.hxx>
+#include <ShapeFix_Edge.hxx>
 
 #include <DBRep.hxx>
 #include <DBRep_DrawableShape.hxx>
@@ -271,6 +274,142 @@ static Standard_Integer wire(Draw_Interpretor& di, Standard_Integer n, const cha
   return 0;
 }
 
+//=======================================================================
+// strongwire
+//=======================================================================
+static Standard_Integer strongwire(Draw_Interpretor& theDI, Standard_Integer theArgC, const char** theArgV)
+{
+  enum StrongWireMode {
+    StrongWireMode_FixTolerance = 1,
+    StrongWireMode_Approximation = 2,
+    StrongWireMode_KeepCurveType = 3
+  };
+
+  BRep_Builder aB;
+
+  TopoDS_Wire aWire;
+  aB.MakeWire(aWire);
+
+  if (theArgC < 3)
+    return 1;
+
+  // Tolerance
+  double aTolerance = Precision::Confusion();
+
+  // mode
+  StrongWireMode aMode = StrongWireMode_KeepCurveType;
+
+  // add edges
+  for (Standard_Integer anArgIter = 2; anArgIter < theArgC; anArgIter++)
+  {
+    TCollection_AsciiString aParam(theArgV[anArgIter]);
+    if (aParam == "-t")
+    {
+      anArgIter++;
+      if (anArgIter < theArgC)
+        aTolerance = Draw::Atof(theArgV[anArgIter]);
+    }
+    else if (aParam == "-m")
+    {
+      anArgIter++;
+      if (anArgIter < theArgC)
+      {
+        if (aParam == "keepType" || Draw::Atoi(theArgV[anArgIter]) == 1)
+        {
+          aMode = StrongWireMode_KeepCurveType;
+          continue;
+        }
+        else if (aParam == "approx" || Draw::Atoi(theArgV[anArgIter]) == 2)
+        {
+          aMode = StrongWireMode_Approximation;
+          continue;
+        }
+        else if (aParam == "fixTol" || Draw::Atoi(theArgV[anArgIter]) == 3)
+        {
+          aMode = StrongWireMode_FixTolerance;
+          continue;
+        }
+      }
+    }
+    else
+    {
+      TopoDS_Shape aShape = DBRep::Get(theArgV[anArgIter]);
+      if (aShape.IsNull())
+        Standard_NullObject::Raise("Shape for wire construction is null");
+      if (aShape.ShapeType() == TopAbs_EDGE || aShape.ShapeType() == TopAbs_WIRE)
+      {
+        TopExp_Explorer anExp(aShape, TopAbs_EDGE);
+        for (; anExp.More(); anExp.Next())
+          aB.Add(aWire, TopoDS::Edge(anExp.Current()));
+      }
+      else
+        Standard_TypeMismatch::Raise("Shape for wire construction is neither an edge nor a wire");
+    }
+  }
+
+  // fix edges order
+  Handle(ShapeFix_Wire) aFW = new ShapeFix_Wire;
+  aFW->Load(aWire);
+  aFW->FixReorder();
+
+  if (aFW->StatusReorder(ShapeExtend_FAIL1)) 
+  {
+    Message::SendFail() << "Error: Wire construction failed: several loops detected";
+    return 1;
+  }
+  else if (aFW->StatusReorder(ShapeExtend_FAIL)) 
+  {
+    Message::SendFail() << "Wire construction failed";
+    return 1;
+  }
+
+  bool isClosed = false;
+  Handle(ShapeAnalysis_Wire) aSaw = aFW->Analyzer();
+  if (aSaw->CheckGap3d(1)) // between last and first edges 
+  {
+    Standard_Real aDist = aSaw->MinDistance3d();
+    if (aDist < aTolerance)
+      isClosed = true;
+  }
+  aFW->ClosedWireMode() = isClosed;
+  aFW->FixConnected(aTolerance);
+  if (aMode == StrongWireMode_KeepCurveType)
+    aFW->FixCurves();
+
+  if (aFW->StatusConnected(ShapeExtend_FAIL))
+  {
+    Message::SendFail() << "Wire construction failed: cannot build connected wire";
+    return 1;
+  }
+
+  if (aFW->StatusConnected(ShapeExtend_DONE3)) 
+  {
+    if (aMode != StrongWireMode_Approximation)
+      aFW->SetPrecision(aTolerance); 
+    aFW->FixGapsByRangesMode() = Standard_True;
+    if (aFW->FixGaps3d())
+    {
+      Handle(ShapeExtend_WireData) sbwd = aFW->WireData();
+      Handle(ShapeFix_Edge) aFe = new ShapeFix_Edge;
+      for (Standard_Integer anIdx = 1; anIdx <= sbwd->NbEdges(); anIdx++) 
+      {
+        TopoDS_Edge aEdge = TopoDS::Edge(sbwd->Edge(anIdx));
+        aFe->FixVertexTolerance(aEdge);
+        aFe->FixSameParameter(aEdge);
+      }
+    }
+    else if (aFW->StatusGaps3d(ShapeExtend_FAIL))
+    {
+      Message::SendFail() << "Wire construction failed: cannot fix 3d gaps";
+      return 1;
+    }
+  }
+  aWire = aFW->WireAPIMake();
+
+  DBRep::Set(theArgV[1], aWire);
+  return 0;
+}
+
 //=======================================================================
 // mkedge
 //=======================================================================
@@ -1997,6 +2136,10 @@ void  BRepTest::CurveCommands(Draw_Interpretor& theCommands)
     "wire wirename [-unsorted] e1/w1 [e2/w2 ...]",__FILE__,
     wire,g);
 
+  theCommands.Add("strongwire",
+    "strongwire wirename e1/w1 [e2/w2 ...] [-t tol] [-m mode(keepType/1 | approx/2 | fixTol/3)]", __FILE__,
+    strongwire, g);
+
   theCommands.Add("profile",
     "profile, no args to get help",__FILE__,
     profile,g);

src/ShapeFix/ShapeFix_Wire.cxx

@@ -52,13 +52,19 @@
 #include <BRepBuilderAPI_MakeEdge.hxx>
 #include <BRepBuilderAPI_MakeVertex.hxx>
 #include <BRepTools.hxx>
+#include <GC_MakeArcOfCircle.hxx>
+#include <GC_MakeArcOfEllipse.hxx>
+#include <GC_MakeLine.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_Line.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
 #include <Geom2dConvert.hxx>
-#include <Geom_BSplineCurve.hxx>
 #include <Geom_Curve.hxx>
+#include <Geom_Circle.hxx>
+#include <Geom_Ellipse.hxx>
+#include <Geom_BezierCurve.hxx>
+#include <Geom_BSplineCurve.hxx>
 #include <Geom_OffsetCurve.hxx>
 #include <Geom_Plane.hxx>
 #include <Geom_SphericalSurface.hxx>
@@ -68,6 +74,7 @@
 #include <GeomAdaptor_Curve.hxx>
 #include <GeomAPI.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
 #include <IntRes2d_IntersectionPoint.hxx>
 #include <IntRes2d_SequenceOfIntersectionPoint.hxx>
@@ -512,6 +519,25 @@ Standard_Boolean ShapeFix_Wire::FixConnected (const Standard_Real prec)
   return StatusConnected ( ShapeExtend_DONE );
 }
 
+//=======================================================================
+//function : FixCurves
+//purpose  : 
+//=======================================================================
+Standard_Boolean ShapeFix_Wire::FixCurves()
+{
+  myStatusConnected = ShapeExtend::EncodeStatus(ShapeExtend_OK);
+  if (!IsLoaded()) return Standard_False;
+
+  Standard_Integer stop = (myClosedMode ? 0 : 1);
+  for (Standard_Integer anIdx = NbEdges(); anIdx > stop; anIdx--) 
+  {
+    FixCurves(anIdx);
+    myStatusConnected |= myLastFixStatus;
+  }
+
+  return StatusConnected(ShapeExtend_DONE);
+}
+
 //=======================================================================
 //function : FixEdgeCurves
 //purpose  : 
@@ -1299,6 +1325,169 @@ Standard_Boolean ShapeFix_Wire::FixConnected (const Standard_Integer num,
   return Standard_True;
 }
 
+//=======================================================================
+//function : FixCurves
+//purpose  : 
+//=======================================================================
+Standard_Boolean ShapeFix_Wire::FixCurves(const Standard_Integer theIdx)
+{
+  // assume fix curves step should be after fix vertices
+  myLastFixStatus = ShapeExtend::EncodeStatus(ShapeExtend_OK);
+  if (!IsLoaded() || NbEdges() <= 0) 
+    return Standard_False;
+
+  // action: replacing curves in edges
+  Handle(ShapeExtend_WireData) anSbwd = WireData();
+  Standard_Integer anIdx = (theIdx > 0 ? theIdx : anSbwd->NbEdges());
+  TopoDS_Edge anEdge = anSbwd->Edge(anIdx);
+  Standard_Real aPrec = BRep_Tool::Tolerance(anEdge);
+  ShapeAnalysis_Edge aSae;
+  ShapeAnalysis_Wire aSaw;
+  Handle(Geom_Curve) aCurve3d;
+  Standard_Real aCurBounds[3];
+  Standard_Boolean IsReversed = Standard_False;
+  aSae.Curve3d(anEdge, aCurve3d, aCurBounds[0], aCurBounds[2], IsReversed);
+  aCurBounds[1] = (aCurBounds[0] + aCurBounds[2]) / 2;
+  gp_Pnt anEnds[3];
+  anEnds[0] = BRep_Tool::Pnt(aSae.FirstVertex(anEdge));
+  aCurve3d->D0(aCurBounds[1], anEnds[1]);
+  anEnds[2] = BRep_Tool::Pnt(aSae.LastVertex(anEdge));
+
+  gp_Pnt aGeomEnds[2];
+  aCurve3d->D0(aCurBounds[0], aGeomEnds[0]);
+  aCurve3d->D0(aCurBounds[2], aGeomEnds[1]);
+
+  // TODO: precise, if IsReversed flag should be considered here
+  Standard_Real aGap0 = Min(anEnds[0].Distance(aGeomEnds[0]), anEnds[0].Distance(aGeomEnds[1]));
+  Standard_Real aGap2 = Min(anEnds[2].Distance(aGeomEnds[0]), anEnds[2].Distance(aGeomEnds[1]));
+  if (Max (aGap0, aGap2) < aPrec) // nothing to do
+    return true;
+
+  if (aCurve3d->IsKind(STANDARD_TYPE(Geom_Circle))) 
+  {
+    Standard_Real anOldR = Handle(Geom_Circle)::DownCast(aCurve3d)->Circ().Radius();
+    gp_Vec anArcNorm = gp_Vec(anEnds[2], anEnds[0]) / 2;
+    gp_Pnt aCenter(anEnds[0].XYZ() - anArcNorm.XYZ());
+    int aSign = anEnds[1].Distance(aCenter) > anOldR ? 1 : -1; // for arc length > PI
+    Standard_Real aD = anOldR*anOldR - anArcNorm.SquareMagnitude();
+    aD = abs(aD) < Precision::Confusion() ? 0. : aD;
+    Standard_Real anArcR = anOldR + aSign * sqrt(aD);
+    gp_Ax2 aNormal(aCenter, anArcNorm);
+    Handle(Geom_Circle) anArc = new Geom_Circle(aNormal, anArcR);
+    GeomAPI_ProjectPointOnCurve projector(anEnds[1], anArc);
+    anEnds[1] = projector.NearestPoint();
+    GC_MakeArcOfCircle arc(anEnds[0], anEnds[1], anEnds[2]);
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(arc.Value()).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_Ellipse)))
+  {
+    Handle(Geom_Ellipse) anOld = Handle(Geom_Ellipse)::DownCast(aCurve3d);
+    Handle(Geom_Plane) aPln = new Geom_Plane(anEnds[0], gp_Vec(anEnds[2], anEnds[0]).Crossed(gp_Vec(anEnds[1], anEnds[0])));
+    GeomAPI_ProjectPointOnSurf aProjector(anOld->Elips().Location(), aPln);
+    gp_Pnt anOrigin = aProjector.NearestPoint();
+    aProjector.Init(anOld->Elips().Location().XYZ() + anOld->Elips().XAxis().Direction().XYZ(), aPln);
+    gp_Ax2 anAx(anOrigin, aPln->Axis().Direction(), aProjector.NearestPoint().XYZ() - anOrigin.XYZ());
+    
+    // compute angle
+    Standard_Real aRec = DBL_MAX;
+    Standard_Real anAngle = 0.;
+    Standard_Integer aSplNum = 10;
+    for (Standard_Integer anIdxI = -aSplNum; anIdxI < aSplNum; ++anIdxI)
+    {
+      Handle(Geom_Ellipse) anEll = new Geom_Ellipse(anAx, anOld->MajorRadius(), anOld->MinorRadius());
+      anEll->Rotate(anAx.Axis(), aPrec*anIdxI / anEll->MajorRadius() / aSplNum);
+      GeomAPI_ProjectPointOnCurve aProjector1(anEnds[0], anEll);
+      Standard_Real aDist = 0.;
+      for (Standard_Integer anIdxJ = 0; anIdxJ < 2; ++anIdxJ)
+      {
+        aProjector1.Perform(anEnds[anIdxJ *2]);
+        aDist += std::fmin (aProjector1.Distance(1), aProjector1.Distance(2));
+      }
+      if (aDist < aRec)
+      {
+        aRec = aDist;
+        anAngle = aPrec*anIdxI / anEll->MajorRadius() / aSplNum;
+      }
+    }
+    gp_Elips aTemp(anAx, anOld->MajorRadius(), anOld->MinorRadius());
+    aTemp.Rotate(anAx.Axis(), anAngle);
+
+    // compute shift
+    gp_Vec aX = aTemp.XAxis().Direction();
+    gp_Vec aY = aTemp.YAxis().Direction();
+    gp_Pnt2d aP1((anEnds[0].XYZ() - anOrigin.XYZ()).Dot(aX.XYZ()), (anEnds[0].XYZ() - anOrigin.XYZ()).Dot(aY.XYZ()));
+    gp_Pnt2d aP2((anEnds[2].XYZ() - anOrigin.XYZ()).Dot(aX.XYZ()), (anEnds[2].XYZ() - anOrigin.XYZ()).Dot(aY.XYZ()));
+    
+    // x = ky + p   linear equation
+    // where (x, y) shift point, 
+    // k, p constant coefficients
+    Standard_Real k = 1, p = 0;
+    Standard_Real R = anOld->MajorRadius();
+    Standard_Real r = anOld->MinorRadius();
+    k = (R / r) * (R / r) *
+      (aP1.Y() - aP2.Y()) / (aP2.X() - aP1.X());
+    p = -(1. / 2) * (R / r) * (R / r) *
+      (aP1.Y()*aP1.Y() - aP2.Y()*aP2.Y()) / (aP2.X() - aP1.X()) + aP1.X() / 2 + aP2.X() / 2;
+    // ax^2 + bx + c = 0  square equation 
+    // a, b, c constant coefficients
+    Standard_Real a = 0., b = 0., c = 0.;
+    a = R*R + k*k*r*r;
+    b = 2 * (k*p*r*r - k*aP1.X()*r*r - aP1.Y()*R*R);
+    c = aP1.X()*aP1.X()*r*r +
+      aP1.Y()*aP1.Y()*R*R -
+      r*r*R*R +
+      p*p*r*r - 2 * aP1.X()*p*r*r;
+    Standard_Real y1 = (-b - sqrt(b*b - 4 * a*c)) / 2 / a;
+    Standard_Real y2 = (-b + sqrt(b*b - 4 * a*c)) / 2 / a;
+    Standard_Real x1 = k*y1 + p;
+    Standard_Real x2 = k*y2 + p;
+    
+    gp_Pnt anOri = anOld->Location();
+    if (x1*x1 + y1*y1 < x2*x2 + y2*y2) 
+      anOri = anOri.XYZ() + aX.XYZ()*x1 + aY.XYZ()*y1;
+    else
+      anOri = anOri.XYZ() + aX.XYZ()*x2 + aY.XYZ()*y2;
+    aTemp.SetLocation(anOri);
+    
+    GC_MakeArcOfEllipse anArc(aTemp, anEnds[2], anEnds[0], true);
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(anArc.Value()).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_Line)))
+  {
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(anEnds[2], anEnds[0]).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) 
+  {
+    Handle(Geom_BSplineCurve) anOld = Handle(Geom_BSplineCurve)::DownCast(aCurve3d);
+    Standard_Integer f, l;
+    int p = anEnds[0].Distance(aGeomEnds[0]) < anEnds[1].Distance(aGeomEnds[0]) ? 0 : 2;
+    anOld->MovePoint(aCurBounds[0], anEnds[p], 1, 1, f, l);
+    anOld->MovePoint(aCurBounds[2], anEnds[2-p], anOld->NbPoles(), anOld->NbPoles(), f, l);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_BezierCurve))) 
+  {
+    Handle(Geom_BezierCurve) anOld = Handle(Geom_BezierCurve)::DownCast(aCurve3d);
+    int p = anEnds[0].Distance(aGeomEnds[0]) < anEnds[1].Distance(aGeomEnds[0]) ? 0 : 2;
+    anOld->SetPole(1, anEnds[p]);
+    anOld->SetPole(anOld->NbPoles(), anEnds[2-p]);
+    return true;
+  }
+
+  // TODO: question: the below code works only for other curve types (not line/arc/circle/ellipse/bspline)
+  //                 Is it really needed?
+  myAnalyzer->CheckConnected(theIdx, aPrec);
+  if (myAnalyzer->LastCheckStatus(ShapeExtend_FAIL))
+    myLastFixStatus |= ShapeExtend::EncodeStatus(ShapeExtend_FAIL1);
+
+  return true;
+}
 
 //=======================================================================
 //function : FixSeam

src/ShapeFix/ShapeFix_Wire.hxx

@@ -277,6 +277,7 @@ public:
   //! flag ClosedMode is True
   //! If <prec> is -1 then MaxTolerance() is taken.
   Standard_EXPORT Standard_Boolean FixConnected (const Standard_Real prec = -1.0);
+  Standard_EXPORT Standard_Boolean FixCurves();
   
   //! Groups the fixes dealing with 3d and pcurves of the edges.
   //! The order of the fixes and the default behaviour are:
@@ -346,6 +347,7 @@ public:
   //! Tests with starting preci or, if given greater, <prec>
   //! If <prec> is -1 then MaxTolerance() is taken.
   Standard_EXPORT Standard_Boolean FixConnected (const Standard_Integer num, const Standard_Real prec);
+  Standard_EXPORT Standard_Boolean FixCurves(const Standard_Integer num);
   
   //! Fixes a seam edge
   //! A Seam edge has two pcurves, one for forward. one for reversed

tests/heal/grids.list

@@ -25,4 +25,5 @@
 025 update_tolerance_locked
 026 checkshape
 027 split_number
-028 split_two_numbers
+028 split_two_numbers
+029 wire_fix_curves

tests/heal/wire_fix_curves/strongwire_circle

@@ -0,0 +1,16 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+circle c1 0 0 0 15
+circle c2 0 0 0 15
+trim c1 c1 0 pi
+trim c2 c2 pi 2*pi
+translate c1 0 0 0.00000005
+translate c2 0 0 -0.00000005
+mkedge c1 c1
+mkedge c2 c2
+strongwire w c1 c2
+checkshape w
+whatis w

tests/heal/wire_fix_curves/strongwire_ellipse

@@ -0,0 +1,16 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+ellipse e1 0 0 0 25 15
+ellipse e2 0 0 0 25 15
+trim e1 e1 0 pi
+trim e2 e2 pi 2*pi
+translate e1 0 0 0.00000005
+translate e2 0 0 -0.00000005
+mkedge e1 e1
+mkedge e2 e2
+strongwire w e1 e2 
+checkshape w
+whatis w

tests/heal/wire_fix_curves/strongwire_halfcircle1

@@ -0,0 +1,18 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+point px 100   0   0
+point py   0 100   0
+point pz   0   0 100
+
+vertex vx px
+vertex vy py
+
+edge exy vx vy
+gcarc arc cir py pz px
+mkedge earc arc
+
+strongwire w exy earc
+whatis w

tests/heal/wire_fix_curves/strongwire_halfcircle2

@@ -0,0 +1,23 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+point p1   0   0  0
+point p2 100   0  0
+point p4 100 110  0
+
+
+gcarc arc cir p1 p2 p4
+
+
+mkedge arc1 arc
+mkedge arc2 arc
+
+tmirror arc2 0 0 0 -110 100 0 -copy
+ttranslate arc2 0 0 0.003 -copy
+
+strongwire w1 arc1 arc2 -t 0.01 -m keepType
+
+checkshape w1
+whatis w1

tests/heal/wire_fix_curves/strongwire_nurb

@@ -0,0 +1,14 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+bsplinecurve c1 3 2  -1.0 4   1.0 4  0 0 0 1  1 4 0 1  2 4 0 1  3 0 0 1
+bsplinecurve c2 3 2  -1.0 4   1.0 4  0 0 0 1  1 -4 0 1  2 -4 0 1  3 0 0 1
+translate c1 0 0 0.00000005
+translate c2 0 0 -0.00000005
+mkedge c1 c1
+mkedge c2 c2
+strongwire w c1 c2
+checkshape w
+whatis w

tests/heal/wire_fix_curves/strongwire_triangle

@@ -0,0 +1,17 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+vertex v1   0   0 0
+vertex v2 100   0 0
+vertex v3  50 100 0
+
+edge e1 v1 v2
+edge e2 v2 v3
+edge e3 v3 v1
+
+strongwire w1 e1 e2 e3
+
+checkshape w1
+whatis w1