From 92d22d7d626fc83b9a42b780de8cc0359ca7d097 Mon Sep 17 00:00:00 2001
From: jgv <jgv@opencascade.com>
Date: Thu, 28 Jul 2022 16:11:45 +0300
Subject: [PATCH] 0033080: Wrong projection point from ShapeAnalysis_Surface

Correct Draw command "projface" so that it really finds only projection points belonging to the face.
---
 src/SWDRAW/SWDRAW_ShapeAnalysis.cxx | 63 +++++++++++++++++++++--------
 tests/bugs/modalg_8/bug33080        | 30 ++++++++++++++
 2 files changed, 76 insertions(+), 17 deletions(-)
 create mode 100644 tests/bugs/modalg_8/bug33080

diff --git a/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx b/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx
index 26b00cbe47..d1045f3e0a 100644
--- a/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx
+++ b/src/SWDRAW/SWDRAW_ShapeAnalysis.cxx
@@ -141,13 +141,18 @@ static Standard_Integer tolerance
 static Standard_Integer projface
   (Draw_Interpretor& di, Standard_Integer argc, const char** argv)
 {
-  if (argc < 4) { di<<"Give FACE name and X Y [Z]\n"; return 1 /* Error */; }
+  if (argc < 4)
+  {
+    di << "Give FACE name and X Y [Z]\n";
+    return 1;
+  }
   Standard_CString arg1 = argv[1];
   TopoDS_Shape Shape = DBRep::Get(arg1);
   if (Shape.IsNull()) { di<<"Shape unknown : "<<arg1<<"\n"; return 1 /* Error */; }
   if (Shape.ShapeType() != TopAbs_FACE) { di<<"Not a face\n"; return 1 /* Error */; }
   TopoDS_Face F = TopoDS::Face (Shape);
   Handle(Geom_Surface) thesurf = BRep_Tool::Surface (F);  // pas locface
+  BRepTopAdaptor_FClass2d aClassifier (F, Precision::Confusion());
 //  On y va
   Standard_Real X,Y,Z,U,V;
   X = U = Draw::Atof (argv[2]);
@@ -167,28 +172,52 @@ static Standard_Integer projface
 
     GeomAPI_ProjectPointOnSurf proj(P3D, thesurf, uf-du, ul+du, vf-dv, vl+dv);
     Standard_Integer sol, nPSurf = proj.NbPoints();
-    di<<" Found "<<nPSurf<<" Points\n";
+    Standard_Integer anIndSol = 0, anIndMin = 0;
+    Standard_Real aMinDist = RealLast();
     for (sol = 1; sol <= nPSurf; sol ++) {
-      di<<"n0 "<<sol<<" Distance "<<proj.Distance(sol);
+      
       proj.Parameters(sol, U,V);
-      di<<"  U = "<<U<<"  V = "<<V<<"\n";
+      TopAbs_State aStatus = aClassifier.Perform (gp_Pnt2d (U,V));
+      if (aStatus == TopAbs_OUT)
+        continue;
+
+      anIndSol++;
+      Standard_Real aDist = proj.Distance(sol);
+      di << "n0 " << anIndSol << " Distance " << aDist;
+      di << "  U = " << U << "  V = " << V << "\n";
+
+      if (aDist < aMinDist)
+      {
+        aMinDist = aDist;
+        anIndMin = sol;
+      }
+      
 //  reprojection
       P3D = thesurf->Value (U,V);
       di<<"  => reproj  X = "<<P3D.X()<<"  Y = "<<P3D.Y()<<"  Z = "<<P3D.Z()<<"\n";
     }
-//    Que donne ShapeTool ?
-    P3D.SetCoord (X,Y,Z);
-    Handle(ShapeAnalysis_Surface) su = new ShapeAnalysis_Surface(thesurf);
-    gp_Pnt2d suval = su->ValueOfUV (P3D,BRep_Tool::Tolerance(F));
-    suval.Coord(U,V);
-    di<<"**  ShapeAnalysis_Surface gives  U = "<<U<<"  V = "<<V<<"\n";
-    P3D = thesurf->Value(U,V);
-      di<<"  => reproj  X = "<<P3D.X()<<"  Y = "<<P3D.Y()<<"  Z = "<<P3D.Z()<<"\n";
+    di<<" Found "<<anIndSol<<" Points\n";
 
-  } else {
-    di<<" Point UV  U = "<<U<<"  V = "<<V<<"\n";
-    gp_Pnt P3D = thesurf->Value(U,V);
-      di<<"  =>   proj  X = "<<P3D.X()<<"  Y = "<<P3D.Y()<<"  Z = "<<P3D.Z()<<"\n";
+    if (anIndMin != 0) //there is at least one suitable solution
+    {
+      di << "** Minimal distance to face = " << aMinDist << "\n";
+      proj.Parameters(anIndMin, U,V);
+      di << "**  Solution of minimal distance:  U = " << U << "  V = " << V << "\n";
+      P3D = thesurf->Value(U,V);
+      di<<"  => reproj  X = "<<P3D.X()<<"  Y = "<<P3D.Y()<<"  Z = "<<P3D.Z()<<"\n";
+    }
+  }
+  else //Check 2D point
+  {
+    di << " Point UV  U = " << U << "  V = " << V << "\n";
+    TopAbs_State aStatus = aClassifier.Perform (gp_Pnt2d (U,V));
+    if (aStatus == TopAbs_OUT)
+      di << "does not belong to the face" << "\n";
+    else
+    {
+      gp_Pnt P3D = thesurf->Value(U,V);
+      di << " => proj  X = " << P3D.X() << "  Y = " << P3D.Y() << "  Z = " << P3D.Z() << "\n";
+    }
   }
   return 0;
 }
@@ -1199,7 +1228,7 @@ Standard_Integer getanacurve(Draw_Interpretor& di,
   Standard_CString g = SWDRAW::GroupName();
 
   theCommands.Add ("tolerance","shape [tolmin tolmax:real]", __FILE__,tolerance,g);
-  theCommands.Add ("projface","nom_face X Y [Z]", __FILE__,projface,g);
+  theCommands.Add ("projface","nom_face X Y [Z] - returns the closest orthogonal projection if exists", __FILE__,projface,g);
   theCommands.Add ("projcurve","nom_edge | curve3d | curve3d first last + X Y Z",
 		   __FILE__,projcurve,g);
   theCommands.Add("projpcurve", "edge face tol x y z [start_param]",
diff --git a/tests/bugs/modalg_8/bug33080 b/tests/bugs/modalg_8/bug33080
new file mode 100644
index 0000000000..32e52a95c6
--- /dev/null
+++ b/tests/bugs/modalg_8/bug33080
@@ -0,0 +1,30 @@
+puts "================================"
+puts "OCC33080: Wrong projection point"
+puts "================================"
+puts ""
+
+restore [locate_data_file bug33080.brep] a
+
+set log [projface a -0.21115 1.17515 1.4504]
+
+regexp {n0 1 Distance ([0-9+-.eE]*)} $log full dist
+
+if { $dist > 0.1240364 } {
+  puts "Error distance: projection is wrong"
+}
+
+regexp {U = ([0-9+-.eE]*)} $log full uparam
+
+if { $uparam < 869.174321 || $uparam > 869.174322 } {
+  puts "Error point: projection is wrong"
+}
+
+regexp {V = ([0-9+-.eE]*)} $log full vparam
+
+if { $vparam < 732.602489 || $vparam > 732.602490 } {
+  puts "Error point: projection is wrong"
+}
+
+if {![regexp "Found 1 Points" $log]} {
+  puts "Error number of points: projection is wrong"
+}