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