[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-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and / or modify it
// under the terms of the GNU Lesser General Public version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#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 Handle(Geom_Surface)& theSurface,
					     const TopLoc_Location& theLocation)
{
  Standard_Real f,l;
  //TopLoc_Location aLoc;
  Handle(Geom2d_Curve) aCurve;
  gp_Pnt2d p1, p2;

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

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