0027105: Make code ISO-compliant [-Wpedantic fixes]
[occt.git] / src / GeomLib / GeomLib_CheckCurveOnSurface.cxx
CommitLineData
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>
36
37class GeomLib_CheckCurveOnSurface_TargetFunc;
38
39static
40Standard_Boolean MinComputing(
41 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
42 const Standard_Real theEpsilon, //1.0e-3
43 const Standard_Integer theNbParticles,
44 Standard_Real& theBestValue,
45 Standard_Real& theBestParameter);
46
47static Standard_Integer FillSubIntervals( const Handle(Geom_Curve)& theCurve3d,
48 const Handle(Geom2d_Curve)& theCurve2d,
49 const Standard_Real theFirst,
50 const Standard_Real theLast,
51 Standard_Integer &theNbParticles,
52 TColStd_Array1OfReal* const theSubIntervals = 0);
53
54//=======================================================================
55//class : GeomLib_CheckCurveOnSurface_TargetFunc
56//purpose : Target function (to be minimized)
57//=======================================================================
58class GeomLib_CheckCurveOnSurface_TargetFunc :
59 public math_MultipleVarFunctionWithHessian
60{
61 public:
62 GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D,
63 const Adaptor3d_Curve& theAdCS,
64 const Standard_Real theFirst,
65 const Standard_Real theLast):
66 myCurve1(theC3D),
67 myCurve2(theAdCS),
68 myFirst(theFirst),
69 myLast(theLast)
70 {
71 }
72
73 //returns the number of parameters of the function
74 //(the function is one-dimension).
75 virtual Standard_Integer NbVariables() const {
76 return 1;
77 }
78
79 //returns value of the function when parameters are equal to theX
80 virtual Standard_Boolean Value(const math_Vector& theX,
81 Standard_Real& theFVal)
82 {
83 return Value(theX(1), theFVal);
84 }
85
86 //returns value of the one-dimension-function when parameter
87 //is equal to theX
88 Standard_Boolean Value( const Standard_Real theX,
89 Standard_Real& theFVal) const
90 {
91 try
92 {
93 OCC_CATCH_SIGNALS
94 if (!CheckParameter(theX))
95 return Standard_False;
96
97 const gp_Pnt aP1(myCurve1.Value(theX)),
98 aP2(myCurve2.Value(theX));
99
100 theFVal = -1.0*aP1.SquareDistance(aP2);
101 }
102 catch(Standard_Failure) {
103 return Standard_False;
104 }
105 //
106 return Standard_True;
107 }
108
109 //see analogical method for abstract owner class math_MultipleVarFunction
110 virtual Standard_Integer GetStateNumber()
111 {
112 return 0;
113 }
114
115 //returns the gradient of the function when parameters are
116 //equal to theX
117 virtual Standard_Boolean Gradient(const math_Vector& theX,
118 math_Vector& theGrad)
119 {
120 return Derive(theX(1), theGrad(1));
121 }
122
123 //returns 1st derivative of the the one-dimension-function when
124 //parameter is equal to theX
125 Standard_Boolean Derive(const Standard_Real theX, Standard_Real& theDeriv) const
126 {
127 try
128 {
129 OCC_CATCH_SIGNALS
130 if (!CheckParameter(theX))
131 {
132 return Standard_False;
133 }
134 //
135 gp_Pnt aP1, aP2;
136 gp_Vec aDC1, aDC2;
137 //
138 myCurve1.D1(theX, aP1, aDC1);
139 myCurve2.D1(theX, aP2, aDC2);
140
141 const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
142 //
143 theDeriv = -2.0*aVec1.Dot(aVec2);
144 }
145 catch(Standard_Failure)
146 {
147 return Standard_False;
148 }
149
150 return Standard_True;
151 }
152
153 //returns value and gradient
154 virtual Standard_Boolean Values(const math_Vector& theX,
155 Standard_Real& theVal,
156 math_Vector& theGrad)
157 {
158 if (!Value(theX, theVal))
159 {
160 return Standard_False;
161 }
162 //
163 if (!Gradient(theX, theGrad)) {
164 return Standard_False;
165 }
166 //
167 return Standard_True;
168 }
169
170 //returns value, gradient and hessian
171 virtual Standard_Boolean Values(const math_Vector& theX,
172 Standard_Real& theVal,
173 math_Vector& theGrad,
174 math_Matrix& theHessian)
175 {
176 if (!Value(theX, theVal))
177 {
178 return Standard_False;
179 }
180 //
181 if (!Gradient(theX, theGrad))
182 {
183 return Standard_False;
184 }
185 //
186 theHessian(1,1) = theGrad(1);
187 //
188 return Standard_True;
189 }
190 //
191 Standard_Real FirstParameter() const
192 {
193 return myFirst;
194 }
195
196 //
197 Standard_Real LastParameter() const
198 {
199 return myLast;
200 }
201
202 private:
203 GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
204
205 //checks if the function can be computed when its parameter is
206 //equal to theParam
207 Standard_Boolean CheckParameter(const Standard_Real theParam) const
208 {
209 return ((myFirst <= theParam) && (theParam <= myLast));
210 }
211
212 const Adaptor3d_Curve& myCurve1;
213 const Adaptor3d_Curve& myCurve2;
214 const Standard_Real myFirst;
215 const Standard_Real myLast;
216};
217
218//=======================================================================
219//class : GeomLib_CheckCurveOnSurface_Local
220//purpose : Created for parallelization possibility only
221//=======================================================================
222class GeomLib_CheckCurveOnSurface_Local
223{
224public:
225 GeomLib_CheckCurveOnSurface_Local(
226 const Handle(Geom_Curve)& theCurve3D,
227 const Handle(Geom2d_Curve)& theCurve2D,
228 const Handle(Geom_Surface)& theSurface,
229 const TColStd_Array1OfReal& theIntervalsArr,
230 const Standard_Real theEpsilonRange,
231 const Standard_Integer theNbParticles):
232 myCurve3D(theCurve3D),
233 myCurve2D(theCurve2D),
234 mySurface(theSurface),
235 mySubIntervals(theIntervalsArr),
236 myEpsilonRange(theEpsilonRange),
237 myNbParticles(theNbParticles),
238 myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1),
239 myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1)
240 {
241 }
242
243 void operator()(const Standard_Integer& theIndex) const
244 {
245 //For every sub-interval (which is set by mySubIntervals array) this method
246 //computes optimal value of GeomLib_CheckCurveOnSurface_TargetFunc function.
247 //This optimal value will be put in corresponding (depending on theIndex - the
248 //identificator of the current interval in mySubIntervals array) cell of
249 //myArrOfDist and myArrOfParam arrays.
250 const GeomAdaptor_Curve anAC(myCurve3D);
251 const Handle(Adaptor2d_HCurve2d) anAd2dC = new Geom2dAdaptor_GHCurve(myCurve2D);
252 const Handle(Adaptor3d_HSurface) anAdS = new GeomAdaptor_HSurface(mySurface);
253
254 const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
255
256 GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
257 mySubIntervals.Value(theIndex),
258 mySubIntervals.Value(theIndex+1));
259
260 Standard_Real aMinDist = RealLast(), aPar = 0.0;
261 if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
262 {
263 myArrOfDist(theIndex) = RealLast();
264 myArrOfParam(theIndex) = aFunc.FirstParameter();
265 return;
266 }
267
268 myArrOfDist(theIndex) = aMinDist;
269 myArrOfParam(theIndex) = aPar;
270 }
271
272 //Returns optimal value (inverse of square of maximal distance)
273 void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
274 {
275 //This method looks for the minimal value of myArrOfDist.
276
277 const Standard_Integer aStartInd = myArrOfDist.Lower();
278 theMinimalValue = myArrOfDist(aStartInd);
279 theParameter = myArrOfParam(aStartInd);
280 for(Standard_Integer i = aStartInd + 1; i <= myArrOfDist.Upper(); i++)
281 {
282 if(myArrOfDist(i) < theMinimalValue)
283 {
284 theMinimalValue = myArrOfDist(i);
285 theParameter = myArrOfParam(i);
286 }
287 }
288 }
289
290private:
291 GeomLib_CheckCurveOnSurface_Local operator=(GeomLib_CheckCurveOnSurface_Local&);
292 const Handle(Geom_Curve)& myCurve3D;
293 const Handle(Geom2d_Curve)& myCurve2D;
294 const Handle(Geom_Surface)& mySurface;
295
296 const TColStd_Array1OfReal& mySubIntervals;
297 const Standard_Real myEpsilonRange;
298 const Standard_Integer myNbParticles;
299 mutable NCollection_Array1<Standard_Real> myArrOfDist;
300 mutable NCollection_Array1<Standard_Real> myArrOfParam;
301};
302
303//=======================================================================
304//function : GeomLib_CheckCurveOnSurface
305//purpose :
306//=======================================================================
307GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
308:
309 myFirst(0.),
310 myLast(0.),
311 myErrorStatus(0),
312 myMaxDistance(RealLast()),
313 myMaxParameter(0.)
314{
315}
316
317//=======================================================================
318//function : GeomLib_CheckCurveOnSurface
319//purpose :
320//=======================================================================
321GeomLib_CheckCurveOnSurface::
322 GeomLib_CheckCurveOnSurface(const Handle(Geom_Curve)& theCurve,
323 const Handle(Geom_Surface)& theSurface,
324 const Standard_Real theFirst,
325 const Standard_Real theLast):
326 myCurve(theCurve),
327 mySurface(theSurface),
328 myFirst(theFirst),
329 myLast(theLast),
330 myErrorStatus(0),
331 myMaxDistance(RealLast()),
332 myMaxParameter(0.)
333{
334}
335
336//=======================================================================
337//function : Init
338//purpose :
339//=======================================================================
340void GeomLib_CheckCurveOnSurface::Init()
341{
342 myCurve.Nullify();
343 mySurface.Nullify();
344 myFirst = 0.0;
345 myLast = 0.0;
346 myErrorStatus = 0;
347 myMaxDistance = RealLast();
348 myMaxParameter = 0.0;
349}
350
351//=======================================================================
352//function : Init
353//purpose :
354//=======================================================================
355void GeomLib_CheckCurveOnSurface::Init( const Handle(Geom_Curve)& theCurve,
356 const Handle(Geom_Surface)& theSurface,
357 const Standard_Real theFirst,
358 const Standard_Real theLast)
359{
360 myCurve = theCurve;
361 mySurface = theSurface;
362 myFirst = theFirst;
363 myLast = theLast;
364 myErrorStatus = 0;
365 myMaxDistance = RealLast();
366 myMaxParameter = 0.0;
367}
368
369//=======================================================================
370//function : Perform
371//purpose :
372//=======================================================================
373
374#ifndef HAVE_TBB
375//After fixing bug # 26365, this fragment should be deleted
376//(together the text "#ifdef HAVE_TBB")
377
378void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
379 const Standard_Boolean)
380{
381 const Standard_Boolean isTheMTDisabled = Standard_True;
382#else
383void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
384 const Standard_Boolean isTheMTDisabled)
385{
386#endif
387 if( myCurve.IsNull() ||
388 mySurface.IsNull() ||
389 thePCurve.IsNull())
390 {
391 myErrorStatus = 1;
392 return;
393 }
394
395 if( (myCurve->FirstParameter() > myFirst) ||
396 (myCurve->LastParameter() < myLast) ||
397 (thePCurve->FirstParameter() > myFirst) ||
398 (thePCurve->LastParameter() < myLast))
399 {
400 myErrorStatus = 2;
401 return;
402 }
403
404 const Standard_Real anEpsilonRange = 1.e-3;
405
406 Standard_Integer aNbParticles = 3;
407
408 //Polynomial function with degree n has not more than n-1 maxima and
409 //minima (degree of 1st derivative is equal to n-1 => 1st derivative has
410 //no greater than n-1 roots). Consequently, this function has
411 //maximum n monotonicity intervals. That is a good idea to try to put
412 //at least one particle in every monotonicity interval. Therefore,
413 //number of particles should be equal to n.
414
415 const Standard_Integer aNbSubIntervals =
416 FillSubIntervals( myCurve, thePCurve,
417 myFirst, myLast, aNbParticles);
418
419 if(!aNbSubIntervals)
420 {
421 myErrorStatus = 3;
422 return;
423 }
424
425 try {
426 OCC_CATCH_SIGNALS
427
428 TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
429 FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
430
431 GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
432 mySurface, anIntervals, anEpsilonRange, aNbParticles);
433
434 OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
435
436 aComp.OptimalValues(myMaxDistance, myMaxParameter);
437
438 myMaxDistance = sqrt(Abs(myMaxDistance));
439 }
440 catch (Standard_Failure) {
441 myErrorStatus = 3;
442 }
443}
444
445//=======================================================================
446// Function : FillSubIntervals
447// purpose : Divides [theFirst, theLast] interval on parts
448// in order to make searching-algorithm more precisely
449// (fills theSubIntervals array).
450// Returns number of subintervals.
451//=======================================================================
452Standard_Integer FillSubIntervals(const Handle(Geom_Curve)& theCurve3d,
453 const Handle(Geom2d_Curve)& theCurve2d,
454 const Standard_Real theFirst,
455 const Standard_Real theLast,
456 Standard_Integer &theNbParticles,
457 TColStd_Array1OfReal* const theSubIntervals)
458{
459 const Standard_Real anArrTempC[2] = {theFirst, theLast};
460 const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
461
462 theNbParticles = 3;
463 Handle(Geom2d_BSplineCurve) aBS2DCurv;
464 Handle(Geom_BSplineCurve) aBS3DCurv;
465
466 //
467 if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
468 {
469 aBS3DCurv = Handle(Geom_BSplineCurve)::
470 DownCast(Handle(Geom_TrimmedCurve)::
471 DownCast(theCurve3d)->BasisCurve());
472 }
473 else
474 {
475 aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
476 }
477
478 if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
479 {
480 aBS2DCurv = Handle(Geom2d_BSplineCurve)::
481 DownCast(Handle(Geom2d_TrimmedCurve)::
482 DownCast(theCurve2d)->BasisCurve());
483 }
484 else
485 {
486 aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
487 }
488
489 const TColStd_Array1OfReal &anArrKnots3D = !aBS3DCurv.IsNull() ?
490 aBS3DCurv->Knots() :
491 anArrTemp;
492 const TColStd_Array1OfReal &anArrKnots2D = !aBS2DCurv.IsNull() ?
493 aBS2DCurv->Knots() :
494 anArrTemp;
495
496 Standard_Integer aNbSubIntervals = 1;
497
498 try
499 {
500 OCC_CATCH_SIGNALS
501 const Standard_Integer anIndMax3D = anArrKnots3D.Upper(),
502 anIndMax2D = anArrKnots2D.Upper();
503
504 Standard_Integer anIndex3D = anArrKnots3D.Lower(),
505 anIndex2D = anArrKnots2D.Lower();
506
507 if(theSubIntervals)
508 theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
509
510 while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
511 {
512 const Standard_Real aVal3D = anArrKnots3D.Value(anIndex3D),
513 aVal2D = anArrKnots2D.Value(anIndex2D);
514 const Standard_Real aDelta = aVal3D - aVal2D;
515
516 if(aDelta < Precision::PConfusion())
517 {//aVal3D <= aVal2D
518 if((aVal3D > theFirst) && (aVal3D < theLast))
519 {
520 aNbSubIntervals++;
521
522 if(theSubIntervals)
523 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
524 }
525
526 anIndex3D++;
527
528 if(-aDelta < Precision::PConfusion())
529 {//aVal3D == aVal2D
530 anIndex2D++;
531 }
532 }
533 else
534 {//aVal2D < aVal3D
535 if((aVal2D > theFirst) && (aVal2D < theLast))
536 {
537 aNbSubIntervals++;
538
539 if(theSubIntervals)
540 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
541 }
542
543 anIndex2D++;
544 }
545 }
546
547 if(theSubIntervals)
548 theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
549
550 if(!aBS3DCurv.IsNull())
551 {
552 theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
553 }
554
555 if(!aBS2DCurv.IsNull())
556 {
557 theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
558 }
559 }
560 catch(Standard_Failure)
561 {
562#ifdef OCCT_DEBUG
563 cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
564 "FillSubIntervals(): Incorrect filling!" << endl;
565#endif
566
567 aNbSubIntervals = 0;
568 }
569
570 return aNbSubIntervals;
571}
572
573//=======================================================================
574//class : MinComputing
575//purpose : Performs computing minimal value
576//=======================================================================
577Standard_Boolean MinComputing (
578 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
579 const Standard_Real theEpsilon, //1.0e-3
580 const Standard_Integer theNbParticles,
581 Standard_Real& theBestValue,
582 Standard_Real& theBestParameter)
583{
584 try
585 {
586 OCC_CATCH_SIGNALS
587
588 //They are used for finding a position of theNbParticles worst places
589 const Standard_Integer aNbControlPoints = 3*theNbParticles;
590 //
591 math_Vector aParInf(1, 1), aParSup(1, 1), anOutputParam(1, 1), aStepPar(1,1);
592 aParInf(1) = theFunction.FirstParameter();
593 aParSup(1) = theFunction.LastParameter();
594 theBestParameter = aParInf(1);
595 theBestValue = RealLast();
596
597 const Standard_Real aDeltaParam = aParSup(1) - aParInf(1);
598 if(aDeltaParam < Precision::PConfusion())
599 return Standard_False;
600
601 aStepPar(1) = theEpsilon*aDeltaParam;
602
603 math_PSOParticlesPool aParticles(theNbParticles, 1);
604
605 const Standard_Real aStep = aDeltaParam/(aNbControlPoints-1);
606 Standard_Integer aCount = 1;
607 for(Standard_Real aPrm = aParInf(1); aCount <= aNbControlPoints; aCount++,
608 aPrm = (aCount == aNbControlPoints)? aParSup(1) : aPrm+aStep)
609 {
610 Standard_Real aVal = RealLast();
611 theFunction.Value(aPrm, aVal);
612
613 PSO_Particle* aParticle = aParticles.GetWorstParticle();
614
615 if(aVal > aParticle->BestDistance)
616 continue;
617
618 aParticle->Position[0] = aPrm;
619 aParticle->BestPosition[0] = aPrm;
620 aParticle->Distance = aVal;
621 aParticle->BestDistance = aVal;
622 }
623
624 math_PSO aPSO(&theFunction, aParInf, aParSup, aStepPar);
625 aPSO.Perform(aParticles, theNbParticles, theBestValue, anOutputParam);
626
627 //Here, anOutputParam contains parameter, which is near to optimal.
628 //It needs to be more precise. Precision is made by math_NewtonMinimum.
629 math_NewtonMinimum anA(theFunction);
630 anA.Perform(theFunction, anOutputParam);
631
632 if(!anA.IsDone())
633 {
634#ifdef OCCT_DEBUG
635 cout << "BRepLib_CheckCurveOnSurface::Compute(): No solution found!" << endl;
636#endif
637 return Standard_False;
638 }
639
640 anA.Location(anOutputParam);
641 theBestParameter = anOutputParam(1);
642 theBestValue = anA.Minimum();
643 }
644 catch(Standard_Failure)
645 {
646#ifdef OCCT_DEBUG
647 cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl;
648#endif
649 return Standard_False;
650 }
651
652 return Standard_True;
653}