projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031642: Visualization - crash in Graphic3d_Structure::SetVisual() on redisplaying...
[occt.git] / src / Extrema / Extrema_GenExtSS.cxx
CommitLineData
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.
30class Extrema_FuncDistSS : public math_MultipleVarFunctionWithGradient
31{
32public:
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
91protected:
92
93private:
94
95 const Adaptor3d_Surface *myS1;
96 const Adaptor3d_Surface *myS2;
97};
98
7fd59977 99//=======================================================================
100//function : Extrema_GenExtSS
101//purpose :
102//=======================================================================
b311480e 103Extrema_GenExtSS::Extrema_GenExtSS()
7fd59977 104{
105 myDone = Standard_False;
106 myInit = Standard_False;
107}
108
109//=======================================================================
110//function : Extrema_GenExtSS
111//purpose :
112//=======================================================================
113
114Extrema_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
130Extrema_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
154void 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
171void 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
220void 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
235void 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/*
283b- 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
393Standard_Boolean Extrema_GenExtSS::IsDone() const
394{
395 return myDone;
396}
397
398//=======================================================================
399//function : NbExt
400//purpose :
401//=======================================================================
402
403Standard_Integer Extrema_GenExtSS::NbExt() const
404{
9775fa61 405 if (!IsDone()) { throw StdFail_NotDone(); }
7fd59977 406 return myF.NbExt();
407
408}
409
410//=======================================================================
638ad7f3 411//function : SquareDistance
7fd59977 412//purpose :
413//=======================================================================
414
415Standard_Real Extrema_GenExtSS::SquareDistance(const Standard_Integer N) const
638ad7f3 416{
417 if (N < 1 || N > NbExt())
418 {
419 throw Standard_OutOfRange();
420 }
421
7fd59977 422 return myF.SquareDistance(N);
423}
424
425//=======================================================================
426//function : PointOnS1
427//purpose :
428//=======================================================================
429
5d99f2c8 430const Extrema_POnSurf& Extrema_GenExtSS::PointOnS1(const Standard_Integer N) const
7fd59977 431{
638ad7f3 432 if (N < 1 || N > NbExt())
433 {
434 throw Standard_OutOfRange();
435 }
436
7fd59977 437 return myF.PointOnS1(N);
438}
439
440//=======================================================================
441//function : PointOnS2
442//purpose :
443//=======================================================================
444
5d99f2c8 445const Extrema_POnSurf& Extrema_GenExtSS::PointOnS2(const Standard_Integer N) const
7fd59977 446{
638ad7f3 447 if (N < 1 || N > NbExt())
448 {
449 throw Standard_OutOfRange();
450 }
451
7fd59977 452 return myF.PointOnS2(N);
453}
454
455//=======================================================================
456//function : Bidon
457//purpose :
458//=======================================================================
459
460Adaptor3d_SurfacePtr Extrema_GenExtSS::Bidon() const
461{
462 return (Adaptor3d_SurfacePtr)0L;
463}
464
Open CASCADE Technology repository