[occt.git] / src / ShapeCustom / ShapeCustom_Surface.cxx

// 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 License 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.

//abv 06.01.99 fix of misprint
//:p6 abv 26.02.99: make ConvertToPeriodic() return Null if nothing done

#include <ElSLib.hxx>
#include <Geom_BezierSurface.hxx>
#include <Geom_BSplineSurface.hxx>
#include <Geom_ConicalSurface.hxx>
#include <Geom_Curve.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_Plane.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_Surface.hxx>
#include <Geom_ToroidalSurface.hxx>
#include <GeomAbs_SurfaceType.hxx>
#include <GeomAdaptor_HSurface.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <gp_Ax3.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Pln.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <ShapeAnalysis_Geom.hxx>
#include <ShapeAnalysis_Surface.hxx>
#include <ShapeCustom_Surface.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array2OfPnt.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_Array2OfReal.hxx>

//=======================================================================
//function : ShapeCustom_Surface
//purpose  : 
//=======================================================================
ShapeCustom_Surface::ShapeCustom_Surface() : myGap (0)
{
}

//=======================================================================
//function : ShapeCustom_Surface
//purpose  : 
//=======================================================================

ShapeCustom_Surface::ShapeCustom_Surface (const Handle(Geom_Surface)& S)
     : myGap (0)
{
  Init ( S );
}

//=======================================================================
//function : Init
//purpose  : 
//=======================================================================

void ShapeCustom_Surface::Init (const Handle(Geom_Surface)& S) 
{
  mySurf = S;
}

//=======================================================================
//function : ConvertToAnalytical
//purpose  : 
//=======================================================================

Handle(Geom_Surface) ShapeCustom_Surface::ConvertToAnalytical (const Standard_Real tol,
							       const Standard_Boolean substitute) 
{
  Handle(Geom_Surface) newSurf;

  Standard_Integer nUP, nVP, nCP, i, j , UDeg, VDeg;
  Standard_Real U1, U2, V1, V2, C1, C2, DU, DV, U=0, V=0;
  Handle(Geom_Curve) iso;
  Standard_Boolean uClosed = Standard_True;

  // seuls cas traites : BSpline et Bezier
  Handle(Geom_BSplineSurface) theBSplneS = 
    Handle(Geom_BSplineSurface)::DownCast(mySurf);
  if (theBSplneS.IsNull()) {
    Handle(Geom_BezierSurface) theBezierS = 
      Handle(Geom_BezierSurface)::DownCast(mySurf);
    if (!theBezierS.IsNull()) {    // Bezier :
      nUP  = theBezierS->NbUPoles();   
      nVP  = theBezierS->NbVPoles();   
      UDeg = theBezierS->UDegree();
      VDeg = theBezierS->VDegree();
    }
    else return newSurf;           // non reconnu : terminus
  }
  else {                           // BSpline :
    nUP  = theBSplneS->NbUPoles();
    nVP  = theBSplneS->NbVPoles();
    UDeg = theBSplneS->UDegree();
    VDeg = theBSplneS->VDegree();
  }


  mySurf->Bounds(U1, U2, V1, V2);
//  mySurf->Bounds(U1, U2, V1, V2);
  TColgp_Array1OfPnt p1(1, 3), p2(1, 3), p3(1, 3);
  TColStd_Array1OfReal R(1,3);
  gp_Pnt origPnt, resPnt;  
  gp_Vec origD1U, resD1U, resD1V;
    
  Standard_Boolean aCySpCo = Standard_False;
  Standard_Boolean aToroid = Standard_False;
  Standard_Boolean aPlanar = Standard_False;

  if (nUP == 2 && nVP == 2) {
    if (UDeg == 1 && VDeg == 1) aPlanar = Standard_True;
  } else if (mySurf->IsUClosed()) {    // VRAI IsUClosed
    if (mySurf->IsVClosed())   aToroid = Standard_True;
    else                        aCySpCo = Standard_True;
  } else {
    if(mySurf->IsVClosed()) {    // VRAI IsVClosed
      aCySpCo = Standard_True;
      uClosed = Standard_False;
    }
  }

  if (aPlanar) {
//    NearestPlane ...
    TColgp_Array1OfPnt Pnts(1,4);
    Pnts.SetValue(1,mySurf->Value(U1,V1));
    Pnts.SetValue(2,mySurf->Value(U2,V1));
    Pnts.SetValue(3,mySurf->Value(U1,V2));
    Pnts.SetValue(4,mySurf->Value(U2,V2));
    gp_Pln aPln;// Standard_Real Dmax;
    Standard_Integer It = ShapeAnalysis_Geom::NearestPlane (Pnts,aPln,myGap/*Dmax*/);

//  ICI, on fabrique le plan, et zou
    if (It == 0 || myGap/*Dmax*/ > tol) return newSurf;    //  pas un plan

//    IL RESTE a verifier l orientation ...
//    On regarde sur chaque surface les vecteurs P(U0->U1),P(V0->V1)
//    On prend la normale : les deux normales doivent etre dans le meme sens
//    Sinon, inverser la normale (pas le Pln entier !) et refaire la Plane
    newSurf = new Geom_Plane (aPln);
    gp_Vec uold (Pnts(1),Pnts(2));
    gp_Vec vold (Pnts(1),Pnts(3));
    gp_Vec nold = uold.Crossed (vold);
    gp_Vec unew (newSurf->Value(U1,V1), newSurf->Value(U2,V1));
    gp_Vec vnew (newSurf->Value(U1,V1), newSurf->Value(U1,V2));
    gp_Vec nnew = unew.Crossed (vnew);
    if (nold.Dot (nnew) < 0.0) {
      gp_Ax3 ax3 = aPln.Position();
      ax3.ZReverse();
      ax3.XReverse();
      aPln = gp_Pln (ax3);
      newSurf = new Geom_Plane (aPln);
    }

    if (substitute) {
      Init (newSurf);
    }
    return newSurf;


  } else if (aCySpCo) {
    if (!uClosed) {
      C1 = U1; C2 = U2;
      U1 = V1; U2 = V2;
      V1 = C1; V2 = C2;
      nCP = nUP; nUP = nVP; nVP = nCP;
    }
  
    for (i=1; i<=3; i++) {
      if      (i==1) V = V1;
      else if (i==2) V = V2;
      else if (i==3) V = 0.5*(V1+V2);

      if(uClosed) iso = mySurf->VIso(V);
      else        iso = mySurf->UIso(V);

      iso->D0(U1, p1(i)); 
      iso->D0(0.5*(U1+U2), p2(i));
      p3(i).SetCoord(0.5*(p1(i).X()+p2(i).X()), 
		     0.5*(p1(i).Y()+p2(i).Y()), 
		     0.5*(p1(i).Z()+p2(i).Z()));
      R(i) = p3(i).Distance(p1(i));
//	std::cout<<"sphere, i="<<i<<" V="<<V<<" R="<<R(i)<<" p1="<<p1(i).X()<<","<<p1(i).Y()<<","<<p1(i).Z()<<" p2="<<p2(i).X()<<","<<p2(i).Y()<<","<<p2(i).Z()<<" p3="<<p3(i).X()<<","<<p3(i).Y()<<","<<p3(i).Z()<<std::endl;
    }
      
    iso->D1 (0.,origPnt,origD1U);
    gp_Vec xVec(p3(3), p1(3));
    gp_Vec aVec(p3(1), p3(2));
//      gp_Dir xDir(xVec);  ne sert pas. Null si R3 = 0
    gp_Dir aDir(aVec);
    gp_Ax3 aAx3 (p3(1),aDir,xVec);
    //  CKY  3-FEV-1997 : verification du sens de description
    //gp_Dir AXY = aAx3.YDirection(); // AXY not used (skl)
    if (aAx3.YDirection().Dot (origD1U) < 0) {
#ifdef OCCT_DEBUG
      std::cout<<" Surface Analytique : sens a inverser"<<std::endl;
#endif
      aAx3.YReverse();    //  mais X reste !
    }
      
    if (nVP > 2) {
      if ((Abs(R(1)) < tol) && 
	  (Abs(R(2)) < tol) && 
	  (Abs(R(3)) > tol)) {
// deja fait	  gp_Ax3 aAx3(p3(1), aDir, xVec);
          //gp_Ax3 aAx3(p3(3), aDir);
	Handle(Geom_SphericalSurface) anObject = 
	  new Geom_SphericalSurface(aAx3, R(3));
	if (!uClosed) anObject->UReverse();
	newSurf = anObject;         
      }
    }
    else if (nVP == 2) {

// deja fait	gp_Ax3 aAx3(p3(1), aDir, xVec); 
        //gp_Ax3 aAx3(p3(1), aDir);

      if (Abs(R(2)-R(1)) < tol) {
	Handle(Geom_CylindricalSurface) anObject = 
	  new Geom_CylindricalSurface(aAx3, R(1));
	if (!uClosed) anObject->UReverse();
	newSurf = anObject;
      }
      else {
	gp_Vec aVec2(p1(1), p1(2));
	Standard_Real angle = aVec.Angle(aVec2);
	if (R(1) < R(2)) {
	  Handle(Geom_ConicalSurface) anObject = 
	    new Geom_ConicalSurface(aAx3, angle, R(1));
	  //if (!uClosed) anObject->UReverse();
	  anObject->UReverse();
	  newSurf = anObject;
	}
	else {
	  aDir.Reverse();
	  gp_Vec anotherXVec(p3(2), p1(2));
	  gp_Dir anotherXDir(anotherXVec);
	  gp_Ax3 anotherAx3(p3(2), aDir, anotherXDir);
	  Handle(Geom_ConicalSurface) anObject = 
	    new Geom_ConicalSurface(anotherAx3, angle, R(2));
	  //if (!uClosed) anObject->UReverse();
	  anObject->UReverse();
	  newSurf = anObject;
	}
      }
    }
  }
  else if (aToroid) {
    // test by iso U and isoV
    Standard_Boolean  isFound = Standard_False;
    for (j=1; (j<=2) && !isFound; j++) {
      if (j==2) {
	C1 = U1; C2 = U2;
	U1 = V1; U2 = V2;
	V1 = C1; V2 = C2;
      }
      for (i=1; i<=3; i++) {
	if      (i==1) U = U1;
	else if (i==2) U = 0.5*(U1+U2);
	else if (i==3) U = 0.25*(U1+U2);

	iso = mySurf->UIso(U);

	iso->D0(V1, p1(i)); 
	iso->D0(0.5*(V1+V2), p2(i));
	p3(i).SetCoord(0.5*(p1(i).X()+p2(i).X()), 
		       0.5*(p1(i).Y()+p2(i).Y()), 
		       0.5*(p1(i).Z()+p2(i).Z()));
	R(i) = p3(i).Distance(p1(i));
      }
      if ((Abs(R(1)-R(2))< tol) && 
	  (Abs(R(1)-R(3))< tol)) {
	gp_Pnt p10(0.5*(p3(1).X()+p3(2).X()), 
		   0.5*(p3(1).Y()+p3(2).Y()), 
		   0.5*(p3(1).Z()+p3(2).Z()));
	gp_Vec aVec(p10, p3(1));
	gp_Vec aVec2(p10, p3(3));
	Standard_Real RR1 = R(1), RR2 = R(2), RR3;
	aVec ^= aVec2;

	if (aVec.Magnitude() <= gp::Resolution()) aVec.SetCoord(0., 0., 1.);

	gp_Dir aDir(aVec);

	gp_Ax3 aAx3(p10, aDir);
	RR1 = p10.Distance(p3(1));
//          modif empirique (pourtant NON DEMONTREE) : inverser roles RR1,RR2
//          CKY, 24-JAN-1997
	if (RR1 < RR2)  {  RR3 = RR1; RR1 = RR2; RR2 = RR3;  }
	Handle(Geom_ToroidalSurface) anObject = 
	  new Geom_ToroidalSurface(aAx3, RR1, RR2);
	if (j==2) anObject->UReverse();
	anObject->D1 (0.,0.,resPnt,resD1U,resD1V);
#ifdef OCCT_DEBUG
	if (resD1U.Dot(origD1U) < 0 && j != 2)
	  std::cout<<" Tore a inverser !"<<std::endl;
#endif
	newSurf = anObject;
	isFound = Standard_True;
      }
    }
  }
  if (newSurf.IsNull()) return newSurf;
  
  //---------------------------------------------------------------------
  //                 verification
  //---------------------------------------------------------------------
  
  Handle(GeomAdaptor_HSurface) NHS = new GeomAdaptor_HSurface (newSurf);
  GeomAdaptor_Surface& SurfAdapt = NHS->ChangeSurface();

  const Standard_Integer NP = 21;
  Standard_Real S = 0., T = 0.;  // U,V deja fait
  gp_Pnt P3d, P3d2;
  Standard_Boolean onSurface = Standard_True;

  Standard_Real dis;  myGap = 0.;

  DU = (U2-U1)/(NP-1);
  DV = (V2-V1)/(NP-1);
  for (j=1; (j<=NP) && onSurface; j++) {
    V = V1 + DV*(j-1);

    if(uClosed) iso = mySurf->VIso(V);
    else        iso = mySurf->UIso(V);

    for (i=1; i<=NP; i++) {
      U = U1 + DU*(i-1);
      iso->D0(U, P3d);
      switch (SurfAdapt.GetType()){    
	  
	case GeomAbs_Cylinder :
	  {
	    gp_Cylinder Cylinder = SurfAdapt.Cylinder();
	    ElSLib::Parameters( Cylinder, P3d, S, T);
	    break;
	  }
	case GeomAbs_Cone :
	  {
	    gp_Cone Cone = SurfAdapt.Cone();
	    ElSLib::Parameters( Cone, P3d, S, T);
	    break;
	  }
	case GeomAbs_Sphere :
	  {
	    gp_Sphere Sphere = SurfAdapt.Sphere();
	    ElSLib::Parameters( Sphere, P3d, S, T);
	    break;
	  }
	case GeomAbs_Torus :
	  {
	    gp_Torus Torus = SurfAdapt.Torus();
	    ElSLib::Parameters( Torus, P3d, S, T);
	    break;
	  }
	default:
	  break;
      }

      newSurf->D0(S, T, P3d2);

      dis = P3d.Distance(P3d2);
      if (dis > myGap) myGap = dis;

      if (dis > tol) {
	onSurface = Standard_False;
	newSurf.Nullify();
          // The presumption is rejected
	break;
      }
    }
  }
  if (substitute && !NHS.IsNull()) {
    Init (newSurf);
  }
  return newSurf;
}
 
//%pdn 30 Nov 98: converting bspline surfaces with degree+1 at ends to periodic
// UKI60591, entity 48720
Handle(Geom_Surface) ShapeCustom_Surface::ConvertToPeriodic (const Standard_Boolean substitute,
                                                             const Standard_Real preci)
{
  Handle(Geom_Surface) newSurf;
  Handle(Geom_BSplineSurface) BSpl = Handle(Geom_BSplineSurface)::DownCast(mySurf);
  if (BSpl.IsNull()) return newSurf;
  
  ShapeAnalysis_Surface sas(mySurf);
  Standard_Boolean uclosed = sas.IsUClosed(preci);
  Standard_Boolean vclosed = sas.IsVClosed(preci);
  
  if ( ! uclosed && ! vclosed ) return newSurf;
  
  Standard_Boolean converted = Standard_False; //:p6

  if ( uclosed && ! BSpl->IsUPeriodic() && BSpl->NbUPoles() >3 ) {
    Standard_Boolean set = Standard_True;
    // if degree+1 at ends, first change it to 1 by rearranging knots
    if ( BSpl->UMultiplicity(1) == BSpl->UDegree() + 1 &&
         BSpl->UMultiplicity(BSpl->NbUKnots()) == BSpl->UDegree() + 1 ) {
      Standard_Integer nbUPoles = BSpl->NbUPoles();
      Standard_Integer nbVPoles = BSpl->NbVPoles();
      TColgp_Array2OfPnt oldPoles(1,nbUPoles,1,nbVPoles);
      TColStd_Array2OfReal oldWeights(1,nbUPoles,1,nbVPoles);
      Standard_Integer nbUKnots = BSpl->NbUKnots();
      Standard_Integer nbVKnots = BSpl->NbVKnots();
      TColStd_Array1OfReal oldUKnots(1,nbUKnots);
      TColStd_Array1OfReal oldVKnots(1,nbVKnots);
      TColStd_Array1OfInteger oldUMults(1,nbUKnots);
      TColStd_Array1OfInteger oldVMults(1,nbVKnots);
      
      BSpl->Poles(oldPoles);
      BSpl->Weights(oldWeights);
      BSpl->UKnots(oldUKnots);
      BSpl->VKnots(oldVKnots);
      BSpl->UMultiplicities(oldUMults);
      BSpl->VMultiplicities(oldVMults);
      
      TColStd_Array1OfReal newUKnots (1,nbUKnots+2);
      TColStd_Array1OfInteger newUMults(1,nbUKnots+2);
      Standard_Real a = 0.5 * ( BSpl->UKnot(2) - BSpl->UKnot(1) + 
			        BSpl->UKnot(nbUKnots) - BSpl->UKnot(nbUKnots-1) );
      
      newUKnots(1) = oldUKnots(1) - a;
      newUKnots(nbUKnots+2) = oldUKnots(nbUKnots) + a;
      newUMults(1) = newUMults(nbUKnots+2) = 1;
      for (Standard_Integer i = 2; i<=nbUKnots+1; i++) {
	newUKnots(i) = oldUKnots(i-1);
	newUMults(i) = oldUMults(i-1);
      }
      newUMults(2) = newUMults(nbUKnots+1) = BSpl->UDegree();
      Handle(Geom_BSplineSurface) res = new Geom_BSplineSurface(oldPoles,
								oldWeights,
								newUKnots,oldVKnots,
								newUMults,oldVMults,
								BSpl->UDegree(),BSpl->VDegree(),
								BSpl->IsUPeriodic(),BSpl->IsVPeriodic());
      BSpl = res;
    }
    else if ( BSpl->UMultiplicity(1) > BSpl->UDegree() ||
	      BSpl->UMultiplicity(BSpl->NbUKnots()) > BSpl->UDegree() + 1 ) set = Standard_False;
    if ( set ) {
      BSpl->SetUPeriodic(); // make periodic
      converted = Standard_True;
    }
  }
  
  if ( vclosed && ! BSpl->IsVPeriodic() && BSpl->NbVPoles() >3 ) {	
    Standard_Boolean set = Standard_True;
    // if degree+1 at ends, first change it to 1 by rearranging knots
    if ( BSpl->VMultiplicity(1) == BSpl->VDegree() + 1 &&
	 BSpl->VMultiplicity(BSpl->NbVKnots()) == BSpl->VDegree() + 1 ) {
      Standard_Integer nbUPoles = BSpl->NbUPoles();
      Standard_Integer nbVPoles = BSpl->NbVPoles();
      TColgp_Array2OfPnt oldPoles(1,nbUPoles,1,nbVPoles);
      TColStd_Array2OfReal oldWeights(1,nbUPoles,1,nbVPoles);
      Standard_Integer nbUKnots = BSpl->NbUKnots();
      Standard_Integer nbVKnots = BSpl->NbVKnots();
      TColStd_Array1OfReal oldUKnots(1,nbUKnots);
      TColStd_Array1OfReal oldVKnots(1,nbVKnots);
      TColStd_Array1OfInteger oldUMults(1,nbUKnots);
      TColStd_Array1OfInteger oldVMults(1,nbVKnots);
      
      BSpl->Poles(oldPoles);
      BSpl->Weights(oldWeights);
      BSpl->UKnots(oldUKnots);
      BSpl->VKnots(oldVKnots);
      BSpl->UMultiplicities(oldUMults);
      BSpl->VMultiplicities(oldVMults);
      
      TColStd_Array1OfReal newVKnots (1,nbVKnots+2);
      TColStd_Array1OfInteger newVMults(1,nbVKnots+2);
      Standard_Real a = 0.5 * ( BSpl->VKnot(2) - BSpl->VKnot(1) + 
			        BSpl->VKnot(nbVKnots) - BSpl->VKnot(nbVKnots-1) );
      
      newVKnots(1) = oldVKnots(1) - a;
      newVKnots(nbVKnots+2) = oldVKnots(nbVKnots) + a;
      newVMults(1) = newVMults(nbVKnots+2) = 1;
      for (Standard_Integer i = 2; i<=nbVKnots+1; i++) {
	newVKnots(i) = oldVKnots(i-1);
	newVMults(i) = oldVMults(i-1);
      }
      newVMults(2) = newVMults(nbVKnots+1) = BSpl->VDegree();
      Handle(Geom_BSplineSurface) res = new Geom_BSplineSurface(oldPoles,
								oldWeights,
								oldUKnots,newVKnots,
								oldUMults,newVMults,
								BSpl->UDegree(),BSpl->VDegree(),
								BSpl->IsUPeriodic(),BSpl->IsVPeriodic());
      BSpl = res;
    }
    else if ( BSpl->VMultiplicity(1) > BSpl->VDegree() ||
	      BSpl->VMultiplicity(BSpl->NbVKnots()) > BSpl->VDegree() + 1 ) set = Standard_False;
    if ( set ) {
      BSpl->SetVPeriodic(); // make periodic
      converted = Standard_True;
    }
  }
  
#ifdef OCCT_DEBUG
  std::cout << "Warning: ShapeCustom_Surface: Closed BSplineSurface is caused to be periodic" << std::endl;
#endif
  if ( ! converted ) return newSurf;
  newSurf = BSpl;
  if ( substitute ) mySurf = newSurf;
  return newSurf;
}
Commit	Line	Data
973c2be1	1	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	2	//
973c2be1	3	// This file is part of Open CASCADE Technology software library.
b311480e	4	//
d5f74e42	5	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	6	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	7	// by the Free Software Foundation, with special exception defined in the file
	8	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	9	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	10	//
973c2be1	11	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	12	// commercial license or contractual agreement.
b311480e	13
7fd59977	14	//abv 06.01.99 fix of misprint
7fd59977	15	//:p6 abv 26.02.99: make ConvertToPeriodic() return Null if nothing done
7fd59977	16
7fd59977	17	#include <ElSLib.hxx>
42cf5bc1	18	#include <Geom_BezierSurface.hxx>
	19	#include <Geom_BSplineSurface.hxx>
	20	#include <Geom_ConicalSurface.hxx>
7fd59977	21	#include <Geom_Curve.hxx>
42cf5bc1	22	#include <Geom_CylindricalSurface.hxx>
7fd59977	23	#include <Geom_Plane.hxx>
7fd59977	24	#include <Geom_SphericalSurface.hxx>
42cf5bc1	25	#include <Geom_Surface.hxx>
7fd59977	26	#include <Geom_ToroidalSurface.hxx>
42cf5bc1	27	#include <GeomAbs_SurfaceType.hxx>
7fd59977	28	#include <GeomAdaptor_HSurface.hxx>
7fd59977	29	#include <GeomAdaptor_Surface.hxx>
42cf5bc1	30	#include <gp_Ax3.hxx>
	31	#include <gp_Cylinder.hxx>
	32	#include <gp_Pln.hxx>
	33	#include <gp_Pnt.hxx>
	34	#include <gp_Vec.hxx>
7fd59977	35	#include <ShapeAnalysis_Geom.hxx>
7fd59977	36	#include <ShapeAnalysis_Surface.hxx>
42cf5bc1	37	#include <ShapeCustom_Surface.hxx>
	38	#include <TColgp_Array1OfPnt.hxx>
	39	#include <TColgp_Array2OfPnt.hxx>
	40	#include <TColStd_Array1OfInteger.hxx>
	41	#include <TColStd_Array1OfReal.hxx>
	42	#include <TColStd_Array2OfReal.hxx>
7fd59977	43
	44	//=======================================================================
	45	//function : ShapeCustom_Surface
	46	//purpose :
	47	//=======================================================================
7fd59977	48	ShapeCustom_Surface::ShapeCustom_Surface() : myGap (0)
	49	{
	50	}
	51
	52	//=======================================================================
	53	//function : ShapeCustom_Surface
	54	//purpose :
	55	//=======================================================================
	56
	57	ShapeCustom_Surface::ShapeCustom_Surface (const Handle(Geom_Surface)& S)
	58	: myGap (0)
	59	{
	60	Init ( S );
	61	}
	62
	63	//=======================================================================
	64	//function : Init
	65	//purpose :
	66	//=======================================================================
	67
	68	void ShapeCustom_Surface::Init (const Handle(Geom_Surface)& S)
	69	{
	70	mySurf = S;
	71	}
	72
	73	//=======================================================================
	74	//function : ConvertToAnalytical
	75	//purpose :
	76	//=======================================================================
	77
	78	Handle(Geom_Surface) ShapeCustom_Surface::ConvertToAnalytical (const Standard_Real tol,
	79	const Standard_Boolean substitute)
	80	{
	81	Handle(Geom_Surface) newSurf;
	82
	83	Standard_Integer nUP, nVP, nCP, i, j , UDeg, VDeg;
	84	Standard_Real U1, U2, V1, V2, C1, C2, DU, DV, U=0, V=0;
	85	Handle(Geom_Curve) iso;
	86	Standard_Boolean uClosed = Standard_True;
	87
	88	// seuls cas traites : BSpline et Bezier
	89	Handle(Geom_BSplineSurface) theBSplneS =
	90	Handle(Geom_BSplineSurface)::DownCast(mySurf);
	91	if (theBSplneS.IsNull()) {
	92	Handle(Geom_BezierSurface) theBezierS =
	93	Handle(Geom_BezierSurface)::DownCast(mySurf);
	94	if (!theBezierS.IsNull()) { // Bezier :
	95	nUP = theBezierS->NbUPoles();
	96	nVP = theBezierS->NbVPoles();
	97	UDeg = theBezierS->UDegree();
	98	VDeg = theBezierS->VDegree();
	99	}
	100	else return newSurf; // non reconnu : terminus
	101	}
	102	else { // BSpline :
	103	nUP = theBSplneS->NbUPoles();
	104	nVP = theBSplneS->NbVPoles();
	105	UDeg = theBSplneS->UDegree();
	106	VDeg = theBSplneS->VDegree();
	107	}
	108
	109
	110	mySurf->Bounds(U1, U2, V1, V2);
	111	// mySurf->Bounds(U1, U2, V1, V2);
112	TColgp_Array1OfPnt p1(1, 3), p2(1, 3), p3(1, 3);
113	TColStd_Array1OfReal R(1,3);
114	gp_Pnt origPnt, resPnt;
115	gp_Vec origD1U, resD1U, resD1V;
116
117	Standard_Boolean aCySpCo = Standard_False;
118	Standard_Boolean aToroid = Standard_False;
119	Standard_Boolean aPlanar = Standard_False;
120
121	if (nUP == 2 && nVP == 2) {
122	if (UDeg == 1 && VDeg == 1) aPlanar = Standard_True;
123	} else if (mySurf->IsUClosed()) { // VRAI IsUClosed
124	if (mySurf->IsVClosed()) aToroid = Standard_True;
125	else aCySpCo = Standard_True;
126	} else {
127	if(mySurf->IsVClosed()) { // VRAI IsVClosed
128	aCySpCo = Standard_True;
129	uClosed = Standard_False;
130	}
131	}
132
133	if (aPlanar) {
134	// NearestPlane ...
135	TColgp_Array1OfPnt Pnts(1,4);
136	Pnts.SetValue(1,mySurf->Value(U1,V1));
137	Pnts.SetValue(2,mySurf->Value(U2,V1));
138	Pnts.SetValue(3,mySurf->Value(U1,V2));
139	Pnts.SetValue(4,mySurf->Value(U2,V2));
140	gp_Pln aPln;// Standard_Real Dmax;
141	Standard_Integer It = ShapeAnalysis_Geom::NearestPlane (Pnts,aPln,myGap/Dmax/);
142
143	// ICI, on fabrique le plan, et zou
144	if (It == 0 \|\| myGap/Dmax/ > tol) return newSurf; // pas un plan
145
146	// IL RESTE a verifier l orientation ...
147	// On regarde sur chaque surface les vecteurs P(U0->U1),P(V0->V1)
148	// On prend la normale : les deux normales doivent etre dans le meme sens
149	// Sinon, inverser la normale (pas le Pln entier !) et refaire la Plane
150	newSurf = new Geom_Plane (aPln);
151	gp_Vec uold (Pnts(1),Pnts(2));
152	gp_Vec vold (Pnts(1),Pnts(3));
153	gp_Vec nold = uold.Crossed (vold);
154	gp_Vec unew (newSurf->Value(U1,V1), newSurf->Value(U2,V1));
155	gp_Vec vnew (newSurf->Value(U1,V1), newSurf->Value(U1,V2));
156	gp_Vec nnew = unew.Crossed (vnew);
157	if (nold.Dot (nnew) < 0.0) {
158	gp_Ax3 ax3 = aPln.Position();
159	ax3.ZReverse();
160	ax3.XReverse();
161	aPln = gp_Pln (ax3);
162	newSurf = new Geom_Plane (aPln);
163	}
164
165	if (substitute) {
166	Init (newSurf);
167	}
168	return newSurf;
169
170
171	} else if (aCySpCo) {
172	if (!uClosed) {
173	C1 = U1; C2 = U2;
174	U1 = V1; U2 = V2;
175	V1 = C1; V2 = C2;
176	nCP = nUP; nUP = nVP; nVP = nCP;
177	}
178
179	for (i=1; i<=3; i++) {
180	if (i==1) V = V1;
181	else if (i==2) V = V2;
182	else if (i==3) V = 0.5*(V1+V2);
183
184	if(uClosed) iso = mySurf->VIso(V);
185	else iso = mySurf->UIso(V);
186
187	iso->D0(U1, p1(i));
188	iso->D0(0.5*(U1+U2), p2(i));
189	p3(i).SetCoord(0.5*(p1(i).X()+p2(i).X()),
190	0.5*(p1(i).Y()+p2(i).Y()),
191	0.5*(p1(i).Z()+p2(i).Z()));
192	R(i) = p3(i).Distance(p1(i));
04232180	193	// std::cout<<"sphere, i="<<i<<" V="<<V<<" R="<<R(i)<<" p1="<<p1(i).X()<<","<<p1(i).Y()<<","<<p1(i).Z()<<" p2="<<p2(i).X()<<","<<p2(i).Y()<<","<<p2(i).Z()<<" p3="<<p3(i).X()<<","<<p3(i).Y()<<","<<p3(i).Z()<<std::endl;
7fd59977	194	}
	195
	196	iso->D1 (0.,origPnt,origD1U);
	197	gp_Vec xVec(p3(3), p1(3));
	198	gp_Vec aVec(p3(1), p3(2));
	199	// gp_Dir xDir(xVec); ne sert pas. Null si R3 = 0
	200	gp_Dir aDir(aVec);
	201	gp_Ax3 aAx3 (p3(1),aDir,xVec);
	202	// CKY 3-FEV-1997 : verification du sens de description
	203	//gp_Dir AXY = aAx3.YDirection(); // AXY not used (skl)
	204	if (aAx3.YDirection().Dot (origD1U) < 0) {
0797d9d3	205	#ifdef OCCT_DEBUG
04232180	206	std::cout<<" Surface Analytique : sens a inverser"<<std::endl;
7fd59977	207	#endif
	208	aAx3.YReverse(); // mais X reste !
	209	}
	210
	211	if (nVP > 2) {
	212	if ((Abs(R(1)) < tol) &&
	213	(Abs(R(2)) < tol) &&
	214	(Abs(R(3)) > tol)) {
	215	// deja fait gp_Ax3 aAx3(p3(1), aDir, xVec);
	216	//gp_Ax3 aAx3(p3(3), aDir);
	217	Handle(Geom_SphericalSurface) anObject =
	218	new Geom_SphericalSurface(aAx3, R(3));
	219	if (!uClosed) anObject->UReverse();
	220	newSurf = anObject;
	221	}
	222	}
	223	else if (nVP == 2) {
	224
	225	// deja fait gp_Ax3 aAx3(p3(1), aDir, xVec);
	226	//gp_Ax3 aAx3(p3(1), aDir);
	227
	228	if (Abs(R(2)-R(1)) < tol) {
	229	Handle(Geom_CylindricalSurface) anObject =
	230	new Geom_CylindricalSurface(aAx3, R(1));
	231	if (!uClosed) anObject->UReverse();
	232	newSurf = anObject;
	233	}
	234	else {
	235	gp_Vec aVec2(p1(1), p1(2));
	236	Standard_Real angle = aVec.Angle(aVec2);
	237	if (R(1) < R(2)) {
	238	Handle(Geom_ConicalSurface) anObject =
	239	new Geom_ConicalSurface(aAx3, angle, R(1));
	240	//if (!uClosed) anObject->UReverse();
	241	anObject->UReverse();
	242	newSurf = anObject;
	243	}
	244	else {
	245	aDir.Reverse();
	246	gp_Vec anotherXVec(p3(2), p1(2));
	247	gp_Dir anotherXDir(anotherXVec);
	248	gp_Ax3 anotherAx3(p3(2), aDir, anotherXDir);
	249	Handle(Geom_ConicalSurface) anObject =
	250	new Geom_ConicalSurface(anotherAx3, angle, R(2));
	251	//if (!uClosed) anObject->UReverse();
	252	anObject->UReverse();
	253	newSurf = anObject;
	254	}
	255	}
	256	}
	257	}
	258	else if (aToroid) {
	259	// test by iso U and isoV
	260	Standard_Boolean isFound = Standard_False;
	261	for (j=1; (j<=2) && !isFound; j++) {
	262	if (j==2) {
	263	C1 = U1; C2 = U2;
	264	U1 = V1; U2 = V2;
	265	V1 = C1; V2 = C2;
	266	}
	267	for (i=1; i<=3; i++) {
	268	if (i==1) U = U1;
	269	else if (i==2) U = 0.5*(U1+U2);
	270	else if (i==3) U = 0.25*(U1+U2);
271
272	iso = mySurf->UIso(U);
273
274	iso->D0(V1, p1(i));
275	iso->D0(0.5*(V1+V2), p2(i));
276	p3(i).SetCoord(0.5*(p1(i).X()+p2(i).X()),
277	0.5*(p1(i).Y()+p2(i).Y()),
278	0.5*(p1(i).Z()+p2(i).Z()));
279	R(i) = p3(i).Distance(p1(i));
280	}
281	if ((Abs(R(1)-R(2))< tol) &&
282	(Abs(R(1)-R(3))< tol)) {
283	gp_Pnt p10(0.5*(p3(1).X()+p3(2).X()),
284	0.5*(p3(1).Y()+p3(2).Y()),
285	0.5*(p3(1).Z()+p3(2).Z()));
286	gp_Vec aVec(p10, p3(1));
287	gp_Vec aVec2(p10, p3(3));
288	Standard_Real RR1 = R(1), RR2 = R(2), RR3;
289	aVec ^= aVec2;
290
291	if (aVec.Magnitude() <= gp::Resolution()) aVec.SetCoord(0., 0., 1.);
292
293	gp_Dir aDir(aVec);
294
295	gp_Ax3 aAx3(p10, aDir);
296	RR1 = p10.Distance(p3(1));
297	// modif empirique (pourtant NON DEMONTREE) : inverser roles RR1,RR2
298	// CKY, 24-JAN-1997
299	if (RR1 < RR2) { RR3 = RR1; RR1 = RR2; RR2 = RR3; }
300	Handle(Geom_ToroidalSurface) anObject =
301	new Geom_ToroidalSurface(aAx3, RR1, RR2);
302	if (j==2) anObject->UReverse();
303	anObject->D1 (0.,0.,resPnt,resD1U,resD1V);
0797d9d3	304	#ifdef OCCT_DEBUG
7fd59977	305	if (resD1U.Dot(origD1U) < 0 && j != 2)
04232180	306	std::cout<<" Tore a inverser !"<<std::endl;
7fd59977	307	#endif
	308	newSurf = anObject;
	309	isFound = Standard_True;
	310	}
	311	}
	312	}
	313	if (newSurf.IsNull()) return newSurf;
	314
	315	//---------------------------------------------------------------------
	316	// verification
	317	//---------------------------------------------------------------------
	318
	319	Handle(GeomAdaptor_HSurface) NHS = new GeomAdaptor_HSurface (newSurf);
	320	GeomAdaptor_Surface& SurfAdapt = NHS->ChangeSurface();
	321
	322	const Standard_Integer NP = 21;
1d47d8d0	323	Standard_Real S = 0., T = 0.; // U,V deja fait
7fd59977	324	gp_Pnt P3d, P3d2;
	325	Standard_Boolean onSurface = Standard_True;
	326
	327	Standard_Real dis; myGap = 0.;
	328
	329	DU = (U2-U1)/(NP-1);
	330	DV = (V2-V1)/(NP-1);
	331	for (j=1; (j<=NP) && onSurface; j++) {
	332	V = V1 + DV*(j-1);
	333
	334	if(uClosed) iso = mySurf->VIso(V);
	335	else iso = mySurf->UIso(V);
	336
	337	for (i=1; i<=NP; i++) {
	338	U = U1 + DU*(i-1);
	339	iso->D0(U, P3d);
	340	switch (SurfAdapt.GetType()){
	341
	342	case GeomAbs_Cylinder :
	343	{
	344	gp_Cylinder Cylinder = SurfAdapt.Cylinder();
	345	ElSLib::Parameters( Cylinder, P3d, S, T);
	346	break;
	347	}
	348	case GeomAbs_Cone :
	349	{
	350	gp_Cone Cone = SurfAdapt.Cone();
	351	ElSLib::Parameters( Cone, P3d, S, T);
	352	break;
	353	}
	354	case GeomAbs_Sphere :
	355	{
	356	gp_Sphere Sphere = SurfAdapt.Sphere();
	357	ElSLib::Parameters( Sphere, P3d, S, T);
	358	break;
	359	}
	360	case GeomAbs_Torus :
	361	{
	362	gp_Torus Torus = SurfAdapt.Torus();
	363	ElSLib::Parameters( Torus, P3d, S, T);
	364	break;
	365	}
	366	default:
	367	break;
	368	}
	369
	370	newSurf->D0(S, T, P3d2);
	371
	372	dis = P3d.Distance(P3d2);
	373	if (dis > myGap) myGap = dis;
	374
	375	if (dis > tol) {
	376	onSurface = Standard_False;
	377	newSurf.Nullify();
	378	// The presumption is rejected
	379	break;
	380	}
	381	}
	382	}
	383	if (substitute && !NHS.IsNull()) {
	384	Init (newSurf);
	385	}
	386	return newSurf;
	387	}
388
389	//%pdn 30 Nov 98: converting bspline surfaces with degree+1 at ends to periodic
390	// UKI60591, entity 48720
391	Handle(Geom_Surface) ShapeCustom_Surface::ConvertToPeriodic (const Standard_Boolean substitute,
392	const Standard_Real preci)
393	{
394	Handle(Geom_Surface) newSurf;
395	Handle(Geom_BSplineSurface) BSpl = Handle(Geom_BSplineSurface)::DownCast(mySurf);
396	if (BSpl.IsNull()) return newSurf;
397
398	ShapeAnalysis_Surface sas(mySurf);
399	Standard_Boolean uclosed = sas.IsUClosed(preci);
400	Standard_Boolean vclosed = sas.IsVClosed(preci);
401
402	if ( ! uclosed && ! vclosed ) return newSurf;
403
404	Standard_Boolean converted = Standard_False; //:p6
405
406	if ( uclosed && ! BSpl->IsUPeriodic() && BSpl->NbUPoles() >3 ) {
407	Standard_Boolean set = Standard_True;
408	// if degree+1 at ends, first change it to 1 by rearranging knots
409	if ( BSpl->UMultiplicity(1) == BSpl->UDegree() + 1 &&
410	BSpl->UMultiplicity(BSpl->NbUKnots()) == BSpl->UDegree() + 1 ) {
411	Standard_Integer nbUPoles = BSpl->NbUPoles();
412	Standard_Integer nbVPoles = BSpl->NbVPoles();
413	TColgp_Array2OfPnt oldPoles(1,nbUPoles,1,nbVPoles);
414	TColStd_Array2OfReal oldWeights(1,nbUPoles,1,nbVPoles);
415	Standard_Integer nbUKnots = BSpl->NbUKnots();
416	Standard_Integer nbVKnots = BSpl->NbVKnots();
417	TColStd_Array1OfReal oldUKnots(1,nbUKnots);
418	TColStd_Array1OfReal oldVKnots(1,nbVKnots);
419	TColStd_Array1OfInteger oldUMults(1,nbUKnots);
420	TColStd_Array1OfInteger oldVMults(1,nbVKnots);
421
422	BSpl->Poles(oldPoles);
423	BSpl->Weights(oldWeights);
424	BSpl->UKnots(oldUKnots);
425	BSpl->VKnots(oldVKnots);
426	BSpl->UMultiplicities(oldUMults);
427	BSpl->VMultiplicities(oldVMults);
428
429	TColStd_Array1OfReal newUKnots (1,nbUKnots+2);
430	TColStd_Array1OfInteger newUMults(1,nbUKnots+2);
431	Standard_Real a = 0.5 * ( BSpl->UKnot(2) - BSpl->UKnot(1) +
432	BSpl->UKnot(nbUKnots) - BSpl->UKnot(nbUKnots-1) );
433
434	newUKnots(1) = oldUKnots(1) - a;
435	newUKnots(nbUKnots+2) = oldUKnots(nbUKnots) + a;
436	newUMults(1) = newUMults(nbUKnots+2) = 1;
437	for (Standard_Integer i = 2; i<=nbUKnots+1; i++) {
438	newUKnots(i) = oldUKnots(i-1);
439	newUMults(i) = oldUMults(i-1);
440	}
441	newUMults(2) = newUMults(nbUKnots+1) = BSpl->UDegree();
442	Handle(Geom_BSplineSurface) res = new Geom_BSplineSurface(oldPoles,
443	oldWeights,
444	newUKnots,oldVKnots,
445	newUMults,oldVMults,
446	BSpl->UDegree(),BSpl->VDegree(),
447	BSpl->IsUPeriodic(),BSpl->IsVPeriodic());
448	BSpl = res;
449	}
450	else if ( BSpl->UMultiplicity(1) > BSpl->UDegree() \|\|
451	BSpl->UMultiplicity(BSpl->NbUKnots()) > BSpl->UDegree() + 1 ) set = Standard_False;
452	if ( set ) {
453	BSpl->SetUPeriodic(); // make periodic
454	converted = Standard_True;
455	}
456	}
457
458	if ( vclosed && ! BSpl->IsVPeriodic() && BSpl->NbVPoles() >3 ) {
459	Standard_Boolean set = Standard_True;
460	// if degree+1 at ends, first change it to 1 by rearranging knots
461	if ( BSpl->VMultiplicity(1) == BSpl->VDegree() + 1 &&
462	BSpl->VMultiplicity(BSpl->NbVKnots()) == BSpl->VDegree() + 1 ) {
463	Standard_Integer nbUPoles = BSpl->NbUPoles();
464	Standard_Integer nbVPoles = BSpl->NbVPoles();
465	TColgp_Array2OfPnt oldPoles(1,nbUPoles,1,nbVPoles);
466	TColStd_Array2OfReal oldWeights(1,nbUPoles,1,nbVPoles);
467	Standard_Integer nbUKnots = BSpl->NbUKnots();
468	Standard_Integer nbVKnots = BSpl->NbVKnots();
469	TColStd_Array1OfReal oldUKnots(1,nbUKnots);
470	TColStd_Array1OfReal oldVKnots(1,nbVKnots);
471	TColStd_Array1OfInteger oldUMults(1,nbUKnots);
472	TColStd_Array1OfInteger oldVMults(1,nbVKnots);
473
474	BSpl->Poles(oldPoles);
475	BSpl->Weights(oldWeights);
476	BSpl->UKnots(oldUKnots);
477	BSpl->VKnots(oldVKnots);
478	BSpl->UMultiplicities(oldUMults);
479	BSpl->VMultiplicities(oldVMults);
480
481	TColStd_Array1OfReal newVKnots (1,nbVKnots+2);
482	TColStd_Array1OfInteger newVMults(1,nbVKnots+2);
483	Standard_Real a = 0.5 * ( BSpl->VKnot(2) - BSpl->VKnot(1) +
484	BSpl->VKnot(nbVKnots) - BSpl->VKnot(nbVKnots-1) );
485
486	newVKnots(1) = oldVKnots(1) - a;
487	newVKnots(nbVKnots+2) = oldVKnots(nbVKnots) + a;
488	newVMults(1) = newVMults(nbVKnots+2) = 1;
489	for (Standard_Integer i = 2; i<=nbVKnots+1; i++) {
490	newVKnots(i) = oldVKnots(i-1);
491	newVMults(i) = oldVMults(i-1);
492	}
493	newVMults(2) = newVMults(nbVKnots+1) = BSpl->VDegree();
494	Handle(Geom_BSplineSurface) res = new Geom_BSplineSurface(oldPoles,
495	oldWeights,
496	oldUKnots,newVKnots,
497	oldUMults,newVMults,
498	BSpl->UDegree(),BSpl->VDegree(),
499	BSpl->IsUPeriodic(),BSpl->IsVPeriodic());
500	BSpl = res;
501	}
502	else if ( BSpl->VMultiplicity(1) > BSpl->VDegree() \|\|
503	BSpl->VMultiplicity(BSpl->NbVKnots()) > BSpl->VDegree() + 1 ) set = Standard_False;
504	if ( set ) {
505	BSpl->SetVPeriodic(); // make periodic
506	converted = Standard_True;
507	}
508	}
509
0797d9d3	510	#ifdef OCCT_DEBUG
04232180	511	std::cout << "Warning: ShapeCustom_Surface: Closed BSplineSurface is caused to be periodic" << std::endl;
7fd59977	512	#endif
	513	if ( ! converted ) return newSurf;
	514	newSurf = BSpl;
	515	if ( substitute ) mySurf = newSurf;
	516	return newSurf;
	517	}