OCCT

// Created on: 1997-01-24
// Created by: Laurent BOURESCHE
// Copyright (c) 1997-1999 Matra Datavision
// Copyright (c) 1999-2014 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 <Adaptor2d_Curve2d.hxx>
#include <Adaptor3d_TopolTool.hxx>
#include <Blend_FuncInv.hxx>
#include <Blend_SurfCurvFuncInv.hxx>
#include <Blend_SurfPointFuncInv.hxx>
#include <Blend_SurfRstFunction.hxx>
#include <BRepBlend_BlendTool.hxx>
#include <BRepBlend_Line.hxx>
#include <BRepBlend_SurfRstLineBuilder.hxx>
#include <gp_Pnt.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec.hxx>
#include <gp_Vec2d.hxx>
#include <IntSurf.hxx>
#include <IntSurf_Transition.hxx>
#include <math_FunctionSetRoot.hxx>
#include <TopAbs.hxx>

//=======================================================================
// function :  ArcToRecadre
// purpose  : Find a suitable arc
//           PrevIndex is used to reject an already tested arc
//=======================================================================

int BRepBlend_SurfRstLineBuilder::ArcToRecadre(const math_Vector& theSol,
                                               const int          PrevIndex,
                                               gp_Pnt2d&          lastpt2d,
                                               gp_Pnt2d&          pt2d,
                                               double&            ponarc)
{
  int    IndexSol = 0, nbarc = 0;
  bool   ok      = false;
  bool   byinter = (line->NbPoints() != 0), okinter = false;
  double distmin = RealLast();
  double uprev = 0., vprev = 0., prm = 0., dist = 0.;

  if (byinter)
  {
    previousP.ParametersOnS(uprev, vprev);
  }
  pt2d.SetCoord(theSol(1), theSol(2));
  lastpt2d.SetCoord(uprev, vprev);
  domain1->Init();

  while (domain1->More())
  {
    nbarc++;
    ok = false;
    if (byinter)
    {
      ok = okinter =
        BRepBlend_BlendTool::Inters(pt2d, lastpt2d, surf1, domain1->Value(), prm, dist);
    }
    if (!ok)
    {
      ok = BRepBlend_BlendTool::Project(pt2d, surf1, domain1->Value(), prm, dist);
    }

    if (ok && (nbarc != PrevIndex))
    {
      if (dist < distmin || okinter)
      {
        distmin  = dist;
        ponarc   = prm;
        IndexSol = nbarc;
        if (okinter && (PrevIndex == 0))
        {
          break;
        }
      }
    }
    domain1->Next();
  }
  return IndexSol;
}

//=================================================================================================

BRepBlend_SurfRstLineBuilder::BRepBlend_SurfRstLineBuilder(
  const occ::handle<Adaptor3d_Surface>&   Surf1,
  const occ::handle<Adaptor3d_TopolTool>& Domain1,
  const occ::handle<Adaptor3d_Surface>&   Surf2,
  const occ::handle<Adaptor2d_Curve2d>&   Rst,
  const occ::handle<Adaptor3d_TopolTool>& Domain2)
    : done(false),
      sol(1, 3),
      surf1(Surf1),
      domain1(Domain1),
      surf2(Surf2),
      rst(Rst),
      domain2(Domain2),
      tolpoint3d(0.0),
      tolpoint2d(0.0),
      tolgui(0.0),
      pasmax(0.0),
      fleche(0.0),
      param(0.0),
      rebrou(false),
      iscomplete(false),
      comptra(false),
      sens(0.0),
      decrochdeb(false),
      decrochfin(false)
{
}

//=================================================================================================

void BRepBlend_SurfRstLineBuilder::Perform(Blend_SurfRstFunction&  Func,
                                           Blend_FuncInv&          Finv,
                                           Blend_SurfPointFuncInv& FinvP,
                                           Blend_SurfCurvFuncInv&  FinvC,
                                           const double            Pdep,
                                           const double            Pmax,
                                           const double            MaxStep,
                                           const double            Tol3d,
                                           const double            Tol2d,
                                           const double            TolGuide,
                                           const math_Vector&      ParDep,
                                           const double            Fleche,
                                           const bool              Appro)
{
  done       = false;
  iscomplete = false;
  comptra    = false;
  line       = new BRepBlend_Line();
  tolpoint3d = Tol3d;
  tolpoint2d = Tol2d;
  tolgui     = std::abs(TolGuide);
  fleche     = std::abs(Fleche);
  rebrou     = false;
  pasmax     = std::abs(MaxStep);

  if (Pmax - Pdep >= 0.)
  {
    sens = 1.;
  }
  else
  {
    sens = -1.;
  }

  Blend_Status State;

  param = Pdep;
  Func.Set(param);

  if (Appro)
  {
    TopAbs_State siturst, situs;
    bool         decroch;
    math_Vector  tolerance(1, 3), infbound(1, 3), supbound(1, 3);
    Func.GetTolerance(tolerance, tolpoint3d);
    Func.GetBounds(infbound, supbound);
    math_FunctionSetRoot rsnld(Func, tolerance, 30);

    rsnld.Perform(Func, ParDep, infbound, supbound);

    if (!rsnld.IsDone())
    {
      return;
    }
    rsnld.Root(sol);
    if (!CheckInside(Func, siturst, situs, decroch))
    {
      return;
    }
  }
  else
  {
    sol = ParDep;
  }

  State = TestArret(Func, false, Blend_OK);
  if (State != Blend_OK)
  {
    return;
  }
  // Update the line.
  line->Append(previousP);
  double U, V;
  previousP.ParametersOnS(U, V);
  //  W = previousP.ParameterOnC();

  BRepBlend_Extremity ptf1(previousP.PointOnS(), U, V, previousP.Parameter(), tolpoint3d);
  BRepBlend_Extremity ptf2(previousP.PointOnC(), U, V, previousP.Parameter(), tolpoint3d);
  if (!previousP.IsTangencyPoint())
  {
    ptf1.SetTangent(previousP.TangentOnS());
    ptf2.SetTangent(previousP.TangentOnC());
  }
  if (sens > 0.)
  {
    line->SetStartPoints(ptf1, ptf2);
  }
  else
  {
    line->SetEndPoints(ptf1, ptf2);
  }
  InternalPerform(Func, Finv, FinvP, FinvC, Pmax);
  done = true;
}

//=================================================================================================

bool BRepBlend_SurfRstLineBuilder::PerformFirstSection(Blend_SurfRstFunction&  Func,
                                                       Blend_FuncInv&          Finv,
                                                       Blend_SurfPointFuncInv& FinvP,
                                                       Blend_SurfCurvFuncInv&  FinvC,
                                                       const double            Pdep,
                                                       const double            Pmax,
                                                       const math_Vector&      ParDep,
                                                       const double            Tol3d,
                                                       const double            Tol2d,
                                                       const double            TolGuide,
                                                       const bool              RecRst,
                                                       const bool              RecP,
                                                       const bool              RecS,
                                                       double&                 Psol,
                                                       math_Vector&            ParSol)
{
  done       = false;
  iscomplete = false;
  comptra    = false;
  line       = new BRepBlend_Line();
  tolpoint3d = Tol3d;
  tolpoint2d = Tol2d;
  tolgui     = std::abs(TolGuide);
  rebrou     = false;

  if (Pmax - Pdep >= 0.)
  {
    sens = 1.;
  }
  else
  {
    sens = -1.;
  }
  Blend_Status                   State = Blend_OnRst12;
  double                         trst  = 0.;
  bool                           recadp, recadrst, recads;
  double                         wp, wrst, ws;
  double                         U = 0., V = 0.;
  math_Vector                    infbound(1, 3), supbound(1, 3), tolerance(1, 3);
  math_Vector                    solinvp(1, 3), solinvrst(1, 4), solinvs(1, 3);
  occ::handle<Adaptor3d_HVertex> Vtxp, Vtxrst, Vtxs, Vtxc;
  bool                           IsVtxp = false, IsVtxrst = false, IsVtxs = false;
  occ::handle<Adaptor2d_Curve2d> Arc;
  wp = wrst = ws = Pmax;
  param          = Pdep;
  Func.Set(param);
  Func.GetTolerance(tolerance, tolpoint3d);
  Func.GetBounds(infbound, supbound);

  math_FunctionSetRoot rsnld(Func, tolerance, 30);
  rsnld.Perform(Func, ParDep, infbound, supbound);
  if (!rsnld.IsDone())
  {
    return false;
  }
  rsnld.Root(sol);

  recads = RecS && Recadre(FinvC, solinvs, Arc, IsVtxs, Vtxs);
  if (recads)
  {
    ws = solinvs(1);
  }
  recadp = RecP && Recadre(FinvP, solinvp, IsVtxp, Vtxp);
  if (recadp)
  {
    wp = solinvp(1);
  }
  recadrst = RecRst && Recadre(Func, Finv, solinvrst, IsVtxrst, Vtxrst);
  if (recadrst)
  {
    wrst = solinvrst(2);
  }
  if (!recads && !recadp && !recadrst)
  {
    return false;
  }
  if (recadp && recadrst)
  {
    if (sens * (wrst - wp) > tolgui)
    { // first one leaves the domain
      wrst     = wp;
      U        = solinvp(2);
      V        = solinvp(3);
      trst     = BRepBlend_BlendTool::Parameter(Vtxp, rst);
      IsVtxrst = IsVtxp;
      Vtxrst   = Vtxp;
    }
    else
    {
      U    = solinvrst(3);
      V    = solinvrst(4);
      trst = solinvrst(1);
    }
  }
  else if (recadp)
  {
    wrst     = wp;
    U        = solinvp(2);
    V        = solinvp(3);
    trst     = BRepBlend_BlendTool::Parameter(Vtxp, rst);
    IsVtxrst = IsVtxp;
    Vtxrst   = Vtxp;
    recadrst = true;
  }
  else if (recadrst)
  {
    U    = solinvrst(3);
    V    = solinvrst(4);
    trst = solinvrst(1);
  }
  if (recads && recadrst)
  {
    if (std::abs(ws - wrst) < tolgui)
    {
      State  = Blend_OnRst12;
      param  = 0.5 * (ws + wrst);
      sol(1) = U;
      sol(2) = V;
      sol(3) = solinvs(2);
    }
    else if (sens * (ws - wrst) < 0)
    {
      // ground on surf
      State = Blend_OnRst1;
      param = ws;
      Arc->Value(solinvs(3)).Coord(U, V);
      sol(1) = U;
      sol(2) = V;
      sol(3) = solinvs(2);
    }
    else
    {
      // ground on rst
      State  = Blend_OnRst2;
      param  = wrst;
      sol(1) = U;
      sol(2) = V;
      sol(3) = trst;
    }
    Func.Set(param);
  }
  else if (recads)
  {
    // ground on surf
    State = Blend_OnRst1;
    param = ws;
    Arc->Value(solinvs(3)).Coord(U, V);
    sol(1) = U;
    sol(2) = V;
    sol(3) = solinvs(2);
    Func.Set(param);
  }
  else if (recadrst)
  {
    // ground on rst
    State  = Blend_OnRst2;
    param  = wrst;
    sol(1) = U;
    sol(2) = V;
    sol(3) = trst;
    Func.Set(param);
  }
  State  = TestArret(Func, false, State);
  Psol   = param;
  ParSol = sol;
  return true;
}

//=================================================================================================

bool BRepBlend_SurfRstLineBuilder::Complete(Blend_SurfRstFunction&  Func,
                                            Blend_FuncInv&          Finv,
                                            Blend_SurfPointFuncInv& FinvP,
                                            Blend_SurfCurvFuncInv&  FinvC,
                                            const double            Pmin)
{
  if (!done)
  {
    throw StdFail_NotDone();
  }
  if (iscomplete)
  {
    return true;
  }
  if (sens > 0.)
  {
    previousP = line->Point(1);
  }
  else
  {
    previousP = line->Point(line->NbPoints());
  }
  sens  = -sens;
  param = previousP.Parameter();
  previousP.ParametersOnS(sol(1), sol(2));
  sol(3) = previousP.ParameterOnC();

  InternalPerform(Func, Finv, FinvP, FinvC, Pmin);
  iscomplete = true;
  return true;
}

//=================================================================================================

void BRepBlend_SurfRstLineBuilder::InternalPerform(Blend_SurfRstFunction&  Func,
                                                   Blend_FuncInv&          Finv,
                                                   Blend_SurfPointFuncInv& FinvP,
                                                   Blend_SurfCurvFuncInv&  FinvC,
                                                   const double            Bound)
{
  double stepw = pasmax;
  int    nbp   = line->NbPoints();
  if (nbp >= 2)
  { // The last step is reproduced if it is not too small.
    if (sens < 0.)
    {
      stepw = (line->Point(2).Parameter() - line->Point(1).Parameter());
    }
    else
    {
      stepw = (line->Point(nbp).Parameter() - line->Point(nbp - 1).Parameter());
    }
    stepw = std::max(stepw, 100. * tolgui);
  }
  double parprec = param;
  if (sens * (parprec - Bound) >= -tolgui)
  {
    return;
  }
  Blend_Status                   State   = Blend_OnRst12;
  TopAbs_State                   situonc = TopAbs_UNKNOWN, situons = TopAbs_UNKNOWN;
  bool                           decroch = false;
  bool                           Arrive, recadp, recadrst, recads, echecrecad;
  double                         wp, wrst, ws;
  double                         U = 0., V = 0.;
  double                         trst = 0.;
  math_Vector                    infbound(1, 3), supbound(1, 3);
  math_Vector                    parinit(1, 3), tolerance(1, 3);
  math_Vector                    solinvp(1, 3), solinvrst(1, 4), solinvs(1, 3);
  occ::handle<Adaptor3d_HVertex> Vtxp, Vtxrst, Vtxs, Vtxc;
  bool                           IsVtxp = false, IsVtxrst = false, IsVtxs = false;
  BRepBlend_Extremity            Extrst, Exts;
  occ::handle<Adaptor2d_Curve2d> Arc;

  // IntSurf_Transition Tline,Tarc;

  Func.GetTolerance(tolerance, tolpoint3d);
  Func.GetBounds(infbound, supbound);

  math_FunctionSetRoot rsnld(Func, tolerance, 30);
  parinit = sol;

  Arrive = false;
  param  = parprec + sens * stepw;
  if (sens * (param - Bound) > 0.)
  {
    stepw = sens * (Bound - parprec) * 0.5;
    param = parprec + sens * stepw;
  }

  while (!Arrive)
  {
    bool bonpoint = true;
    Func.Set(param);
    rsnld.Perform(Func, parinit, infbound, supbound);

    if (rsnld.IsDone())
    {
      rsnld.Root(sol);
      if (!CheckInside(Func, situonc, situons, decroch) && line->NbPoints() == 1)
      {
        State    = Blend_StepTooLarge;
        bonpoint = false;
      }
    }
    else
    {
      State    = Blend_StepTooLarge;
      bonpoint = false;
    }
    if (bonpoint)
    {
      wp = wrst = ws = Bound;
      recadp = recadrst = recads = false;
      echecrecad                 = false;
      if (situons == TopAbs_OUT || situons == TopAbs_ON)
      {
        // pb inverse rst/rst
        recads = Recadre(FinvC, solinvs, Arc, IsVtxs, Vtxs);
        if (recads)
        {
          ws = solinvs(1);
          // It is necessary to reevaluate the deviation (BUC60360)
          gp_Vec t, n;
          Func.Set(ws);
          Arc->Value(solinvs(3)).Coord(U, V);
          sol(1)  = U;
          sol(2)  = V;
          sol(3)  = solinvs(2);
          decroch = Func.Decroch(sol, t, n);
        }
        else
        {
          echecrecad = true;
        }
      }
      if (situonc == TopAbs_OUT || situonc == TopAbs_ON)
      {
        // pb inverse point/surf
        recadp = Recadre(FinvP, solinvp, IsVtxp, Vtxp);
        if (recadp)
        {
          wp = solinvp(1);
        }
        else
        {
          echecrecad = true;
        }
      }
      if (decroch)
      {
        // pb inverse rst/surf
        recadrst = Recadre(Func, Finv, solinvrst, IsVtxrst, Vtxrst);
        if (recadrst)
        {
          wrst = solinvrst(2);
        }
        else
        {
          echecrecad = true;
        }
      }
      decroch = false;
      if (recadp || recads || recadrst)
      {
        echecrecad = false;
      }
      if (!echecrecad)
      {
        if (recadp && recadrst)
        {
          if (sens * (wrst - wp) > tolgui)
          { // first one leaves the domain
            wrst     = wp;
            U        = solinvp(2);
            V        = solinvp(3);
            trst     = BRepBlend_BlendTool::Parameter(Vtxp, rst);
            IsVtxrst = IsVtxp;
            Vtxrst   = Vtxp;
          }
          else
          {
            decroch = true;
            U       = solinvrst(3);
            V       = solinvrst(4);
            trst    = solinvrst(1);
          }
        }
        else if (recadp)
        {
          wrst     = wp;
          U        = solinvp(2);
          V        = solinvp(3);
          trst     = BRepBlend_BlendTool::Parameter(Vtxp, rst);
          IsVtxrst = IsVtxp;
          Vtxrst   = Vtxp;
          recadrst = true;
        }
        else if (recadrst)
        {
          decroch = true;
          U       = solinvrst(3);
          V       = solinvrst(4);
          trst    = solinvrst(1);
        }
        if (recads && recadrst)
        {
          if (std::abs(ws - wrst) < tolgui)
          {
            State  = Blend_OnRst12;
            param  = 0.5 * (ws + wrst);
            sol(1) = U;
            sol(2) = V;
            sol(3) = solinvs(3);
          }
          else if (sens * (ws - wrst) < 0)
          {
            // ground on surf
            decroch = false;
            State   = Blend_OnRst1;
            param   = ws;
            Arc->Value(solinvs(3)).Coord(U, V);
            sol(1) = U;
            sol(2) = V;
            sol(3) = solinvs(2);
          }
          else
          {
            // ground on rst
            State  = Blend_OnRst2;
            param  = wrst;
            sol(1) = U;
            sol(2) = V;
            sol(3) = trst;
          }
          Func.Set(param);
        }
        else if (recads)
        {
          // ground on surf
          State = Blend_OnRst1;
          param = ws;
          Arc->Value(solinvs(3)).Coord(U, V);
          sol(1) = U;
          sol(2) = V;
          sol(3) = solinvs(2);
          Func.Set(param);
        }
        else if (recadrst)
        {
          // ground on rst
          State  = Blend_OnRst2;
          param  = wrst;
          sol(1) = U;
          sol(2) = V;
          sol(3) = trst;
          Func.Set(param);
        }
        else
        {
          State = Blend_OK;
        }
        State = TestArret(Func, true, State);
      }
      else
      {
        // Failed reframing. Leave with PointsConfondus
#ifdef OCCT_DEBUG
        std::cout << "SurfRstLineBuilder : failed reframing" << std::endl;
#endif
        State = Blend_SamePoints;
      }
    }

    switch (State)
    {
      case Blend_OK: {
        // Update the line.
        if (sens > 0.)
        {
          line->Append(previousP);
        }
        else
        {
          line->Prepend(previousP);
        }
        parinit = sol;
        parprec = param;

        if (param == Bound)
        {
          Arrive = true;
          Exts.SetValue(previousP.PointOnS(), sol(1), sol(2), previousP.Parameter(), tolpoint3d);
          MakeExtremity(Extrst, false, rst, sol(3), IsVtxrst, Vtxrst);
          // Indicate end on Bound.
        }
        else
        {
          param = param + sens * stepw;
          if (sens * (param - Bound) > -tolgui)
          {
            param = Bound;
          }
        }
      }
      break;

      case Blend_StepTooLarge: {
        stepw = stepw / 2.;
        if (std::abs(stepw) < tolgui)
        {
          previousP.ParametersOnS(U, V);
          Exts.SetValue(previousP.PointOnS(), U, V, previousP.Parameter(), tolpoint3d);
          Extrst.SetValue(previousP.PointOnC(),
                          previousP.ParameterOnC(),
                          previousP.Parameter(),
                          tolpoint3d);
          Arrive = true;
          if (line->NbPoints() >= 2)
          {
            // Indicate that one stops during the processing
#ifdef OCCT_DEBUG
            std::cout << "SurfRstLineBuilder : No advancement in the processing" << std::endl;
#endif
          }
        }
        else
        {
          param = parprec + sens * stepw; // no risk to exceed Bound.
        }
      }
      break;

      case Blend_StepTooSmall: {
        // Update the line.
        if (sens > 0.)
        {
          line->Append(previousP);
        }
        else
        {
          line->Prepend(previousP);
        }
        parinit = sol;
        parprec = param;

        stepw = std::min(1.5 * stepw, pasmax);
        if (param == Bound)
        {
          Arrive = true;
          Exts.SetValue(previousP.PointOnS(), sol(1), sol(2), previousP.Parameter(), tolpoint3d);
          MakeExtremity(Extrst, false, rst, sol(3), IsVtxrst, Vtxrst);
          // Indicate end on Bound.
        }
        else
        {
          param = param + sens * stepw;
          if (sens * (param - Bound) > -tolgui)
          {
            param = Bound;
          }
        }
      }
      break;

      case Blend_OnRst1: {
        if (sens > 0.)
        {
          line->Append(previousP);
        }
        else
        {
          line->Prepend(previousP);
        }
        MakeExtremity(Exts, true, Arc, solinvs(3), IsVtxs, Vtxs);
        MakeExtremity(Extrst, false, rst, sol(3), IsVtxrst, Vtxrst);
        Arrive = true;
      }
      break;

      case Blend_OnRst2: {
        if (sens > 0.)
        {
          line->Append(previousP);
        }
        else
        {
          line->Prepend(previousP);
        }
        Exts.SetValue(previousP.PointOnS(), sol(1), sol(2), previousP.Parameter(), tolpoint3d);
        MakeExtremity(Extrst, false, rst, sol(3), IsVtxrst, Vtxrst);
        Arrive = true;
      }
      break;

      case Blend_OnRst12: {
        if (sens > 0.)
        {
          line->Append(previousP);
        }
        else
        {
          line->Prepend(previousP);
        }
        MakeExtremity(Exts, true, Arc, solinvs(1), IsVtxs, Vtxs);
        MakeExtremity(Extrst, false, rst, sol(3), IsVtxrst, Vtxrst);
        Arrive = true;
      }
      break;

      case Blend_SamePoints: {
        // Stop
#ifdef OCCT_DEBUG
        std::cout << "SurfRstLineBuilder Points mixed in the processing" << std::endl;
#endif
        previousP.ParametersOnS(U, V);
        Exts.SetValue(previousP.PointOnS(), U, V, previousP.Parameter(), tolpoint3d);
        Extrst.SetValue(previousP.PointOnC(),
                        previousP.ParameterOnC(),
                        previousP.Parameter(),
                        tolpoint3d);
        Arrive = true;
      }
      break;
      default:
        break;
    }
    if (Arrive)
    {
      if (sens > 0.)
      {
        line->SetEndPoints(Exts, Extrst);
        decrochfin = decroch;
      }
      else
      {
        line->SetStartPoints(Exts, Extrst);
        decrochdeb = decroch;
      }
    }
  }
}

//=======================================================================
// function : Recadre
// purpose  : Reframe section Surface / Restriction
//=======================================================================

bool BRepBlend_SurfRstLineBuilder::Recadre(Blend_SurfCurvFuncInv&          FinvC,
                                           math_Vector&                    Solinv,
                                           occ::handle<Adaptor2d_Curve2d>& Arc,
                                           bool&                           IsVtx,
                                           occ::handle<Adaptor3d_HVertex>& Vtx)
{
  bool recadre = false;

  gp_Pnt2d pt2d, lastpt2d;
  int      IndexSol, nbarc;
  double   pmin;

  IndexSol = ArcToRecadre(sol, 0, lastpt2d, pt2d, pmin);

  IsVtx = false;
  if (IndexSol == 0)
  {
    return false;
  }

  domain1->Init();
  nbarc = 1;
  while (nbarc < IndexSol)
  {
    nbarc++;
    domain1->Next();
  }
  Arc = domain1->Value();

  FinvC.Set(Arc);

  math_Vector toler(1, 3), infb(1, 3), supb(1, 3);
  // use reduced Tol argument value to pass testcase
  // blend complex A6 with scale factor of model 0.1 (base scale = 1000)
  // So, here we using 1.0e-5 rather than 1.0e-4 of tolerance of point in 3d
  FinvC.GetTolerance(toler, 0.1 * tolpoint3d);
  FinvC.GetBounds(infb, supb);
  Solinv(1) = param;
  Solinv(2) = sol(3);
  Solinv(3) = pmin;

  math_FunctionSetRoot rsnld(FinvC, toler, 30);
  rsnld.Perform(FinvC, Solinv, infb, supb);

  if (!rsnld.IsDone())
  {
#ifdef OCCT_DEBUG
    std::cout << "SurfRstLineBuilder : RSNLD not done " << std::endl << std::endl;
#endif
  }
  else
  {
    // It is necessary to check the value of the function
    rsnld.Root(Solinv);
    recadre = FinvC.IsSolution(Solinv, tolpoint3d);
  }

  // In case of fail, it is checked if another arc
  // can be useful (case of output at the proximity of a vertex)
  if (!recadre)
  {

    IndexSol = ArcToRecadre(sol, IndexSol, lastpt2d, pt2d, pmin);
    if (IndexSol == 0)
    {
      return false; // No other solution
    }

    domain1->Init();
    nbarc = 1;
    while (nbarc < IndexSol)
    {
      nbarc++;
      domain1->Next();
    }

    Arc = domain1->Value();
    FinvC.Set(Arc);

    FinvC.GetTolerance(toler, tolpoint3d);
    FinvC.GetBounds(infb, supb);

    Solinv(3) = pmin;

    math_FunctionSetRoot aRsnld(FinvC, toler, 30);
    aRsnld.Perform(FinvC, Solinv, infb, supb);

    if (!aRsnld.IsDone())
    {
#ifdef OCCT_DEBUG
      std::cout << "SurfRstLineBuilder : RSNLD not done " << std::endl << std::endl;
#endif
    }
    else
    {
      // It is necessary to check the value of the function
      aRsnld.Root(Solinv);
      recadre = FinvC.IsSolution(Solinv, tolpoint3d);
    }
  }

  if (recadre)
  {
    double w = Solinv(2);
    if (w < rst->FirstParameter() - toler(2) || w > rst->LastParameter() + toler(2))
    {
      return false;
    }
    domain1->Initialize(Arc);
    domain1->InitVertexIterator();
    IsVtx = !domain1->MoreVertex();
    while (!IsVtx)
    {
      Vtx = domain1->Vertex();
      if (std::abs(BRepBlend_BlendTool::Parameter(Vtx, Arc) - Solinv(3))
          <= BRepBlend_BlendTool::Tolerance(Vtx, Arc))
      {
        IsVtx = true;
      }
      else
      {
        domain1->NextVertex();
        IsVtx = !domain1->MoreVertex();
      }
    }
    if (!domain1->MoreVertex())
    {
      IsVtx = false;
    }
    return true;
  }
  return false;
}

//=================================================================================================

bool BRepBlend_SurfRstLineBuilder::Recadre(Blend_SurfRstFunction&          Func,
                                           Blend_FuncInv&                  Finv,
                                           math_Vector&                    Solinv,
                                           bool&                           IsVtx,
                                           occ::handle<Adaptor3d_HVertex>& Vtx)
{
  math_Vector toler(1, 4), infb(1, 4), supb(1, 4);
  Finv.GetTolerance(toler, tolpoint3d);
  Finv.GetBounds(infb, supb);
  Solinv(1) = sol(3);
  Solinv(2) = param;
  Solinv(3) = sol(1);
  Solinv(4) = sol(2);

  math_FunctionSetRoot rsnld(Finv, toler, 30);
  rsnld.Perform(Finv, Solinv, infb, supb);
  if (!rsnld.IsDone())
  {
#ifdef OCCT_DEBUG
    std::cout << "SurfRstLineBuilder :RSNLD not done " << std::endl;
#endif
    return false;
  }
  rsnld.Root(Solinv);

  if (Finv.IsSolution(Solinv, tolpoint3d))
  {
    gp_Pnt2d     p2d(Solinv(3), Solinv(4));
    TopAbs_State situ = domain1->Classify(p2d, tolpoint2d, false);
    if ((situ != TopAbs_IN) && (situ != TopAbs_ON))
    {
      return false;
    }
    domain2->Initialize(rst);
    domain2->InitVertexIterator();
    IsVtx = !domain2->MoreVertex();
    while (!IsVtx)
    {
      Vtx = domain2->Vertex();
      if (std::abs(BRepBlend_BlendTool::Parameter(Vtx, rst) - Solinv(1))
          <= BRepBlend_BlendTool::Tolerance(Vtx, rst))
      {
        IsVtx = true;
      }
      else
      {
        domain2->NextVertex();
        IsVtx = !domain2->MoreVertex();
      }
    }
    if (!domain2->MoreVertex())
    {
      IsVtx = false;
    }
    // The section is recalculated by direct resolution, otherwise
    // incoherences between the parameter and the ground caused by yawn are returned.

    math_Vector infbound(1, 3), supbound(1, 3);
    math_Vector parinit(1, 3), tolerance(1, 3);
    Func.GetTolerance(tolerance, tolpoint3d);
    Func.GetBounds(infbound, supbound);

    math_FunctionSetRoot rsnld2(Func, tolerance, 30);
    parinit(1) = Solinv(3);
    parinit(2) = Solinv(4);
    parinit(3) = Solinv(1);
    Func.Set(Solinv(2));
    rsnld2.Perform(Func, parinit, infbound, supbound);
    if (!rsnld2.IsDone())
    {
      return false;
    }
    rsnld2.Root(parinit);
    Solinv(3) = parinit(1);
    Solinv(4) = parinit(2);
    Solinv(1) = parinit(3);
    return true;
  }
  return false;
}

//=================================================================================================

bool BRepBlend_SurfRstLineBuilder::Recadre(Blend_SurfPointFuncInv&         FinvP,
                                           math_Vector&                    Solinv,
                                           bool&                           IsVtx,
                                           occ::handle<Adaptor3d_HVertex>& Vtx)
{
  // Le point.
  gp_Pnt2d p2drst;
  double   firstrst = rst->FirstParameter();
  double   lastrst  = rst->LastParameter();
  double   wpoint   = firstrst;
  if ((sol(3) - firstrst) > (lastrst - sol(3)))
  {
    wpoint = lastrst;
  }
  p2drst          = rst->Value(wpoint);
  gp_Pnt thepoint = surf2->Value(p2drst.X(), p2drst.Y());

  FinvP.Set(thepoint);
  math_Vector toler(1, 3), infb(1, 3), supb(1, 3);
  FinvP.GetTolerance(toler, tolpoint3d);
  FinvP.GetBounds(infb, supb);
  Solinv(1) = param;
  Solinv(2) = sol(1);
  Solinv(3) = sol(2);

  math_FunctionSetRoot rsnld(FinvP, toler, 30);
  rsnld.Perform(FinvP, Solinv, infb, supb);
  if (!rsnld.IsDone())
  {
#ifdef OCCT_DEBUG
    std::cout << "SurfRstLineBuilder :RSNLD not done " << std::endl;
#endif
    return false;
  }
  rsnld.Root(Solinv);

  if (FinvP.IsSolution(Solinv, tolpoint3d))
  {
    gp_Pnt2d     p2d(Solinv(2), Solinv(3));
    TopAbs_State situ = domain1->Classify(p2d, tolpoint2d, false);
    if ((situ != TopAbs_IN) && (situ != TopAbs_ON))
    {
      return false;
    }
    domain2->Initialize(rst);
    domain2->InitVertexIterator();
    IsVtx = !domain2->MoreVertex();
    while (!IsVtx)
    {
      Vtx = domain2->Vertex();
      if (std::abs(BRepBlend_BlendTool::Parameter(Vtx, rst) - wpoint)
          <= BRepBlend_BlendTool::Tolerance(Vtx, rst))
      {
        IsVtx = true;
      }
      else
      {
        domain2->NextVertex();
        IsVtx = !domain2->MoreVertex();
      }
    }
    if (!domain2->MoreVertex())
    {
      IsVtx = false;
    }
    return true;
  }
  return false;
}

//=================================================================================================

void BRepBlend_SurfRstLineBuilder::Transition(const bool                            OnFirst,
                                              const occ::handle<Adaptor2d_Curve2d>& Arc,
                                              const double                          Param,
                                              IntSurf_Transition&                   TLine,
                                              IntSurf_Transition&                   TArc)
{
  bool        computetranstionaveclacorde = false;
  gp_Vec      tgline;
  Blend_Point prevprev;

  if (previousP.IsTangencyPoint())
  {
    if (line->NbPoints() < 2)
    {
      return;
    }
    computetranstionaveclacorde = true;
    if (sens < 0)
    {
      prevprev = line->Point(2);
    }
    else
    {
      prevprev = line->Point(line->NbPoints() - 1);
    }
  }
  gp_Pnt2d p2d;
  gp_Vec2d dp2d;

  gp_Pnt pbid;
  gp_Vec d1u, d1v, normale, tgrst;

  Arc->D1(Param, p2d, dp2d);
  if (OnFirst)
  {
    surf1->D1(p2d.X(), p2d.Y(), pbid, d1u, d1v);
    if (!computetranstionaveclacorde)
    {
      tgline = previousP.TangentOnS1();
    }
    else
    {
      tgline = gp_Vec(prevprev.PointOnS(), previousP.PointOnS());
    }
  }
  else
  {
    surf2->D1(p2d.X(), p2d.Y(), pbid, d1u, d1v);
    if (!computetranstionaveclacorde)
    {
      tgline = previousP.TangentOnS2();
    }
    else
    {
      tgline = gp_Vec(prevprev.PointOnC(), previousP.PointOnC());
    }
  }

  tgrst.SetLinearForm(dp2d.X(), d1u, dp2d.Y(), d1v);
  normale = d1u.Crossed(d1v);

  IntSurf::MakeTransition(tgline, tgrst, normale, TLine, TArc);
}

//=================================================================================================

void BRepBlend_SurfRstLineBuilder::MakeExtremity(BRepBlend_Extremity&                  Extrem,
                                                 const bool                            OnFirst,
                                                 const occ::handle<Adaptor2d_Curve2d>& Arc,
                                                 const double                          Param,
                                                 const bool                            IsVtx,
                                                 const occ::handle<Adaptor3d_HVertex>& Vtx)
{
  IntSurf_Transition               Tline, Tarc;
  double                           prm;
  occ::handle<Adaptor3d_TopolTool> Iter;
  if (OnFirst)
  {
    Extrem.SetValue(previousP.PointOnS(), sol(1), sol(2), previousP.Parameter(), tolpoint3d);
    if (!previousP.IsTangencyPoint())
    {
      Extrem.SetTangent(previousP.TangentOnS());
    }
    Iter = domain1;
  }
  else
  {
    Extrem.SetValue(previousP.PointOnC(), sol(3), previousP.Parameter(), tolpoint3d);
    if (!previousP.IsTangencyPoint())
    {
      Extrem.SetTangent(previousP.TangentOnC());
    }
    Iter = domain2;
  }

  Iter->Init();
  if (!IsVtx)
  {
    Transition(OnFirst, Arc, Param, Tline, Tarc);
    Extrem.AddArc(Arc, Param, Tline, Tarc);
  }
  else
  {
    Extrem.SetVertex(Vtx);
    while (Iter->More())
    {
      occ::handle<Adaptor2d_Curve2d> arc = Iter->Value();
      if (arc != Arc)
      {
        Iter->Initialize(arc);
        Iter->InitVertexIterator();
        while (Iter->MoreVertex())
        {
          if (Iter->Identical(Vtx, Iter->Vertex()))
          {
            prm = BRepBlend_BlendTool::Parameter(Vtx, arc);
            Transition(OnFirst, arc, prm, Tline, Tarc);
            Extrem.AddArc(arc, prm, Tline, Tarc);
          }
          Iter->NextVertex();
        }
      }
      else
      {
        Transition(OnFirst, arc, Param, Tline, Tarc);
        Extrem.AddArc(arc, Param, Tline, Tarc);
      }
      Iter->Next();
    }
  }
}

//=================================================================================================

Blend_Status BRepBlend_SurfRstLineBuilder::CheckDeflectionOnSurf(const Blend_Point& CurPoint)
{
  // Controls 3d of Blend_CSWalking.

  // rule by tests in U4 corresponds to 11.478 d
  const double CosRef3D = 0.98;
  double       Cosi = 0, Cosi2 = 0;
  bool         curpointistangent  = CurPoint.IsTangencyPoint();
  bool         prevpointistangent = previousP.IsTangencyPoint();

  gp_Pnt Psurf = CurPoint.PointOnS();
  gp_Vec Tgsurf;
  if (!curpointistangent)
  {
    Tgsurf = CurPoint.TangentOnS();
  }
  gp_Pnt prevP = previousP.PointOnS();
  gp_Vec prevTg;
  if (!prevpointistangent)
  {
    prevTg = previousP.TangentOnS();
  }
  double Norme, prevNorme = 0.;
  gp_Vec Corde(prevP, Psurf);
  Norme = Corde.SquareMagnitude();
  //  if(!curpointistangent) curNorme = Tgsurf.SquareMagnitude();
  if (!prevpointistangent)
  {
    prevNorme = prevTg.SquareMagnitude();
  }

  const double toler3d = 0.01 * tolpoint3d;
  if (Norme <= toler3d * toler3d)
  {
    // it can be necessary to force same point
    return Blend_SamePoints;
  }
  if (!prevpointistangent)
  {
    if (prevNorme <= toler3d * toler3d)
    {
      return Blend_SamePoints;
    }
    Cosi = sens * Corde * prevTg;
    if (Cosi < 0.)
    { // angle 3d>pi/2. --> return back
      return Blend_Backward;
    }

    Cosi2 = Cosi * Cosi / prevNorme / Norme;
    if (Cosi2 < CosRef3D)
    {
      return Blend_StepTooLarge;
    }
  }

  if (!curpointistangent)
  {
    // Check if it is necessary to control the sign of prevtg*Tgsurf
    Cosi  = sens * Corde * Tgsurf;
    Cosi2 = Cosi * Cosi / Tgsurf.SquareMagnitude() / Norme;
    if (Cosi2 < CosRef3D || Cosi < 0.)
    {
      return Blend_StepTooLarge;
    }
  }

  if (!curpointistangent && !prevpointistangent)
  {
    // Estimation of the current arrow
    double FlecheCourante =
      (prevTg.Normalized().XYZ() - Tgsurf.Normalized().XYZ()).SquareModulus() * Norme / 64.;

    if (FlecheCourante <= 0.25 * fleche * fleche)
    {
      return Blend_StepTooSmall;
    }
    if (FlecheCourante > fleche * fleche)
    {
      // not too great :
      return Blend_StepTooLarge;
    }
  }
  return Blend_OK;
}

//=================================================================================================

Blend_Status BRepBlend_SurfRstLineBuilder::CheckDeflectionOnRst(const Blend_Point& CurPoint)
{
  // Controls 3D of Blend_CSWalking.

  // rule by tests in U4 corresponds to 11.478 d
  const double CosRef3D = 0.98;
  double       Cosi, Cosi2;
  bool         curpointistangent  = CurPoint.IsTangencyPoint();
  bool         prevpointistangent = previousP.IsTangencyPoint();

  gp_Pnt Psurf = CurPoint.PointOnC();
  gp_Vec Tgsurf;
  if (!curpointistangent)
  {
    Tgsurf = CurPoint.TangentOnC();
  }
  gp_Pnt prevP = previousP.PointOnC();
  gp_Vec prevTg;
  if (!prevpointistangent)
  {
    prevTg = previousP.TangentOnC();
  }
  double Norme, prevNorme = 0.;
  gp_Vec Corde(prevP, Psurf);
  Norme = Corde.SquareMagnitude();
  //  if(!curpointistangent) curNorme = Tgsurf.SquareMagnitude();
  if (!prevpointistangent)
  {
    prevNorme = prevTg.SquareMagnitude();
  }

  const double toler3d = 0.01 * tolpoint3d;
  if (Norme <= toler3d * toler3d)
  {
    // it can be necessary to force same point
    return Blend_SamePoints;
  }
  if (!prevpointistangent)
  {
    if (prevNorme <= toler3d * toler3d)
    {
      return Blend_SamePoints;
    }
    Cosi = sens * Corde * prevTg;
    if (Cosi < 0.)
    { // angle 3d>pi/2. --> return back
      return Blend_Backward;
    }

    Cosi2 = Cosi * Cosi / prevNorme / Norme;
    if (Cosi2 < CosRef3D)
    {
      return Blend_StepTooLarge;
    }
  }

  if (!curpointistangent)
  {
    // Check if it is necessary to control the sign of prevtg*Tgsurf
    Cosi  = sens * Corde * Tgsurf;
    Cosi2 = Cosi * Cosi / Tgsurf.SquareMagnitude() / Norme;
    if (Cosi2 < CosRef3D || Cosi < 0.)
    {
      return Blend_StepTooLarge;
    }
  }

  if (!curpointistangent && !prevpointistangent)
  {
    // Estimation of the current arrow
    double FlecheCourante =
      (prevTg.Normalized().XYZ() - Tgsurf.Normalized().XYZ()).SquareModulus() * Norme / 64.;

    if (FlecheCourante <= 0.25 * fleche * fleche)
    {
      return Blend_StepTooSmall;
    }
    if (FlecheCourante > fleche * fleche)
    {
      // not too great
      return Blend_StepTooLarge;
    }
  }
  return Blend_OK;
}

static IntSurf_TypeTrans ConvOrToTra(const TopAbs_Orientation O)
{
  if (O == TopAbs_FORWARD)
  {
    return IntSurf_In;
  }
  return IntSurf_Out;
}

//=================================================================================================

Blend_Status BRepBlend_SurfRstLineBuilder::TestArret(Blend_SurfRstFunction& Func,
                                                     const bool             TestDeflection,
                                                     const Blend_Status     State)
{
  gp_Pnt            pts, ptrst;
  gp_Pnt2d          pt2drst;
  gp_Vec            tgs, tgrst;
  gp_Vec2d          tg2ds, tg2drst;
  Blend_Status      StateS, StateRst;
  IntSurf_TypeTrans tras = IntSurf_Undecided, trarst = IntSurf_Undecided;
  Blend_Point       curpoint;

  if (Func.IsSolution(sol, tolpoint3d))
  {
    bool curpointistangent = Func.IsTangencyPoint();
    pts                    = Func.PointOnS();
    ptrst                  = Func.PointOnRst();
    pt2drst                = Func.Pnt2dOnRst();
    if (curpointistangent)
    {
      curpoint.SetValue(pts, ptrst, param, sol(1), sol(2), pt2drst.X(), pt2drst.Y(), sol(3));
    }
    else
    {
      tgs     = Func.TangentOnS();
      tgrst   = Func.TangentOnRst();
      tg2ds   = Func.Tangent2dOnS();
      tg2drst = Func.Tangent2dOnRst();

      curpoint.SetValue(pts,
                        ptrst,
                        param,
                        sol(1),
                        sol(2),
                        pt2drst.X(),
                        pt2drst.Y(),
                        sol(3),
                        tgs,
                        tgrst,
                        tg2ds,
                        tg2drst);
    }
    if (TestDeflection)
    {
      StateS   = CheckDeflectionOnSurf(curpoint);
      StateRst = CheckDeflectionOnRst(curpoint);
    }
    else
    {
      StateS = StateRst = Blend_OK;
    }
    if (StateS == Blend_Backward)
    {
      StateS = Blend_StepTooLarge;
      rebrou = true;
    }
    if (StateRst == Blend_Backward)
    {
      StateRst = Blend_StepTooLarge;
      rebrou   = true;
    }
    if (StateS == Blend_StepTooLarge || StateRst == Blend_StepTooLarge)
    {
      return Blend_StepTooLarge;
    }

    if (!comptra && !curpointistangent)
    {
      gp_Vec tgsecs, nors;
      Func.Decroch(sol, nors, tgsecs);
      nors.Normalize();
      double testra = tgsecs.Dot(nors.Crossed(tgs));
      if (std::abs(testra) > tolpoint3d)
      {
        if (testra < 0.)
        {
          tras = IntSurf_In;
        }
        else if (testra > 0.)
        {
          tras = IntSurf_Out;
        }
        gp_Pnt2d p2drstref;
        gp_Vec2d tg2drstref;
        rst->D1(sol(3), p2drstref, tg2drstref);
        testra                = tg2drst.Dot(tg2drstref);
        TopAbs_Orientation Or = domain2->Orientation(rst);
        if (std::abs(testra) > 1.e-8)
        {
          if (testra < 0.)
          {
            trarst = ConvOrToTra(TopAbs::Reverse(Or));
          }
          else if (testra > 0.)
          {
            trarst = ConvOrToTra(Or);
          }
          comptra = true;
          line->Set(tras, trarst);
        }
      }
    }
    if (StateS == Blend_OK || StateRst == Blend_OK)
    {
      previousP = curpoint;
      return State;
    }
    if (StateS == Blend_StepTooSmall && StateRst == Blend_StepTooSmall)
    {
      previousP = curpoint;
      if (State == Blend_OK)
      {
        return Blend_StepTooSmall;
      }
      else
      {
        return State;
      }
    }
    if (State == Blend_OK)
    {
      return Blend_SamePoints;
    }
    else
    {
      return State;
    }
  }
  return Blend_StepTooLarge;
}

//=================================================================================================

bool BRepBlend_SurfRstLineBuilder::CheckInside(Blend_SurfRstFunction& Func,
                                               TopAbs_State&          SituOnC,
                                               TopAbs_State&          SituOnS,
                                               bool&                  Decroch)
{
  math_Vector tolerance(1, 3);
  Func.GetTolerance(tolerance, tolpoint3d);
  // face pcurve.
  double w = sol(3);
  if (w < rst->FirstParameter() - tolerance(3) || w > rst->LastParameter() + tolerance(3))
  {
    SituOnC = TopAbs_OUT;
  }
  else if (w > rst->FirstParameter() && w < rst->LastParameter())
  {
    SituOnC = TopAbs_IN;
  }
  else
  {
    SituOnC = TopAbs_ON;
  }

  // face surface
  gp_Pnt2d p2d(sol(1), sol(2));
  SituOnS = domain1->Classify(p2d, tolpoint2d, false);

  // lost contact
  gp_Vec tgs, nors;
  Decroch = Func.Decroch(sol, tgs, nors);

  return (SituOnC == TopAbs_IN && SituOnS == TopAbs_IN && !Decroch);
}