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