0023830: BRepExtrema_DistShapeShape does not find intersection of face with edge
[occt.git] / src / Extrema / Extrema_GenExtCS.cxx
index 33261d1..3b2639c 100755 (executable)
@@ -117,7 +117,6 @@ void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S,
                                  const Standard_Real Tol2)
 {
   myS = (Adaptor3d_SurfacePtr)&S;
-  mypoints2 = new TColgp_HArray2OfPnt(0,NbU+1,0,NbV+1);
   myusample = NbU;
   myvsample = NbV;
   myumin = Umin;
@@ -125,29 +124,6 @@ void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S,
   myvmin = Vmin;
   myvsup = Vsup;
   mytol2 = Tol2;
-
-// Parametrage de l echantillon sur S
-
-  Standard_Real PasU = myusup - myumin;
-  Standard_Real PasV = myvsup - myvmin;
-  Standard_Real U0 = PasU / myusample / 100.;
-  Standard_Real V0 = PasV / myvsample / 100.;
-  gp_Pnt P1;
-  PasU = (PasU - U0) / (myusample - 1);
-  PasV = (PasV - V0) / (myvsample - 1);
-  U0 = myumin + U0/2.;
-  V0 = myvmin + V0/2.;
-
-// Calcul des distances
-
-  Standard_Integer NoU, NoV;
-  Standard_Real U, V;
-  for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
-    for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
-      P1 = myS->Value(U, V);
-      mypoints2->SetValue(NoU,NoV,P1);
-    }
-  }
 }
 
 //=======================================================================
@@ -181,13 +157,7 @@ void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
   mytsup = tsup;
   mytol1 = Tol1;
   mytsample = NbT;
-  mypoints1 = new TColgp_HArray1OfPnt(0,NbT);
-
-  Standard_Real t1, U, V;
-  Standard_Integer NoT, NoU, NoV;
-  gp_Pnt P1, P2;
-
-// Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000
+  // Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000
 
   Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin;
   if (Precision::IsInfinite(trimusup)){
@@ -202,40 +172,7 @@ void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
   if (Precision::IsInfinite(trimvmin)){
     trimvmin = -1.0e+10;
   }
-  
-
-// Parametrage de l echantillon sur C
-
-  Standard_Real PasT = mytsup - mytmin;
-  Standard_Real t0 = PasT / mytsample / 100.;
-  PasT = (PasT - t0) / (mytsample - 1);
-  t0 = mytmin + t0/2.;
-
-// Parametrage de l echantillon sur S
-
-  Standard_Real PasU = trimusup - trimumin;
-  Standard_Real PasV = trimvsup - trimvmin;
-  Standard_Real U0 = PasU / myusample / 100.;
-  Standard_Real V0 = PasV / myvsample / 100.;
-  PasU = (PasU - U0) / (myusample - 1);
-  PasV = (PasV - V0) / (myvsample - 1);
-  U0 = trimumin + U0/2.;
-  V0 = trimvmin + V0/2.;
-
-// Calcul des distances
-
-  for ( NoT = 1, t1 = t0; NoT <= mytsample; NoT++, t1 += PasT) {
-    P1 = C.Value(t1);
-    mypoints1->SetValue(NoT,P1);
-  }
-  
-
-/*
-b- Calcul des minima:
-   -----------------
-   b.a) Initialisations:
-*/
-
+  //
   math_Vector Tol(1, 3);
   Tol(1) = mytol1;
   Tol(2) = mytol2;
@@ -248,7 +185,6 @@ b- Calcul des minima:
   UVsup(1) = mytsup;
   UVsup(2) = trimusup;
   UVsup(3) = trimvsup;
-
   // 18/02/02 akm vvv : (OCC163) bad extremas - extrusion surface versus the line.
   //                    Try to preset the initial solution as extrema between
   //                    extrusion direction and the curve.
@@ -268,80 +204,168 @@ b- Calcul des minima:
       Extrema_POnCurv aP1, aP2;
       for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
       {
-       aLocator.Points (iExt, aP1, aP2);
-       // Parameter on curve
-       UV(1) = aP1.Parameter();
-       // To find parameters on surf, try ExtPS
-       Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2);
-       if (aPreciser.IsDone())
-       {
-         // Managed to find extremas between point and surface
-         Standard_Integer iPExt;
-         for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++)
-         {
-           aPreciser.Point(iPExt).Parameter(UV(2),UV(3));
-           math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
-         }
-       }
-       else
-       {
-         // Failed... try the point on iso line
-         UV(2) = dfUFirst;
-         UV(3) = aP2.Parameter();
-         math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
-       }
+  aLocator.Points (iExt, aP1, aP2);
+  // Parameter on curve
+  UV(1) = aP1.Parameter();
+  // To find parameters on surf, try ExtPS
+  Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2);
+  if (aPreciser.IsDone())
+  {
+    // Managed to find extremas between point and surface
+    Standard_Integer iPExt;
+    for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++)
+    {
+      aPreciser.Point(iPExt).Parameter(UV(2),UV(3));
+      math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
+    }
+  }
+  else
+  {
+    // Failed... try the point on iso line
+    UV(2) = dfUFirst;
+    UV(3) = aP2.Parameter();
+    math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
+  }
       } // for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
     } // if (aLocator.IsDone() && aLocator.NbExt()>0)
   } // if (myS.Type() == GeomAbs_ExtrusionSurface)
-  else 
-  // 18/02/02 akm ^^^
+  else
   {
-    Standard_Real distmin = RealLast(), distmax = 0.0, TheDist;
-
-    Standard_Integer NTmin=0,NUmin=0,NVmin=0;
-    gp_Pnt PP1min, PP2min;
-    Standard_Integer NTmax=0,NUmax=0,NVmax=0;
-    gp_Pnt PP1max, PP2max;
-    for ( NoT = 1, t1 = t0; NoT <= mytsample; NoT++, t1 += PasT) {
-      P1 = mypoints1->Value(NoT);
-      for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
-       for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
-         P2 = mypoints2->Value(NoU, NoV);
-         TheDist = P1.SquareDistance(P2);
-         if (TheDist < distmin) {
-           distmin = TheDist;
-           NTmin = NoT;
-           NUmin = NoU;
-           NVmin = NoV;
-           PP1min = P1;
-           PP2min = P2;
-         }
-         if (TheDist > distmax) {
-           distmax = TheDist;
-           NTmax = NoT;
-           NUmax = NoU;
-           NVmax = NoV;
-           PP1max = P1;
-           PP2max = P2;
-         }
-       }
+    Standard_Real aCUAdd = (mytsup - mytmin) / mytsample;
+    Standard_Real aSUAdd = (myusup - myumin) / myusample;
+    Standard_Real aSVAdd = (myvsup - myvmin) / myvsample;
+    TColgp_HArray1OfPnt aCPs(1, mytsample);
+    TColgp_HArray2OfPnt aSPs(1, myusample, 1, myvsample);
+    Standard_Integer aRestIterCount = 3;
+      // The value is calculated by the bug CR23830.
+    Standard_Integer aCUDen = 2, aSUDen = 2, aSVDen = 2;
+    Standard_Boolean anAreAvSqsInited = Standard_False;
+    Standard_Real aCUSq = 0, aSUSq = 0, aSVSq = 0;
+    while (aRestIterCount--)
+    {
+      Standard_Real aMinCU, aMinSU, aMinSV, aMaxCU, aMaxSU, aMaxSV;
+      Standard_Real aMinSqDist = DBL_MAX, aMaxSqDist = -DBL_MAX;
+      for (Standard_Integer aSUNom = 1; aSUNom < aSUDen; aSUNom += 2)
+      {
+        Standard_Real aSU0 = myumin + (aSUNom * aSUAdd) / aSUDen;
+        for (Standard_Integer aSVNom = 1; aSVNom < aSVDen; aSVNom += 2)
+        {
+          Standard_Real aSV0 = myvmin + (aSVNom * aSVAdd) / aSVDen;
+          for (Standard_Integer aCUNom = 1; aCUNom < aCUDen; aCUNom += 2)
+          {
+            Standard_Real aCU0 = mytmin + (aCUNom * aCUAdd) / aCUDen;
+            Standard_Real aCU = aCU0;
+            for (Standard_Integer aCUI = 1; aCUI <= mytsample;
+              aCUI++, aCU += aCUAdd)
+            {
+              aCPs.SetValue(aCUI, C.Value(aCU));
+            }
+            //
+            aCU = aCU0;
+            Standard_Real aSU = aSU0;
+            for (Standard_Integer aSUI = 1; aSUI <= myusample;
+              aSUI++, aSU += aSUAdd)
+            {
+              Standard_Real aSV = aSV0;
+              for (Standard_Integer aSVI = 1; aSVI <= myvsample;
+                aSVI++, aSV += aSVAdd)
+              {
+                gp_Pnt aSP = myS->Value(aSU, aSV);
+                aSPs.ChangeValue(aSUI, aSVI) = aSP;
+                Standard_Real aCU = aCU0;
+                for (Standard_Integer aCUI = 1; aCUI <= mytsample;
+                  aCUI++, aCU += aCUAdd)
+                {
+                  Standard_Real aSqDist = aSP.SquareDistance(aCPs.Value(aCUI));
+                  if (aSqDist < aMinSqDist)
+                  {
+                    aMinSqDist = aSqDist;
+                    aMinCU = aCU;
+                    aMinSU = aSU;
+                    aMinSV = aSV;
+                  }
+                  if (aMaxSqDist < aSqDist)
+                  {
+                    aMaxSqDist = aSqDist;
+                    aMaxCU = aCU;
+                    aMaxSU = aSU;
+                    aMaxSV = aSV;
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+      // Find min approximation.
+      UV(1) = aMinCU;
+      UV(2) = aMinSU;
+      UV(3) = aMinSV;
+      math_FunctionSetRoot anA(myF, UV, Tol, UVinf, UVsup, 100, aRestIterCount);
+      // Find max approximation.
+      if (!anA.IsDivergent() || !aRestIterCount)
+      {
+        UV(1) = aMaxCU;
+        UV(2) = aMaxSU;
+        UV(3) = aMaxSV;
+        math_FunctionSetRoot(myF, UV, Tol, UVinf, UVsup);
+        break;
+      }
+      //
+      if (!anAreAvSqsInited)
+      {
+        anAreAvSqsInited = Standard_True;
+        //
+        const gp_Pnt * aCP1 = &aCPs.Value(1);
+        for (Standard_Integer aCUI = 2; aCUI <= mytsample; aCUI++)
+        {
+          const gp_Pnt & aCP2 = aCPs.Value(aCUI);
+          aCUSq += aCP1->SquareDistance(aCP2);
+          aCP1 = &aCP2;
+        }
+        aCUSq /= mytsample - 1;
+        //
+        for (Standard_Integer aSVI = 1; aSVI <= myvsample; aSVI++)
+        {
+          const gp_Pnt * aSP1 = &aSPs.Value(1, aSVI);
+          for (Standard_Integer aSUI = 2; aSUI <= myusample; aSUI++)
+          {
+            const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
+            aSUSq += aSP1->SquareDistance(aSP2);
+            aSP1 = &aSP2;
+          }
+        }
+        aSUSq /= myvsample * (myusample - 1);
+        //
+        for (Standard_Integer aSUI = 1; aSUI <= myusample; aSUI++)
+        {
+          const gp_Pnt * aSP1 = &aSPs.Value(aSUI, 1);
+          for (Standard_Integer aSVI = 2; aSVI <= myvsample; aSVI++)
+          {
+            const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
+            aSVSq += aSP1->SquareDistance(aSP2);
+            aSP1 = &aSP2;
+          }
+        }
+        aSVSq /= myusample * (myvsample - 1);
+      }
+      //
+      if ((aSUSq <= aCUSq) && (aSVSq <= aCUSq))
+      {
+        aCUDen += aCUDen;
+        aCUSq /= 4;
+      }
+      else if ((aCUSq <= aSUSq) && (aSVSq <= aSUSq))
+      {
+        aSUDen += aSUDen;
+        aSUSq /= 4;
+      }
+      else
+      {
+        aSVDen += aSVDen;
+        aSVSq /= 4;
       }
     }
-  
-    // minimum:
-    UV(1) = t0 + (NTmin - 1) * PasT;
-    UV(2) = U0 + (NUmin - 1) * PasU;
-    UV(3) = V0 + (NVmin - 1) * PasV;
-
-    math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
-
-    // maximum:
-    UV(1) = t0 + (NTmax - 1) * PasT;
-    UV(2) = U0 + (NUmax - 1) * PasU;
-    UV(3) = V0 + (NVmax - 1) * PasV;
-
-    math_FunctionSetRoot S2 (myF,UV,Tol,UVinf,UVsup);
-
   }
 
   myDone = Standard_True;