OCC22247 BRepMesh pure performance
[occt.git] / src / BRepMesh / BRepMesh_FastDiscretFace.cxx
index 8f96927..f79dfc0 100755 (executable)
@@ -648,6 +648,80 @@ Standard_Boolean BRepMesh_FastDiscretFace::Update(const TopoDS_Edge&  edge,
 //function : InternalVertices
 //purpose  : 
 //=======================================================================
+
+
+static void filterParameters(const TColStd_IndexedMapOfReal& theParams,
+                             const Standard_Real theMinDist,
+                             const Standard_Real theFilterDist,
+                             TColStd_SequenceOfReal& theResult)
+{
+  // Sort sequence of parameters
+  TColStd_SequenceOfReal aParamTmp;
+  Standard_Integer aParamLength = 1;
+  const Standard_Integer anInitLen = theParams.Extent();
+    
+  TColStd_Array1OfReal aParamArray(1, anInitLen);
+  Standard_Integer j;
+  for (j = 1; j <= anInitLen; j++)
+    aParamArray(j) = theParams(j);
+
+  TCollection_CompareOfReal aCompare;
+  SortTools_ShellSortOfReal::Sort(aParamArray, aCompare);
+
+  // mandadory pre filtering using the first (minimal) filter value
+  Standard_Real aP1, aP2;
+  aP1 = aParamArray(1);
+  aParamTmp.Append(aP1);
+  for (j = 2; j <= anInitLen; j++) 
+  {
+    aP2 = aParamArray(j);
+    if ((aP2-aP1) > theMinDist)
+    {
+      aParamTmp.Append(aP2);
+      aP1 = aP2;
+      aParamLength++;
+    }
+  }
+
+  //add last point if required
+  if(aParamArray(anInitLen)-theParams(aParamLength) > theMinDist)
+  {
+    aParamTmp.Append(aParamArray(anInitLen));        
+    aParamLength++;
+  }
+  
+  //perform filtering on series
+  Standard_Real aLastAdded, aLastCandidate;
+  Standard_Boolean isCandidateDefined = Standard_False;
+  aLastAdded = aParamTmp.First();
+  aLastCandidate = aLastAdded;
+  theResult.Append(aParamTmp.First());
+  
+  for(j=2;j<aParamTmp.Length();j++) 
+  {
+    Standard_Real aVal = aParamTmp.Value(j);
+    if(aVal-aLastAdded > theFilterDist) 
+    {
+      //adds the parameter
+      if(isCandidateDefined) {
+        aLastAdded = aLastCandidate;
+        isCandidateDefined = Standard_False;
+        j--;
+      }
+      else 
+      {
+        aLastAdded = aVal;
+      }
+      theResult.Append(aLastAdded);
+      continue;
+    }
+    
+    aLastCandidate = aVal;
+    isCandidateDefined = Standard_True;
+  }
+  theResult.Append(aParamTmp.Last());
+}
+
 void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurface)& caro,
                                                   BRepMesh_ListOfVertex&              InternalV,
                                                   const Standard_Real                 defface,
@@ -938,101 +1012,49 @@ void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurf
   }
   else if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
   {
-    Standard_Integer i, j;
+    //define resolutions
+    Standard_Real uRes = BS.UResolution(defface);
+    Standard_Real vRes = BS.VResolution(defface);
+
+    // Sort and filter sequence of parameters
+    Standard_Real aMinDu = Precision::PConfusion();
+    if(deltaX < 1.)
+      aMinDu/=deltaX;
+
+    Standard_Real aDuMaxLim = 0.1*(umax-umin);
+    Standard_Real ddu = Min(aDuMaxLim,Max(0.005*(umax-umin),2.*uRes));
+    TColStd_SequenceOfReal ParamU; 
+    filterParameters(myUParam,aMinDu,ddu,ParamU);
+    Standard_Integer ParamULength = ParamU.Length();
+   
+    Standard_Real aMinDv = Precision::PConfusion();
+    if(deltaY < 1)
+      aMinDv/=deltaY;
+
+    Standard_Real aDvMaxLim = 0.1*(vmax-vmin);
+    Standard_Real ddv = Min(aDvMaxLim,Max(0.005*(vmax-vmin),2.*vRes));
+    TColStd_SequenceOfReal ParamV; 
+    filterParameters(myVParam,aMinDv,ddv,ParamV);
+    Standard_Integer ParamVLength = ParamV.Length();
 
-    Handle(TColStd_HArray1OfReal) anUKnots, anVKnots;
-    Standard_Integer aNbUNkots,aNbVNkots;
+    // check intermediate isolines
     Handle(Geom_Surface) B;
     if (thetype == GeomAbs_BezierSurface) {
-      anUKnots = new TColStd_HArray1OfReal(1,2);
-      anVKnots = new TColStd_HArray1OfReal(1,2);
-      anUKnots->SetValue(1,0.);
-      anUKnots->SetValue(2,1.);
-      anVKnots->SetValue(1,0.);
-      anVKnots->SetValue(2,1.);
-      aNbUNkots = 2;
-      aNbVNkots = 2;
       B = BS.Bezier();
     }
     else {
-      Handle(Geom_BSplineSurface) aSurf = BS.BSpline();
-      B = aSurf;
-      aNbUNkots = aSurf->NbUKnots();
-      aNbVNkots = aSurf->NbVKnots();
-      anUKnots = new TColStd_HArray1OfReal(1,aNbUNkots);
-      anVKnots = new TColStd_HArray1OfReal(1,aNbVNkots);
-
-      aSurf->UKnots(anUKnots->ChangeArray1());
-      aSurf->VKnots(anVKnots->ChangeArray1());
-    }
-
-    Standard_Real ddu = caro->UResolution(defface)*5.;
-
-    Standard_Real aDUmin = (umax-umin)/5.;
-    if(ddu > aDUmin)
-      ddu = aDUmin;
-    // Sort sequence of U parameters
-    TColStd_SequenceOfReal ParamU;
-    Standard_Integer anUdegree = caro->UDegree();
-     
-    Standard_Real aU1 = anUKnots->Value(1);
-    for( i = 1; i < aNbUNkots; i++)
-    {
-      Standard_Real aStart = anUKnots->Value(i);
-      Standard_Real aEnd = anUKnots->Value(i+1);
-      Standard_Integer aNbPoints = anUdegree-1;
-      Standard_Real aStep = (aEnd-aStart)/(aNbPoints+1);
-      for( Standard_Integer anInd = 0; anInd<= aNbPoints; anInd++)
-      {
-        Standard_Real aU2 = aStart+anInd*aStep;
-        if((aU2 - aU1) > ddu)
-        {
-          ParamU.Append(aU2);
-          aU1 = aU2;
-        }
-      }
+      B = BS.BSpline();
     }
-    ParamU.Append(anUKnots->Value(aNbUNkots));
-    Standard_Integer ParamULength = ParamU.Length();
-
-    Standard_Real ddv = caro->VResolution(defface)*5.;
-
-    Standard_Real aDVmin = (vmax-vmin)/5.;
-    if(ddv > aDVmin)
-      ddv = aDVmin;
-    // Sort sequence of V parameters
-    TColStd_SequenceOfReal ParamV; 
-    Standard_Integer anVdegree = caro->VDegree();
-      
-    Standard_Real aV1 = anVKnots->Value(1);
-    for( i = 1; i < aNbVNkots; i++)
-    {
-      Standard_Real aStart = anVKnots->Value(i);
-      Standard_Real aEnd = anVKnots->Value(i+1);
-      Standard_Integer aNbPoints = anVdegree-1;
-      Standard_Real aStep = (aEnd-aStart)/(aNbPoints+1);
-      for( Standard_Integer anInd = 0; anInd<= aNbPoints; anInd++)
-      {
-        Standard_Real aV2 = aStart+anInd*aStep;
-        if ((aV2 - aV1) > ddv)
-        { 
-          ParamV.Append(aV2);
-          aV1 = aV2;
-        }
-      }
-    }
-    ParamV.Append(anVKnots->Value(aNbVNkots));
-    Standard_Integer ParamVLength = ParamV.Length();
-
-    // controle des isos U et insertion eventuelle:
 
     gp_Pnt P1, P2, PControl;
-    Standard_Real u, v, dist, V1, V2, U1, U2;
+    Standard_Real u, v, dist;
 
     // precision for compare square distances
     double dPreci = Precision::Confusion()*Precision::Confusion();
 
     // Insert V parameters by deflection criterion
+    Standard_Integer i,j;
+    Standard_Real V1, V2, U1, U2;
     for (i = 1; i <= ParamULength; i++) {
       Handle(Geom_Curve) IsoU = B->UIso(ParamU.Value(i));
       V1 = ParamV.Value(1);
@@ -1056,9 +1078,22 @@ void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurf
           ParamVLength++;
         }
         else {
-          V1 = V2;
-          P1 = P2;
-          j++;
+          //put regular grig for normals
+          gp_Dir N1(0,0,1),N2(0,0,1);
+          Standard_Boolean aSt1 = GeomLib::NormEstim(B, gp_Pnt2d(ParamU.Value(i),V1), Precision::Confusion(), N1);
+          Standard_Boolean aSt2 = GeomLib::NormEstim(B, gp_Pnt2d(ParamU.Value(i),v), Precision::Confusion(), N2);
+          
+          Standard_Real anAngle1 = N2.Angle(N1);
+               if(aSt1 < 1 && aSt2 < 1 && anAngle1 > angle ) {
+            // insertion 
+            ParamV.InsertBefore(j, v);
+            ParamVLength++;
+          }
+          else {
+            V1 = V2;
+            P1 = P2;
+            j++;
+          }
              }
       }
     }
@@ -1088,22 +1123,35 @@ void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurf
         }
         else {
           //check normal
-          U1 = U2;
-          P1 = P2;
-          if (j < ParamULength) {
-            // Classify intersection point
-            if (classifier->Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
-            {
-              // Record 3d point
-              nbLocat++;
-              Location3d.Bind(nbLocat, P1);
-              // Record 2d point
-              p2d.SetCoord((U1-umin)/deltaX, (v-vmin)/deltaY);
-              newV.Initialize(p2d.XY(), nbLocat, MeshDS_Free);
-              InternalV.Append(newV);
+          //put regular grig for normals
+          gp_Dir N1(0,0,1),N2(0,0,1);
+          Standard_Boolean aSt1 = GeomLib::NormEstim(B, gp_Pnt2d(U1,v), Precision::Confusion(), N1);
+          Standard_Boolean aSt2 = GeomLib::NormEstim(B, gp_Pnt2d(u,v), Precision::Confusion(), N2);
+          
+          Standard_Real anAngle1 = N2.Angle(N1);
+               if(aSt1 < 1 && aSt2 < 1 && anAngle1 > angle) {
+            // insertion 
+            ParamU.InsertBefore(j, u);
+            ParamULength++;
+          }
+          else {
+            U1 = U2;
+            P1 = P2;
+            if (j < ParamULength) {
+              // Classify intersection point
+              if (classifier->Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
+              {
+                // Record 3d point
+                nbLocat++;
+                Location3d.Bind(nbLocat, P1);
+                // Record 2d point
+                p2d.SetCoord((U1-umin)/deltaX, (v-vmin)/deltaY);
+                newV.Initialize(p2d.XY(), nbLocat, MeshDS_Free);
+                InternalV.Append(newV);
+              }
             }
+            j++;
           }
-          j++;
              }
       }
     }
@@ -1112,7 +1160,7 @@ void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurf
     const Standard_Real theangle = 0.35;
 
     Standard_Integer i, j, nbpointsU = 10, nbpointsV = 10;
-    Adaptor3d_IsoCurve tabu[10], tabv[10];
+    Adaptor3d_IsoCurve tabu[11], tabv[11];
 
     TColStd_SequenceOfReal ParamU, ParamV;
     Standard_Real u, v, du, dv;
@@ -1121,23 +1169,23 @@ void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurf
 
     du = (umax-umin) / (nbpointsU+1); dv = (vmax-vmin) / (nbpointsV+1);
 
-    for (iu = 1; iu <= nbpointsU; iu++) {
+    for (iu = 0; iu <= nbpointsU; iu++) {
       u = umin + iu*du;
-      tabu[iu-1].Load(caro);
-      tabu[iu-1].Load(GeomAbs_IsoU, u);
+      tabu[iu].Load(caro);
+      tabu[iu].Load(GeomAbs_IsoU, u);
     }
 
-    for (iv = 1; iv <= nbpointsV; iv++) {
+    for (iv = 0; iv <= nbpointsV; iv++) {
       v = vmin + iv*dv;
-      tabv[iv-1].Load(caro);
-      tabv[iv-1].Load(GeomAbs_IsoV, v);
+      tabv[iv].Load(caro);
+      tabv[iv].Load(GeomAbs_IsoV, v);
     }
 
     Standard_Integer imax = 1, MaxV = 0;
     
-    GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU];
+    GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU+1];
 
-    for (i = 0; i <= nbpointsU-1; i++) {
+    for (i = 0; i <= nbpointsU; i++) {
       f = Max(vmin, tabu[i].FirstParameter());
       l = Min(vmax, tabu[i].LastParameter());
       GCPnts_TangentialDeflection theDeflection(tabu[i], f, l, theangle, 0.7*defface, 2);
@@ -1158,9 +1206,9 @@ void BRepMesh_FastDiscretFace::InternalVertices ( const Handle(BRepAdaptor_HSurf
     imax = 1;
     Standard_Integer MaxU = 0;
 
-    GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV];
+    GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV+1];
 
-    for (i = 0; i <= nbpointsV-1; i++) {
+    for (i = 0; i <= nbpointsV; i++) {
       f = Max(umin, tabv[i].FirstParameter());
       l = Min(umax, tabv[i].LastParameter());
       GCPnts_TangentialDeflection thedeflection2(tabv[i], f, l, theangle, 0.7*defface, 2);
@@ -1286,9 +1334,6 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
   gp_Pnt pDef;
   Standard_Real dv = 0, defl = 0, maxdef = -1;
   Standard_Integer pass = 1, nf = 0, nl = 0;
-  Standard_Integer v1, v2, v3, e1, e2, e3;
-  Standard_Boolean o1, o2, o3;
-  Standard_Boolean m1, m2, m3;
   BRepMesh_Vertex InsVertex;
   Standard_Boolean caninsert;
 
@@ -1414,7 +1459,11 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
       // Check deflection on triangle
       mi2d = (xy1+xy2+xy3)/3.0;
       caro->D0(mi2d.X(), mi2d.Y(), pDef);
-      defl = pDef.SquareDistance((p1+p2+p3)/3.);
+      defl = Abs(normal*(pDef.XYZ()-p1));
+      defl = defl*defl;     
+      /*mi2d = (xy1+xy2+xy3)/3.0;
+      caro->D0(mi2d.X(), mi2d.Y(), pDef);
+      defl = pDef.SquareDistance((p1+p2+p3)/3.);*/
       if (defl > maxdef) maxdef = defl;
       if (defl > sqdefface)
       {
@@ -1441,7 +1490,8 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
           // Check deflection on edge 1
           mi2d = (xy2+xy3)*0.5;
           caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          defl = pDef.SquareDistance((p2+p3)/2.);
+          gp_Lin L (p2, gp_Vec(p2, p3));
+          defl = L.SquareDistance(pDef);
           if (defl > maxdef) maxdef = defl;
           if (defl > sqdefface)
           {
@@ -1470,7 +1520,8 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
           // Check deflection on edge 2
           mi2d = (xy3+xy1)*0.5;
           caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          defl = pDef.SquareDistance((p1+p3)/2.);
+          gp_Lin L (p1, gp_Vec(p1, p3));
+          defl = L.SquareDistance(pDef);
           if (defl > maxdef) maxdef = defl;
           if (defl > sqdefface)
           {
@@ -1499,7 +1550,8 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
           // Check deflection on edge 3
           mi2d = (xy1+xy2)*0.5;
           caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          defl = pDef.SquareDistance((p1+p2)/2.);
+          gp_Lin L (p1, gp_Vec(p1, p2));
+          defl = L.SquareDistance(pDef);
           if (defl > maxdef) maxdef = defl;
           if (defl > sqdefface)
           {
@@ -1519,8 +1571,8 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
         }
       }
       
-      //check normal
-      if(isSpline && !BSpl.IsNull() && (badTriangles.IsEmpty() || badTriangles.Last() != TriId))      
+      //check normal on bsplines
+      if(isfirst && isSpline && !BSpl.IsNull() )      
       {
              gp_Dir N1(0,0,1), N2(0,0,1), N3(0,0,1);
         Standard_Integer aSt1, aSt2, aSt3;
@@ -1559,71 +1611,8 @@ Standard_Real BRepMesh_FastDiscretFace::Control (const Handle(BRepAdaptor_HSurfa
         Standard_Real anAngle3 = N1.Angle(N3);
              if(aSt1 < 1 && aSt2 < 1 && aSt3 < 1 && 
           (anAngle1 > angle || anAngle2 > angle || anAngle3 > angle)) {
-           if (isfirst) {
             maxdef = -1;
             break;
-          }
-          // Record new vertex
-          aNbPnt++;
-          nbLocat++;
-          mi2d = (xy1+xy2+xy3)/3.;
-          caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          Location3d.Bind(nbLocat,pDef);
-          mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);              
-          InsVertex.Initialize(mi2d,nbLocat,MeshDS_Free);
-          InternalV.Append(InsVertex);
-          badTriangles.Append(TriId);
-             }   
-      
-        //In case if triangle is considerd as OK, we have to check three intermediate 
-        //points to be sure that we free from wave effect. If it is OK triangle passed if not split in middle point
-        if(badTriangles.IsEmpty() || badTriangles.Last() != TriId)
-        {
-          Bnd_Box2d aB2d;
-          aB2d.Add(gp_Pnt2d(xy1)); aB2d.Add(gp_Pnt2d(xy2)); aB2d.Add(gp_Pnt2d(xy3));
-
-          Standard_Real aXMin1, aXMax1, aYMin1, aYMax1;
-          aB2d.Get(aXMin1, aYMin1, aXMax1, aYMax1);
-
-          if(aXMax1 - aXMin1 < ddu && aYMax1 - aYMin1 < ddv)
-            continue;
-          
-          mi2d = (xy1+xy2+xy3)/3.;
-          gp_Pnt2d aP[3] = {mi2d+(xy1-mi2d)*2/3., mi2d+(xy2-mi2d)*2/3, mi2d+(xy3-mi2d)*2/3.};
-          gp_Dir midDir(0,0,1);
-          Standard_Integer aSt[4];
-          aSt[0] = GeomLib::NormEstim(BSpl, gp_Pnt2d(mi2d), Precision::Confusion(), midDir);
-          Standard_Integer i = 0;
-          gp_Dir dir[3] = {gp_Dir(0,0,1), gp_Dir(0,0,1), gp_Dir(0,0,1)};
-          Standard_Real anAngle[3];
-          for(; i < 3; i++)
-          {
-            gp_Dir dir(0,0,1);
-            aSt[i+1] = GeomLib::NormEstim(BSpl, aP[i], Precision::Confusion(), dir);
-            anAngle[i] = dir.Angle(midDir);
-          }
-          
-          caro->D0(mi2d.X(), mi2d.Y(), pDef);
-
-          aB2d.Add(gp_Pnt2d(xy1)); aB2d.Add(gp_Pnt2d(xy2)); aB2d.Add(gp_Pnt2d(xy3));
-          
-          if(aSt[0] < 1 && aSt[1] < 1 && aSt[2] < 1 && aSt[3] < 1 &&
-            (anAngle[0] > angle || anAngle[1] > angle || anAngle[2] > angle) &&
-            (aXMax1 - aXMin1 > ddu || aYMax1 - aYMin1 > ddv))
-          {
-            if (isfirst) {
-              maxdef = -1;
-              break;
-            }
-            aNbPnt++;
-            nbLocat++;
-            caro->D0(mi2d.X(), mi2d.Y(), pDef);
-            Location3d.Bind(nbLocat,pDef);
-            mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);              
-            InsVertex.Initialize(mi2d,nbLocat,MeshDS_Free);
-            InternalV.Append(InsVertex);            
-            badTriangles.Append(TriId);
-          }
         }
       }
     }