5adae760 |
1 | // Created by: Nikolai BUKHALOV |
2 | // Copyright (c) 2015 OPEN CASCADE SAS |
3 | // |
4 | // This file is part of Open CASCADE Technology software library. |
5 | // |
6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
11 | // |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
14 | |
15 | #include <Adaptor2d_HCurve2d.hxx> |
16 | #include <Adaptor3d_Curve.hxx> |
17 | #include <Adaptor3d_CurveOnSurface.hxx> |
18 | #include <Adaptor3d_HSurface.hxx> |
19 | #include <Geom_BSplineCurve.hxx> |
20 | #include <Geom_TrimmedCurve.hxx> |
21 | #include <Geom2d_BSplineCurve.hxx> |
22 | #include <Geom2d_TrimmedCurve.hxx> |
23 | #include <Geom2dAdaptor_GHCurve.hxx> |
24 | #include <GeomAdaptor_Curve.hxx> |
25 | #include <GeomAdaptor_HSurface.hxx> |
26 | #include <GeomLib_CheckCurveOnSurface.hxx> |
27 | #include <gp_Pnt.hxx> |
28 | #include <math_Matrix.hxx> |
29 | #include <math_MultipleVarFunctionWithHessian.hxx> |
30 | #include <math_NewtonMinimum.hxx> |
31 | #include <math_PSO.hxx> |
32 | #include <math_PSOParticlesPool.hxx> |
33 | #include <OSD_Parallel.hxx> |
34 | #include <Standard_ErrorHandler.hxx> |
35 | #include <TColStd_Array1OfReal.hxx> |
e841c38c |
36 | #include <TColStd_HArray1OfReal.hxx> |
5adae760 |
37 | |
38 | class GeomLib_CheckCurveOnSurface_TargetFunc; |
39 | |
40 | static |
41 | Standard_Boolean MinComputing( |
42 | GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, |
43 | const Standard_Real theEpsilon, //1.0e-3 |
44 | const Standard_Integer theNbParticles, |
45 | Standard_Real& theBestValue, |
46 | Standard_Real& theBestParameter); |
47 | |
48 | static Standard_Integer FillSubIntervals( const Handle(Geom_Curve)& theCurve3d, |
49 | const Handle(Geom2d_Curve)& theCurve2d, |
50 | const Standard_Real theFirst, |
51 | const Standard_Real theLast, |
52 | Standard_Integer &theNbParticles, |
53 | TColStd_Array1OfReal* const theSubIntervals = 0); |
54 | |
55 | //======================================================================= |
56 | //class : GeomLib_CheckCurveOnSurface_TargetFunc |
57 | //purpose : Target function (to be minimized) |
58 | //======================================================================= |
59 | class GeomLib_CheckCurveOnSurface_TargetFunc : |
60 | public math_MultipleVarFunctionWithHessian |
61 | { |
62 | public: |
63 | GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D, |
64 | const Adaptor3d_Curve& theAdCS, |
65 | const Standard_Real theFirst, |
66 | const Standard_Real theLast): |
67 | myCurve1(theC3D), |
68 | myCurve2(theAdCS), |
69 | myFirst(theFirst), |
70 | myLast(theLast) |
71 | { |
72 | } |
73 | |
74 | //returns the number of parameters of the function |
75 | //(the function is one-dimension). |
76 | virtual Standard_Integer NbVariables() const { |
77 | return 1; |
78 | } |
79 | |
80 | //returns value of the function when parameters are equal to theX |
81 | virtual Standard_Boolean Value(const math_Vector& theX, |
82 | Standard_Real& theFVal) |
83 | { |
84 | return Value(theX(1), theFVal); |
85 | } |
86 | |
87 | //returns value of the one-dimension-function when parameter |
88 | //is equal to theX |
89 | Standard_Boolean Value( const Standard_Real theX, |
90 | Standard_Real& theFVal) const |
91 | { |
92 | try |
93 | { |
94 | OCC_CATCH_SIGNALS |
95 | if (!CheckParameter(theX)) |
96 | return Standard_False; |
97 | |
98 | const gp_Pnt aP1(myCurve1.Value(theX)), |
99 | aP2(myCurve2.Value(theX)); |
100 | |
101 | theFVal = -1.0*aP1.SquareDistance(aP2); |
102 | } |
a738b534 |
103 | catch(Standard_Failure const&) { |
5adae760 |
104 | return Standard_False; |
105 | } |
106 | // |
107 | return Standard_True; |
108 | } |
109 | |
110 | //see analogical method for abstract owner class math_MultipleVarFunction |
111 | virtual Standard_Integer GetStateNumber() |
112 | { |
113 | return 0; |
114 | } |
115 | |
116 | //returns the gradient of the function when parameters are |
117 | //equal to theX |
118 | virtual Standard_Boolean Gradient(const math_Vector& theX, |
119 | math_Vector& theGrad) |
120 | { |
121 | return Derive(theX(1), theGrad(1)); |
122 | } |
123 | |
124 | //returns 1st derivative of the the one-dimension-function when |
125 | //parameter is equal to theX |
8e45500e |
126 | Standard_Boolean Derive(const Standard_Real theX, Standard_Real& theDeriv1, Standard_Real* const theDeriv2 = 0) const |
5adae760 |
127 | { |
128 | try |
129 | { |
130 | OCC_CATCH_SIGNALS |
131 | if (!CheckParameter(theX)) |
132 | { |
133 | return Standard_False; |
134 | } |
135 | // |
136 | gp_Pnt aP1, aP2; |
8e45500e |
137 | gp_Vec aDC1, aDC2, aDCC1, aDCC2; |
5adae760 |
138 | // |
8e45500e |
139 | if (!theDeriv2) |
140 | { |
141 | myCurve1.D1(theX, aP1, aDC1); |
142 | myCurve2.D1(theX, aP2, aDC2); |
143 | } |
144 | else |
145 | { |
146 | myCurve1.D2(theX, aP1, aDC1, aDCC1); |
147 | myCurve2.D2(theX, aP2, aDC2, aDCC2); |
148 | } |
5adae760 |
149 | |
150 | const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1); |
151 | // |
8e45500e |
152 | theDeriv1 = -2.0*aVec1.Dot(aVec2); |
153 | |
154 | if (theDeriv2) |
155 | { |
156 | const gp_Vec aVec3(aDCC2 - aDCC1); |
157 | *theDeriv2 = -2.0*(aVec2.SquareMagnitude() + aVec1.Dot(aVec3)); |
158 | } |
5adae760 |
159 | } |
a738b534 |
160 | catch(Standard_Failure const&) |
5adae760 |
161 | { |
162 | return Standard_False; |
163 | } |
164 | |
165 | return Standard_True; |
166 | } |
167 | |
168 | //returns value and gradient |
169 | virtual Standard_Boolean Values(const math_Vector& theX, |
170 | Standard_Real& theVal, |
171 | math_Vector& theGrad) |
172 | { |
173 | if (!Value(theX, theVal)) |
174 | { |
175 | return Standard_False; |
176 | } |
177 | // |
178 | if (!Gradient(theX, theGrad)) { |
179 | return Standard_False; |
180 | } |
181 | // |
182 | return Standard_True; |
183 | } |
184 | |
185 | //returns value, gradient and hessian |
186 | virtual Standard_Boolean Values(const math_Vector& theX, |
187 | Standard_Real& theVal, |
188 | math_Vector& theGrad, |
189 | math_Matrix& theHessian) |
190 | { |
191 | if (!Value(theX, theVal)) |
192 | { |
193 | return Standard_False; |
194 | } |
195 | // |
8e45500e |
196 | if (!Derive(theX(1), theGrad(1), &theHessian(1, 1))) |
5adae760 |
197 | { |
198 | return Standard_False; |
199 | } |
200 | // |
5adae760 |
201 | return Standard_True; |
202 | } |
203 | // |
204 | Standard_Real FirstParameter() const |
205 | { |
206 | return myFirst; |
207 | } |
208 | |
209 | // |
210 | Standard_Real LastParameter() const |
211 | { |
212 | return myLast; |
213 | } |
214 | |
215 | private: |
216 | GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&); |
217 | |
218 | //checks if the function can be computed when its parameter is |
219 | //equal to theParam |
220 | Standard_Boolean CheckParameter(const Standard_Real theParam) const |
221 | { |
222 | return ((myFirst <= theParam) && (theParam <= myLast)); |
223 | } |
224 | |
225 | const Adaptor3d_Curve& myCurve1; |
226 | const Adaptor3d_Curve& myCurve2; |
227 | const Standard_Real myFirst; |
228 | const Standard_Real myLast; |
229 | }; |
230 | |
231 | //======================================================================= |
232 | //class : GeomLib_CheckCurveOnSurface_Local |
233 | //purpose : Created for parallelization possibility only |
234 | //======================================================================= |
235 | class GeomLib_CheckCurveOnSurface_Local |
236 | { |
237 | public: |
238 | GeomLib_CheckCurveOnSurface_Local( |
239 | const Handle(Geom_Curve)& theCurve3D, |
240 | const Handle(Geom2d_Curve)& theCurve2D, |
241 | const Handle(Geom_Surface)& theSurface, |
242 | const TColStd_Array1OfReal& theIntervalsArr, |
243 | const Standard_Real theEpsilonRange, |
244 | const Standard_Integer theNbParticles): |
245 | myCurve3D(theCurve3D), |
246 | myCurve2D(theCurve2D), |
247 | mySurface(theSurface), |
248 | mySubIntervals(theIntervalsArr), |
249 | myEpsilonRange(theEpsilonRange), |
250 | myNbParticles(theNbParticles), |
251 | myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1), |
252 | myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1) |
253 | { |
254 | } |
255 | |
256 | void operator()(const Standard_Integer& theIndex) const |
257 | { |
258 | //For every sub-interval (which is set by mySubIntervals array) this method |
259 | //computes optimal value of GeomLib_CheckCurveOnSurface_TargetFunc function. |
260 | //This optimal value will be put in corresponding (depending on theIndex - the |
261 | //identificator of the current interval in mySubIntervals array) cell of |
262 | //myArrOfDist and myArrOfParam arrays. |
263 | const GeomAdaptor_Curve anAC(myCurve3D); |
264 | const Handle(Adaptor2d_HCurve2d) anAd2dC = new Geom2dAdaptor_GHCurve(myCurve2D); |
265 | const Handle(Adaptor3d_HSurface) anAdS = new GeomAdaptor_HSurface(mySurface); |
266 | |
267 | const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS); |
268 | |
269 | GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS, |
270 | mySubIntervals.Value(theIndex), |
271 | mySubIntervals.Value(theIndex+1)); |
272 | |
273 | Standard_Real aMinDist = RealLast(), aPar = 0.0; |
274 | if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar)) |
275 | { |
276 | myArrOfDist(theIndex) = RealLast(); |
277 | myArrOfParam(theIndex) = aFunc.FirstParameter(); |
278 | return; |
279 | } |
280 | |
281 | myArrOfDist(theIndex) = aMinDist; |
282 | myArrOfParam(theIndex) = aPar; |
283 | } |
284 | |
285 | //Returns optimal value (inverse of square of maximal distance) |
286 | void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const |
287 | { |
288 | //This method looks for the minimal value of myArrOfDist. |
289 | |
290 | const Standard_Integer aStartInd = myArrOfDist.Lower(); |
291 | theMinimalValue = myArrOfDist(aStartInd); |
292 | theParameter = myArrOfParam(aStartInd); |
293 | for(Standard_Integer i = aStartInd + 1; i <= myArrOfDist.Upper(); i++) |
294 | { |
295 | if(myArrOfDist(i) < theMinimalValue) |
296 | { |
297 | theMinimalValue = myArrOfDist(i); |
298 | theParameter = myArrOfParam(i); |
299 | } |
300 | } |
301 | } |
302 | |
303 | private: |
304 | GeomLib_CheckCurveOnSurface_Local operator=(GeomLib_CheckCurveOnSurface_Local&); |
305 | const Handle(Geom_Curve)& myCurve3D; |
306 | const Handle(Geom2d_Curve)& myCurve2D; |
307 | const Handle(Geom_Surface)& mySurface; |
308 | |
309 | const TColStd_Array1OfReal& mySubIntervals; |
310 | const Standard_Real myEpsilonRange; |
311 | const Standard_Integer myNbParticles; |
312 | mutable NCollection_Array1<Standard_Real> myArrOfDist; |
313 | mutable NCollection_Array1<Standard_Real> myArrOfParam; |
314 | }; |
315 | |
316 | //======================================================================= |
317 | //function : GeomLib_CheckCurveOnSurface |
318 | //purpose : |
319 | //======================================================================= |
320 | GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface() |
321 | : |
322 | myFirst(0.), |
323 | myLast(0.), |
324 | myErrorStatus(0), |
325 | myMaxDistance(RealLast()), |
6ca1c746 |
326 | myMaxParameter(0.), |
327 | myTolRange(Precision::PConfusion()) |
5adae760 |
328 | { |
329 | } |
330 | |
331 | //======================================================================= |
332 | //function : GeomLib_CheckCurveOnSurface |
333 | //purpose : |
334 | //======================================================================= |
335 | GeomLib_CheckCurveOnSurface:: |
336 | GeomLib_CheckCurveOnSurface(const Handle(Geom_Curve)& theCurve, |
337 | const Handle(Geom_Surface)& theSurface, |
338 | const Standard_Real theFirst, |
6ca1c746 |
339 | const Standard_Real theLast, |
340 | const Standard_Real theTolRange): |
5adae760 |
341 | myCurve(theCurve), |
342 | mySurface(theSurface), |
343 | myFirst(theFirst), |
344 | myLast(theLast), |
345 | myErrorStatus(0), |
346 | myMaxDistance(RealLast()), |
6ca1c746 |
347 | myMaxParameter(0.), |
348 | myTolRange(theTolRange) |
5adae760 |
349 | { |
350 | } |
351 | |
352 | //======================================================================= |
353 | //function : Init |
354 | //purpose : |
355 | //======================================================================= |
356 | void GeomLib_CheckCurveOnSurface::Init() |
357 | { |
358 | myCurve.Nullify(); |
359 | mySurface.Nullify(); |
360 | myFirst = 0.0; |
361 | myLast = 0.0; |
362 | myErrorStatus = 0; |
363 | myMaxDistance = RealLast(); |
364 | myMaxParameter = 0.0; |
6ca1c746 |
365 | myTolRange = Precision::PConfusion(); |
5adae760 |
366 | } |
367 | |
368 | //======================================================================= |
369 | //function : Init |
370 | //purpose : |
371 | //======================================================================= |
372 | void GeomLib_CheckCurveOnSurface::Init( const Handle(Geom_Curve)& theCurve, |
373 | const Handle(Geom_Surface)& theSurface, |
374 | const Standard_Real theFirst, |
6ca1c746 |
375 | const Standard_Real theLast, |
376 | const Standard_Real theTolRange) |
5adae760 |
377 | { |
378 | myCurve = theCurve; |
379 | mySurface = theSurface; |
380 | myFirst = theFirst; |
381 | myLast = theLast; |
382 | myErrorStatus = 0; |
383 | myMaxDistance = RealLast(); |
384 | myMaxParameter = 0.0; |
6ca1c746 |
385 | myTolRange = theTolRange; |
5adae760 |
386 | } |
387 | |
388 | //======================================================================= |
389 | //function : Perform |
390 | //purpose : |
391 | //======================================================================= |
5adae760 |
392 | #ifndef HAVE_TBB |
393 | //After fixing bug # 26365, this fragment should be deleted |
394 | //(together the text "#ifdef HAVE_TBB") |
395 | |
396 | void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve, |
397 | const Standard_Boolean) |
398 | { |
399 | const Standard_Boolean isTheMTDisabled = Standard_True; |
400 | #else |
401 | void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve, |
402 | const Standard_Boolean isTheMTDisabled) |
403 | { |
404 | #endif |
405 | if( myCurve.IsNull() || |
406 | mySurface.IsNull() || |
407 | thePCurve.IsNull()) |
408 | { |
409 | myErrorStatus = 1; |
410 | return; |
411 | } |
412 | |
6ca1c746 |
413 | if(((myCurve->FirstParameter() - myFirst) > myTolRange) || |
414 | ((myCurve->LastParameter() - myLast) < -myTolRange) || |
415 | ((thePCurve->FirstParameter() - myFirst) > myTolRange) || |
416 | ((thePCurve->LastParameter() - myLast) < -myTolRange)) |
5adae760 |
417 | { |
418 | myErrorStatus = 2; |
419 | return; |
420 | } |
421 | |
422 | const Standard_Real anEpsilonRange = 1.e-3; |
423 | |
424 | Standard_Integer aNbParticles = 3; |
425 | |
426 | //Polynomial function with degree n has not more than n-1 maxima and |
427 | //minima (degree of 1st derivative is equal to n-1 => 1st derivative has |
428 | //no greater than n-1 roots). Consequently, this function has |
429 | //maximum n monotonicity intervals. That is a good idea to try to put |
430 | //at least one particle in every monotonicity interval. Therefore, |
431 | //number of particles should be equal to n. |
432 | |
433 | const Standard_Integer aNbSubIntervals = |
434 | FillSubIntervals( myCurve, thePCurve, |
435 | myFirst, myLast, aNbParticles); |
436 | |
437 | if(!aNbSubIntervals) |
438 | { |
439 | myErrorStatus = 3; |
440 | return; |
441 | } |
442 | |
443 | try { |
444 | OCC_CATCH_SIGNALS |
445 | |
446 | TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1); |
447 | FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals); |
448 | |
449 | GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve, |
450 | mySurface, anIntervals, anEpsilonRange, aNbParticles); |
451 | |
452 | OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled); |
453 | |
454 | aComp.OptimalValues(myMaxDistance, myMaxParameter); |
455 | |
456 | myMaxDistance = sqrt(Abs(myMaxDistance)); |
457 | } |
a738b534 |
458 | catch (Standard_Failure const&) { |
5adae760 |
459 | myErrorStatus = 3; |
460 | } |
461 | } |
462 | |
463 | //======================================================================= |
464 | // Function : FillSubIntervals |
465 | // purpose : Divides [theFirst, theLast] interval on parts |
466 | // in order to make searching-algorithm more precisely |
467 | // (fills theSubIntervals array). |
468 | // Returns number of subintervals. |
469 | //======================================================================= |
470 | Standard_Integer FillSubIntervals(const Handle(Geom_Curve)& theCurve3d, |
471 | const Handle(Geom2d_Curve)& theCurve2d, |
472 | const Standard_Real theFirst, |
473 | const Standard_Real theLast, |
474 | Standard_Integer &theNbParticles, |
475 | TColStd_Array1OfReal* const theSubIntervals) |
476 | { |
e841c38c |
477 | const Standard_Integer aMaxKnots = 101; |
5adae760 |
478 | const Standard_Real anArrTempC[2] = {theFirst, theLast}; |
479 | const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2); |
480 | |
481 | theNbParticles = 3; |
482 | Handle(Geom2d_BSplineCurve) aBS2DCurv; |
483 | Handle(Geom_BSplineCurve) aBS3DCurv; |
e841c38c |
484 | Standard_Boolean isTrimmed3D = Standard_False, isTrimmed2D = Standard_False; |
5adae760 |
485 | |
486 | // |
487 | if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) |
488 | { |
489 | aBS3DCurv = Handle(Geom_BSplineCurve):: |
490 | DownCast(Handle(Geom_TrimmedCurve):: |
491 | DownCast(theCurve3d)->BasisCurve()); |
e841c38c |
492 | isTrimmed3D = Standard_True; |
5adae760 |
493 | } |
494 | else |
495 | { |
496 | aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d); |
497 | } |
498 | |
499 | if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) |
500 | { |
501 | aBS2DCurv = Handle(Geom2d_BSplineCurve):: |
502 | DownCast(Handle(Geom2d_TrimmedCurve):: |
503 | DownCast(theCurve2d)->BasisCurve()); |
e841c38c |
504 | isTrimmed2D = Standard_True; |
5adae760 |
505 | } |
506 | else |
507 | { |
508 | aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d); |
509 | } |
510 | |
e841c38c |
511 | Handle(TColStd_HArray1OfReal) anArrKnots3D, anArrKnots2D; |
512 | |
513 | if(!aBS3DCurv.IsNull()) |
514 | { |
515 | if(aBS3DCurv->NbKnots() <= aMaxKnots) |
516 | { |
517 | anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots()); |
518 | } |
519 | else |
520 | { |
521 | Standard_Integer KnotCount; |
522 | if(isTrimmed3D) |
523 | { |
524 | Standard_Integer i; |
525 | KnotCount = 0; |
526 | const TColStd_Array1OfReal& aKnots = aBS3DCurv->Knots(); |
527 | for(i = aBS3DCurv->FirstUKnotIndex(); i <= aBS3DCurv->LastUKnotIndex(); ++i) |
528 | { |
529 | if(aKnots(i) > theFirst && aKnots(i) < theLast) |
530 | { |
531 | ++KnotCount; |
532 | } |
533 | } |
534 | KnotCount += 2; |
535 | } |
536 | else |
537 | { |
538 | KnotCount = aBS3DCurv->LastUKnotIndex() - aBS3DCurv->FirstUKnotIndex() + 1; |
539 | } |
540 | if(KnotCount <= aMaxKnots) |
541 | { |
542 | anArrKnots3D = new TColStd_HArray1OfReal(aBS3DCurv->Knots()); |
543 | } |
544 | else |
545 | { |
546 | anArrKnots3D = new TColStd_HArray1OfReal(1, aMaxKnots); |
547 | anArrKnots3D->SetValue(1, theFirst); |
548 | anArrKnots3D->SetValue(aMaxKnots, theLast); |
549 | Standard_Integer i; |
550 | Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1); |
551 | Standard_Real t = theFirst + dt; |
552 | for(i = 2; i < aMaxKnots; ++i, t += dt) |
553 | { |
554 | anArrKnots3D->SetValue(i, t); |
555 | } |
556 | } |
557 | } |
558 | } |
559 | else |
560 | { |
561 | anArrKnots3D = new TColStd_HArray1OfReal(anArrTemp); |
562 | } |
563 | if(!aBS2DCurv.IsNull()) |
564 | { |
565 | if(aBS2DCurv->NbKnots() <= aMaxKnots) |
566 | { |
567 | anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots()); |
568 | } |
569 | else |
570 | { |
571 | Standard_Integer KnotCount; |
572 | if(isTrimmed2D) |
573 | { |
574 | Standard_Integer i; |
575 | KnotCount = 0; |
576 | const TColStd_Array1OfReal& aKnots = aBS2DCurv->Knots(); |
577 | for(i = aBS2DCurv->FirstUKnotIndex(); i <= aBS2DCurv->LastUKnotIndex(); ++i) |
578 | { |
579 | if(aKnots(i) > theFirst && aKnots(i) < theLast) |
580 | { |
581 | ++KnotCount; |
582 | } |
583 | } |
584 | KnotCount += 2; |
585 | } |
586 | else |
587 | { |
588 | KnotCount = aBS2DCurv->LastUKnotIndex() - aBS2DCurv->FirstUKnotIndex() + 1; |
589 | } |
590 | if(KnotCount <= aMaxKnots) |
591 | { |
592 | anArrKnots2D = new TColStd_HArray1OfReal(aBS2DCurv->Knots()); |
593 | } |
594 | else |
595 | { |
596 | anArrKnots2D = new TColStd_HArray1OfReal(1, aMaxKnots); |
597 | anArrKnots2D->SetValue(1, theFirst); |
598 | anArrKnots2D->SetValue(aMaxKnots, theLast); |
599 | Standard_Integer i; |
600 | Standard_Real dt = (theLast - theFirst) / (aMaxKnots - 1); |
601 | Standard_Real t = theFirst + dt; |
602 | for(i = 2; i < aMaxKnots; ++i, t += dt) |
603 | { |
604 | anArrKnots2D->SetValue(i, t); |
605 | } |
606 | } |
607 | } |
608 | } |
609 | else |
610 | { |
611 | anArrKnots2D = new TColStd_HArray1OfReal(anArrTemp); |
612 | } |
613 | |
5adae760 |
614 | |
615 | Standard_Integer aNbSubIntervals = 1; |
616 | |
617 | try |
618 | { |
619 | OCC_CATCH_SIGNALS |
e841c38c |
620 | const Standard_Integer anIndMax3D = anArrKnots3D->Upper(), |
621 | anIndMax2D = anArrKnots2D->Upper(); |
5adae760 |
622 | |
e841c38c |
623 | Standard_Integer anIndex3D = anArrKnots3D->Lower(), |
624 | anIndex2D = anArrKnots2D->Lower(); |
5adae760 |
625 | |
626 | if(theSubIntervals) |
627 | theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst; |
628 | |
629 | while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D)) |
630 | { |
e841c38c |
631 | const Standard_Real aVal3D = anArrKnots3D->Value(anIndex3D), |
632 | aVal2D = anArrKnots2D->Value(anIndex2D); |
5adae760 |
633 | const Standard_Real aDelta = aVal3D - aVal2D; |
634 | |
635 | if(aDelta < Precision::PConfusion()) |
636 | {//aVal3D <= aVal2D |
637 | if((aVal3D > theFirst) && (aVal3D < theLast)) |
638 | { |
639 | aNbSubIntervals++; |
640 | |
641 | if(theSubIntervals) |
642 | theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D; |
643 | } |
644 | |
645 | anIndex3D++; |
646 | |
647 | if(-aDelta < Precision::PConfusion()) |
648 | {//aVal3D == aVal2D |
649 | anIndex2D++; |
650 | } |
651 | } |
652 | else |
653 | {//aVal2D < aVal3D |
654 | if((aVal2D > theFirst) && (aVal2D < theLast)) |
655 | { |
656 | aNbSubIntervals++; |
657 | |
658 | if(theSubIntervals) |
659 | theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D; |
660 | } |
661 | |
662 | anIndex2D++; |
663 | } |
664 | } |
665 | |
666 | if(theSubIntervals) |
667 | theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast; |
668 | |
669 | if(!aBS3DCurv.IsNull()) |
670 | { |
671 | theNbParticles = Max(theNbParticles, aBS3DCurv->Degree()); |
672 | } |
673 | |
674 | if(!aBS2DCurv.IsNull()) |
675 | { |
676 | theNbParticles = Max(theNbParticles, aBS2DCurv->Degree()); |
677 | } |
678 | } |
a738b534 |
679 | catch(Standard_Failure const&) |
5adae760 |
680 | { |
681 | #ifdef OCCT_DEBUG |
682 | cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, " |
683 | "FillSubIntervals(): Incorrect filling!" << endl; |
684 | #endif |
685 | |
686 | aNbSubIntervals = 0; |
687 | } |
688 | |
689 | return aNbSubIntervals; |
690 | } |
691 | |
243505b8 |
692 | //======================================================================= |
693 | //class : PSO_Perform |
694 | //purpose : Searches minimal distance with math_PSO class |
695 | //======================================================================= |
696 | Standard_Boolean PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, |
697 | const math_Vector &theParInf, |
698 | const math_Vector &theParSup, |
699 | const Standard_Real theEpsilon, |
700 | const Standard_Integer theNbParticles, |
701 | Standard_Real& theBestValue, |
702 | math_Vector &theOutputParam) |
703 | { |
704 | const Standard_Real aDeltaParam = theParSup(1) - theParInf(1); |
705 | if(aDeltaParam < Precision::PConfusion()) |
706 | return Standard_False; |
707 | |
708 | math_Vector aStepPar(1, 1); |
709 | aStepPar(1) = theEpsilon*aDeltaParam; |
710 | |
711 | math_PSOParticlesPool aParticles(theNbParticles, 1); |
712 | |
713 | //They are used for finding a position of theNbParticles worst places |
714 | const Standard_Integer aNbControlPoints = 3*theNbParticles; |
715 | |
716 | const Standard_Real aStep = aDeltaParam/(aNbControlPoints-1); |
717 | Standard_Integer aCount = 1; |
718 | for(Standard_Real aPrm = theParInf(1); aCount <= aNbControlPoints; aCount++, |
719 | aPrm = (aCount == aNbControlPoints)? theParSup(1) : aPrm+aStep) |
720 | { |
721 | Standard_Real aVal = RealLast(); |
722 | if(!theFunction.Value(aPrm, aVal)) |
723 | continue; |
724 | |
725 | PSO_Particle* aParticle = aParticles.GetWorstParticle(); |
726 | |
727 | if(aVal > aParticle->BestDistance) |
728 | continue; |
729 | |
730 | aParticle->Position[0] = aPrm; |
731 | aParticle->BestPosition[0] = aPrm; |
732 | aParticle->Distance = aVal; |
733 | aParticle->BestDistance = aVal; |
734 | } |
735 | |
736 | math_PSO aPSO(&theFunction, theParInf, theParSup, aStepPar); |
737 | aPSO.Perform(aParticles, theNbParticles, theBestValue, theOutputParam); |
738 | |
739 | return Standard_True; |
740 | } |
741 | |
5adae760 |
742 | //======================================================================= |
743 | //class : MinComputing |
744 | //purpose : Performs computing minimal value |
745 | //======================================================================= |
746 | Standard_Boolean MinComputing ( |
747 | GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, |
748 | const Standard_Real theEpsilon, //1.0e-3 |
749 | const Standard_Integer theNbParticles, |
750 | Standard_Real& theBestValue, |
751 | Standard_Real& theBestParameter) |
752 | { |
753 | try |
754 | { |
755 | OCC_CATCH_SIGNALS |
756 | |
5adae760 |
757 | // |
243505b8 |
758 | math_Vector aParInf(1, 1), aParSup(1, 1), anOutputParam(1, 1); |
5adae760 |
759 | aParInf(1) = theFunction.FirstParameter(); |
760 | aParSup(1) = theFunction.LastParameter(); |
761 | theBestParameter = aParInf(1); |
762 | theBestValue = RealLast(); |
763 | |
243505b8 |
764 | if(!PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles, |
765 | theBestValue, anOutputParam)) |
5adae760 |
766 | { |
243505b8 |
767 | #ifdef OCCT_DEBUG |
768 | cout << "BRepLib_CheckCurveOnSurface::Compute(): math_PSO is failed!" << endl; |
769 | #endif |
770 | return Standard_False; |
5adae760 |
771 | } |
772 | |
243505b8 |
773 | theBestParameter = anOutputParam(1); |
5adae760 |
774 | |
775 | //Here, anOutputParam contains parameter, which is near to optimal. |
776 | //It needs to be more precise. Precision is made by math_NewtonMinimum. |
243505b8 |
777 | math_NewtonMinimum aMinSol(theFunction); |
778 | aMinSol.Perform(theFunction, anOutputParam); |
779 | |
780 | if(aMinSol.IsDone() && (aMinSol.GetStatus() == math_OK)) |
781 | {//math_NewtonMinimum has precised the value. We take it. |
782 | aMinSol.Location(anOutputParam); |
783 | theBestParameter = anOutputParam(1); |
784 | theBestValue = aMinSol.Minimum(); |
785 | } |
786 | else |
787 | {//Use math_PSO again but on smaller range. |
788 | const Standard_Real aStep = theEpsilon*(aParSup(1) - aParInf(1)); |
789 | aParInf(1) = theBestParameter - 0.5*aStep; |
790 | aParSup(1) = theBestParameter + 0.5*aStep; |
791 | |
792 | Standard_Real aValue = RealLast(); |
793 | if(PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles, |
794 | aValue, anOutputParam)) |
795 | { |
796 | if(aValue < theBestValue) |
797 | { |
798 | theBestValue = aValue; |
799 | theBestParameter = anOutputParam(1); |
800 | } |
801 | } |
5adae760 |
802 | } |
5adae760 |
803 | } |
a738b534 |
804 | catch(Standard_Failure const&) |
5adae760 |
805 | { |
806 | #ifdef OCCT_DEBUG |
807 | cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl; |
808 | #endif |
809 | return Standard_False; |
810 | } |
811 | |
812 | return Standard_True; |
813 | } |