Using gradient descent algorithm to find the lowest point on a surface.

jianbaoxia

CAD master
Hello, I‘m JianBaoxia.
I'm trying to find the lowest point (the point with the minimum z-coordinate) on a surface. It seems that directly traversing all (u, v) coordinates uniformly may not guarantee accuracy and efficiency. So, I'm attempting to use the gradient descent algorithm to find the lowest point.
The idea behind this algorithm is to move along the negative direction of the u and v derivatives on the surface to determine the next (u, v) values. However, the derivative type returned by Geom_Surface->D1() is gp_Vec, which represents (x, y, z) coordinates. This is where I'm stuck. I can't simply calculate the next point on the surface based on the (x, y, z) values. Although I can easily update the (x, y, z) coordinates of a point and use GeomAPI_ProjectPointOnSurf to obtain the point on the surface closest to it, this approach doesn't perform well in terms of accuracy and efficiency.
Can you please provide some suggestions? I would greatly appreciate it.
 

Quaoar

Administrator
Staff member
If your optimization is done in the UV space of a surface, then I do not quite get why would you need point projection. It looks like a two-dimensional problem where your objective function F = Min_Z( S(u,v) ) is defined by two unknowns: u and v.

I quickly tried the following scheme:
  1. Apply PSO (particle swarm optimizer) to find a course solution.
  2. Apply gradient descent with a gradient computed as finite difference (that was the easiest way but I bet you can do something more intelligent to compute a gradient). The result of PSO is used as the initial guess in the gradient descent optimizer.
1684338185131.png

I used two functions defined as follows (note that it's all based on Analysis Situs framework but it still should be understandable enough):

Code:
//! Objective function.
class LowestZ : public asiAlgo_Function<Prim_UV>
{
public:

  //! Ctor.
  LowestZ(const Handle(Geom_Surface)& surf) : m_surf(surf) {}

public:

  inline virtual double Value(const Prim_UV& uv)
  {
    gp_Pnt P = m_surf->Value( uv.U(), uv.V() );

    return P.Z();
  }

protected:

  Handle(Geom_Surface) m_surf;

};

//! Objective function with gradient.
class LowestZ_Grad : public asiAlgo_FunctionWithGradient<Prim_UV>
{
public:

  LowestZ_Grad(const Handle(Geom_Surface)& surf) : m_surf(surf)
  {}

public:

  inline virtual double Value(const Prim_UV& uv)
  {
    gp_Pnt P = m_surf->Value( uv.U(), uv.V() );

    return P.Z();
  }

  inline virtual t_coord Gradient(const Prim_UV& uv)
  {
    const double d = 0.01;
    const double z1 = Value(uv);
    const double z2x = Value(uv + Prim_UV(d,0));
    const double z2y = Value(uv + Prim_UV(0,d));
    const Prim_UV grad( (z2x - z1)/d, (z2y - z1)/d );
    return grad;
  }

protected:

  Handle(Geom_Surface) m_surf;

};

The first function (LowestZ) comes without a gradient and the second one (LowestZ_Grad) computes gradient with this naive forward finite difference approach. The Prim_UV is like gp_XY but with zero-based indexation, so that I could use PSO and gradient descent from Analysis Situs (these classes are defined here: https://gitlab.com/ssv/AnalysisSitus/-/tree/master/src/asiAlgo/opt).

To run optimization, I used:

Code:
  Handle(Geom_Surface) surf = BRep_Tool::Surface(F);
  double uMin, uMax, vMin, vMax;
  BRepTools::UVBounds(F, uMin, uMax, vMin, vMax);

  // Objective function.
  LowestZ func(surf);

  asiAlgo_PSO<Prim_UV>::t_search_params pso_params;
  //
  pso_params.num_particles   = 10;                   // In PSO the number of agents is typically small
  pso_params.num_iterations  = 5000;                 // More than enough in typical situations
  pso_params.area.corner_min = Prim_UV(uMin, vMin);  // Min corner of the search area
  pso_params.area.corner_max = Prim_UV(uMax, vMax);  // Max corner of the search area
  pso_params.precision       = 1.0e-6;               // Less values are of no practical sense
  pso_params.pFunc           = &func;                // Objective function
  pso_params.m               = 0.729;                // Retain weight
  pso_params.mu_cognition    = 0.4;                  // Pure social behavior
  pso_params.nu_social       = 1.49445;              // Social determinant of an agent
  //
  asiAlgo_PSO<Prim_UV> pso(pso_params);

  int pso_num_iters = 0;
  pso.Perform(pso_num_iters);
  const asiAlgo_PSO<Prim_UV>::t_measuring& pso_sol = pso.GetBestGlobal();

  interp->GetProgress().SendLogMessage(LogInfo(Normal) << "Iterations done: %1." << pso_num_iters);
  interp->GetPlotter().REDRAW_POINT("sol", surf->Value(pso_sol.p.U(), pso_sol.p.V() ), Color_Red);

  double u_sol, v_sol;
  LowestZ_Grad F_Grad(surf);

  // Gradient descent parameters
  asiAlgo_GradientDescent<Prim_UV>::t_search_params grad_params;
  grad_params.num_iterations = 10000;
  grad_params.pFunc          = &F_Grad;
  grad_params.precision      = 1.0e-9;
  grad_params.start          = pso_sol.p;
  grad_params.is_adaptive    = true;
  grad_params.step           = 1.0e-1; // To be corrected by Armijo rule

  // Run local optimization
  int grad_num_iters = 0;
  asiAlgo_GradientDescent<Prim_UV> Descent(grad_params);
  if ( !Descent.Perform(grad_num_iters) )
  {
    interp->GetProgress().SendLogMessage(LogWarn(Normal) << "Gradient descent did not converge.");
  }
  const Prim_UV& grad_sol = Descent.Solution();

  interp->GetPlotter().REDRAW_POINT("grad_sol", surf->Value(grad_sol.U(), grad_sol.V() ), Color_Green);

The funny thing is that gradient descent did not really converge so I simply used what it returned at the end. Not sure if all this makes sense, and I did not really test it, because it would be too time consuming ;) But it's kinda how I would start looking at this problem myself probably. Also note that neither PSO parameters nor finite difference step are calibrated, and I used kinda random numbers to give it a try.

The test surface is attached (test-pso.stp).
 

Attachments

  • test-pso.stp
    19.7 KB · Views: 1
Last edited:

jianbaoxia

CAD master
If your optimization is done in the UV space of a surface, then I do not quite get why would you need point projection.
Thank you very much for your professional response and suggestions.

As you mentioned, I indeed want to implement the gradient descent algorithm within the two-dimensional UV parameter domain of the surface. This requires obtaining the partial derivatives of the surface along the u and v directions at the current point. However, the derivatives computed by Geom_Surface->D1() in OCCT are in the three-dimensional Cartesian coordinate system (x, y, z), which has been a challenge for me.

Your suggestion of using finite differences to calculate the gradients has given me inspiration and proved to be very helpful. I sincerely appreciate your advice. 😃 😃
 
Top