Integrate fixes for issues 32995, 33104

src/GeomConvert/FILES

@@ -18,3 +18,14 @@ GeomConvert_CompBezierSurfacesToBSplineSurface.hxx
 GeomConvert_CompBezierSurfacesToBSplineSurface.lxx
 GeomConvert_CompCurveToBSplineCurve.cxx
 GeomConvert_CompCurveToBSplineCurve.hxx
+GeomConvert_CurveToAnaCurve.cxx
+GeomConvert_CurveToAnaCurve.hxx
+GeomConvert_SurfToAnaSurf.cxx
+GeomConvert_SurfToAnaSurf.hxx
+GeomConvert_ConvType.hxx
+GeomConvert_FuncSphereLSDist.cxx
+GeomConvert_FuncSphereLSDist.hxx
+GeomConvert_FuncCylinderLSDist.cxx
+GeomConvert_FuncCylinderLSDist.hxx
+GeomConvert_FuncConeLSDist.cxx
+GeomConvert_FuncConeLSDist.hxx

src/GeomConvert/GeomConvert_ConvType.hxx

@@ -0,0 +1,23 @@
+// Copyright (c) 2022 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.
+
+#ifndef _GeomConvert_ConvType_HeaderFile
+#define _GeomConvert_ConvType_HeaderFile
+enum GeomConvert_ConvType
+{
+  GeomConvert_Target,
+  GeomConvert_Simplest,
+  GeomConvert_MinGap
+};
+
+#endif // _GeomConvert_ConvType_HeaderFile

src/GeomConvert/GeomConvert_CurveToAnaCurve.cxx

@@ -0,0 +1,768 @@
+// Created: 2001-05-21 
+// 
+// Copyright (c) 2001-2013 OPEN CASCADE SAS 
+// 
+// This file is part of commercial software by OPEN CASCADE SAS, 
+// furnished in accordance with the terms and conditions of the contract 
+// and with the inclusion of this copyright notice. 
+// This file or any part thereof may not be provided or otherwise 
+// made available to any third party. 
+// 
+// No ownership title to the software is transferred hereby. 
+// 
+// OPEN CASCADE SAS makes no representation or warranties with respect to the 
+// performance of this software, and specifically disclaims any responsibility 
+// for any damages, special or consequential, connected with its use. 
+
+
+#include <ElCLib.hxx>
+#include <gce_MakeCirc.hxx>
+#include <Geom_BezierCurve.hxx>
+#include <Geom_BSplineCurve.hxx>
+#include <Geom_Circle.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Ellipse.hxx>
+#include <Geom_Line.hxx>
+#include <Geom_TrimmedCurve.hxx>
+#include <gp_Ax2.hxx>
+#include <gp_Ax3.hxx>
+#include <gp_Circ.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <Precision.hxx>
+#include <GeomConvert_CurveToAnaCurve.hxx>
+#include <TColgp_Array1OfPnt.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <TColStd_Array2OfReal.hxx>
+#include <GeomAbs_CurveType.hxx>
+#include <math_Vector.hxx>
+#include <math_Matrix.hxx>
+#include <math_Gauss.hxx>
+
+
+GeomConvert_CurveToAnaCurve::GeomConvert_CurveToAnaCurve():
+  myGap(Precision::Infinite()),
+  myConvType(GeomConvert_MinGap),
+  myTarget(GeomAbs_Line)
+{
+}
+
+GeomConvert_CurveToAnaCurve::GeomConvert_CurveToAnaCurve(const Handle(Geom_Curve)& C) :
+  myGap(Precision::Infinite()),
+  myConvType(GeomConvert_MinGap),
+  myTarget(GeomAbs_Line)
+{
+  myCurve = C;
+}
+
+void GeomConvert_CurveToAnaCurve::Init(const Handle(Geom_Curve)& C) 
+{
+  myCurve = C;
+  myGap = Precision::Infinite();
+}
+
+//=======================================================================
+//function : ConvertToAnalytical
+//purpose  : 
+//=======================================================================
+
+Standard_Boolean GeomConvert_CurveToAnaCurve::ConvertToAnalytical(const Standard_Real tol, 
+                                                         Handle(Geom_Curve)& theResultCurve, 
+                                                         const Standard_Real F, const Standard_Real L,
+                                                         Standard_Real& NewF, Standard_Real& NewL)
+{
+  if(myCurve.IsNull())
+    return Standard_False;
+  
+  Handle(Geom_Curve) aCurve = myCurve;
+  while (aCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
+    Handle(Geom_TrimmedCurve) aTrimmed = Handle(Geom_TrimmedCurve)::
+      DownCast(aCurve);
+    aCurve = aTrimmed->BasisCurve();
+  }
+  
+  Handle(Geom_Curve) C = ComputeCurve(aCurve,tol,F, L, NewF, NewL, myGap, myConvType, myTarget);
+
+  if(C.IsNull()) return Standard_False;
+  theResultCurve = C; 
+  return Standard_True;
+}
+
+//=======================================================================
+//function : IsLinear
+//purpose  : 
+//=======================================================================
+
+Standard_Boolean GeomConvert_CurveToAnaCurve::IsLinear(const TColgp_Array1OfPnt& aPoles,
+                                              const Standard_Real tolerance,
+                                              Standard_Real& Deviation)
+{
+  Standard_Integer nbPoles = aPoles.Length();
+  if(nbPoles < 2)
+    return Standard_False;
+  
+  Standard_Real dMax = 0;
+  Standard_Integer iMax1=0,iMax2=0;
+  
+  Standard_Integer i;
+  for(i = 1; i < nbPoles; i++)
+    for(Standard_Integer j = i+1; j <= nbPoles; j++) {
+      Standard_Real dist = aPoles(i).SquareDistance(aPoles(j));
+      if(dist > dMax) {
+        dMax = dist;
+        iMax1 = i;
+        iMax2 = j;
+      }
+    }
+  
+  if (dMax < Precision::SquareConfusion())
+    return Standard_False;
+  
+  Standard_Real tol2 = tolerance*tolerance;
+  gp_Vec avec (aPoles(iMax1),aPoles(iMax2));  gp_Dir adir (avec);  gp_Lin alin (aPoles(iMax1),adir);
+  
+  Standard_Real aMax = 0.;
+  for(i = 1; i <= nbPoles; i++) {
+    Standard_Real dist = alin.SquareDistance(aPoles(i));
+    if(dist > tol2)
+      return Standard_False;
+    if(dist > aMax)
+      aMax = dist;
+  }
+  Deviation = sqrt(aMax);
+  
+  return Standard_True;
+}
+
+//=======================================================================
+//function : GetLine
+//purpose  : 
+//=======================================================================
+
+gp_Lin GeomConvert_CurveToAnaCurve::GetLine(const gp_Pnt& P1, const gp_Pnt& P2,
+                                             Standard_Real& cf, Standard_Real& cl)
+{
+  gp_Vec avec(P1, P2);  gp_Dir adir(avec);  gp_Lin alin(P1, adir);
+  cf = ElCLib::Parameter(alin, P1);
+  cl = ElCLib::Parameter(alin, P2);
+  return alin;
+}
+//=======================================================================
+//function : ComputeLine
+//purpose  : 
+//=======================================================================
+
+Handle(Geom_Line) GeomConvert_CurveToAnaCurve::ComputeLine (const Handle(Geom_Curve)& curve,
+                                                   const Standard_Real tolerance,
+                                                   const Standard_Real c1, const Standard_Real c2,
+                                                   Standard_Real& cf, Standard_Real& cl,
+                                                   Standard_Real& Deviation)
+{
+  Handle(Geom_Line) line;
+  if (curve.IsNull()) return line;
+  line = Handle(Geom_Line)::DownCast(curve);  // qui sait
+  if (!line.IsNull()) {
+    cf = c1;
+    cl = c2;
+    Deviation = 0.;
+    return line;
+  }
+
+  gp_Pnt P1 = curve->Value (c1);
+  gp_Pnt P2 = curve->Value (c2);
+  if(P1.SquareDistance(P2) < Precision::SquareConfusion())
+    return line;
+  cf = c1;  cl = c2;
+  
+  Handle(TColgp_HArray1OfPnt) Poles;
+  Standard_Integer nbPoles;
+  Handle(Geom_BSplineCurve) bsc = Handle(Geom_BSplineCurve)::DownCast(curve);
+  if (!bsc.IsNull()) {
+    nbPoles = bsc->NbPoles();
+    Poles = new TColgp_HArray1OfPnt(1, nbPoles);
+    bsc->Poles(Poles->ChangeArray1());
+  }
+  else
+  {
+    Handle(Geom_BezierCurve) bzc = Handle(Geom_BezierCurve)::DownCast(curve);
+    if (!bzc.IsNull()) {
+      nbPoles = bzc->NbPoles();
+      Poles = new TColgp_HArray1OfPnt(1, nbPoles);
+      bzc->Poles(Poles->ChangeArray1());
+    }
+    else
+    {
+      nbPoles = 23;
+      Poles = new TColgp_HArray1OfPnt(1, nbPoles);
+      Standard_Real dt = (c2 - c1) / (nbPoles - 1);
+      Poles->SetValue(1, P1);
+      Poles->SetValue(nbPoles, P2);
+      Standard_Integer i;
+      for (i = 2; i < nbPoles; ++i)
+      {
+        Poles->SetValue(i, curve->Value(c1 + (i - 1) * dt));
+      }
+    }
+  }
+  if(!IsLinear(Poles->Array1(),tolerance,Deviation)) return line;  // non
+  gp_Lin alin = GetLine (P1, P2, cf, cl);
+  line = new Geom_Line (alin);
+  return line;
+}
+
+//=======================================================================
+//function : GetCircle
+//purpose  : 
+//=======================================================================
+
+Standard_Boolean GeomConvert_CurveToAnaCurve::GetCircle (gp_Circ& crc,   
+                                                const gp_Pnt& P0,const gp_Pnt& P1, const gp_Pnt& P2) 
+{
+//  Control if points are not aligned (should be done by MakeCirc 
+  Standard_Real aMaxCoord = Sqrt(Precision::Infinite());
+  if (Abs(P0.X()) > aMaxCoord || Abs(P0.Y()) > aMaxCoord || Abs(P0.Z()) > aMaxCoord)
+    return Standard_False;
+  if (Abs(P1.X()) > aMaxCoord || Abs(P1.Y()) > aMaxCoord || Abs(P1.Z()) > aMaxCoord)
+    return Standard_False;
+  if (Abs(P2.X()) > aMaxCoord || Abs(P2.Y()) > aMaxCoord || Abs(P2.Z()) > aMaxCoord)
+    return Standard_False;
+
+//  Building the circle
+  gce_MakeCirc mkc (P0,P1,P2);
+  if (!mkc.IsDone()) return Standard_False;
+  crc = mkc.Value();
+  if (crc.Radius() < gp::Resolution()) return Standard_False;
+  //  Recalage sur P0
+  gp_Pnt PC = crc.Location();
+  gp_Ax2 axe = crc.Position();
+  gp_Vec VX (PC,P0);
+  axe.SetXDirection (VX);
+  crc.SetPosition (axe);
+  return Standard_True;
+}
+
+//=======================================================================
+//function : ComputeCircle
+//purpose  : 
+//=======================================================================
+
+Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeCircle (const Handle(Geom_Curve)& c3d, 
+                                                      const Standard_Real tol,
+                                                      const Standard_Real c1, const Standard_Real c2,
+                                                      Standard_Real& cf, Standard_Real& cl,
+                                                      Standard_Real& Deviation)
+{
+  if (c3d->IsKind (STANDARD_TYPE(Geom_Circle))) {
+    cf = c1;
+    cl = c2;
+    Deviation = 0.;
+    Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(c3d);
+    return aCirc;
+  }
+  
+  Handle(Geom_Circle) circ;
+  gp_Pnt P0,P1,P2;
+  Standard_Real ca = (c1+c1+c2) / 3;  Standard_Real cb = (c1+c2+c2) / 3;
+  P0 = c3d->Value(c1);
+  P1 = c3d->Value(ca);
+  P2 = c3d->Value(cb);
+  
+  gp_Circ crc;
+  if (!GetCircle (crc,P0,P1,P2)) return circ;
+  
+  //  Reste a controler que c est bien un cercle : prendre 20 points
+  Standard_Real du = (c2-c1)/20;
+  Standard_Integer i;
+  Standard_Real aMax = 0.;
+  for (i = 0; i <= 20; i ++) {
+    Standard_Real u = c1+(du*i);
+    gp_Pnt PP = c3d->Value(u);
+    Standard_Real dist = crc.Distance(PP);
+    if (dist > tol) return circ;  // not done
+    if (dist > aMax)
+      aMax = dist;
+  }
+  Deviation = aMax;
+  
+  // defining the parameters
+  Standard_Real PI2 = 2 * M_PI;
+  
+  cf = ElCLib::Parameter (crc,c3d->Value (c1));
+  cf = ElCLib::InPeriod(cf, 0., PI2);
+
+  //first parameter should be closed to zero
+  
+  if(Abs(cf) < Precision::PConfusion() || Abs(PI2-cf) < Precision::PConfusion())
+    cf = 0.;
+
+  Standard_Real cm = ElCLib::Parameter (crc,c3d->Value ((c1+c2)/2.));
+  cm = ElCLib::InPeriod(cm, cf, cf + PI2);
+
+  cl = ElCLib::Parameter (crc,c3d->Value (c2));
+  cl = ElCLib::InPeriod(cl, cm, cm + PI2);
+
+  circ = new Geom_Circle (crc);
+  return circ;
+}
+
+//=======================================================================
+//              Compute Ellipse
+//=======================================================================
+
+//=======================================================================
+//function : IsArrayPntPlanar
+//purpose  : 
+//=======================================================================
+
+static Standard_Boolean IsArrayPntPlanar(const Handle(TColgp_HArray1OfPnt)& HAP,
+                                         gp_Dir& Norm, const Standard_Real prec)
+{
+  Standard_Integer size = HAP->Length();
+  if(size<3)
+    return Standard_False;
+  gp_Pnt P1 = HAP->Value(1);
+  gp_Pnt P2 = HAP->Value(2);
+  gp_Pnt P3 = HAP->Value(3);
+  Standard_Real dist1 = P1.Distance(P2);
+  Standard_Real dist2 = P1.Distance(P3);
+  if( dist1<prec || dist2<prec )
+    return Standard_False;
+  gp_Vec V1(P1,P2);
+  gp_Vec V2(P1,P3);
+  if(V1.IsParallel(V2,prec))
+    return Standard_False;
+  gp_Vec NV = V1.Crossed(V2);
+
+  Standard_Integer i;
+  for (i = 1; i <= 3; ++i)
+  {
+    if (Precision::IsInfinite(NV.Coord(i)))
+      return Standard_False;
+  }
+
+  if(NV.Magnitude() < gp::Resolution())
+    return Standard_False;
+
+  if(size>3) {
+    for(i=4; i<=size; i++) {
+      gp_Pnt PN = HAP->Value(i);
+      dist1 = P1.Distance(PN);
+      if (dist1 < prec || Precision::IsInfinite(dist1))
+      {
+        return Standard_False;
+      }
+      gp_Vec VN(P1,PN);
+      if(!NV.IsNormal(VN,prec))
+        return Standard_False;
+    }
+  }
+  Norm = NV;
+  return Standard_True;
+}
+
+//=======================================================================
+//function : ConicdDefinition
+//purpose  : 
+//=======================================================================
+
+static Standard_Boolean ConicDefinition
+       ( const Standard_Real a, const Standard_Real b1, const Standard_Real c,
+         const Standard_Real d1, const Standard_Real e1, const Standard_Real f,
+         const Standard_Boolean IsParab, const Standard_Boolean IsEllip,
+         gp_Pnt& Center, gp_Dir& MainAxis, Standard_Real& Rmin, Standard_Real& Rmax )
+{
+  Standard_Real Xcen = 0.,Ycen = 0., Xax = 0.,Yax = 0.;
+  Standard_Real b,d,e;
+  //  conic : a*x2 + 2*b*x*y + c*y2 + 2*d*x + 2*e*y + f = 0.
+  //Equation (a,b,c,d,e,f);
+  b = b1/2.;  d = d1/2.;  e = e1/2.;    // chgt de variable
+
+  Standard_Real eps = 1.E-08;  // ?? comme ComputedForm
+
+  if (IsParab) {
+
+  } 
+  else {
+    //   -> Conique a centre, cas general
+    //  On utilise les Determinants des matrices :
+    //               | a b d |
+    //  gdet (3x3) = | b c e |  et pdet (2X2) = | a b |
+    //               | d e f |                  | b c |
+    
+    Standard_Real gdet = a*c*f + 2*b*d*e - c*d*d - a*e*e - b*b*f;
+    Standard_Real pdet = a*c - b*b;
+    
+    Xcen = (b*e - c*d) / pdet;
+    Ycen = (b*d - a*e) / pdet;
+
+    Standard_Real term1 = a-c;
+    Standard_Real term2 = 2*b;
+    Standard_Real cos2t;
+    Standard_Real auxil;
+
+    if (Abs(term2) <= eps && Abs(term1) <= eps) {
+      cos2t = 1.;
+      auxil = 0.;
+    }
+    else {
+      if (Abs(term1) < eps)
+      {
+        return Standard_False;
+      }
+      Standard_Real t2d = term2/term1; //skl 21.11.2001
+      cos2t = 1./sqrt(1+t2d*t2d);
+      auxil = sqrt (term1*term1 + term2*term2);
+    }
+    
+    Standard_Real cost = sqrt ( (1+cos2t)/2. );
+    Standard_Real sint = sqrt ( (1-cos2t)/2. );
+
+    Standard_Real aprim = (a+c+auxil)/2.;
+    Standard_Real cprim = (a+c-auxil)/2.;
+
+    if (Abs(aprim) < gp::Resolution() || Abs(cprim) < gp::Resolution())
+      return Standard_False;
+    
+    term1 = -gdet/(aprim*pdet);
+    term2 = -gdet/(cprim*pdet);
+      
+    if (IsEllip) {
+      Xax = cost;
+      Yax = sint;
+      Rmin = sqrt ( term1);
+      Rmax = sqrt ( term2);
+      if(Rmax<Rmin){
+        Rmax = sqrt ( term1);
+        Rmin = sqrt ( term2);
+      }
+    }
+    else if (term1 <= eps){
+      Xax  = -sint;
+      Yax  =  cost;
+      Rmin = sqrt (-term1);
+      Rmax = sqrt (term2);
+    } 
+    else {
+      Xax  =  cost;
+      Yax  =  sint;
+      Rmin = sqrt (-term2);
+      Rmax = sqrt (term1);
+    }
+  }
+  Center.SetCoord (Xcen,Ycen,0.);
+  MainAxis.SetCoord (Xax,Yax,0.);
+  return Standard_True;
+}
+
+//=======================================================================
+//function : ComputeEllipse
+//purpose  : 
+//=======================================================================
+Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeEllipse(const Handle(Geom_Curve)& c3d,
+                                                       const Standard_Real tol,
+                                                       const Standard_Real c1, const Standard_Real c2,
+                                                       Standard_Real& cf, Standard_Real& cl,
+                                                       Standard_Real& Deviation)
+{
+  if (c3d->IsKind (STANDARD_TYPE(Geom_Ellipse))) {
+    cf = c1;
+    cl = c2;
+    Deviation = 0.;
+    Handle(Geom_Ellipse) anElips = Handle(Geom_Ellipse)::DownCast(c3d);
+    return anElips;
+  }
+
+  Handle(Geom_Curve) res;
+  Standard_Real prec = Precision::PConfusion();
+  
+  Standard_Real AF,BF,CF,DF,EF,Q1,Q2,Q3,c2n;
+  Standard_Integer i;
+  
+  gp_Pnt PStart = c3d->Value(c1);
+  gp_Pnt PEnd = c3d->Value(c2);
+
+  const Standard_Boolean IsClos = PStart.Distance(PEnd) < prec;
+  if (IsClos)
+  {
+    c2n=c2-(c2-c1)/5;
+  }
+  else
+    c2n=c2;
+  //
+  gp_XYZ aBC;
+  Handle(TColgp_HArray1OfPnt) AP = new TColgp_HArray1OfPnt(1,5);
+  AP->SetValue(1,PStart);
+  aBC += PStart.XYZ();
+  Standard_Real dc=(c2n-c1)/4;
+  for (i = 1; i < 5; i++)
+  {
+    gp_Pnt aP = c3d->Value(c1 + dc*i);
+    AP->SetValue(i + 1, aP);
+    aBC += aP.XYZ();
+  }
+  aBC /= 5;
+  aBC *= -1;
+  gp_Vec aTrans(aBC);
+  for (i = 1; i <= 5; ++i)
+  {
+    AP->ChangeValue(i).Translate(aTrans);
+  }
+  gp_Dir ndir;
+  if(!IsArrayPntPlanar(AP,ndir,prec))
+    return res;
+
+  if (Abs(ndir.X()) < gp::Resolution() && Abs(ndir.Y()) < gp::Resolution()
+    && Abs(ndir.Z()) < gp::Resolution())
+    return res;
+
+  gp_Ax3 AX(gp_Pnt(0,0,0),ndir);
+  gp_Trsf Tr;
+  Tr.SetTransformation(AX);
+  gp_Trsf Tr2 = Tr.Inverted();
+  
+  math_Matrix Dt(1, 5, 1, 5);
+  math_Vector F(1, 5), Sl(1, 5);
+
+  Standard_Real XN,YN,ZN = 0.;
+  gp_Pnt PT,PP;
+  for(i=1; i<=5; i++) {
+    PT = AP->Value(i).Transformed(Tr);
+    PT.Coord(XN,YN,ZN);
+    Dt(i, 1) = XN*XN;
+    Dt(i, 2) = XN*YN;
+    Dt(i, 3) = YN*YN;
+    Dt(i, 4) = XN;
+    Dt(i, 5) = YN;
+    F(i) = -1.;
+  }
+
+  math_Gauss aSolver(Dt);
+  if (!aSolver.IsDone())
+    return res;
+
+  aSolver.Solve(F, Sl);
+  
+  AF=Sl(1);
+  BF=Sl(2);
+  CF=Sl(3);
+  DF=Sl(4);
+  EF=Sl(5);
+
+  Q1=AF*CF+BF*EF*DF/4-CF*DF*DF/4-BF*BF/4-AF*EF*EF/4;
+  Q2=AF*CF-BF*BF/4;
+  Q3=AF+CF;
+
+  Standard_Real Rmax, Rmin;
+  gp_Pnt Center;
+  gp_Dir MainAxis;
+  Standard_Boolean IsParab = Standard_False, IsEllip = Standard_False;
+
+if (Q2 > 0 && Q1*Q3 < 0) {
+  // ellipse
+  IsEllip = Standard_True;
+  if (ConicDefinition(AF, BF, CF, DF, EF, 1., IsParab, IsEllip,
+    Center, MainAxis, Rmin, Rmax)) {
+    // create ellipse
+    if (Rmax - Rmin < Precision::Confusion())
+    {
+      return res; //really it is circle, which must be recognized in other method
+    }
+    aTrans *= -1;
+    Center.SetZ(ZN);
+    gp_Pnt NewCenter = Center.Transformed(Tr2);
+    gp_Pnt Ptmp(Center.X() + MainAxis.X() * 10,
+      Center.Y() + MainAxis.Y() * 10,
+      Center.Z() + MainAxis.Z() * 10);
+    gp_Pnt NewPtmp = Ptmp.Transformed(Tr2);
+    gp_Dir NewMainAxis(NewPtmp.X() - NewCenter.X(),
+      NewPtmp.Y() - NewCenter.Y(),
+      NewPtmp.Z() - NewCenter.Z());
+    gp_Ax2 ax2(NewCenter, ndir, NewMainAxis);
+
+    gp_Elips anEllipse(ax2, Rmax, Rmin);
+    anEllipse.Translate(aTrans);
+    Handle(Geom_Ellipse) gell = new Geom_Ellipse(anEllipse);
+
+    // test for 20 points
+    Standard_Real param2 = 0;
+    dc = (c2 - c1) / 20;
+    for (i = 1; i <= 20; i++) {
+      PP = c3d->Value(c1 + i*dc);
+      Standard_Real aPar = ElCLib::Parameter(anEllipse, PP);
+      Standard_Real dist = gell->Value(aPar).Distance(PP);
+      if (dist > tol) return res;  // not done
+      if (dist > param2)
+        param2 = dist;
+    }
+
+
+    Deviation = param2;
+
+    Standard_Real PI2 = 2 * M_PI;
+    cf = ElCLib::Parameter(anEllipse, c3d->Value(c1));
+    cf = ElCLib::InPeriod(cf, 0., PI2);
+
+    //first parameter should be closed to zero
+
+    if (Abs(cf) < Precision::PConfusion() || Abs(PI2 - cf) < Precision::PConfusion())
+      cf = 0.;
+
+    Standard_Real cm = ElCLib::Parameter(anEllipse, c3d->Value((c1 + c2) / 2.));
+    cm = ElCLib::InPeriod(cm, cf, cf + PI2);
+
+    cl = ElCLib::Parameter(anEllipse, c3d->Value(c2));
+    cl = ElCLib::InPeriod(cl, cm, cm + PI2);
+
+    res = gell;
+  }
+}
+/*
+if (Q2 < 0 && Q1 != 0) {
+  // hyberbola
+}
+
+if (Q2 == 0 && Q1 != 0) {
+  // parabola
+}
+*/
+return res;
+}
+
+
+//=======================================================================
+//function : ComputeCurve
+//purpose  : 
+//=======================================================================
+
+Handle(Geom_Curve) GeomConvert_CurveToAnaCurve::ComputeCurve(const Handle(Geom_Curve)& theC3d,
+  const Standard_Real tolerance,
+  const Standard_Real c1, const Standard_Real c2,
+  Standard_Real& cf, Standard_Real& cl,
+  Standard_Real& theGap,
+  const GeomConvert_ConvType theConvType, const GeomAbs_CurveType theTarget)
+{
+  cf = c1;  cl = c2;
+  Handle(Geom_Curve) c3d, newc3d[3];
+  Standard_Integer i, imin = -1;
+  c3d = theC3d;
+  if (c3d.IsNull()) return newc3d[imin];
+  gp_Pnt P1 = c3d->Value(c1);
+  gp_Pnt P2 = c3d->Value(c2);
+  gp_Pnt P3 = c3d->Value(c1 + (c2 - c1) / 2);
+  Standard_Real d[3] = { RealLast(), RealLast(), RealLast() };
+  Standard_Real fp[3], lp[3];
+
+  if (c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
+    Handle(Geom_TrimmedCurve) aTc = Handle(Geom_TrimmedCurve)::DownCast(c3d);
+    c3d = aTc->BasisCurve();
+  }
+
+  if (theConvType == GeomConvert_Target)
+  {
+    theGap = RealLast();
+    if (theTarget == GeomAbs_Line)
+    {
+      newc3d[0] = ComputeLine(c3d, tolerance, c1, c2, fp[0], lp[0], theGap);
+      cf = fp[0];
+      cl = lp[0];
+      return  newc3d[0];
+    }
+    if (theTarget == GeomAbs_Circle)
+    {
+      newc3d[1] = ComputeCircle(c3d, tolerance, c1, c2, fp[1], lp[1], theGap);
+      cf = fp[1];
+      cl = lp[1];
+      return  newc3d[1];
+    }
+    if (theTarget == GeomAbs_Ellipse)
+    {
+      newc3d[2] = ComputeEllipse(c3d, tolerance, c1, c2, fp[2], lp[2], theGap);
+      cf = fp[2];
+      cl = lp[2];
+      return  newc3d[2];
+    }
+  }
+  //
+  if (theConvType == GeomConvert_Simplest)
+  {
+    theGap = RealLast();
+    newc3d[0] = ComputeLine(c3d, tolerance, c1, c2, fp[0], lp[0], theGap);
+    if (!newc3d[0].IsNull())
+    {
+      cf = fp[0];
+      cl = lp[0];
+      return  newc3d[0];
+    }
+    theGap = RealLast();
+    newc3d[1] = ComputeCircle(c3d, tolerance, c1, c2, fp[1], lp[1], theGap);
+    if (!newc3d[1].IsNull())
+    {
+      cf = fp[1];
+      cl = lp[1];
+      return  newc3d[1];
+    }
+    theGap = RealLast();
+    newc3d[2] = ComputeEllipse(c3d, tolerance, c1, c2, fp[2], lp[2], theGap);
+    if (!newc3d[2].IsNull())
+    {
+      cf = fp[2];
+      cl = lp[2];
+      return  newc3d[2];
+    }
+    // Conversion failed, returns null curve
+    return newc3d[0];
+  }
+
+  //  theConvType == GeomConvert_MinGap
+  // recognition in case of small curve
+  imin = -1;
+  if((P1.Distance(P2) < 2*tolerance) && (P1.Distance(P3) < 2*tolerance)) {
+    newc3d[1] = ComputeCircle(c3d, tolerance, c1, c2, fp[1], lp[1], d[1]);  
+    newc3d[0] = ComputeLine(c3d, tolerance, c1, c2, fp[0], lp[0], d[0]);
+    imin = 1;
+    if (newc3d[1].IsNull() || d[0] < d[1])
+    {
+      imin = 0;
+    }
+  }
+  else {
+    d[0] = RealLast();
+    newc3d[0] = ComputeLine (c3d,tolerance,c1,c2,fp[0],lp[0],d[0]);  
+    Standard_Real tol = Min(tolerance, d[0]);
+    if (!Precision::IsInfinite(c1) && !Precision::IsInfinite(c2))
+    {
+      d[1] = RealLast();
+      newc3d[1] = ComputeCircle(c3d, tol, c1, c2, fp[1], lp[1], d[1]);
+      tol = Min(tol, d[1]);
+      d[2] = RealLast();
+      newc3d[2] = ComputeEllipse(c3d, tol, c1, c2, fp[2], lp[2], d[2]);
+    }
+    Standard_Real dd = RealLast();
+    for (i = 0; i < 3; ++i)
+    {
+      if (newc3d[i].IsNull()) continue;
+      if (d[i] < dd)
+      {
+        dd = d[i];
+        imin = i;
+      }
+    }
+  }
+
+  if (imin >= 0)
+  {
+    cf = fp[imin];
+    cl = lp[imin];
+    theGap = d[imin];
+    return newc3d[imin];
+  }
+  else
+  {
+    cf = c1;
+    cl = c2;
+    theGap = -1.;
+    return newc3d[0]; // must be null curve; 
+  }
+}

src/GeomConvert/GeomConvert_CurveToAnaCurve.hxx

@@ -0,0 +1,133 @@
+// Created: 2001-05-21 
+// 
+// Copyright (c) 2001-2013 OPEN CASCADE SAS 
+// 
+// This file is part of commercial software by OPEN CASCADE SAS, 
+// furnished in accordance with the terms and conditions of the contract 
+// and with the inclusion of this copyright notice. 
+// This file or any part thereof may not be provided or otherwise 
+// made available to any third party. 
+// 
+// No ownership title to the software is transferred hereby. 
+// 
+// OPEN CASCADE SAS makes no representation or warranties with respect to the 
+// performance of this software, and specifically disclaims any responsibility 
+// for any damages, special or consequential, connected with its use. 
+
+#ifndef _GeomConvert_CurveToAnaCurve_HeaderFile
+#define _GeomConvert_CurveToAnaCurve_HeaderFile
+
+#include <Standard.hxx>
+#include <Standard_DefineAlloc.hxx>
+#include <Standard_Handle.hxx>
+#include <Standard_Boolean.hxx>
+#include <Standard_Real.hxx>
+#include <TColgp_Array1OfPnt.hxx>
+#include <GeomConvert_ConvType.hxx>
+#include <GeomAbs_CurveType.hxx>
+
+class Geom_Curve;
+class Geom_Line;
+class gp_Lin;
+class gp_Pnt;
+class gp_Circ;
+
+
+
+class GeomConvert_CurveToAnaCurve 
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+ 
+  Standard_EXPORT GeomConvert_CurveToAnaCurve();
+  
+  Standard_EXPORT GeomConvert_CurveToAnaCurve(const Handle(Geom_Curve)& C);
+  
+  Standard_EXPORT void Init (const Handle(Geom_Curve)& C);
+  
+  //! Converts me to analytical if possible with given
+  //! tolerance. The new first and last parameters are
+  //! returned to newF, newL
+  Standard_EXPORT Standard_Boolean ConvertToAnalytical (const Standard_Real theTol, Handle(Geom_Curve)& theResultCurve, const Standard_Real F, const Standard_Real L, Standard_Real& newF, Standard_Real& newL);
+  
+  Standard_EXPORT static Handle(Geom_Curve) ComputeCurve (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, 
+    const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, 
+    Standard_Real& theGap, const GeomConvert_ConvType theCurvType = GeomConvert_MinGap, const GeomAbs_CurveType theTarget = GeomAbs_Line);
+  
+  //! Tries to convert the given curve to circle with given
+  //! tolerance. Returns NULL curve if conversion is
+  //! not possible.
+  Standard_EXPORT static Handle(Geom_Curve) ComputeCircle (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, Standard_Real& Deviation);
+  
+  //! Tries to convert the given curve to ellipse with given
+  //! tolerance. Returns NULL curve if conversion is
+  //! not possible.
+  Standard_EXPORT static Handle(Geom_Curve) ComputeEllipse (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, Standard_Real& Deviation);
+  
+  //! Tries to convert the given curve to line with given
+  //! tolerance. Returns NULL curve if conversion is
+  //! not possible.
+  Standard_EXPORT static Handle(Geom_Line) ComputeLine (const Handle(Geom_Curve)& curve, const Standard_Real tolerance, const Standard_Real c1, const Standard_Real c2, Standard_Real& cf, Standard_Real& cl, Standard_Real& Deviation);
+  
+  //! Returns true if the set of points is linear with given
+  //! tolerance
+  Standard_EXPORT static Standard_Boolean IsLinear (const TColgp_Array1OfPnt& aPoints, const Standard_Real tolerance, Standard_Real& Deviation);
+   
+  //! Creates line on two points.
+  //! Resulting parameters returned
+  Standard_EXPORT static gp_Lin GetLine(const gp_Pnt& P1, const gp_Pnt& P2, Standard_Real& cf, Standard_Real& cl);
+
+  //! Creates circle on points. Returns true if OK.
+  Standard_EXPORT static Standard_Boolean GetCircle(gp_Circ& Circ, const gp_Pnt& P0, const gp_Pnt& P1, const gp_Pnt& P2); 
+
+  //! Returns maximal deviation of converted surface from the original
+  //! one computed by last call to ConvertToAnalytical
+  Standard_Real Gap() const
+  {
+    return myGap;
+  }
+
+  //! Returns conversion type 
+  GeomConvert_ConvType GetConvType() const
+  {
+    return myConvType;
+  }
+
+  //! Sets type of convertion
+  void SetConvType(const GeomConvert_ConvType theConvType)
+  {
+    myConvType = theConvType;
+  }
+
+  //! Returns target curve type 
+  GeomAbs_CurveType GetTarget() const
+  {
+    return myTarget;
+  }
+
+  //! Sets target curve type
+  void SetTarget(const GeomAbs_CurveType theTarget)
+  {
+    myTarget = theTarget;
+  }
+
+
+protected:
+
+
+
+
+
+private:
+
+  Handle(Geom_Curve) myCurve;
+  Standard_Real myGap;
+  GeomConvert_ConvType myConvType;
+  GeomAbs_CurveType myTarget;
+
+};
+
+
+
+#endif // _GeomConvert_CurveToAnaCurve_HeaderFile

src/GeomConvert/GeomConvert_FuncConeLSDist.cxx

@@ -0,0 +1,68 @@
+// Copyright (c) 1995-1999 Matra Datavision
+// Copyright (c) 1999-2022 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.
+
+
+#include <GeomConvert_FuncConeLSDist.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <gp_Ax3.hxx>
+#include <math_Vector.hxx>
+#include <ElSLib.hxx>
+
+//=======================================================================
+//function : GeomConvert_FuncConeLSDist
+//purpose  : 
+//=======================================================================
+GeomConvert_FuncConeLSDist::GeomConvert_FuncConeLSDist(
+                                  const Handle(TColgp_HArray1OfXYZ)& thePoints,
+                                  const gp_Dir& theDir):
+  myPoints(thePoints), myDir(theDir)
+{
+}
+
+//=======================================================================
+//function : NbVariables
+//purpose  : 
+//=======================================================================
+Standard_Integer GeomConvert_FuncConeLSDist::NbVariables () const
+{
+  return 5;
+}
+
+//=======================================================================
+//function : Value
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncConeLSDist::Value(const math_Vector& X, Standard_Real& F)
+{
+  gp_Pnt aLoc(X(1), X(2), X(3));
+  Standard_Real aSemiAngle = X(4), anR = X(5);
+  gp_Ax3 aPos(aLoc, myDir);
+
+  F = 0.;
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    Standard_Real u, v;
+    gp_Pnt aPi(myPoints->Value(i));
+    ElSLib::ConeParameters(aPos, anR, aSemiAngle, aPi, u, v);
+    gp_Pnt aPp;
+    ElSLib::ConeD0(u, v, aPos, anR, aSemiAngle, aPp);
+    F += aPi.SquareDistance(aPp);
+  }
+
+  return Standard_True;
+}
+
+

src/GeomConvert/GeomConvert_FuncConeLSDist.hxx

@@ -0,0 +1,66 @@
+// Copyright (c) 1991-1999 Matra Datavision
+// Copyright (c) 1999-2022 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.
+
+#ifndef _GeomConvert_FuncConeLSDist_HeaderFile
+#define _GeomConvert_FuncConeLSDist_HeaderFile
+
+#include <Standard.hxx>
+#include <Standard_DefineAlloc.hxx>
+
+#include <math_MultipleVarFunction.hxx>
+#include <TColgp_HArray1OfXYZ.hxx>
+#include <math_Vector.hxx>
+#include <gp_Dir.hxx>
+
+//! Function for search of Cone canonic parameters: coordinates of center local coordinate system, 
+//! direction of axis, radius and semi-angle from set of points
+//! by least square method.
+//! 
+//!
+class GeomConvert_FuncConeLSDist : public math_MultipleVarFunction
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  //! Constructor.
+  Standard_EXPORT GeomConvert_FuncConeLSDist() {};
+  
+  Standard_EXPORT GeomConvert_FuncConeLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints,
+                                                   const gp_Dir& theDir);
+
+  void SetPoints(const Handle(TColgp_HArray1OfXYZ)& thePoints)
+  {
+    myPoints = thePoints;
+  }
+
+  void SetDir(const gp_Dir& theDir)
+  {
+    myDir = theDir;
+  }
+
+  //! Number of variables.
+  Standard_EXPORT Standard_Integer NbVariables() const Standard_OVERRIDE;
+
+  //! Value.
+  Standard_EXPORT Standard_Boolean Value(const math_Vector& X,Standard_Real& F) Standard_OVERRIDE;
+
+
+private:
+
+  Handle(TColgp_HArray1OfXYZ) myPoints;
+  gp_Dir myDir;
+  
+};
+#endif // _GeomConvert_FuncConeLSDist_HeaderFile

src/GeomConvert/GeomConvert_FuncCylinderLSDist.cxx

@@ -0,0 +1,140 @@
+// Copyright (c) 1995-1999 Matra Datavision
+// Copyright (c) 1999-2022 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.
+
+
+#include <GeomConvert_FuncCylinderLSDist.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <math_Vector.hxx>
+
+//=======================================================================
+//function : GeomConvert_FuncCylinderLSDist
+//purpose  : 
+//=======================================================================
+GeomConvert_FuncCylinderLSDist::GeomConvert_FuncCylinderLSDist(
+                                  const Handle(TColgp_HArray1OfXYZ)& thePoints,
+                                  const gp_Dir& theDir):
+  myPoints(thePoints), myDir(theDir)
+{
+}
+
+//=======================================================================
+//function : NbVariables
+//purpose  : 
+//=======================================================================
+Standard_Integer GeomConvert_FuncCylinderLSDist::NbVariables () const
+{
+  return 4;
+}
+
+//=======================================================================
+//function : Value
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncCylinderLSDist::Value(const math_Vector& X,Standard_Real& F)
+{
+  gp_XYZ aLoc(X(1), X(2), X(3));
+  Standard_Real anR2 = X(4)*X(4);
+
+  F = 0.;
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    gp_Vec aV(myPoints->Value(i) - aLoc);
+    Standard_Real aD2 = aV.CrossSquareMagnitude(myDir);
+    Standard_Real d = aD2 - anR2;
+    F += d * d;
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : Gradient
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncCylinderLSDist::Gradient(const math_Vector& X,math_Vector& G)
+
+{
+  gp_XYZ aLoc(X(1), X(2), X(3));
+  Standard_Real anR = X(4), anR2 = anR * anR;
+  Standard_Real x = myDir.X(), y = myDir.Y(), z = myDir.Z();
+  G.Init(0.);
+
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    gp_Vec aV(myPoints->Value(i) - aLoc);
+    Standard_Real aD2 = aV.CrossSquareMagnitude(myDir);
+    Standard_Real d = aD2 - anR2;
+    Standard_Real Dx0 =  2.*(aV.Z()*x - aV.X()*z)*z
+                        -2.*(aV.X()*y - aV.Y()*x)*y;
+    Standard_Real Dy0 = -2.*(aV.Y()*z - aV.Z()*y)*z
+                        +2.*(aV.X()*y - aV.Y()*x)*x;
+    Standard_Real Dz0 =  2.*(aV.Y()*z - aV.Z()*y)*y
+                        -2.*(aV.Z()*x - aV.X()*z)*x;
+
+    G(1) += d * Dx0;
+    G(2) += d * Dy0;
+    G(3) += d * Dz0;
+    //
+    G(4) += d;
+  }
+
+  G *= 2;
+  G(6) *= -2.*anR;
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : Values
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncCylinderLSDist::Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
+{
+  gp_XYZ aLoc(X(1), X(2), X(3));
+  Standard_Real anR = X(4), anR2 = anR * anR;
+  Standard_Real x = myDir.X(), y = myDir.Y(), z = myDir.Z();
+
+  F = 0.;
+  G.Init(0.);
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    gp_Vec aV(myPoints->Value(i) - aLoc);
+    Standard_Real aD2 = aV.CrossSquareMagnitude(myDir);
+    Standard_Real d = aD2 - anR2;
+    Standard_Real Dx0 = 2.*(aV.Z()*x - aV.X()*z)*z
+      - 2.*(aV.X()*y - aV.Y()*x)*y;
+    Standard_Real Dy0 = -2.*(aV.Y()*z - aV.Z()*y)*z
+      + 2.*(aV.X()*y - aV.Y()*x)*x;
+    Standard_Real Dz0 = 2.*(aV.Y()*z - aV.Z()*y)*y
+      - 2.*(aV.Z()*x - aV.X()*z)*x;
+
+    G(1) += d * Dx0;
+    G(2) += d * Dy0;
+    G(3) += d * Dz0;
+    //
+    G(4) += d;
+    //
+    F += d * d;
+  }
+
+  G *= 2;
+  G(4) *= -2.*anR;
+
+  return true;
+}
+

src/GeomConvert/GeomConvert_FuncCylinderLSDist.hxx

@@ -0,0 +1,100 @@
+// Copyright (c) 1991-1999 Matra Datavision
+// Copyright (c) 1999-2022 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.
+
+#ifndef _GeomConvert_FuncCylinderLSDist_HeaderFile
+#define _GeomConvert_FuncCylinderLSDist_HeaderFile
+
+#include <Standard.hxx>
+#include <Standard_DefineAlloc.hxx>
+
+#include <math_MultipleVarFunctionWithGradient.hxx>
+#include <TColgp_HArray1OfXYZ.hxx>
+#include <math_Vector.hxx>
+#include <gp_Dir.hxx>
+
+//! Function for search of cylinder canonic parameters: coordinates of center local coordinate system, 
+//! direction of axis and radius from set of points
+//! by least square method.
+//! 
+//! The class inherits math_MultipleVarFunctionWithGradient and thus is intended
+//! for use in math_BFGS algorithm.
+//!
+//! Parametrisation:
+//! Cylinder is defined by its axis and radius. Axis is defined by 3 cartesian coordinats it location x0, y0, z0 
+//! and direction, which is constant and set by user: 
+//! dir.x, dir.y, dir.z
+//! The criteria is:
+//! F(x0, y0, z0, theta, phi, R) = Sum[|(P(i) - Loc)^dir|^2 - R^2]^2 => min
+//! P(i) is i-th sample point, Loc, dir - axis location and direction, R - radius
+//!
+//! The square vector product |(P(i) - Loc)^dir|^2 is:
+//!
+//! [(y - y0)*dir.z - (z - z0)*dir.y]^2 +
+//! [(z - z0)*dir.x - (x - x0)*dir.z]^2 +
+//! [(x - x0)*dir.y - (y - y0)*dir.x]^2
+//!
+//! First derivative of square vector product are:
+//! Dx0 =  2*[(z - z0)*dir.x - (x - x0)*dir.z]*dir.z 
+//!       -2*[(x - x0)*dir.y - (y - y0)*dir.x]*dir.y
+//! Dy0 = -2*[(y - y0)*dir.z - (z - z0)*dir.y]*dir.z
+//!       +2*[(x - x0)*dir.y - (y - y0)*dir.x]*dir.x
+//! Dz0 =  2*[(y - y0)*dir.z - (z - z0)*dir.y]*dir.y
+//!       -2*[(z - z0)*dir.x - (x - x0)*dir.z]*dir.x
+//!
+//! dF/dx0 : G1(...) = 2*Sum{[...]*Dx0}
+//! dF/dy0 : G2(...) = 2*Sum{[...]*Dy0}
+//! dF/dz0 : G3(...) = 2*Sum{[...]*Dz0}
+//! dF/dR : G4(...) = -4*R*Sum[...]
+//! [...] = [|(P(i) - Loc)^dir|^2 - R^2]
+class GeomConvert_FuncCylinderLSDist : public math_MultipleVarFunctionWithGradient
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  //! Constructor.
+  Standard_EXPORT GeomConvert_FuncCylinderLSDist() {};
+  
+  Standard_EXPORT GeomConvert_FuncCylinderLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints,
+                                                   const gp_Dir& theDir);
+
+  void SetPoints(const Handle(TColgp_HArray1OfXYZ)& thePoints)
+  {
+    myPoints = thePoints;
+  }
+
+  void SetDir(const gp_Dir& theDir)
+  {
+    myDir = theDir;
+  }
+
+  //! Number of variables.
+  Standard_EXPORT Standard_Integer NbVariables() const Standard_OVERRIDE;
+
+  //! Value.
+  Standard_EXPORT Standard_Boolean Value(const math_Vector& X,Standard_Real& F) Standard_OVERRIDE;
+
+  //! Gradient.
+  Standard_EXPORT Standard_Boolean Gradient(const math_Vector& X,math_Vector& G) Standard_OVERRIDE;
+
+  //! Value and gradient.
+  Standard_EXPORT Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G) Standard_OVERRIDE;
+
+private:
+
+  Handle(TColgp_HArray1OfXYZ) myPoints;
+  gp_Dir myDir;
+  
+};
+#endif // _GeomConvert_FuncCylinderLSDist_HeaderFile

src/GeomConvert/GeomConvert_FuncSphereLSDist.cxx

@@ -0,0 +1,115 @@
+// Created on: 2016-05-10
+// Created by: Alexander MALYSHEV
+// Copyright (c) 1995-1999 Matra Datavision
+// Copyright (c) 1999-2016 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.
+
+
+#include <GeomConvert_FuncSphereLSDist.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <math_Vector.hxx>
+
+//=======================================================================
+//function : GeomConvert_FuncSphereLSDist
+//purpose  : 
+//=======================================================================
+GeomConvert_FuncSphereLSDist::GeomConvert_FuncSphereLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints):
+  myPoints(thePoints)
+{
+}
+
+//=======================================================================
+//function : NbVariables
+//purpose  : 
+//=======================================================================
+Standard_Integer GeomConvert_FuncSphereLSDist::NbVariables () const
+{
+  return 4;
+}
+
+//=======================================================================
+//function : Value
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncSphereLSDist::Value(const math_Vector& X,Standard_Real& F)
+{
+  gp_XYZ aLoc(X(1), X(2), X(3));
+  Standard_Real anR2 = X(4)*X(4);
+
+  F = 0.;
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    Standard_Real d = (myPoints->Value(i) - aLoc).SquareModulus() - anR2;
+    F += d * d;
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : Gradient
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncSphereLSDist::Gradient(const math_Vector& X,math_Vector& G)
+
+{
+  gp_XYZ aLoc(X(1), X(2), X(3));
+  Standard_Real anR = X(4), anR2 = anR * anR;
+
+  G.Init(0.);
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    gp_XYZ dLoc = myPoints->Value(i) - aLoc;
+    Standard_Real d = dLoc.SquareModulus() - anR2;
+    G(1) += d * dLoc.X();
+    G(2) += d * dLoc.Y();
+    G(3) += d * dLoc.Z();
+    G(4) += d;
+  }
+  G *= -4;
+  G(4) *= anR;
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : Values
+//purpose  : 
+//=======================================================================
+Standard_Boolean GeomConvert_FuncSphereLSDist::Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
+{
+  gp_XYZ aLoc(X(1), X(2), X(3));
+  Standard_Real anR = X(4), anR2 = anR * anR;
+
+  G.Init(0.);
+  F = 0.;
+  Standard_Integer i;
+  for (i = myPoints->Lower(); i <= myPoints->Upper(); ++i)
+  {
+    gp_XYZ dLoc = myPoints->Value(i) - aLoc;
+    Standard_Real d = dLoc.SquareModulus() - anR2;
+    G(1) += d * dLoc.X();
+    G(2) += d * dLoc.Y();
+    G(3) += d * dLoc.Z();
+    G(4) += d;
+    F += d * d;
+  }
+  G *= -4;
+  G(4) *= anR;
+
+  return true;
+}
+

src/GeomConvert/GeomConvert_FuncSphereLSDist.hxx

@@ -0,0 +1,77 @@
+
+// Copyright (c) 1991-1999 Matra Datavision
+// Copyright (c) 1999-2022 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.
+
+#ifndef _GeomConvert_FuncSphereLSDist_HeaderFile
+#define _GeomConvert_FuncSphereLSDist_HeaderFile
+
+#include <Standard.hxx>
+#include <Standard_DefineAlloc.hxx>
+
+#include <math_MultipleVarFunctionWithGradient.hxx>
+#include <TColgp_HArray1OfXYZ.hxx>
+#include <math_Vector.hxx>
+
+//! Function for search of sphere canonic parameters: coordinates of center and radius from set of moints
+//! by least square method.
+//! //!
+//! The class inherits math_MultipleVarFunctionWithGradient and thus is intended
+//! for use in math_BFGS algorithm.
+//!
+//! The criteria is:
+//! F(x0, y0, z0, R) = Sum[(x(i) - x0)^2 + (y(i) - y0)^2 + (z(i) - z0)^2 - R^2]^2 => min,
+//! x(i), y(i), z(i) - coordinates of sample points, x0, y0, z0, R - coordinates of center and radius of sphere,
+//! which must be defined
+//!
+//! The first derivative are:
+//! dF/dx0 : G1(x0, y0, z0, R) = -4*Sum{[...]*(x(i) - x0)} 
+//! dF/dy0 : G2(x0, y0, z0, R) = -4*Sum{[...]*(y(i) - y0)}
+//! dF/dz0 : G3(x0, y0, z0, R) = -4*Sum{[...]*(z(i) - z0)}
+//! dF/dR : G4(x0, y0, z0, R) = -4*R*Sum[...]
+//! [...] = [(x(i) - x0)^2 + (y(i) - y0)^2 + (z(i) - z0)^2 - R^2]
+//!
+class GeomConvert_FuncSphereLSDist : public math_MultipleVarFunctionWithGradient
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  //! Constructor.
+  Standard_EXPORT GeomConvert_FuncSphereLSDist() {};
+  
+  Standard_EXPORT GeomConvert_FuncSphereLSDist(const Handle(TColgp_HArray1OfXYZ)& thePoints);
+
+  void SetPoints(const Handle(TColgp_HArray1OfXYZ)& thePoints)
+  {
+    myPoints = thePoints;
+  }
+
+  //! Number of variables.
+  Standard_EXPORT Standard_Integer NbVariables() const Standard_OVERRIDE;
+
+  //! Value.
+  Standard_EXPORT Standard_Boolean Value(const math_Vector& X,Standard_Real& F) Standard_OVERRIDE;
+
+  //! Gradient.
+  Standard_EXPORT Standard_Boolean Gradient(const math_Vector& X,math_Vector& G) Standard_OVERRIDE;
+
+  //! Value and gradient.
+  Standard_EXPORT Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G) Standard_OVERRIDE;
+
+private:
+
+  Handle(TColgp_HArray1OfXYZ) myPoints;
+  
+};
+#endif // _GeomConvert_FuncSphereLSDist_HeaderFile

src/GeomConvert/GeomConvert_SurfToAnaSurf.hxx

@@ -0,0 +1,135 @@
+// Created: 1998-06-03 
+// 
+// Copyright (c) 1999-2013 OPEN CASCADE SAS 
+// 
+// This file is part of commercial software by OPEN CASCADE SAS, 
+// furnished in accordance with the terms and conditions of the contract 
+// and with the inclusion of this copyright notice. 
+// This file or any part thereof may not be provided or otherwise 
+// made available to any third party. 
+// 
+// No ownership title to the software is transferred hereby. 
+// 
+// OPEN CASCADE SAS makes no representation or warranties with respect to the 
+// performance of this software, and specifically disclaims any responsibility 
+// for any damages, special or consequential, connected with its use. 
+
+#ifndef _GeomConvert_SurfToAnaSurf_HeaderFile
+#define _GeomConvert_SurfToAnaSurf_HeaderFile
+
+#include <Standard.hxx>
+#include <Standard_DefineAlloc.hxx>
+#include <Standard_Handle.hxx>
+
+#include <Standard_Real.hxx>
+#include <Standard_Boolean.hxx>
+#include <GeomConvert_ConvType.hxx>
+#include <GeomAbs_SurfaceType.hxx>
+#include <TColgp_HArray1OfXYZ.hxx>
+class Geom_Surface;
+class Geom_SurfaceOfRevolution;
+class Geom_Circle;
+
+//! Converts a surface to the analitical form with given
+//! precision. Conversion is done only the surface is bspline
+//! of bezier and this can be approximed by some analytical
+//! surface with that precision.
+class GeomConvert_SurfToAnaSurf 
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+  
+  Standard_EXPORT GeomConvert_SurfToAnaSurf();
+  
+  Standard_EXPORT GeomConvert_SurfToAnaSurf(const Handle(Geom_Surface)& S);
+  
+  Standard_EXPORT void Init (const Handle(Geom_Surface)& S);
+
+  void SetConvType(const GeomConvert_ConvType theConvType = GeomConvert_Simplest)
+  {
+    myConvType = theConvType;
+  }
+  void SetTarget(const GeomAbs_SurfaceType theSurfType = GeomAbs_Plane)
+  {
+    myTarget = theSurfType;
+  }
+
+  //! Returns maximal deviation of converted surface from the original
+  //! one computed by last call to ConvertToAnalytical
+  Standard_Real Gap() const
+  {
+    return myGap;
+  }
+  
+  //! Tries to convert the Surface to an Analytic form
+  //! Returns the result
+  //! In case of failure, returns a Null Handle
+  //!
+  Standard_EXPORT Handle(Geom_Surface) ConvertToAnalytical (const Standard_Real InitialToler);
+  Standard_EXPORT Handle(Geom_Surface) ConvertToAnalytical (const Standard_Real InitialToler,
+                                                            const Standard_Real Umin, const Standard_Real Umax,
+                                                            const Standard_Real Vmin, const Standard_Real Vmax);
+ 
+  //! Returns true if surfaces is same with the given tolerance
+  Standard_EXPORT static Standard_Boolean IsSame (const Handle(Geom_Surface)& S1, const Handle(Geom_Surface)& S2, const Standard_Real tol);
+  
+  //! Returns true, if surface is canonical
+  Standard_EXPORT static Standard_Boolean IsCanonical (const Handle(Geom_Surface)& S);
+
+private:
+  //!static method for checking surface of revolution
+  //!To avoid two-parts cone-like surface
+  static void CheckVTrimForRevSurf(const Handle(Geom_SurfaceOfRevolution)& aRevSurf,
+    Standard_Real& V1, Standard_Real& V2);
+
+  //!static method to try create cylindrical or conical surface
+  static Handle(Geom_Surface) TryCylinerCone(const Handle(Geom_Surface)& theSurf, const Standard_Boolean theVCase,
+    const Handle(Geom_Curve)& theUmidiso, const Handle(Geom_Curve)& theVmidiso,
+    const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, const Standard_Real theV2,
+    const Standard_Real theToler);
+
+  //!static method to try create cylinrical surface using least square method
+  static Standard_Boolean GetCylByLS(const Handle(TColgp_HArray1OfXYZ)& thePoints,
+    const Standard_Real theTol,
+    gp_Ax3& thePos, Standard_Real& theR,
+    Standard_Real& theGap);
+
+  //!static method to try create cylinrical surface based on its Gauss field
+  static Handle(Geom_Surface) TryCylinderByGaussField(const Handle(Geom_Surface)& theSurf,
+    const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, const Standard_Real theV2,
+    const Standard_Real theToler, const Standard_Integer theNbU = 20, const Standard_Integer theNbV = 20,
+    const Standard_Boolean theLeastSquare = Standard_False);
+
+  //! static method to try create toroidal surface.
+  //! In case <isTryUMajor> = Standard_True try to use V isoline radius as minor radaius.
+  static Handle(Geom_Surface) TryTorusSphere(const Handle(Geom_Surface)& theSurf,
+    const Handle(Geom_Circle)& circle,
+    const Handle(Geom_Circle)& otherCircle,
+    const Standard_Real Param1,
+    const Standard_Real Param2,
+    const Standard_Real aParam1ToCrv,
+    const Standard_Real aParam2ToCrv,
+    const Standard_Real toler,
+    const Standard_Boolean isTryUMajor);
+
+  static Standard_Real ComputeGap(const Handle(Geom_Surface)& theSurf,
+    const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, const Standard_Real theV2,
+    const Handle(Geom_Surface) theNewSurf, const Standard_Real theTol = RealLast());
+
+
+
+
+protected:
+
+private:
+
+  Handle(Geom_Surface) mySurf;
+  Standard_Real myGap;
+  GeomConvert_ConvType myConvType;
+  GeomAbs_SurfaceType myTarget;
+
+};
+
+
+#endif // _GeomConvert_SurfToAnaSurf_HeaderFile

src/GeomliteTest/GeomliteTest_SurfaceCommands.cxx

@@ -72,7 +72,9 @@
 #include <Geom2dConvert.hxx>
 #include <Geom2dConvert_BSplineCurveToBezierCurve.hxx>
 #include <GeomLProp_SLProps.hxx>
-
+#include <GeomConvert_SurfToAnaSurf.hxx>
+#include <GeomConvert_CurveToAnaCurve.hxx>
+#include <GeomConvert_ConvType.hxx>
 
 #include <DrawTrSurf_BezierSurface.hxx>
 #include <DrawTrSurf_BSplineSurface.hxx>
@@ -517,6 +519,111 @@ static Standard_Integer converting(Draw_Interpretor& , Standard_Integer n, const
   return 0;
 }
 
+//=======================================================================
+//function : converting to canonical
+//purpose  : 
+//=======================================================================
+
+static Standard_Integer tocanon(Draw_Interpretor& di, Standard_Integer n, const char ** a)
+{
+  if (n < 3) return 1;
+
+  GeomConvert_ConvType aConvType = GeomConvert_Simplest;
+  GeomAbs_CurveType aCurv = GeomAbs_Line;
+  GeomAbs_SurfaceType aSurf = GeomAbs_Plane;
+  if (n > 4)
+  {
+    if (strcmp(a[4], "sim") == 0) {
+      aConvType = GeomConvert_Simplest;
+    }
+    else if (strcmp(a[4], "gap") == 0) {
+      aConvType = GeomConvert_MinGap;
+    }
+    else if (strcmp(a[4], "lin") == 0) {
+      aConvType = GeomConvert_Target;
+      aCurv = GeomAbs_Line;
+    }
+    else if (strcmp(a[4], "cir") == 0) {
+      aConvType = GeomConvert_Target;
+      aCurv = GeomAbs_Circle;
+    }
+    else if (strcmp(a[4], "ell") == 0) {
+      aConvType = GeomConvert_Target;
+      aCurv = GeomAbs_Ellipse;
+    }
+    else if (strcmp(a[4], "pln") == 0) {
+      aConvType = GeomConvert_Target;
+      aSurf = GeomAbs_Plane;
+    }
+    else if (strcmp(a[4], "cyl") == 0) {
+      aConvType = GeomConvert_Target;
+      aSurf = GeomAbs_Cylinder;
+    }
+    else if (strcmp(a[4], "con") == 0) {
+      aConvType = GeomConvert_Target;
+      aSurf = GeomAbs_Cone;
+    }
+    else if (strcmp(a[4], "sph") == 0) {
+      aConvType = GeomConvert_Target;
+      aSurf = GeomAbs_Sphere;
+    }
+    else if (strcmp(a[4], "tor") == 0) {
+      aConvType = GeomConvert_Target;
+      aSurf = GeomAbs_Torus;
+    }
+  }
+
+  Standard_Real tol = Precision::Confusion();
+  if (n > 3)
+  {
+    tol = Draw::Atof(a[3]);
+  }
+
+  Handle(Geom_Curve) GC = DrawTrSurf::GetCurve(a[2]);
+  if (GC.IsNull()) {
+    Handle(Geom_Surface) GS = DrawTrSurf::GetSurface(a[2]);
+    if (GS.IsNull()) {
+      return 1;
+    }
+    else {
+      GeomConvert_SurfToAnaSurf aSurfToAna(GS);
+      aSurfToAna.SetConvType(aConvType);
+      if(aConvType == GeomConvert_Target)
+        aSurfToAna.SetTarget(aSurf);
+      Handle(Geom_Surface) anAnaSurf = aSurfToAna.ConvertToAnalytical(tol);
+      if (!anAnaSurf.IsNull())
+      {
+        DrawTrSurf::Set(a[1], anAnaSurf);
+        Standard_Real aGap = aSurfToAna.Gap();
+        di << "Gap = " << aGap << "\n";
+      }
+      else
+        di << "Conversion failed" << "\n";
+    }
+  }
+  else {
+    GeomConvert_CurveToAnaCurve aCurvToAna(GC);
+    aCurvToAna.SetConvType(aConvType);
+    if (aConvType == GeomConvert_Target)
+      aCurvToAna.SetTarget(aCurv);
+
+    Handle(Geom_Curve) anAnaCurv;
+    Standard_Real tf = GC->FirstParameter(), tl = GC->LastParameter(), ntf, ntl;
+    Standard_Boolean isdone = aCurvToAna.ConvertToAnalytical(tol, anAnaCurv, tf, tl, ntf, ntl);
+    if (isdone)
+    {
+      anAnaCurv = new Geom_TrimmedCurve(anAnaCurv, ntf, ntl);
+      DrawTrSurf::Set(a[1], anAnaCurv);
+      Standard_Real aGap = aCurvToAna.Gap();
+      di << "Gap = " << aGap << "\n";
+    }
+    else
+      di << "Conversion failed" << "\n";
+  }
+
+  return 0;
+}
+
 
 //=======================================================================
 //function : tobezier
@@ -1661,6 +1768,11 @@ void  GeomliteTest::SurfaceCommands(Draw_Interpretor& theCommands)
 		  __FILE__,
 		  converting,g);
 
+  theCommands.Add("tocanon",
+    "tocanon result c3d/surf [tol [sim gap lin cir ell pln cyl con sph tor]]",
+    __FILE__,
+    tocanon, g);
+
   theCommands.Add("tobezier",
 		  "tobezier result c2d/c3d/surf [ufirst, ulast / ufirst, ulast, vfirst, vlast]",
 		  __FILE__,

src/SWDRAW/SWDRAW_ShapeAnalysis.cxx

@@ -31,6 +31,7 @@
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <gp_Pnt.hxx>
 #include <gp_Pnt2d.hxx>
+#include <gp_Elips.hxx>
 #include <GProp_GProps.hxx>
 #include <Precision.hxx>
 #include <ShapeAnalysis.hxx>
@@ -63,6 +64,14 @@
 #include <Adaptor3d_CurveOnSurface.hxx>
 #include <BRepAdaptor_HCurve2d.hxx>
 #include <BRepAdaptor_HSurface.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_CylindricalSurface.hxx>
+#include <Geom_ConicalSurface.hxx>
+#include <Geom_SphericalSurface.hxx>
+#include <Geom_Line.hxx>
+#include <Geom_Circle.hxx>
+#include <Geom_Ellipse.hxx>
+#include <ShapeAnalysis_CanonicalRecognition.hxx>
 
 #include <stdio.h>
 static Standard_Integer tolerance
@@ -978,6 +987,203 @@ static Standard_Integer checkedge(Draw_Interpretor& di, Standard_Integer argc, c
 
   return 0;
 }
+//=======================================================================
+// getanasurf
+//=======================================================================
+static Standard_Integer getanasurf(Draw_Interpretor& di,
+  Standard_Integer n, const char** a)
+
+{
+  if (n < 3) {
+    di << "Usage: \n";
+    di << "getanasurf res shape [target [tol [sample]]] \n";
+    di << "target is reqired type of surface and can be: pln, cyl, con sph \n";
+    di << "sample is surface of required type, which parameters are used as starting \n";
+    di << "point for seaching parametrs of surface by Least Square method when input shape \n";
+    di << "is edge or wire \n";
+      return 1;
+  }
+  TopoDS_Shape sh = DBRep::Get(a[2]);
+  if (sh.IsNull()) return 1;
+  TopAbs_ShapeEnum aShType = sh.ShapeType();
+  if (!(aShType == TopAbs_SHELL || aShType == TopAbs_FACE || aShType == TopAbs_EDGE || aShType == TopAbs_WIRE))
+  {
+    di << "Wrong shape type, shape can be shell or face or edge or wire\n";
+    return 1;
+  }
+
+  GeomAbs_SurfaceType aTargets[] = { GeomAbs_Plane, GeomAbs_Cylinder, GeomAbs_Cone, GeomAbs_Sphere };
+  Standard_Integer isurf = 0;
+  if (n > 3)
+  {
+    if (strcmp(a[3], "pln") == 0)
+      isurf = 0;
+    else if (strcmp(a[3], "cyl") == 0)
+      isurf = 1;
+    else if (strcmp(a[3], "con") == 0)
+      isurf = 2;
+    else if (strcmp(a[3], "sph") == 0)
+      isurf = 3;
+  }
+
+  Standard_Real tol = 1.e-7;
+  if (n > 4)
+    tol = Draw::Atof(a[4]);
+
+  // get sample  target for edge and wire
+  GeomAdaptor_Surface aSampleSurf;
+  if (n > 5 && (sh.ShapeType() == TopAbs_EDGE || sh.ShapeType() == TopAbs_WIRE ))
+  {
+    Handle(Geom_Surface) aGSurf = DrawTrSurf::GetSurface(a[5]);
+    if (aGSurf.IsNull())
+    {
+      di << "Sample surface is null" << "\n";
+      return 1;
+    }
+    aSampleSurf.Load(aGSurf);
+    GeomAbs_SurfaceType aSType = aSampleSurf.GetType();
+    if (aSType != aTargets[isurf])
+    {
+      di << "Sample surface has wrong type" << "\n";
+      return 1;
+    }
+  }
+
+  ShapeAnalysis_CanonicalRecognition aCanonRec(sh);
+  Handle(Geom_Surface) aRes;
+  switch (aTargets[isurf])
+  {
+  case GeomAbs_Plane:
+  {
+    gp_Pln aPln;
+    if (aSampleSurf.GetType() == GeomAbs_Plane)
+      aPln = aSampleSurf.Plane();
+    if (aCanonRec.IsPlane(tol, aPln))
+      aRes = new Geom_Plane(aPln);
+    break;
+  }
+  case GeomAbs_Cylinder:
+  {
+    gp_Cylinder aCyl;
+    if (aSampleSurf.GetType() == GeomAbs_Cylinder)
+      aCyl = aSampleSurf.Cylinder();
+    if (aCanonRec.IsCylinder(tol, aCyl))
+      aRes = new Geom_CylindricalSurface(aCyl);
+    break;
+  }
+  case GeomAbs_Cone:
+  {
+    gp_Cone aCon;
+    if (aSampleSurf.GetType() == GeomAbs_Cone)
+      aCon = aSampleSurf.Cone();
+    if (aCanonRec.IsCone(tol, aCon))
+      aRes = new Geom_ConicalSurface(aCon);
+    break;
+  }
+  case GeomAbs_Sphere:
+  {
+    gp_Sphere aSph;
+    if (aSampleSurf.GetType() == GeomAbs_Sphere)
+      aSph = aSampleSurf.Sphere();
+    if (aCanonRec.IsSphere(tol, aSph))
+      aRes = new Geom_SphericalSurface(aSph);
+    break;
+  }
+  default:
+    break;
+  }
+
+  if (!aRes.IsNull())
+  {
+    DrawTrSurf::Set(a[1], aRes);
+    di << "Gap = " << aCanonRec.GetGap() << "\n";
+  }
+  else
+  {
+    di << "Cannot get required surface" << "\n";
+  }
+  return 0;
+}
+//=======================================================================
+//function : getanacurve
+//purpose  : 
+//=======================================================================
+
+Standard_Integer getanacurve(Draw_Interpretor& di,
+  Standard_Integer n, const char** a)
+{
+  if (n < 3) {
+    di << "Usage: \n";
+    di << "getanacurve res shape [target [tol]] \n";
+    di << "target is reqired type of curve and can be: lin, cir, ell  \n";
+    return 1;
+  }
+  TopoDS_Shape sh = DBRep::Get(a[2]);
+  if (sh.IsNull()) return 1;
+  TopAbs_ShapeEnum aShType = sh.ShapeType();
+  if (!(aShType == TopAbs_WIRE || aShType == TopAbs_EDGE))
+  {
+    di << "Wrong shape type, shape can be wire or or edge \n";
+    return 1;
+  }
+
+  GeomAbs_CurveType aTargets[] = { GeomAbs_Line, GeomAbs_Circle, GeomAbs_Ellipse };
+  Standard_Integer icurv = 0;
+  if (n > 3)
+  {
+    if (strcmp(a[3],"lin") == 0)
+      icurv = 0;
+    else if (strcmp(a[3], "cir") == 0)
+      icurv = 1;
+    else if (strcmp(a[3], "ell") == 0)
+      icurv = 2;
+  }
+
+  Standard_Real tol = Precision::Confusion();
+  if (n > 4)
+      tol = Draw::Atof(a[4]);
+
+  ShapeAnalysis_CanonicalRecognition aCanonRec(sh);
+  Handle(Geom_Curve) aRes;
+  switch (aTargets[icurv])
+  {
+  case GeomAbs_Line:
+  {
+    gp_Lin aLin;
+    if (aCanonRec.IsLine(tol, aLin))
+      aRes = new Geom_Line(aLin);
+    break;
+  }
+  case GeomAbs_Circle:
+  {
+    gp_Circ aCirc;
+    if (aCanonRec.IsCircle(tol, aCirc))
+      aRes = new Geom_Circle(aCirc);
+    break;
+  }
+  case GeomAbs_Ellipse:
+  {
+    gp_Elips anElips;
+    if (aCanonRec.IsEllipse(tol, anElips))
+      aRes = new Geom_Ellipse(anElips);
+    break;
+  }
+  default:
+    break;
+  }
+
+  if (!aRes.IsNull())
+  {
+    DrawTrSurf::Set(a[1], aRes);
+    di << "Gap = " << aCanonRec.GetGap() << "\n";
+  }
+  else
+  {
+    di << "Cannot get required curve" << "\n";
+  }
+  return 0;
+
+}
 
 //=======================================================================
 //function : InitCommands
@@ -1020,4 +1226,8 @@ static Standard_Integer checkedge(Draw_Interpretor& di, Standard_Integer argc, c
   theCommands.Add("getareacontour","wire ",__FILE__, getareacontour, groupold);
   theCommands.Add ("checkselfintersection","wire [face]", __FILE__,checkselfintersection,g);
   theCommands.Add ("checkedge","edge [face]", __FILE__,checkedge,g);
+  theCommands.Add("getanasurf", "getanasurf res shape [target [tol [sample]]] ", __FILE__, getanasurf, g);
+  theCommands.Add("getanacurve", "getanacurve res shape [target [tol]]", __FILE__, getanacurve, g);
+
+
 }

src/ShapeAnalysis/FILES

@@ -45,3 +45,5 @@ ShapeAnalysis_WireOrder.cxx
 ShapeAnalysis_WireOrder.hxx
 ShapeAnalysis_WireVertex.cxx
 ShapeAnalysis_WireVertex.hxx
+ShapeAnalysis_CanonicalRecognition.cxx
+ShapeAnalysis_CanonicalRecognition.hxx

src/ShapeAnalysis/ShapeAnalysis_CanonicalRecognition.hxx

@@ -0,0 +1,169 @@
+// Copyright (c) 2022 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.
+
+#ifndef _ShapeAnalysis_CanonicalRecognition_HeaderFile
+#define _ShapeAnalysis_CanonicalRecognition_HeaderFile
+
+#include <Standard.hxx>
+#include <Standard_DefineAlloc.hxx>
+#include <Standard_Handle.hxx>
+
+#include <Standard_Boolean.hxx>
+#include <Standard_Integer.hxx>
+#include <TopAbs_ShapeEnum.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Shell.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Wire.hxx>
+#include <TopoDS_Edge.hxx>
+#include <GeomAbs_SurfaceType.hxx>
+#include <GeomAbs_CurveType.hxx>
+#include <GeomConvert_ConvType.hxx>
+#include <TColStd_Array1OfReal.hxx>
+
+class gp_Pln;
+class gp_Cone;
+class gp_Cylinder;
+class gp_Sphere;
+class gp_Lin;
+class gp_Circ;
+class gp_Elips;
+class Geom_Curve;
+class Geom_Surface;
+
+//! This class provides operators for analysis surfaces and curves of shapes
+//! in order to find out more simple geometry entities, which could replace
+//! existing complex (for exampe, BSpline) geometry objects with given tolerance.
+class ShapeAnalysis_CanonicalRecognition  
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+
+  //! Empty constructor
+  Standard_EXPORT ShapeAnalysis_CanonicalRecognition();
+  
+  //! constructor with shape initialisation
+  Standard_EXPORT ShapeAnalysis_CanonicalRecognition(const TopoDS_Shape& theShape);
+
+  //! Sets shape
+  Standard_EXPORT void SetShape(const TopoDS_Shape& theShape);
+  
+  //! Returns input shape
+  const TopoDS_Shape& GetShape() const
+  {
+    return myShape;
+  }
+
+  //! Returns deviation between input geometry entity and analytical entity
+  Standard_Real GetGap() const
+  {
+    return myGap;
+  }
+
+  //! Returns status of operation.
+  //! Current meaning of possible values of status:
+  //! -1 - algorithm is not initalazed by shape
+  //!  0 - no errors
+  //!  1 - error during any operation (usually - because of wrong input data)
+  //! Any operation (calling any methods like IsPlane(...), ...) can be performed 
+  //! when current staue is equal 0. 
+  //! If after any operation status != 0, it is necessary to set it 0 by method ClearStatus()
+  //! before calling other operation.
+  Standard_Integer GetStatus() const
+  {
+    return myStatus;
+  }
+
+  //! Returns status to be equal 0.
+  void ClearStatus() 
+  {
+    myStatus = 0;
+  }
+
+  //! Returns true if the underlined surface can be represent by plane with tolerance theTol
+  //! and sets in thePln the result plane.  
+  Standard_EXPORT Standard_Boolean IsPlane(const Standard_Real theTol, gp_Pln& thePln);
+
+  //! Returns true if the underlined surface can be represent by cylindrical one with tolerance theTol
+  //! and sets in theCyl the result cylinrical surface.  
+  Standard_EXPORT Standard_Boolean IsCylinder(const Standard_Real theTol, gp_Cylinder& theCyl);
+
+  //! Returns true if the underlined surface can be represent by conical one with tolerance theTol
+  //! and sets in theCone the result conical surface.  
+  Standard_EXPORT Standard_Boolean IsCone(const Standard_Real theTol, gp_Cone& theCone);
+
+  //! Returns true if the underlined surface can be represent by spherical one with tolerance theTol
+  //! and sets in theSphere the result spherical surface.  
+  Standard_EXPORT Standard_Boolean IsSphere(const Standard_Real theTol, gp_Sphere& theSphere);
+
+  //! Returns true if the underlined curve can be represent by line with tolerance theTol
+  //! and sets in theLin the result line.  
+  Standard_EXPORT Standard_Boolean IsLine(const Standard_Real theTol, gp_Lin& theLin);
+
+  //! Returns true if the underlined curve can be represent by circle with tolerance theTol
+  //! and sets in theCirc the result circle.  
+  Standard_EXPORT Standard_Boolean IsCircle(const Standard_Real theTol, gp_Circ& theCirc);
+
+  //! Returns true if the underlined curve can be represent by ellipse with tolerance theTol
+  //! and sets in theCirc the result ellipse.  
+  Standard_EXPORT Standard_Boolean IsEllipse(const Standard_Real theTol, gp_Elips& theElips);
+
+
+private:
+  Standard_Boolean IsElementarySurf(const GeomAbs_SurfaceType theTarget, const Standard_Real theTol,
+    gp_Ax3& thePos, TColStd_Array1OfReal& theParams);
+
+  Standard_Boolean IsConic(const GeomAbs_CurveType theTarget, const Standard_Real theTol,
+    gp_Ax2& thePos, TColStd_Array1OfReal& theParams);
+
+  static Handle(Geom_Surface) GetSurface(const TopoDS_Face& theFace, const Standard_Real theTol,
+    const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget,
+    Standard_Real& theGap, Standard_Integer& theStatus);
+
+  static Handle(Geom_Surface) GetSurface(const TopoDS_Shell& theShell, const Standard_Real theTol,
+    const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget,
+    Standard_Real& theGap, Standard_Integer& theStatus);
+
+  static Handle(Geom_Surface) GetSurface(const TopoDS_Edge& theEdge, const Standard_Real theTol,
+    const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget, 
+    gp_Ax3& thePos, TColStd_Array1OfReal& theParams, 
+    Standard_Real& theGap, Standard_Integer& theStatus);
+
+  static Handle(Geom_Surface) GetSurface(const TopoDS_Wire& theWire, const Standard_Real theTol,
+    const GeomConvert_ConvType theType, const GeomAbs_SurfaceType theTarget,
+    gp_Ax3& thePos, TColStd_Array1OfReal& theParams,
+    Standard_Real& theGap, Standard_Integer& theStatus);
+
+  static Handle(Geom_Curve) GetCurve(const TopoDS_Edge& theEdge, const Standard_Real theTol,
+    const GeomConvert_ConvType theType, const GeomAbs_CurveType theTarget,
+    Standard_Real& theGap, Standard_Integer& theStatus);
+
+  static Standard_Boolean GetSurfaceByLS(const TopoDS_Wire& theWire, const Standard_Real theTol,
+    const GeomAbs_SurfaceType theTarget,
+    gp_Ax3& thePos, TColStd_Array1OfReal& theParams,
+    Standard_Real& theGap, Standard_Integer& theStatus);
+
+  void Init(const TopoDS_Shape& theShape);
+
+private:
+
+  TopoDS_Shape myShape;
+  TopAbs_ShapeEnum mySType;
+  Standard_Real myGap;
+  Standard_Integer myStatus;
+
+};
+
+#endif // _ShapeAnalysis_CanonicalRecognition_HeaderFile

src/gce/gce_MakeCirc.cxx

@@ -54,7 +54,7 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt&  P1 ,
 //=========================================================================
   Standard_Real dist1, dist2, dist3, aResolution;
   //
-  aResolution=gp::Resolution();
+  aResolution = gp::Resolution();
   //
   dist1 = P1.Distance(P2);
   dist2 = P1.Distance(P3);
@@ -77,7 +77,7 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt&  P1 ,
   P2.Coord(x2,y2,z2);
   P3.Coord(x3,y3,z3);
   gp_Dir Dir1(x2-x1,y2-y1,z2-z1);
-  gp_Dir Dir2(x3-x2,y3-y2,z3-z2);
+  gp_Vec VDir2(x3-x2,y3-y2,z3-z2);
   //
   gp_Ax1 anAx1(P1, Dir1);
   gp_Lin aL12 (anAx1);
@@ -86,14 +86,21 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt&  P1 ,
     return;
   }
   //
-  gp_Dir Dir3 = Dir1.Crossed(Dir2);
+  gp_Vec VDir1(Dir1);
+  gp_Vec VDir3 = VDir1.Crossed(VDir2);
+  if(VDir3.SquareMagnitude() < aResolution)
+  {
+    TheError = gce_ColinearPoints;
+    return;
+  }
   //
+  gp_Dir Dir3(VDir3);
   gp_Dir dir = Dir1.Crossed(Dir3);
   gp_Lin L1(gp_Pnt((P1.XYZ()+P2.XYZ())/2.),dir);
-  dir = Dir2.Crossed(Dir3);
+  dir = VDir2.Crossed(Dir3);
   gp_Lin L2(gp_Pnt((P3.XYZ()+P2.XYZ())/2.),dir);
 
-  Standard_Real Tol = 0.000000001;
+  Standard_Real Tol = Precision::PConfusion();
   Extrema_ExtElC distmin(L1,L2,Tol);
   
   if (!distmin.IsDone()) { 
@@ -140,7 +147,7 @@ gce_MakeCirc::gce_MakeCirc(const gp_Pnt&  P1 ,
       //Dir2 = gp_Dir(x2-x3,y2-y3,z2-z3);
       //modified by NIZNHY-PKV Thu Mar  3 11:31:13 2005t
       //
-      TheCirc = gp_Circ(gp_Ax2(pInt, Dir3, Dir1),(dist1+dist2+dist3)/3.);
+      TheCirc = gp_Circ(gp_Ax2(pInt, gp_Dir(VDir3), Dir1),(dist1+dist2+dist3)/3.);
       TheError = gce_Done;
     }
   }