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