b311480e |
1 | // Created on: 1995-07-18 |
2 | // Created by: Modelistation |
3 | // Copyright (c) 1995-1999 Matra Datavision |
973c2be1 |
4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
5 | // |
973c2be1 |
6 | // This file is part of Open CASCADE Technology software library. |
b311480e |
7 | // |
d5f74e42 |
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 |
973c2be1 |
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. |
b311480e |
13 | // |
973c2be1 |
14 | // Alternatively, this file may be used under the terms of Open CASCADE |
15 | // commercial license or contractual agreement. |
7fd59977 |
16 | |
42cf5bc1 |
17 | |
18 | #include <Adaptor3d_Surface.hxx> |
19 | #include <Extrema_GenExtSS.hxx> |
20 | #include <Extrema_POnSurf.hxx> |
6220ba10 |
21 | #include <math_BFGS.hxx> |
7fd59977 |
22 | #include <math_FunctionSetRoot.hxx> |
6220ba10 |
23 | #include <math_MultipleVarFunctionWithGradient.hxx> |
42cf5bc1 |
24 | #include <math_Vector.hxx> |
25 | #include <Standard_OutOfRange.hxx> |
26 | #include <Standard_TypeMismatch.hxx> |
27 | #include <StdFail_NotDone.hxx> |
7fd59977 |
28 | |
6220ba10 |
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 | |
7fd59977 |
99 | //======================================================================= |
100 | //function : Extrema_GenExtSS |
101 | //purpose : |
102 | //======================================================================= |
b311480e |
103 | Extrema_GenExtSS::Extrema_GenExtSS() |
7fd59977 |
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 | |
6220ba10 |
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); |
7fd59977 |
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 | |
6220ba10 |
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. |
859a47c3 |
382 | math_FunctionSetRoot SR2(myF, Tol); |
383 | SR2.Perform(myF, UV, UVinf, UVsup); |
7fd59977 |
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 | { |
9775fa61 |
405 | if (!IsDone()) { throw StdFail_NotDone(); } |
7fd59977 |
406 | return myF.NbExt(); |
407 | |
408 | } |
409 | |
410 | //======================================================================= |
411 | //function : Value |
412 | //purpose : |
413 | //======================================================================= |
414 | |
415 | Standard_Real Extrema_GenExtSS::SquareDistance(const Standard_Integer N) const |
416 | { |
9775fa61 |
417 | if (!IsDone()) { throw StdFail_NotDone(); } |
7fd59977 |
418 | return myF.SquareDistance(N); |
419 | } |
420 | |
421 | //======================================================================= |
422 | //function : PointOnS1 |
423 | //purpose : |
424 | //======================================================================= |
425 | |
5d99f2c8 |
426 | const Extrema_POnSurf& Extrema_GenExtSS::PointOnS1(const Standard_Integer N) const |
7fd59977 |
427 | { |
9775fa61 |
428 | if (!IsDone()) { throw StdFail_NotDone(); } |
7fd59977 |
429 | return myF.PointOnS1(N); |
430 | } |
431 | |
432 | //======================================================================= |
433 | //function : PointOnS2 |
434 | //purpose : |
435 | //======================================================================= |
436 | |
5d99f2c8 |
437 | const Extrema_POnSurf& Extrema_GenExtSS::PointOnS2(const Standard_Integer N) const |
7fd59977 |
438 | { |
9775fa61 |
439 | if (!IsDone()) { throw StdFail_NotDone(); } |
7fd59977 |
440 | return myF.PointOnS2(N); |
441 | } |
442 | |
443 | //======================================================================= |
444 | //function : Bidon |
445 | //purpose : |
446 | //======================================================================= |
447 | |
448 | Adaptor3d_SurfacePtr Extrema_GenExtSS::Bidon() const |
449 | { |
450 | return (Adaptor3d_SurfacePtr)0L; |
451 | } |
452 | |