0025880: fuzzy booleans with multiple tools
[occt.git] / src / GeomLib / GeomLib.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 7c51d69..a38fda2
@@ -1,10 +1,21 @@
-// File:       GeomLib.cxx
-// Created:    Wed Jul  7 15:32:09 1993
-// Author:     Jean Claude VAUTHIER
+// Created on: 1993-07-07
+// Created by: Jean Claude VAUTHIER
+// Copyright (c) 1993-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 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.
 
-// Copyright:    Matra Datavision      
 // Version:    
-// History:    pmn 24/09/96 Ajout du prolongement de courbe.
+//pmn 24/09/96 Ajout du prolongement de courbe.
 //              jct 15/04/97 Ajout du prolongement de surface.
 //              jct 24/04/97 simplification ou suppression de calculs
 //                           inutiles dans ExtendSurfByLength
 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
 #include <GeomConvert_ApproxSurface.hxx>
 
+#include <CSLib.hxx>
+#include <CSLib_NormalStatus.hxx>
+
 
 #include <Standard_ConstructionError.hxx>
 
@@ -938,18 +952,9 @@ void GeomLib::SameRange(const Standard_Real         Tolerance,
   else { // On segmente le resultat
     Handle(Geom2d_TrimmedCurve) TC =
       new Geom2d_TrimmedCurve( CurvePtr, FirstOnCurve, LastOnCurve );
-
-    Standard_Real newFirstOnCurve = TC->FirstParameter(), newLastOnCurve = TC->LastParameter();
-    
+    //
     Handle(Geom2d_BSplineCurve) BS =
       Geom2dConvert::CurveToBSplineCurve(TC);
-
-    if (BS->IsPeriodic()) 
-      BS->Segment( newFirstOnCurve, newLastOnCurve) ;
-    else 
-      BS->Segment( Max(newFirstOnCurve, BS->FirstParameter()),
-                  Min(newLastOnCurve,  BS->LastParameter()) );
-
     TColStd_Array1OfReal Knots(1,BS->NbKnots());
     BS->Knots(Knots);
     
@@ -1039,7 +1044,7 @@ void GeomLib::BuildCurve3d(const Standard_Real           Tolerance,
                           Adaptor3d_CurveOnSurface&       Curve, 
                           const Standard_Real           FirstParameter,
                           const Standard_Real           LastParameter,
-                          Handle_Geom_Curve&            NewCurvePtr, 
+                          Handle(Geom_Curve)&            NewCurvePtr, 
                           Standard_Real&                MaxDeviation,
                           Standard_Real&                AverageDeviation,
                           const GeomAbs_Shape           Continuity,
@@ -1455,8 +1460,10 @@ void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
   }     
 
 
-  Standard_Boolean rational = ( InU && BS->IsURational() ) 
-                                  || ( !InU && BS->IsVRational() ) ;
+// IFV Fix OCC bug 0022694 - wrong result extrapolating rational surfaces
+//   Standard_Boolean rational = ( InU && BS->IsURational() ) 
+//                                   || ( !InU && BS->IsVRational() ) ;
+  Standard_Boolean rational = (BS->IsURational() ||  BS->IsVRational());
   Standard_Boolean NullWeight;
    Standard_Real EpsW = 10*Precision::PConfusion();
   Standard_Integer gap = 3;
@@ -1464,10 +1471,10 @@ void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
 
 
         
-  Standard_Integer Cdeg, Cdim, NbP, Ksize, Psize ;
+  Standard_Integer Cdeg = 0, Cdim = 0, NbP = 0, Ksize = 0, Psize = 1;
   Standard_Integer ii, jj, ipole, Kount;  
   Standard_Real Tbord, lambmin=Length;
-  Standard_Real * Padr;
+  Standard_Real * Padr = NULL;
   Standard_Boolean Ok;
   Handle(TColStd_HArray1OfReal)  FKnots, Point, lambda, Tgte, Poles;
 
@@ -1650,7 +1657,7 @@ void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
   }
   
 //  tableaux necessaires pour l'extension
-  Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots;
+  Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots = 0;
   TColStd_Array1OfReal  FK(1, Ksize2) ; 
   Standard_Real * FKRadr = &FK(1);
 
@@ -1732,7 +1739,7 @@ void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
     }
 
     if (NullWeight) {
-#if DEB
+#ifdef OCCT_DEBUG
       cout << "Echec de l'Extension rationnelle" << endl;    
 #endif
       lambmin /= 3.;
@@ -1832,7 +1839,7 @@ void GeomLib::Inertia(const TColgp_Array1OfPnt& Points,
 
   math_Jacobi J(M);
   if (!J.IsDone()) {
-#if DEB
+#ifdef OCCT_DEBUG
     cout << "Erreur dans Jacobbi" << endl;
     M.Dump(cout);
 #endif
@@ -2323,12 +2330,6 @@ Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S,
 
   Standard_Real MDU = DU.SquareMagnitude(), MDV = DV.SquareMagnitude();
 
-  Standard_Real h, sign;
-  Standard_Boolean AlongV;
-  Handle(Geom_Curve) Iso;
-  Standard_Real t, first, last, bid1, bid2;
-  gp_Vec Tang;
-
   if(MDU >= aTol2 && MDV >= aTol2) {
     gp_Vec Norm = DU^DV;
     Standard_Real Magn = Norm.SquareMagnitude();
@@ -2339,125 +2340,94 @@ Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S,
 
     return 0;
   }
-  else if(MDU < aTol2 && MDV >= aTol2) {
-    AlongV = Standard_True;
-    Iso = S->UIso(UV.X());
-    t = UV.Y();
-    S->Bounds(bid1, bid2, first, last);
-  }
-  else if(MDU >= aTol2 && MDV < aTol2) {
-    AlongV = Standard_False;
-    Iso = S->VIso(UV.Y());
-    t = UV.X();
-    S->Bounds(first, last, bid1, bid2);
-  }
   else {
-    return 3;
-  }
-
-  Standard_Real L = .001;
-
-  if(Precision::IsInfinite(Abs(first))) first = t - 1.;
-  if(Precision::IsInfinite(Abs(last)))  last  = t + 1.;
-  
-  if(last - t >= t - first) {
-    sign = 1.;
-  }
-  else {
-    sign = -1.;
-  }
-
-  Standard_Real hmax = .01*(last - first);
-  if(AlongV) {
-    h = Min(L/sqrt(MDV), hmax);
-    S->D1(UV.X(), UV.Y() + sign*h, DummyPnt, DU, DV);
-  }
-  else {
-    h = Min(L/sqrt(MDU), hmax);
-    S->D1(UV.X() + sign*h, UV.Y(), DummyPnt, DU, DV);
-  }
-
-  gp_Vec DD;
-
-  gp_Vec NAux = DU^DV;
-  Standard_Real h1 = h;
-  while(NAux.SquareMagnitude() < aTol2) {
-    h1 += h;
-    if(AlongV) {
-      Standard_Real v = UV.Y() + sign*h1;
-
-      if(v < first || v > last) return 3;
-
-      S->D1(UV.X(), v, DummyPnt, DU, DV);
-    }
-    else {
-      Standard_Real v = UV.X() + sign*h1;
-
-      if(v < first || v > last) return 3;
-
-      S->D1(v, UV.Y(), DummyPnt, DU, DV);
-
-    }
-    NAux = DU^DV;
-  }
-
-
-  Iso->D2(t, DummyPnt, Tang, DD);
-
-  if(DD.SquareMagnitude() >= aTol2) {
-    gp_Vec NV = DD * (Tang * Tang) - Tang * (Tang * DD);
-    Standard_Real MagnNV = NV.SquareMagnitude();
-    if(MagnNV < aTol2) return 3;
-
-    MagnNV = sqrt(MagnNV);
-    N.SetXYZ(NV.XYZ()/MagnNV);
-
-    Standard_Real par = .5*(bid2+bid1);
-
-    if(AlongV) {
-      Iso = S->UIso(par);
-    }
-    else {
-      Iso = S->VIso(par);
-    }
-
-    Iso->D2(t, DummyPnt, Tang, DD);
-
-    gp_Vec N1V = DD * (Tang * Tang) - Tang * (Tang * DD);
-    Standard_Real MagnN1V = N1V.SquareMagnitude();
-    if(MagnN1V < aTol2) return 3;
-
-    MagnN1V = sqrt(MagnN1V);
-    gp_Dir N1(N1V.XYZ()/MagnN1V);
-
-    Standard_Integer res = 1;
+    gp_Vec D2U, D2V, D2UV;
+    Standard_Boolean isDone;
+    CSLib_NormalStatus aStatus;
+    gp_Dir aNormal;
+
+    S->D2(UV.X(), UV.Y(), DummyPnt, DU, DV, D2U, D2V, D2UV);
+    CSLib::Normal(DU, DV, D2U, D2V, D2UV, Tol, isDone, aStatus, aNormal);
+
+    if (isDone) {
+     Standard_Real Umin, Umax, Vmin, Vmax;
+     Standard_Real step = 1.0e-5;
+     Standard_Real eps = 1.0e-16;
+     Standard_Real sign = -1.0;
+     
+     S->Bounds(Umin, Umax, Vmin, Vmax);
+
+     // check for cone apex singularity point
+     if ((UV.Y() > Vmin + step) && (UV.Y() < Vmax - step))
+     {
+       gp_Dir aNormal1, aNormal2;
+       Standard_Real aConeSingularityAngleEps = 1.0e-4;
+       S->D1(UV.X(), UV.Y() - sign * step, DummyPnt, DU, DV);
+       if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
+         aNormal1 = DU^DV;
+         S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
+         if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
+           aNormal2 = DU^DV;
+           if (aNormal1.IsOpposite(aNormal2, aConeSingularityAngleEps))
+             return 2;
+         }
+       }
+     }
 
-    if(N*N1 < 1. - Tol) res = 2;
+     // Along V
+     if(MDU < aTol2 && MDV >= aTol2) {
+       if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
+         sign = 1.0;
+       S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
+       gp_Vec Norm = DU^DV;
+       if (Norm.SquareMagnitude() < eps) {
+         Standard_Real sign1 = -1.0;
+         if ((Umax - UV.X()) > (UV.X() - Umin))
+           sign1 = 1.0;
+         S->D1(UV.X() + sign1 * step, UV.Y() + sign * step, DummyPnt, DU, DV);
+         Norm = DU^DV;
+       }
+       if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
+         aNormal.Reverse();
+     }
 
-    if(N*NAux <= 0.) N.Reverse();
+     // Along U
+     if(MDV < aTol2 && MDU >= aTol2) {
+       if ((Umax - UV.X()) > (UV.X() - Umin))
+         sign = 1.0;
+       S->D1(UV.X() + sign * step, UV.Y(), DummyPnt, DU, DV);
+       gp_Vec Norm = DU^DV;
+       if (Norm.SquareMagnitude() < eps) {
+         Standard_Real sign1 = -1.0;
+         if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
+           sign1 = 1.0;
+         S->D1(UV.X() + sign * step, UV.Y() + sign1 * step, DummyPnt, DU, DV);
+         Norm = DU^DV;
+       }
+       if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
+         aNormal.Reverse();
+     }
 
-    return res;
-  }
-  else {
-    //Seems to be conical singular point
-    
-    if(AlongV) {
-      NAux = DU^Tang;
+      // quasysingular
+      if ((aStatus == CSLib_D1NuIsNull) || (aStatus == CSLib_D1NvIsNull) || 
+          (aStatus == CSLib_D1NuIsParallelD1Nv)) {
+            N.SetXYZ(aNormal.XYZ());
+            return 1;
+      }
+      // conical
+      if (aStatus == CSLib_InfinityOfSolutions)
+          return 2;
     }
+    // computation is impossible
     else {
-      NAux = Tang^DV;
+      // conical
+      if (aStatus == CSLib_D1NIsNull) {
+        return 2;
+      }
+      return 3;
     }
-
-    sign = NAux.Magnitude();
-
-    if(sign < Tol) return 3;
-
-    N = NAux;
-
-    return 2;
-   
   }
-
+  return 3;
 }