[occt.git] / src / BRepLib / BRepLib_FindSurface.cxx

// Created on: 1994-07-22
// Created by: Remi LEQUETTE
// Copyright (c) 1994-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.


#include <BRepLib_FindSurface.ixx>

#include <Precision.hxx>
#include <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <math_Gauss.hxx>

#include <gp_Lin.hxx>
#include <gp_Circ.hxx>
#include <gp_Elips.hxx>
#include <gp_Hypr.hxx>
#include <gp_Parab.hxx>
#include <gp_Ax2.hxx>
#include <gp_Ax3.hxx>
#include <gp_Vec.hxx>
#include <TColgp_SequenceOfPnt.hxx>
#include <TColStd_SequenceOfReal.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_HArray1OfPnt.hxx>
#include <Geom_Plane.hxx>

#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepLib_MakeFace.hxx>
#include <BRepTools_WireExplorer.hxx>
#include <BRep_Tool.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS_Vertex.hxx>
#include <TopoDS_Wire.hxx>
#include <TopoDS.hxx>

#include <GeomLib.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom_BezierCurve.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_RectangularTrimmedSurface.hxx>
#include <Standard_ErrorHandler.hxx>

//=======================================================================
//function : Controle
//purpose  : 
//=======================================================================
static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
			      const Handle(Geom_Plane)& thePlane)
{
  Standard_Real dfMaxDist=0.;
  Standard_Real a,b,c,d, dist;
  Standard_Integer ii;
  thePlane->Coefficients(a,b,c,d);
  for (ii=1; ii<=thePoints.Length(); ii++) {
    const gp_XYZ& xyz = thePoints(ii).XYZ();
    dist = Abs(a*xyz.X() + b*xyz.Y() + c*xyz.Z() + d);
    if (dist > dfMaxDist)
      dfMaxDist = dist;
  }

  return dfMaxDist;
}

//=======================================================================
//function : Is2DConnected
//purpose  : Return true if the last vertex of theEdge1 coincides with
//           the first vertex of theEdge2 in parametric space of theFace
//=======================================================================

inline static Standard_Boolean Is2DConnected( const TopoDS_Edge& theEdge1,
                                              const TopoDS_Edge& theEdge2,
                                              const TopoDS_Face& theFace )
{
  Standard_Real f,l;
  Handle(Geom2d_Curve) aCurve;
  gp_Pnt2d p1, p2;

  // get 2D points
  aCurve   = BRep_Tool::CurveOnSurface( theEdge1, theFace,f,l );
  p1       = aCurve->Value( theEdge1.Orientation() == TopAbs_FORWARD ? l : f );
  aCurve   = BRep_Tool::CurveOnSurface( theEdge2,  theFace,f,l );
  p2       = aCurve->Value(  theEdge2.Orientation() == TopAbs_FORWARD ? f : l );

  // compare 2D points
  BRepAdaptor_Surface aSurface( theFace );
  TopoDS_Vertex    aV = TopExp::FirstVertex( theEdge2, /*CumOri=*/Standard_True );
  Standard_Real tol3D = BRep_Tool::Tolerance( aV );
  Standard_Real tol2D = aSurface.UResolution( tol3D ) + aSurface.VResolution( tol3D );
  Standard_Real dist2 = p1.SquareDistance( p2 );
  return dist2 < tol2D * tol2D;
}

//=======================================================================
//function : Is2DClosed
//purpose  : Return true if edges of theShape form a closed wire in
//           parametric space of theSurface
//=======================================================================

static Standard_Boolean Is2DClosed( const TopoDS_Shape&         theShape,
                                    const Handle(Geom_Surface)& theSurface)
{
  try
  {
    // get a wire theShape 
    TopExp_Explorer aWireExp( theShape, TopAbs_WIRE );
    if ( !aWireExp.More() )
      return Standard_False;
    TopoDS_Wire aWire = TopoDS::Wire( aWireExp.Current() );
    // a tmp face
    TopoDS_Face aTmpFace = BRepLib_MakeFace( theSurface, Precision::PConfusion() );

    // check topological closeness using wire explorer, if the wire is not closed
    // the 1st and the last vertices of wire are different
    BRepTools_WireExplorer aWireExplorer( aWire, aTmpFace );
    if ( !aWireExplorer.More())
      return Standard_False;
    // remember the 1st and the last edges of aWire
    TopoDS_Edge aFisrtEdge = aWireExplorer.Current(), aLastEdge = aFisrtEdge;
    // check if edges connected topologically (that is assured by BRepTools_WireExplorer)
    // are connected in 2D
    TopoDS_Edge aPrevEdge = aFisrtEdge;
    for ( aWireExplorer.Next(); aWireExplorer.More(); aWireExplorer.Next() )
    {
      aLastEdge = aWireExplorer.Current();
      if ( !Is2DConnected( aPrevEdge, aLastEdge, aTmpFace ))
        return false;
      aPrevEdge = aLastEdge;
    }
    // wire is closed if ( 1st vertex of aFisrtEdge ) ==
    // ( last vertex of aLastEdge ) in 2D
    TopoDS_Vertex aV1 = TopExp::FirstVertex( aFisrtEdge, /*CumOri=*/Standard_True );
    TopoDS_Vertex aV2 = TopExp::LastVertex( aLastEdge, /*CumOri=*/Standard_True );
    return ( aV1.IsSame( aV2 ) && Is2DConnected( aLastEdge, aFisrtEdge, aTmpFace ));
  }
  catch ( Standard_Failure )
  {
    return Standard_False;
  }
}

//=======================================================================
//function : BRepLib_FindSurface
//purpose  : 
//=======================================================================
BRepLib_FindSurface::BRepLib_FindSurface() 
{
}
//=======================================================================
//function : BRepLib_FindSurface
//purpose  : 
//=======================================================================
BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape&    S, 
					 const Standard_Real    Tol,
					 const Standard_Boolean OnlyPlane,
                                         const Standard_Boolean OnlyClosed)
{
  Init(S,Tol,OnlyPlane,OnlyClosed);
}
//=======================================================================
//function : Init
//purpose  : 
//=======================================================================
void BRepLib_FindSurface::Init(const TopoDS_Shape&    S, 
			       const Standard_Real    Tol,
			       const Standard_Boolean OnlyPlane,
                               const Standard_Boolean OnlyClosed)
{
  myTolerance = Tol;
  myTolReached = 0.;
  isExisted = Standard_False;
  myLocation.Identity();
  mySurface.Nullify();

  // compute the tolerance
  TopExp_Explorer ex;

  for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
    Standard_Real t = BRep_Tool::Tolerance(TopoDS::Edge(ex.Current()));
    if (t > myTolerance) myTolerance = t;
  }

  // search an existing surface
  ex.Init(S,TopAbs_EDGE);
  if (!ex.More()) return;    // no edges ....

  TopoDS_Edge E = TopoDS::Edge(ex.Current());
  Standard_Real f,l,ff,ll;
  Handle(Geom2d_Curve) PC,PPC;
  Handle(Geom_Surface) SS;
  TopLoc_Location L;
  Standard_Integer i = 0,j;

  // iterate on the surfaces of the first edge
  while ( Standard_True) {
    i++;
    BRep_Tool::CurveOnSurface(E,PC,mySurface,myLocation,f,l,i);
    if (mySurface.IsNull()) {
      break;
    }
    // check the other edges
    for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
      if (!E.IsSame(ex.Current())) {
	j = 0;
	while (Standard_True) {
	  j++;
	  BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
				    PPC,SS,L,ff,ll,j);
	  if (SS.IsNull()) {
	    break;
	  }
	  if (SS == mySurface) {
	    break;
	  }
	  SS.Nullify();
	}

	if (SS.IsNull()) {
	  mySurface.Nullify();
	  break;
	}
      }
    }

    // if OnlyPlane, eval if mySurface is a plane.
    if ( OnlyPlane && !mySurface.IsNull() ) 
    {
      if ( mySurface->IsKind( STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
        mySurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurface)->BasisSurface();
      mySurface = Handle(Geom_Plane)::DownCast(mySurface);
    }

    if (!mySurface.IsNull())
      // if S is e.g. the bottom face of a cylinder, mySurface can be the
      // lateral (cylindrical) face of the cylinder; reject an improper mySurface
      if ( !OnlyClosed || Is2DClosed( S, mySurface ))
        break;
  }

  if (!mySurface.IsNull()) {
    isExisted = Standard_True;
    return;
  }
  //
  // no existing surface, search a plane
  // 07/02/02 akm vvv : (OCC157) changed algorithm
  //                    1. Collect the points along all edges of the shape
  //                       For each point calculate the WEIGHT = sum of
  //                       distances from neighboring points (_only_ same edge)
  //                    2. Minimizing the weighed sum of squared deviations
  //                       compute coefficients of the sought plane.
  
  TColgp_SequenceOfPnt aPoints;
  TColStd_SequenceOfReal aWeight;

  // ======================= Step #1
  for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) 
  {
    BRepAdaptor_Curve c(TopoDS::Edge(ex.Current()));

    Standard_Real dfUf = c.FirstParameter();
    Standard_Real dfUl = c.LastParameter();
    if (IsEqual(dfUf,dfUl)) {
      // Degenerate
      continue;
    }
    Standard_Integer iNbPoints=0;

    // Add the points with weights to the sequences
    switch (c.GetType()) 
    {
    case GeomAbs_BezierCurve:
      {
	// Put all poles for bezier
	Handle(Geom_BezierCurve) GC = c.Bezier();
	Standard_Integer iNbPol = GC->NbPoles();
	if ( iNbPol < 2)
	  // Degenerate
	  continue;
	else
	{
	  Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
	  GC->Poles(aPoles->ChangeArray1());
	  gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
	  Standard_Real dfDistPrev = 0., dfDistNext;
	  for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
	  {
	    if (iPol<iNbPol)
	    {
	      aPoleNext = aPoles->Value(iPol+1);
	      dfDistNext = aPolePrev.Distance(aPoleNext);
	    }
	    else
	      dfDistNext = 0.;
	    aPoints.Append (aPolePrev);
	    aWeight.Append (dfDistPrev+dfDistNext);
	    dfDistPrev = dfDistNext;
	    aPolePrev = aPoleNext;
	  }
	}
      }
      break;
    case GeomAbs_BSplineCurve:
      {
	// Put all poles for bspline
	Handle(Geom_BSplineCurve) GC = c.BSpline();
	Standard_Integer iNbPol = GC->NbPoles();
	if ( iNbPol < 2)
	  // Degenerate
	  continue;
	else
	{
	  Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
	  GC->Poles(aPoles->ChangeArray1());
	  gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
	  Standard_Real dfDistPrev = 0., dfDistNext;
	  for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
	  {
	    if (iPol<iNbPol)
	    {
	      aPoleNext = aPoles->Value(iPol+1);
	      dfDistNext = aPolePrev.Distance(aPoleNext);
	    }
	    else
	      dfDistNext = 0.;
	    aPoints.Append (aPolePrev);
	    aWeight.Append (dfDistPrev+dfDistNext);
	    dfDistPrev = dfDistNext;
	    aPolePrev = aPoleNext;
	  }
	}
      }
      break;

    case GeomAbs_Line:
    case GeomAbs_Circle:
    case GeomAbs_Ellipse:
    case GeomAbs_Hyperbola:
    case GeomAbs_Parabola:
      if (c.GetType() == GeomAbs_Line)
	// Two points on straight segment
	iNbPoints=2;
      else
	// Four points on otheranalitical curves
	iNbPoints=4;
    default:
      { 
	// Put some points on other curves
	if (iNbPoints==0)
	  iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
	Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
	Standard_Integer iPoint;
	Standard_Real dfU;
	gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
	Standard_Real dfDistPrev = 0., dfDistNext;
	for (iPoint=1, dfU=dfUf+dfDelta; 
	     iPoint<=iNbPoints; 
	     iPoint++, dfU+=dfDelta) 
	{
	  if (iPoint<iNbPoints)
	  {
	    aPointNext = c.Value(dfU);
	    dfDistNext = aPointPrev.Distance(aPointNext);
	  }
	  else
	    dfDistNext = 0.;
	  aPoints.Append (aPointPrev);
	  aWeight.Append (dfDistPrev+dfDistNext);
	  dfDistPrev = dfDistNext;
	  aPointPrev = aPointNext;
	}
      } // default:
    } // switch (c.GetType()) ...
  } // for (ex.Init(S,TopAbs_EDGE); ex.More() && control; ex.Next()) ...
  
  if (aPoints.Length() < 3) {
    return;
  }

  // ======================= Step #2
  myLocation.Identity();
  Standard_Integer iPoint;
  math_Matrix aMat (1,3,1,3, 0.);
  math_Vector aVec (1,3, 0.);
  // Find the barycenter and normalize weights 
  Standard_Real dfMaxWeight=0.;
  gp_XYZ aBaryCenter(0.,0.,0.);
  Standard_Real dfSumWeight=0.;
  for (iPoint=1; iPoint<=aPoints.Length(); iPoint++)  {
    Standard_Real dfW = aWeight(iPoint);
    aBaryCenter += dfW*aPoints(iPoint).XYZ();
    dfSumWeight += dfW;
    if (dfW > dfMaxWeight) {
      dfMaxWeight = dfW;
    }
  }
  aBaryCenter /= dfSumWeight;

  // Fill the matrix and the right vector
  for (iPoint=1; iPoint<=aPoints.Length(); iPoint++)  {
    gp_XYZ p=aPoints(iPoint).XYZ()-aBaryCenter;
    Standard_Real w=aWeight(iPoint)/dfMaxWeight;
    aMat(1,1)+=w*p.X()*p.X(); 
        aMat(1,2)+=w*p.X()*p.Y(); 
            aMat(1,3)+=w*p.X()*p.Z();
    aMat(2,1)+=w*p.Y()*p.X();  
        aMat(2,2)+=w*p.Y()*p.Y();  
            aMat(2,3)+=w*p.Y()*p.Z();
    aMat(3,1)+=w*p.Z()*p.X();  
        aMat(3,2)+=w*p.Z()*p.Y(); 
            aMat(3,3)+=w*p.Z()*p.Z();
    aVec(1) -= w*p.X();
    aVec(2) -= w*p.Y();
    aVec(3) -= w*p.Z();
  }

  // Solve the system of equations to get plane coefficients
  math_Gauss aSolver(aMat);
  Standard_Boolean isSolved = aSolver.IsDone();
  //
  //  let us be more tolerant (occ415)
  Standard_Real dfDist = RealLast();
  Handle(Geom_Plane) aPlane;
  //
  if (isSolved)  {
    aSolver.Solve(aVec);
    if (aVec.Norm2()<gp::Resolution()) {
      isSolved = Standard_False;
    }
  }
  //
  if (isSolved)  {
    aPlane = new Geom_Plane(aBaryCenter,gp_Dir(aVec(1),aVec(2),aVec(3)));
    dfDist = Controle (aPoints, aPlane);
  }
  //
  if (!isSolved || myTolerance < dfDist)  {
    gp_Pnt aFirstPnt=aPoints(1);
    for (iPoint=2; iPoint<=aPoints.Length(); iPoint++)  {
      gp_Vec aDir(aFirstPnt,aPoints(iPoint));
      Standard_Real dfSide=aDir.Magnitude();
      if (dfSide<myTolerance) {
	continue; // degeneration
      }
      for (Standard_Integer iP1=iPoint+1; iP1<=aPoints.Length(); iP1++) {
	gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
	if (aCross.Magnitude() > dfSide*myTolerance) {
	  Handle(Geom_Plane) aPlane2 = new Geom_Plane(aFirstPnt, aCross);
	  Standard_Real dfDist2 = Controle (aPoints, aPlane2);
	  if (dfDist2 < myTolerance)  {
	    myTolReached = dfDist2;
	    mySurface = aPlane2;
	    return;
	  }
	  if (dfDist2 < dfDist)  {
	    dfDist = dfDist2;
	    aPlane = aPlane2;
	  }
	}
      }
    }
  }
  //
  //XXf
  //static Standard_Real weakness = 5.0;
  Standard_Real weakness = 5.0;
  //XXf
  if(dfDist <= myTolerance || dfDist < myTolerance*weakness && Tol<0) { 
    //XXf 
    //myTolReached = dfDist;
    //XXt
    mySurface = aPlane;
  }
  //XXf
  myTolReached = dfDist;
  //XXt
}
//=======================================================================
//function : Found
//purpose  : 
//=======================================================================
Standard_Boolean BRepLib_FindSurface::Found() const 
{
  return !mySurface.IsNull();
}
//=======================================================================
//function : Surface
//purpose  : 
//=======================================================================
Handle(Geom_Surface) BRepLib_FindSurface::Surface() const 
{
  return mySurface;
}
//=======================================================================
//function : Tolerance
//purpose  : 
//=======================================================================
Standard_Real BRepLib_FindSurface::Tolerance() const 
{
  return myTolerance;
}
//=======================================================================
//function : ToleranceReached
//purpose  : 
//=======================================================================
Standard_Real BRepLib_FindSurface::ToleranceReached() const 
{
  return myTolReached;
}
//=======================================================================
//function : Existed
//purpose  : 
//=======================================================================
Standard_Boolean BRepLib_FindSurface::Existed() const 
{
  return isExisted;
}
//=======================================================================
//function : Location
//purpose  : 
//=======================================================================
TopLoc_Location BRepLib_FindSurface::Location() const
{
  return myLocation;
}
Commit	Line	Data
b311480e	1	// Created on: 1994-07-22
	2	// Created by: Remi LEQUETTE
	3	// Copyright (c) 1994-1999 Matra Datavision
	4	// Copyright (c) 1999-2012 OPEN CASCADE SAS
	5	//
	6	// The content of this file is subject to the Open CASCADE Technology Public
	7	// License Version 6.5 (the "License"). You may not use the content of this file
	8	// except in compliance with the License. Please obtain a copy of the License
	9	// at http://www.opencascade.org and read it completely before using this file.
	10	//
	11	// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
	12	// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
	13	//
	14	// The Original Code and all software distributed under the License is
	15	// distributed on an "AS IS" basis, without warranty of any kind, and the
	16	// Initial Developer hereby disclaims all such warranties, including without
	17	// limitation, any warranties of merchantability, fitness for a particular
	18	// purpose or non-infringement. Please see the License for the specific terms
	19	// and conditions governing the rights and limitations under the License.
	20
7fd59977	21
	22	#include <BRepLib_FindSurface.ixx>
	23
	24	#include <Precision.hxx>
	25	#include <math_Matrix.hxx>
	26	#include <math_Vector.hxx>
	27	#include <math_Gauss.hxx>
	28
	29	#include <gp_Lin.hxx>
	30	#include <gp_Circ.hxx>
	31	#include <gp_Elips.hxx>
	32	#include <gp_Hypr.hxx>
	33	#include <gp_Parab.hxx>
	34	#include <gp_Ax2.hxx>
	35	#include <gp_Ax3.hxx>
	36	#include <gp_Vec.hxx>
	37	#include <TColgp_SequenceOfPnt.hxx>
	38	#include <TColStd_SequenceOfReal.hxx>
	39	#include <TColgp_Array1OfPnt.hxx>
	40	#include <TColgp_HArray1OfPnt.hxx>
	41	#include <Geom_Plane.hxx>
	42
7fd59977	43	#include <BRepAdaptor_Curve.hxx>
0f5cd7d5	44	#include <BRepAdaptor_Surface.hxx>
	45	#include <BRepLib_MakeFace.hxx>
	46	#include <BRepTools_WireExplorer.hxx>
	47	#include <BRep_Tool.hxx>
	48	#include <TopExp.hxx>
	49	#include <TopExp_Explorer.hxx>
	50	#include <TopoDS_Vertex.hxx>
	51	#include <TopoDS_Wire.hxx>
	52	#include <TopoDS.hxx>
7fd59977	53
	54	#include <GeomLib.hxx>
	55	#include <Geom2d_Curve.hxx>
	56	#include <Geom_BezierCurve.hxx>
	57	#include <Geom_BSplineCurve.hxx>
0f5cd7d5	58	#include <Geom_RectangularTrimmedSurface.hxx>
0f5cd7d5	59	#include <Standard_ErrorHandler.hxx>
7fd59977	60
	61	//=======================================================================
	62	//function : Controle
	63	//purpose :
	64	//=======================================================================
	65	static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
	66	const Handle(Geom_Plane)& thePlane)
	67	{
	68	Standard_Real dfMaxDist=0.;
	69	Standard_Real a,b,c,d, dist;
	70	Standard_Integer ii;
	71	thePlane->Coefficients(a,b,c,d);
	72	for (ii=1; ii<=thePoints.Length(); ii++) {
	73	const gp_XYZ& xyz = thePoints(ii).XYZ();
	74	dist = Abs(axyz.X() + bxyz.Y() + c*xyz.Z() + d);
	75	if (dist > dfMaxDist)
	76	dfMaxDist = dist;
	77	}
	78
	79	return dfMaxDist;
	80	}
	81
0f5cd7d5	82	//=======================================================================
	83	//function : Is2DConnected
	84	//purpose : Return true if the last vertex of theEdge1 coincides with
	85	// the first vertex of theEdge2 in parametric space of theFace
	86	//=======================================================================
	87
	88	inline static Standard_Boolean Is2DConnected( const TopoDS_Edge& theEdge1,
	89	const TopoDS_Edge& theEdge2,
	90	const TopoDS_Face& theFace )
	91	{
	92	Standard_Real f,l;
	93	Handle(Geom2d_Curve) aCurve;
	94	gp_Pnt2d p1, p2;
	95
	96	// get 2D points
	97	aCurve = BRep_Tool::CurveOnSurface( theEdge1, theFace,f,l );
	98	p1 = aCurve->Value( theEdge1.Orientation() == TopAbs_FORWARD ? l : f );
	99	aCurve = BRep_Tool::CurveOnSurface( theEdge2, theFace,f,l );
	100	p2 = aCurve->Value( theEdge2.Orientation() == TopAbs_FORWARD ? f : l );
	101
	102	// compare 2D points
	103	BRepAdaptor_Surface aSurface( theFace );
	104	TopoDS_Vertex aV = TopExp::FirstVertex( theEdge2, /CumOri=/Standard_True );
	105	Standard_Real tol3D = BRep_Tool::Tolerance( aV );
	106	Standard_Real tol2D = aSurface.UResolution( tol3D ) + aSurface.VResolution( tol3D );
	107	Standard_Real dist2 = p1.SquareDistance( p2 );
	108	return dist2 < tol2D * tol2D;
	109	}
	110
	111	//=======================================================================
	112	//function : Is2DClosed
	113	//purpose : Return true if edges of theShape form a closed wire in
	114	// parametric space of theSurface
	115	//=======================================================================
	116
	117	static Standard_Boolean Is2DClosed( const TopoDS_Shape& theShape,
	118	const Handle(Geom_Surface)& theSurface)
	119	{
	120	try
	121	{
	122	// get a wire theShape
	123	TopExp_Explorer aWireExp( theShape, TopAbs_WIRE );
	124	if ( !aWireExp.More() )
	125	return Standard_False;
	126	TopoDS_Wire aWire = TopoDS::Wire( aWireExp.Current() );
	127	// a tmp face
	128	TopoDS_Face aTmpFace = BRepLib_MakeFace( theSurface, Precision::PConfusion() );
	129
	130	// check topological closeness using wire explorer, if the wire is not closed
	131	// the 1st and the last vertices of wire are different
	132	BRepTools_WireExplorer aWireExplorer( aWire, aTmpFace );
	133	if ( !aWireExplorer.More())
	134	return Standard_False;
	135	// remember the 1st and the last edges of aWire
	136	TopoDS_Edge aFisrtEdge = aWireExplorer.Current(), aLastEdge = aFisrtEdge;
	137	// check if edges connected topologically (that is assured by BRepTools_WireExplorer)
	138	// are connected in 2D
	139	TopoDS_Edge aPrevEdge = aFisrtEdge;
	140	for ( aWireExplorer.Next(); aWireExplorer.More(); aWireExplorer.Next() )
	141	{
	142	aLastEdge = aWireExplorer.Current();
	143	if ( !Is2DConnected( aPrevEdge, aLastEdge, aTmpFace ))
	144	return false;
	145	aPrevEdge = aLastEdge;
146	}
147	// wire is closed if ( 1st vertex of aFisrtEdge ) ==
148	// ( last vertex of aLastEdge ) in 2D
149	TopoDS_Vertex aV1 = TopExp::FirstVertex( aFisrtEdge, /CumOri=/Standard_True );
150	TopoDS_Vertex aV2 = TopExp::LastVertex( aLastEdge, /CumOri=/Standard_True );
151	return ( aV1.IsSame( aV2 ) && Is2DConnected( aLastEdge, aFisrtEdge, aTmpFace ));
152	}
153	catch ( Standard_Failure )
154	{
155	return Standard_False;
156	}
157	}
158
7fd59977	159	//=======================================================================
	160	//function : BRepLib_FindSurface
	161	//purpose :
	162	//=======================================================================
	163	BRepLib_FindSurface::BRepLib_FindSurface()
	164	{
	165	}
	166	//=======================================================================
	167	//function : BRepLib_FindSurface
	168	//purpose :
	169	//=======================================================================
	170	BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape& S,
	171	const Standard_Real Tol,
0f5cd7d5	172	const Standard_Boolean OnlyPlane,
0f5cd7d5	173	const Standard_Boolean OnlyClosed)
7fd59977	174	{
0f5cd7d5	175	Init(S,Tol,OnlyPlane,OnlyClosed);
7fd59977	176	}
	177	//=======================================================================
	178	//function : Init
	179	//purpose :
	180	//=======================================================================
	181	void BRepLib_FindSurface::Init(const TopoDS_Shape& S,
	182	const Standard_Real Tol,
0f5cd7d5	183	const Standard_Boolean OnlyPlane,
0f5cd7d5	184	const Standard_Boolean OnlyClosed)
7fd59977	185	{
	186	myTolerance = Tol;
	187	myTolReached = 0.;
	188	isExisted = Standard_False;
	189	myLocation.Identity();
	190	mySurface.Nullify();
	191
	192	// compute the tolerance
	193	TopExp_Explorer ex;
	194
	195	for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
	196	Standard_Real t = BRep_Tool::Tolerance(TopoDS::Edge(ex.Current()));
	197	if (t > myTolerance) myTolerance = t;
	198	}
	199
	200	// search an existing surface
	201	ex.Init(S,TopAbs_EDGE);
	202	if (!ex.More()) return; // no edges ....
	203
	204	TopoDS_Edge E = TopoDS::Edge(ex.Current());
	205	Standard_Real f,l,ff,ll;
	206	Handle(Geom2d_Curve) PC,PPC;
	207	Handle(Geom_Surface) SS;
	208	TopLoc_Location L;
	209	Standard_Integer i = 0,j;
	210
	211	// iterate on the surfaces of the first edge
	212	while ( Standard_True) {
	213	i++;
	214	BRep_Tool::CurveOnSurface(E,PC,mySurface,myLocation,f,l,i);
	215	if (mySurface.IsNull()) {
	216	break;
	217	}
	218	// check the other edges
	219	for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
	220	if (!E.IsSame(ex.Current())) {
	221	j = 0;
	222	while (Standard_True) {
	223	j++;
	224	BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
	225	PPC,SS,L,ff,ll,j);
	226	if (SS.IsNull()) {
	227	break;
	228	}
	229	if (SS == mySurface) {
	230	break;
	231	}
	232	SS.Nullify();
	233	}
	234
	235	if (SS.IsNull()) {
	236	mySurface.Nullify();
	237	break;
	238	}
	239	}
	240	}
	241
	242	// if OnlyPlane, eval if mySurface is a plane.
	243	if ( OnlyPlane && !mySurface.IsNull() )
0f5cd7d5	244	{
	245	if ( mySurface->IsKind( STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
	246	mySurface = Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurface)->BasisSurface();
7fd59977	247	mySurface = Handle(Geom_Plane)::DownCast(mySurface);
0f5cd7d5	248	}
7fd59977	249
0f5cd7d5	250	if (!mySurface.IsNull())
	251	// if S is e.g. the bottom face of a cylinder, mySurface can be the
	252	// lateral (cylindrical) face of the cylinder; reject an improper mySurface
	253	if ( !OnlyClosed \|\| Is2DClosed( S, mySurface ))
	254	break;
7fd59977	255	}
	256
	257	if (!mySurface.IsNull()) {
	258	isExisted = Standard_True;
	259	return;
	260	}
	261	//
	262	// no existing surface, search a plane
	263	// 07/02/02 akm vvv : (OCC157) changed algorithm
	264	// 1. Collect the points along all edges of the shape
	265	// For each point calculate the WEIGHT = sum of
	266	// distances from neighboring points (_only_ same edge)
	267	// 2. Minimizing the weighed sum of squared deviations
	268	// compute coefficients of the sought plane.
	269
	270	TColgp_SequenceOfPnt aPoints;
	271	TColStd_SequenceOfReal aWeight;
	272
	273	// ======================= Step #1
	274	for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next())
	275	{
	276	BRepAdaptor_Curve c(TopoDS::Edge(ex.Current()));
	277
	278	Standard_Real dfUf = c.FirstParameter();
	279	Standard_Real dfUl = c.LastParameter();
	280	if (IsEqual(dfUf,dfUl)) {
	281	// Degenerate
	282	continue;
	283	}
	284	Standard_Integer iNbPoints=0;
	285
	286	// Add the points with weights to the sequences
	287	switch (c.GetType())
	288	{
	289	case GeomAbs_BezierCurve:
	290	{
	291	// Put all poles for bezier
	292	Handle(Geom_BezierCurve) GC = c.Bezier();
	293	Standard_Integer iNbPol = GC->NbPoles();
	294	if ( iNbPol < 2)
	295	// Degenerate
	296	continue;
	297	else
	298	{
	299	Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
	300	GC->Poles(aPoles->ChangeArray1());
	301	gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
	302	Standard_Real dfDistPrev = 0., dfDistNext;
	303	for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
	304	{
	305	if (iPol<iNbPol)
	306	{
	307	aPoleNext = aPoles->Value(iPol+1);
	308	dfDistNext = aPolePrev.Distance(aPoleNext);
	309	}
	310	else
	311	dfDistNext = 0.;
	312	aPoints.Append (aPolePrev);
	313	aWeight.Append (dfDistPrev+dfDistNext);
	314	dfDistPrev = dfDistNext;
	315	aPolePrev = aPoleNext;
	316	}
	317	}
	318	}
319	break;
320	case GeomAbs_BSplineCurve:
321	{
322	// Put all poles for bspline
323	Handle(Geom_BSplineCurve) GC = c.BSpline();
324	Standard_Integer iNbPol = GC->NbPoles();
325	if ( iNbPol < 2)
326	// Degenerate
327	continue;
328	else
329	{
330	Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
331	GC->Poles(aPoles->ChangeArray1());
332	gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
333	Standard_Real dfDistPrev = 0., dfDistNext;
334	for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
335	{
336	if (iPol<iNbPol)
337	{
338	aPoleNext = aPoles->Value(iPol+1);
339	dfDistNext = aPolePrev.Distance(aPoleNext);
340	}
341	else
342	dfDistNext = 0.;
343	aPoints.Append (aPolePrev);
344	aWeight.Append (dfDistPrev+dfDistNext);
345	dfDistPrev = dfDistNext;
346	aPolePrev = aPoleNext;
347	}
348	}
349	}
350	break;
351
352	case GeomAbs_Line:
353	case GeomAbs_Circle:
354	case GeomAbs_Ellipse:
355	case GeomAbs_Hyperbola:
356	case GeomAbs_Parabola:
357	if (c.GetType() == GeomAbs_Line)
358	// Two points on straight segment
359	iNbPoints=2;
360	else
361	// Four points on otheranalitical curves
362	iNbPoints=4;
363	default:
364	{
365	// Put some points on other curves
366	if (iNbPoints==0)
367	iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
368	Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
369	Standard_Integer iPoint;
370	Standard_Real dfU;
371	gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
372	Standard_Real dfDistPrev = 0., dfDistNext;
373	for (iPoint=1, dfU=dfUf+dfDelta;
374	iPoint<=iNbPoints;
375	iPoint++, dfU+=dfDelta)
376	{
377	if (iPoint<iNbPoints)
378	{
379	aPointNext = c.Value(dfU);
380	dfDistNext = aPointPrev.Distance(aPointNext);
381	}
382	else
383	dfDistNext = 0.;
384	aPoints.Append (aPointPrev);
385	aWeight.Append (dfDistPrev+dfDistNext);
386	dfDistPrev = dfDistNext;
387	aPointPrev = aPointNext;
388	}
389	} // default:
390	} // switch (c.GetType()) ...
391	} // for (ex.Init(S,TopAbs_EDGE); ex.More() && control; ex.Next()) ...
392
393	if (aPoints.Length() < 3) {
394	return;
395	}
396
397	// ======================= Step #2
398	myLocation.Identity();
399	Standard_Integer iPoint;
400	math_Matrix aMat (1,3,1,3, 0.);
401	math_Vector aVec (1,3, 0.);
402	// Find the barycenter and normalize weights
403	Standard_Real dfMaxWeight=0.;
404	gp_XYZ aBaryCenter(0.,0.,0.);
405	Standard_Real dfSumWeight=0.;
406	for (iPoint=1; iPoint<=aPoints.Length(); iPoint++) {
407	Standard_Real dfW = aWeight(iPoint);
408	aBaryCenter += dfW*aPoints(iPoint).XYZ();
409	dfSumWeight += dfW;
410	if (dfW > dfMaxWeight) {
411	dfMaxWeight = dfW;
412	}
413	}
414	aBaryCenter /= dfSumWeight;
415
416	// Fill the matrix and the right vector
417	for (iPoint=1; iPoint<=aPoints.Length(); iPoint++) {
418	gp_XYZ p=aPoints(iPoint).XYZ()-aBaryCenter;
419	Standard_Real w=aWeight(iPoint)/dfMaxWeight;
420	aMat(1,1)+=wp.X()p.X();
421	aMat(1,2)+=wp.X()p.Y();
422	aMat(1,3)+=wp.X()p.Z();
423	aMat(2,1)+=wp.Y()p.X();
424	aMat(2,2)+=wp.Y()p.Y();
425	aMat(2,3)+=wp.Y()p.Z();
426	aMat(3,1)+=wp.Z()p.X();
427	aMat(3,2)+=wp.Z()p.Y();
428	aMat(3,3)+=wp.Z()p.Z();
429	aVec(1) -= w*p.X();
430	aVec(2) -= w*p.Y();
431	aVec(3) -= w*p.Z();
432	}
433
434	// Solve the system of equations to get plane coefficients
435	math_Gauss aSolver(aMat);
436	Standard_Boolean isSolved = aSolver.IsDone();
437	//
438	// let us be more tolerant (occ415)
439	Standard_Real dfDist = RealLast();
440	Handle(Geom_Plane) aPlane;
441	//
442	if (isSolved) {
443	aSolver.Solve(aVec);
444	if (aVec.Norm2()<gp::Resolution()) {
445	isSolved = Standard_False;
446	}
447	}
448	//
449	if (isSolved) {
450	aPlane = new Geom_Plane(aBaryCenter,gp_Dir(aVec(1),aVec(2),aVec(3)));
451	dfDist = Controle (aPoints, aPlane);
452	}
453	//
454	if (!isSolved \|\| myTolerance < dfDist) {
455	gp_Pnt aFirstPnt=aPoints(1);
456	for (iPoint=2; iPoint<=aPoints.Length(); iPoint++) {
457	gp_Vec aDir(aFirstPnt,aPoints(iPoint));
458	Standard_Real dfSide=aDir.Magnitude();
459	if (dfSide<myTolerance) {
460	continue; // degeneration
461	}
462	for (Standard_Integer iP1=iPoint+1; iP1<=aPoints.Length(); iP1++) {
463	gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
464	if (aCross.Magnitude() > dfSide*myTolerance) {
465	Handle(Geom_Plane) aPlane2 = new Geom_Plane(aFirstPnt, aCross);
466	Standard_Real dfDist2 = Controle (aPoints, aPlane2);
467	if (dfDist2 < myTolerance) {
468	myTolReached = dfDist2;
469	mySurface = aPlane2;
470	return;
471	}
472	if (dfDist2 < dfDist) {
473	dfDist = dfDist2;
474	aPlane = aPlane2;
475	}
476	}
477	}
478	}
479	}
480	//
481	//XXf
482	//static Standard_Real weakness = 5.0;
483	Standard_Real weakness = 5.0;
484	//XXf
485	if(dfDist <= myTolerance \|\| dfDist < myTolerance*weakness && Tol<0) {
486	//XXf
487	//myTolReached = dfDist;
488	//XXt
489	mySurface = aPlane;
490	}
491	//XXf
492	myTolReached = dfDist;
493	//XXt
494	}
495	//=======================================================================
496	//function : Found
497	//purpose :
498	//=======================================================================
499	Standard_Boolean BRepLib_FindSurface::Found() const
500	{
501	return !mySurface.IsNull();
502	}
503	//=======================================================================
504	//function : Surface
505	//purpose :
506	//=======================================================================
507	Handle(Geom_Surface) BRepLib_FindSurface::Surface() const
508	{
509	return mySurface;
510	}
511	//=======================================================================
512	//function : Tolerance
513	//purpose :
514	//=======================================================================
515	Standard_Real BRepLib_FindSurface::Tolerance() const
516	{
517	return myTolerance;
518	}
519	//=======================================================================
520	//function : ToleranceReached
521	//purpose :
522	//=======================================================================
523	Standard_Real BRepLib_FindSurface::ToleranceReached() const
524	{
525	return myTolReached;
526	}
527	//=======================================================================
528	//function : Existed
529	//purpose :
530	//=======================================================================
531	Standard_Boolean BRepLib_FindSurface::Existed() const
532	{
533	return isExisted;
534	}
535	//=======================================================================
536	//function : Location
537	//purpose :
538	//=======================================================================
539	TopLoc_Location BRepLib_FindSurface::Location() const
540	{
541	return myLocation;
542	}