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