src/Extrema/Extrema_GenExtSS.cxx

   1 // Created on: 1995-07-18
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <Adaptor3d_Surface.hxx>
  19 #include <Extrema_GenExtSS.hxx>
  20 #include <Extrema_POnSurf.hxx>
  21 #include <math_BFGS.hxx>
  22 #include <math_FunctionSetRoot.hxx>
  23 #include <math_MultipleVarFunctionWithGradient.hxx>
  24 #include <math_Vector.hxx>
  25 #include <Standard_OutOfRange.hxx>
  26 #include <StdFail_NotDone.hxx>
  27
  28 //! This class represents distance objective function for surface / surface.
  29 class Extrema_FuncDistSS  : public math_MultipleVarFunctionWithGradient
  30 {
  31 public:
  32     DEFINE_STANDARD_ALLOC
  33
  34     Standard_EXPORT Extrema_FuncDistSS(const Adaptor3d_Surface& S1,
  35                                        const Adaptor3d_Surface& S2)
  36     : myS1(&S1),
  37       myS2(&S2)
  38     {
  39     }
  40
  41   Standard_EXPORT Standard_Integer NbVariables() const
  42   {
  43     return 4;
  44   }
  45
  46   Standard_EXPORT virtual Standard_Boolean Value(const math_Vector& X,Standard_Real& F)
  47   {
  48     F = myS1->Value(X(1), X(2)).SquareDistance(myS2->Value(X(3), X(4)));
  49     return true;
  50   }
  51
  52   Standard_EXPORT Standard_Boolean Gradient(const math_Vector& X,math_Vector& G)
  53   {
  54     gp_Pnt P1, P2;
  55     gp_Vec Du1s1, Dv1s1;
  56     gp_Vec Du2s2, Dv2s2;
  57     myS1->D1(X(1),X(2),P1,Du1s1,Dv1s1);
  58     myS2->D1(X(3),X(4),P2,Du2s2,Dv2s2);
  59
  60     gp_Vec P1P2 (P2,P1);
  61
  62     G(1) = P1P2.Dot(Du1s1);
  63     G(2) = P1P2.Dot(Dv1s1);
  64     G(3) = -P1P2.Dot(Du2s2);
  65     G(4) = -P1P2.Dot(Dv2s2);
  66
  67     return true;
  68   }
  69
  70   Standard_EXPORT virtual  Standard_Boolean Values(const math_Vector& X,Standard_Real& F,math_Vector& G)
  71   {
  72     F = myS1->Value(X(1), X(2)).SquareDistance(myS2->Value(X(3), X(4)));
  73
  74     gp_Pnt P1, P2;
  75     gp_Vec Du1s1, Dv1s1;
  76     gp_Vec Du2s2, Dv2s2;
  77     myS1->D1(X(1),X(2),P1,Du1s1,Dv1s1);
  78     myS2->D1(X(3),X(4),P2,Du2s2,Dv2s2);
  79
  80     gp_Vec P1P2 (P2,P1);
  81
  82     G(1) = P1P2.Dot(Du1s1);
  83     G(2) = P1P2.Dot(Dv1s1);
  84     G(3) = -P1P2.Dot(Du2s2);
  85     G(4) = -P1P2.Dot(Dv2s2);
  86
  87     return true;
  88   }
  89
  90 protected:
  91
  92 private:
  93
  94   const Adaptor3d_Surface *myS1;
  95   const Adaptor3d_Surface *myS2;
  96 };
  97
  98 //=======================================================================
  99 //function : Extrema_GenExtSS
 100 //purpose  :
 101 //=======================================================================
 102 Extrema_GenExtSS::Extrema_GenExtSS()
 103 : myu1min(0.0),
 104   myu1sup(0.0),
 105   myv1min(0.0),
 106   myv1sup(0.0),
 107   myu2min(0.0),
 108   myu2sup(0.0),
 109   myv2min(0.0),
 110   myv2sup(0.0),
 111   myusample(0),
 112   myvsample(0),
 113   mytol1(0.0),
 114   mytol2(0.0),
 115   myS2(NULL)
 116 {
 117   myDone = Standard_False;
 118   myInit = Standard_False;
 119 }
 120
 121 //=======================================================================
 122 //function : Extrema_GenExtSS
 123 //purpose  :
 124 //=======================================================================
 125
 126 Extrema_GenExtSS::Extrema_GenExtSS(const Adaptor3d_Surface& S1,
 127                                    const Adaptor3d_Surface& S2,
 128                                    const Standard_Integer NbU,
 129                                    const Standard_Integer NbV,
 130                                    const Standard_Real    Tol1,
 131                                    const Standard_Real    Tol2) : myF(S1,S2)
 132 {
 133   Initialize(S2, NbU, NbV, Tol2);
 134   Perform(S1, Tol1);
 135 }
 136
 137 //=======================================================================
 138 //function : Extrema_GenExtSS
 139 //purpose  :
 140 //=======================================================================
 141
 142 Extrema_GenExtSS::Extrema_GenExtSS(const Adaptor3d_Surface& S1,
 143                                    const Adaptor3d_Surface& S2,
 144                                    const Standard_Integer NbU,
 145                                    const Standard_Integer NbV,
 146                                    const Standard_Real U1min,
 147                                    const Standard_Real U1sup,
 148                                    const Standard_Real V1min,
 149                                    const Standard_Real V1sup,
 150                                    const Standard_Real U2min,
 151                                    const Standard_Real U2sup,
 152                                    const Standard_Real V2min,
 153                                    const Standard_Real V2sup,
 154                                    const Standard_Real Tol1,
 155                                    const Standard_Real Tol2): myF(S1, S2)
 156 {
 157   Initialize(S2, NbU, NbV, U2min,U2sup,V2min,V2sup,Tol2);
 158   Perform(S1, U1min,U1sup,V1min,V1sup,Tol1);
 159 }
 160
 161 //=======================================================================
 162 //function : Initialize
 163 //purpose  :
 164 //=======================================================================
 165
 166 void Extrema_GenExtSS::Initialize(const Adaptor3d_Surface& S2,
 167                                   const Standard_Integer NbU,
 168                                   const Standard_Integer NbV,
 169                                   const Standard_Real Tol2)
 170 {
 171   myu2min = S2.FirstUParameter();
 172   myu2sup = S2.LastUParameter();
 173   myv2min = S2.FirstVParameter();
 174   myv2sup = S2.LastVParameter();
 175   Initialize(S2,NbU,NbV,myu2min,myu2sup,myv2min,myv2sup,Tol2);
 176 }
 177
 178 //=======================================================================
 179 //function : Initialize
 180 //purpose  :
 181 //=======================================================================
 182
 183 void Extrema_GenExtSS::Initialize(const Adaptor3d_Surface& S2,
 184                                   const Standard_Integer NbU,
 185                                   const Standard_Integer NbV,
 186                                   const Standard_Real U2min,
 187                                   const Standard_Real U2sup,
 188                                   const Standard_Real V2min,
 189                                   const Standard_Real V2sup,
 190                                   const Standard_Real Tol2)
 191 {
 192   myS2 = &S2;
 193   mypoints1 = new TColgp_HArray2OfPnt(0,NbU+1,0,NbV+1);
 194   mypoints2 = new TColgp_HArray2OfPnt(0,NbU+1,0,NbV+1);
 195   myusample = NbU;
 196   myvsample = NbV;
 197   myu2min = U2min;
 198   myu2sup = U2sup;
 199   myv2min = V2min;
 200   myv2sup = V2sup;
 201   mytol2 = Tol2;
 202
 203 // Parametrage de l echantillon sur S2
 204
 205   Standard_Real PasU = myu2sup - myu2min;
 206   Standard_Real PasV = myv2sup - myv2min;
 207   Standard_Real U0 = PasU / myusample / 100.;
 208   Standard_Real V0 = PasV / myvsample / 100.;
 209   gp_Pnt P1;
 210   PasU = (PasU - U0) / (myusample - 1);
 211   PasV = (PasV - V0) / (myvsample - 1);
 212   U0 = myu2min + U0/2.;
 213   V0 = myv2min + V0/2.;
 214
 215 // Calcul des distances
 216
 217   Standard_Integer NoU, NoV;
 218   Standard_Real U, V;
 219   for ( NoU = 1, U = U0; NoU <= myusample; NoU++, U += PasU) {
 220     for ( NoV = 1, V = V0; NoV <= myvsample; NoV++, V += PasV) {
 221       P1 = myS2->Value(U, V);
 222       mypoints2->SetValue(NoU,NoV,P1);
 223     }
 224   }
 225 }
 226
 227 //=======================================================================
 228 //function : Perform
 229 //purpose  :
 230 //=======================================================================
 231
 232 void Extrema_GenExtSS::Perform(const Adaptor3d_Surface& S1,
 233                                const Standard_Real    Tol1)
 234 {
 235   myu1min = S1.FirstUParameter();
 236   myu1sup = S1.LastUParameter();
 237   myv1min = S1.FirstVParameter();
 238   myv1sup = S1.LastVParameter();
 239   Perform(S1, myu1min, myu1sup,myv1min,myv1sup,Tol1);
 240 }
 241
 242 //=======================================================================
 243 //function : Perform
 244 //purpose  :
 245 //=======================================================================
 246
 247 void Extrema_GenExtSS::Perform(const Adaptor3d_Surface& S1,
 248                                const Standard_Real U1min,
 249                                const Standard_Real U1sup,
 250                                const Standard_Real V1min,
 251                                const Standard_Real V1sup,
 252                                const Standard_Real Tol1)
 253 {
 254   myF.Initialize(S1,*myS2);
 255   myu1min = U1min;
 256   myu1sup = U1sup;
 257   myv1min = V1min;
 258   myv1sup = V1sup;
 259   mytol1 = Tol1;
 260
 261   Standard_Real U1, V1, U2, V2;
 262   Standard_Integer NoU1, NoV1, NoU2, NoV2;
 263   gp_Pnt P1, P2;
 264
 265 // Parametrage de l echantillon sur S1
 266
 267   Standard_Real PasU1 = myu1sup - myu1min;
 268   Standard_Real PasV1 = myv1sup - myv1min;
 269   Standard_Real U10 = PasU1 / myusample / 100.;
 270   Standard_Real V10 = PasV1 / myvsample / 100.;
 271   PasU1 = (PasU1 - U10) / (myusample - 1);
 272   PasV1 = (PasV1 - V10) / (myvsample - 1);
 273   U10 = myu1min + U10/2.;
 274   V10 = myv1min + V10/2.;
 275
 276   Standard_Real PasU2 = myu2sup - myu2min;
 277   Standard_Real PasV2 = myv2sup - myv2min;
 278   Standard_Real U20 = PasU2 / myusample / 100.;
 279   Standard_Real V20 = PasV2 / myvsample / 100.;
 280   PasU2 = (PasU2 - U20) / (myusample - 1);
 281   PasV2 = (PasV2 - V20) / (myvsample - 1);
 282   U20 = myu2min + U20/2.;
 283   V20 = myv2min + V20/2.;
 284
 285 // Calcul des distances
 286
 287   for ( NoU1 = 1, U1 = U10; NoU1 <= myusample; NoU1++, U1 += PasU1) {
 288     for ( NoV1 = 1, V1 = V10; NoV1 <= myvsample; NoV1++, V1 += PasV1) {
 289       P1 = S1.Value(U1, V1);
 290       mypoints1->SetValue(NoU1,NoV1,P1);
 291     }
 292   }
 293
 294 /*
 295 b- Calcul des minima:
 296    -----------------
 297    b.a) Initialisations:
 298 */
 299
 300   math_Vector Tol(1, 4);
 301   Tol(1) = mytol1;
 302   Tol(2) = mytol1;
 303   Tol(3) = mytol2;
 304   Tol(4) = mytol2;
 305   math_Vector UV(1,4), UVinf(1,4), UVsup(1,4);
 306   UVinf(1) = myu1min;
 307   UVinf(2) = myv1min;
 308   UVinf(3) = myu2min;
 309   UVinf(4) = myv2min;
 310   UVsup(1) = myu1sup;
 311   UVsup(2) = myv1sup;
 312   UVsup(3) = myu2sup;
 313   UVsup(4) = myv2sup;
 314
 315
 316   Standard_Real distmin = RealLast(), distmax = 0.0, TheDist;
 317
 318   Standard_Integer N1Umin=0,N1Vmin=0,N2Umin=0,N2Vmin=0;
 319   gp_Pnt PP1min, PP2min;
 320   Standard_Integer N1Umax=0,N1Vmax=0,N2Umax=0,N2Vmax=0;
 321   gp_Pnt PP1max, PP2max;
 322
 323   for ( NoU1 = 1, U1 = U10; NoU1 <= myusample; NoU1++, U1 += PasU1) {
 324     for ( NoV1 = 1, V1 = V10; NoV1 <= myvsample; NoV1++, V1 += PasV1) {
 325       P1 = mypoints1->Value(NoU1, NoV1);
 326       for ( NoU2 = 1, U2 = U20; NoU2 <= myusample; NoU2++, U2 += PasU2) {
 327         for ( NoV2 = 1, V2 = V20; NoV2 <= myvsample; NoV2++, V2 += PasV2) {
 328           P2 = mypoints2->Value(NoU2, NoV2);
 329           TheDist = P1.SquareDistance(P2);
 330           if (TheDist < distmin) {
 331             distmin = TheDist;
 332             N1Umin = NoU1;
 333             N1Vmin = NoV1;
 334             N2Umin = NoU2;
 335             N2Vmin = NoV2;
 336             PP1min = P1;
 337             PP2min = P2;
 338           }
 339           if (TheDist > distmax) {
 340             distmax = TheDist;
 341             N1Umax = NoU1;
 342             N1Vmax = NoV1;
 343             N2Umax = NoU2;
 344             N2Vmax = NoV2;
 345             PP1max = P1;
 346             PP2max = P2;
 347           }
 348         }
 349       }
 350     }
 351   }
 352
 353   UV(1) = U10 + (N1Umin - 1) * PasU1;
 354   UV(2) = V10 + (N1Vmin - 1) * PasV1;
 355   UV(3) = U20 + (N2Umin - 1) * PasU2;
 356   UV(4) = V20 + (N2Vmin - 1) * PasV2;
 357
 358   Extrema_FuncDistSS aGFSS(S1, *myS2);
 359   math_BFGS aBFGSSolver(4);
 360   aBFGSSolver.Perform(aGFSS, UV);
 361   if (aBFGSSolver.IsDone())
 362   {
 363     aBFGSSolver.Location(UV);
 364
 365     //  Store result in myF.
 366     myF.Value(UV , UV);
 367     myF.GetStateNumber();
 368   }
 369   else
 370   {
 371     // If optimum is not computed successfully then compute by old approach.
 372
 373     // Restore initial point.
 374     UV(1) = U10 + (N1Umin - 1) * PasU1;
 375     UV(2) = V10 + (N1Vmin - 1) * PasV1;
 376     UV(3) = U20 + (N2Umin - 1) * PasU2;
 377     UV(4) = V20 + (N2Vmin - 1) * PasV2;
 378
 379     math_FunctionSetRoot SR1(myF, Tol);
 380     SR1.Perform(myF, UV, UVinf, UVsup);
 381   }
 382
 383   //math_FunctionSetRoot SR1(myF, Tol);
 384   //SR1.Perform(myF, UV, UVinf, UVsup);
 385
 386   UV(1) = U10 + (N1Umax - 1) * PasU1;
 387   UV(2) = V10 + (N1Vmax - 1) * PasV1;
 388   UV(3) = U20 + (N2Umax - 1) * PasU2;
 389   UV(4) = V20 + (N2Vmax - 1) * PasV2;
 390
 391   // It is impossible to compute max distance in the same manner,
 392   // since for the distance functional for max have bad definition.
 393   // So, for max computation old approach is used.
 394   math_FunctionSetRoot SR2(myF, Tol);
 395   SR2.Perform(myF, UV, UVinf, UVsup);
 396
 397   myDone = Standard_True;
 398 }
 399
 400 //=======================================================================
 401 //function : IsDone
 402 //purpose  :
 403 //=======================================================================
 404
 405 Standard_Boolean Extrema_GenExtSS::IsDone() const
 406 {
 407   return myDone;
 408 }
 409
 410 //=======================================================================
 411 //function : NbExt
 412 //purpose  :
 413 //=======================================================================
 414
 415 Standard_Integer Extrema_GenExtSS::NbExt() const
 416 {
 417   if (!IsDone()) { throw StdFail_NotDone(); }
 418   return myF.NbExt();
 419
 420 }
 421
 422 //=======================================================================
 423 //function : SquareDistance
 424 //purpose  :
 425 //=======================================================================
 426
 427 Standard_Real Extrema_GenExtSS::SquareDistance(const Standard_Integer N) const
 428 {
 429   if (N < 1 || N > NbExt())
 430   {
 431     throw Standard_OutOfRange();
 432   }
 433
 434   return myF.SquareDistance(N);
 435 }
 436
 437 //=======================================================================
 438 //function : PointOnS1
 439 //purpose  :
 440 //=======================================================================
 441
 442 const Extrema_POnSurf& Extrema_GenExtSS::PointOnS1(const Standard_Integer N) const
 443 {
 444   if (N < 1 || N > NbExt())
 445   {
 446     throw Standard_OutOfRange();
 447   }
 448
 449   return myF.PointOnS1(N);
 450 }
 451
 452 //=======================================================================
 453 //function : PointOnS2
 454 //purpose  :
 455 //=======================================================================
 456
 457 const Extrema_POnSurf& Extrema_GenExtSS::PointOnS2(const Standard_Integer N) const
 458 {
 459   if (N < 1 || N > NbExt())
 460   {
 461     throw Standard_OutOfRange();
 462   }
 463
 464   return myF.PointOnS2(N);
 465 }