0026936: Drawbacks of inlining in new type system in OCCT 7.0 -- automatic
[occt.git] / src / Bisector / Bisector_BisecCC.cxx
CommitLineData
b311480e 1// Created on: 1994-03-10
2// Created by: Yves FRICAUD
3// Copyright (c) 1994-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
42cf5bc1 17
7fd59977 18#include <Bisector.hxx>
42cf5bc1 19#include <Bisector_BisecCC.hxx>
20#include <Bisector_BisecPC.hxx>
7fd59977 21#include <Bisector_Curve.hxx>
22#include <Bisector_FunctionH.hxx>
23#include <Bisector_PointOnBis.hxx>
42cf5bc1 24#include <Bisector_PolyBis.hxx>
25#include <Geom2d_Circle.hxx>
7fd59977 26#include <Geom2d_Curve.hxx>
42cf5bc1 27#include <Geom2d_Geometry.hxx>
28#include <Geom2d_Line.hxx>
29#include <Geom2d_TrimmedCurve.hxx>
30#include <Geom2dAdaptor_Curve.hxx>
31#include <Geom2dAPI_ProjectPointOnCurve.hxx>
7fd59977 32#include <Geom2dGcc.hxx>
33#include <Geom2dGcc_Circ2d2TanRad.hxx>
34#include <Geom2dGcc_QualifiedCurve.hxx>
7fd59977 35#include <Geom2dInt_GInter.hxx>
42cf5bc1 36#include <Geom2dLProp_CLProps2d.hxx>
7fd59977 37#include <gp.hxx>
42cf5bc1 38#include <gp_Pnt2d.hxx>
39#include <gp_Trsf2d.hxx>
40#include <gp_Vec2d.hxx>
7fd59977 41#include <IntRes2d_IntersectionPoint.hxx>
42cf5bc1 42#include <math_BissecNewton.hxx>
7fd59977 43#include <math_FunctionRoot.hxx>
44#include <math_FunctionRoots.hxx>
42cf5bc1 45#include <Precision.hxx>
7fd59977 46#include <Standard_DivideByZero.hxx>
42cf5bc1 47#include <Standard_DomainError.hxx>
7fd59977 48#include <Standard_NotImplemented.hxx>
42cf5bc1 49#include <Standard_OutOfRange.hxx>
50#include <Standard_RangeError.hxx>
51#include <Standard_Type.hxx>
7fd59977 52
92efcf78 53IMPLEMENT_STANDARD_RTTIEXT(Bisector_BisecCC,Bisector_Curve)
54
91bb31f3 55static Standard_Boolean ProjOnCurve (const gp_Pnt2d& P,
873c119f 56 const Handle(Geom2d_Curve)& C,
57 Standard_Real& theParam);
7fd59977 58
59static Standard_Real Curvature (const Handle(Geom2d_Curve)& C,
873c119f 60 Standard_Real U,
61 Standard_Real Tol) ;
7fd59977 62
63static Standard_Boolean TestExtension (const Handle(Geom2d_Curve)& C1,
873c119f 64 const Handle(Geom2d_Curve)& C2,
65 const Standard_Integer Start_End);
7fd59977 66
67static Standard_Boolean DiscretPar(const Standard_Real DU,
873c119f 68 const Standard_Real EpsMin,
69 const Standard_Real EpsMax,
70 const Standard_Integer NbMin,
71 const Standard_Integer NbMax,
72 Standard_Real& Eps,
73 Standard_Integer& Nb);
7fd59977 74
75//=============================================================================
76//function :
77//purpose :
78//=============================================================================
79Bisector_BisecCC::Bisector_BisecCC()
80{
81 shiftParameter = 0;
82 isEmpty = Standard_False;
83}
84
85//=============================================================================
86//function :
87//purpose :
88//=============================================================================
89Bisector_BisecCC::Bisector_BisecCC(const Handle(Geom2d_Curve)& Cu1,
873c119f 90 const Handle(Geom2d_Curve)& Cu2,
91 const Standard_Real Side1,
92 const Standard_Real Side2,
93 const gp_Pnt2d& Origin,
94 const Standard_Real DistMax)
7fd59977 95{
96 Perform (Cu1,Cu2,Side1,Side2,Origin,DistMax);
97}
98
99//=============================================================================
100//function : Perform
101//purpose :
102//=============================================================================
103void Bisector_BisecCC::Perform(const Handle(Geom2d_Curve)& Cu1,
873c119f 104 const Handle(Geom2d_Curve)& Cu2,
105 const Standard_Real Side1,
106 const Standard_Real Side2,
107 const gp_Pnt2d& Origin,
108 const Standard_Real DistMax)
7fd59977 109{
110 isEmpty = Standard_False;
111 distMax = DistMax;
112
113 curve1 = Handle (Geom2d_Curve)::DownCast(Cu1->Copy());
114 curve2 = Handle (Geom2d_Curve)::DownCast(Cu2->Copy());
115
116 sign1 = Side1;
117 sign2 = Side2;
118
119 isConvex1 = Bisector::IsConvex(curve1,sign1);
120 isConvex2 = Bisector::IsConvex(curve2,sign2);
121
122 Standard_Real U,UC1,UC2,Dist,dU,USol;
123 gp_Pnt2d P;
124 Standard_Integer NbPnts = 21;
125 Standard_Real EpsMin = 10*Precision::Confusion();
126 Standard_Boolean YaPoly = Standard_True;
127 Standard_Boolean OriInPoly = Standard_False;
128 //---------------------------------------------
0d969553 129 // Calculate first point of the polygon.
7fd59977 130 //---------------------------------------------
91bb31f3 131 Standard_Boolean isProjDone = ProjOnCurve (Origin,curve1, U);
91bb31f3 132
133 if(!isProjDone)
134 {
135 isEmpty = Standard_True;
136 return;
137 }
7fd59977 138
873c119f 139 P = ValueByInt (U,UC1,UC2,Dist);
140 if(Dist < Precision::Confusion())
141 {
142 gp_Pnt2d aP1 = curve1->Value(UC1);
143 gp_Pnt2d aP2 = curve2->Value(UC2);
144 Standard_Real dp = (aP1.Distance(P)+aP2.Distance(P));
145 Standard_Real dorig = (aP1.Distance(Origin)+aP2.Distance(Origin));
146 if(dp < dorig)
147 {
148 isEmpty = Standard_True;
149 return;
150 }
151 }
152
7fd59977 153 if (Dist < Precision::Infinite()) {
154 //----------------------------------------------------
0d969553
Y
155 // the parameter of the origin point gives a point
156 // on the polygon.
7fd59977 157 //----------------------------------------------------
158 myPolygon.Append(Bisector_PointOnBis(UC1,UC2,U,Dist,P));
159 startIntervals.Append(U);
160 if (P.IsEqual(Origin,Precision::Confusion())) {
161 //----------------------------------------
0d969553 162 // test if the first point is the origin.
7fd59977 163 //----------------------------------------
164 OriInPoly = Standard_True;
165 }
166 }
167 else {
168 //-------------------------------------------------------
0d969553
Y
169 // The origin point is on the extension.
170 // Find the first point of the polygon by dichotomy.
7fd59977 171 //-------------------------------------------------------
172 dU = (curve1->LastParameter() - U)/(NbPnts - 1);
173 U += dU;
174 for (Standard_Integer i = 1; i <= NbPnts - 1; i++) {
175 P = ValueByInt(U,UC1,UC2,Dist);
176 if (Dist < Precision::Infinite()) {
91bb31f3 177 USol = SearchBound(U - dU,U);
178 P = ValueByInt(USol,UC1,UC2,Dist);
179 startIntervals.Append(USol);
180 myPolygon.Append(Bisector_PointOnBis(UC1,UC2,USol,Dist,P));
181 break;
7fd59977 182 }
183 U += dU;
184 }
185 }
186
c29a9290 187 if ( myPolygon.Length() != 0 ) {
7fd59977 188 SupLastParameter();
189 //----------------------------------------------
0d969553 190 // Construction of the polygon of the bissectrice.
7fd59977 191 //---------------------------------------------
192 U = FirstParameter();
193 Standard_Real DU = LastParameter() - U;
194
195 if (DU < EpsMin) {NbPnts = 3;}
196 dU = DU/(NbPnts - 1);
197
198 U += dU;
91bb31f3 199 // modified by NIZHNY-EAP Fri Jan 21 09:33:20 2000 ___BEGIN___
200 // prevent addition of the same point
7fd59977 201 gp_Pnt2d prevPnt = P;
202 for (Standard_Integer i = 1; i <= NbPnts - 1; i++) {
203 P = ValueByInt(U,UC1,UC2,Dist);
204 if (Dist < Precision::Infinite()) {
91bb31f3 205 if (P.Distance (prevPnt) > Precision::Confusion())
206 myPolygon.Append(Bisector_PointOnBis(UC1,UC2,U,Dist,P));
7fd59977 207 }
208 else {
91bb31f3 209 USol = SearchBound(U - dU,U);
210 P = ValueByInt(USol,UC1,UC2,Dist);
211 endIntervals.SetValue(1,USol);
212 if (P.Distance (prevPnt) > Precision::Confusion())
213 myPolygon.Append(Bisector_PointOnBis(UC1,UC2,USol,Dist,P));
214 break;
7fd59977 215 }
216 U += dU;
217 prevPnt=P;
91bb31f3 218 // modified by NIZHNY-EAP Fri Jan 21 09:33:24 2000 ___END___
7fd59977 219 }
220 }
221 else {
222 //----------------
0d969553 223 // Empty Polygon.
7fd59977 224 //----------------
225 YaPoly = Standard_False;
226 }
91bb31f3 227
7fd59977 228 extensionStart = Standard_False;
229 extensionEnd = Standard_False;
230 pointStart = Origin;
231
232 if (isConvex1 && isConvex2) {
233 if (YaPoly) pointEnd = myPolygon.Last().Point();
234 }
235 else {
236 //-----------------------------------------------------------------------------
0d969553
Y
237 // Extension : The curve is extended at the beginning and/or the end if
238 // - one of two curves is concave.
239 // - the curves have a common point at the beginning and/or the end
240 // - the angle of opening at the common point between two curves
c6541a0c 241 // values M_PI.
0d969553
Y
242 // the extension at the beginning is taken into account if the origin is found above.
243 // ie : the origin is not the in the polygon.
7fd59977 244 //-----------------------------------------------------------------------------
91bb31f3 245
7fd59977 246 //---------------------------------
0d969553 247 // Do the extensions exist ?
7fd59977 248 //---------------------------------
249 if (OriInPoly) {
250 extensionStart = Standard_False;
251 }
252 else {
253 extensionStart = TestExtension(curve1,curve2,1);
254 }
255 extensionEnd = TestExtension(curve1,curve2,2);
256
257 //-----------------
0d969553 258 // Calculate pointEnd.
7fd59977 259 //-----------------
260 if (extensionEnd) {
261 pointEnd = curve1->Value(curve1->LastParameter());
262 }
263 else if (YaPoly) {
264 pointEnd = myPolygon.Last().Point();
265 }
266 else {
267 ComputePointEnd();
268 }
269 //------------------------------------------------------
0d969553 270 // Update the Limits of intervals of definition.
7fd59977 271 //------------------------------------------------------
272 if (YaPoly) {
273 if (extensionStart) {
91bb31f3 274 gp_Pnt2d P1 = myPolygon.First().Point();
275 Standard_Real UFirst = startIntervals.First() - pointStart.Distance(P1);
276 startIntervals.InsertBefore(1,UFirst);
277 endIntervals .InsertBefore(1,startIntervals.Value(2));
7fd59977 278 }
279 if (extensionEnd) {
91bb31f3 280 gp_Pnt2d P1;
281 Standard_Real UFirst,ULast;
282 P1 = myPolygon.Last().Point();
283 UFirst = endIntervals.Last();
284 ULast = UFirst + pointEnd.Distance(P1);
285 startIntervals.Append(UFirst);
286 endIntervals .Append(ULast );
7fd59977 287 }
288 }
289 else {
290 //--------------------------------------------------
0d969553 291 // No polygon => the bissectrice is a segment.
7fd59977 292 //--------------------------------------------------
293 startIntervals.Append(0.);
294 endIntervals .Append(pointEnd.Distance(pointStart));
295 }
296 }
297 if (!YaPoly && !extensionStart && !extensionEnd)
298 isEmpty = Standard_True;
91bb31f3 299 // modified by NIZHNY-EAP Mon Jan 17 17:32:40 2000 ___BEGIN___
7fd59977 300 if (myPolygon.Length() <= 2)
301 isEmpty = Standard_True;
91bb31f3 302 // modified by NIZHNY-EAP Mon Jan 17 17:32:42 2000 ___END___
7fd59977 303}
304
305//=============================================================================
306//function : IsExtendAtStart
307//purpose :
308//=============================================================================
309Standard_Boolean Bisector_BisecCC::IsExtendAtStart() const
310{
311 return extensionStart;
312}
313
314//=============================================================================
315//function : IsExtendAtEnd
316//purpose :
317//=============================================================================
318Standard_Boolean Bisector_BisecCC::IsExtendAtEnd() const
319{
320 return extensionEnd;
321}
322
323//=============================================================================
324//function : IsEmpty
325//purpose :
326//=============================================================================
327Standard_Boolean Bisector_BisecCC::IsEmpty() const
328{
329 return isEmpty;
330}
331
332//=============================================================================
333//function : Reverse
334//purpose :
335//=============================================================================
336void Bisector_BisecCC::Reverse()
337{
338 Standard_NotImplemented::Raise();
339}
340
341//=============================================================================
342//function : ReversedParameter
343// purpose :
344//=============================================================================
345Standard_Real Bisector_BisecCC::ReversedParameter(const Standard_Real U) const
346{
347 return LastParameter() + FirstParameter() - U;
348}
349
350//=============================================================================
351//function : Copy
352//purpose :
353//=============================================================================
354Handle(Geom2d_Geometry) Bisector_BisecCC::Copy() const
355{
356 Handle(Geom2d_Curve) CopyCurve1
357 = Handle(Geom2d_Curve)::DownCast(curve1->Copy());
358 Handle(Geom2d_Curve) CopyCurve2
359 = Handle(Geom2d_Curve)::DownCast(curve2->Copy());
360
361 Handle(Bisector_BisecCC) C = new Bisector_BisecCC();
362
363 C -> Curve (1, CopyCurve1) ; C -> Curve (2, CopyCurve2);
364 C -> Sign (1, sign1 ) ; C -> Sign (2, sign2 );
365 C -> IsConvex (1, isConvex1) ; C -> IsConvex (2, isConvex2);
366 C -> Polygon (myPolygon);
367 C -> IsEmpty (isEmpty) ;
368 C -> DistMax (distMax) ;
369 C -> StartIntervals (startIntervals); C -> EndIntervals (endIntervals);
370 C -> ExtensionStart (extensionStart); C -> ExtensionEnd (extensionEnd);
371 C -> PointStart (pointStart) ; C -> PointEnd (pointEnd) ;
372
373 return C;
374}
375
376//=============================================================================
377//function : ChangeGuide
0d969553
Y
378//purpose : Changet of the guideline for the parameters of the bissectrice
379// ATTENTION : - This can invert the direction of parameterization.
380// - This concerns only the part of the curve
381// corresponding to the polygon.
7fd59977 382//=============================================================================
383Handle(Bisector_BisecCC) Bisector_BisecCC::ChangeGuide() const
384{
385 Handle(Bisector_BisecCC) C = new Bisector_BisecCC();
873c119f 386
7fd59977 387 C -> Curve (1, curve2) ; C -> Curve (2, curve1);
388 C -> Sign (1, sign2 ) ; C -> Sign (2, sign1 );
389 C -> IsConvex (1, isConvex2); C -> IsConvex (2, isConvex1);
873c119f 390
7fd59977 391 //-------------------------------------------------------------------------
0d969553
Y
392 // Construction of the new polygon from the initial one.
393 // inversion of PointOnBis and Calculation of new parameters on the bissectrice.
7fd59977 394 //-------------------------------------------------------------------------
395 Bisector_PolyBis Poly;
396 if (sign1 == sign2 ) {
397 //---------------------------------------------------------------
0d969553 398 // elements of the new polygon are ranked in the other direction.
7fd59977 399 //---------------------------------------------------------------
400 for (Standard_Integer i = myPolygon.Length(); i >=1; i--) {
401 Bisector_PointOnBis P = myPolygon.Value(i);
402 Bisector_PointOnBis NewP (P.ParamOnC2(), P.ParamOnC1(),
873c119f 403 P.ParamOnC2(), P.Distance (),
404 P.Point());
7fd59977 405 Poly.Append (NewP);
406 }
407 }
408 else {
409 for (Standard_Integer i = 1; i <= myPolygon.Length(); i ++) {
410 Bisector_PointOnBis P = myPolygon.Value(i);
411 Bisector_PointOnBis NewP (P.ParamOnC2(), P.ParamOnC1(),
873c119f 412 P.ParamOnC2(), P.Distance (),
413 P.Point());
7fd59977 414 Poly.Append (NewP);
415 }
416 }
417 C -> Polygon (Poly);
418 C -> FirstParameter (Poly.First().ParamOnBis());
419 C -> LastParameter (Poly.Last() .ParamOnBis());
873c119f 420
7fd59977 421 return C;
422}
423
424//=============================================================================
425//function : Transform
426//purpose :
427//=============================================================================
428void Bisector_BisecCC::Transform (const gp_Trsf2d& T)
429{
430 curve1 ->Transform(T);
431 curve2 ->Transform(T);
432 myPolygon . Transform(T);
433 pointStart. Transform(T);
434 pointEnd . Transform(T);
435}
436
437//=============================================================================
438//function : IsCN
439//purpose :
440//=============================================================================
441Standard_Boolean Bisector_BisecCC::IsCN (const Standard_Integer N) const
442{
443 return (curve1->IsCN(N+1) && curve2->IsCN(N+1));
444}
445
446//=============================================================================
447//function : FirstParameter
448//purpose :
449//=============================================================================
450Standard_Real Bisector_BisecCC::FirstParameter() const
451{
873c119f 452 return startIntervals.First();
7fd59977 453}
454
455//=============================================================================
456//function : LastParameter
457//purpose :
458//=============================================================================
459Standard_Real Bisector_BisecCC::LastParameter() const
460{
873c119f 461 return endIntervals.Last();
7fd59977 462}
463
464//=============================================================================
465//function : Continuity
466//purpose :
467//=============================================================================
468GeomAbs_Shape Bisector_BisecCC::Continuity() const
469{
470 GeomAbs_Shape Cont = curve1->Continuity();
471 switch (Cont) {
472 case GeomAbs_C1 : return GeomAbs_C0;
473 case GeomAbs_C2 : return GeomAbs_C1;
474 case GeomAbs_C3 : return GeomAbs_C2;
475 case GeomAbs_CN : return GeomAbs_CN;
7fd59977 476 default: break;
7fd59977 477 }
478 return GeomAbs_C0;
479}
480
481//=============================================================================
482//function : NbIntervals
483//purpose :
484//=============================================================================
485Standard_Integer Bisector_BisecCC::NbIntervals() const
486{
487 return startIntervals.Length();
488}
489
490//=============================================================================
491//function : IntervalFirst
492//purpose :
493//=============================================================================
494Standard_Real Bisector_BisecCC::IntervalFirst(const Standard_Integer Index) const
495{
496 return startIntervals.Value(Index);
497}
873c119f 498
7fd59977 499//=============================================================================
500//function : IntervalLast
501//purpose :
502//=============================================================================
503Standard_Real Bisector_BisecCC::IntervalLast(const Standard_Integer Index) const
504{
505 return endIntervals.Value(Index);
506}
507
508//=============================================================================
509//function : IntervalContinuity
510//purpose :
511//=============================================================================
512GeomAbs_Shape Bisector_BisecCC::IntervalContinuity() const
513{
514 GeomAbs_Shape Cont = curve1->Continuity();
515 switch (Cont) {
516 case GeomAbs_C1 : return GeomAbs_C0;
517 case GeomAbs_C2 : return GeomAbs_C1;
518 case GeomAbs_C3 : return GeomAbs_C2;
519 case GeomAbs_CN : return GeomAbs_CN;
7fd59977 520 default: break;
7fd59977 521 }
522 return GeomAbs_C0;
523}
524
525//=============================================================================
526//function : IsClosed
527//purpose :
528//=============================================================================
529Standard_Boolean Bisector_BisecCC::IsClosed() const
530{
531 if (curve1->IsClosed()) {
532 if (startIntervals.First() == curve1->FirstParameter() &&
873c119f 533 endIntervals .Last () == curve1->LastParameter () )
7fd59977 534 return Standard_True;
535 }
536 return Standard_False;
537}
538
539//=============================================================================
540//function : IsPeriodic
541//purpose :
542//=============================================================================
543Standard_Boolean Bisector_BisecCC::IsPeriodic() const
544{
545 return Standard_False;
546}
547
548
549//=============================================================================
550//function : Curvature
551//purpose :
552//=============================================================================
553static Standard_Real Curvature (const Handle(Geom2d_Curve)& C,
873c119f 554 Standard_Real U,
555 Standard_Real Tol)
7fd59977 556{
557 Standard_Real K1;
558 gp_Vec2d D1,D2;
559 gp_Pnt2d P;
560 C->D2(U,P,D1,D2);
561 Standard_Real Norm2 = D1.SquareMagnitude();;
562 if (Norm2 < Tol) {
563 K1 = 0.0;
564 }
565 else {
566 K1 = (D1^D2)/(Norm2*sqrt(Norm2));
567 }
568 return K1;
569}
570
571//=============================================================================
572//function : Value
0d969553 573//purpose : CALCULATE THE CURRENT POINT BY ITERATIVE METHOD.
7fd59977 574// ----------------------------------------------
0d969553
Y
575// Calculate the current point, the distance from the current point to
576// both curves, the parameters on each curve of the projection
577// of the current point.
7fd59977 578//
0d969553
Y
579//method : - Find start parameter by using <myPolygon>.
580// - Calculate parameter U2 on curve C2 solution of H(U,V)= 0
7fd59977 581// - P(U) = F(U,U2)
582//
0d969553 583// or :
7fd59977 584// ||P2(v0)P1(u)||**2
585// F(u,v) = P1(u) - 1/2 *----------------* N(u)
586// (N(u).P2(v0)P1(u))
587//
588// H(u,v) = (Tu.P1(u)P2(v))**2||Tv||**2 - (Tv.P1(u)P2(v))**2||Tu||**2
589//
590//=============================================================================
591gp_Pnt2d Bisector_BisecCC::ValueAndDist (const Standard_Real U,
873c119f 592 Standard_Real& U1,
593 Standard_Real& U2,
594 Standard_Real& Dist) const
7fd59977 595{
596 gp_Vec2d T;
597
598 //-----------------------------------------------
0d969553 599 // is the polygon reduced to a point or empty?
7fd59977 600 //-----------------------------------------------
601 if (myPolygon.Length() <= 1) {
602 return Extension(U,U1,U2,Dist,T);
603 }
604
605 //-----------------------------------------------
0d969553 606 // test U out of the limits of the polygon.
7fd59977 607 //-----------------------------------------------
608 if (U < myPolygon.First().ParamOnBis()) {
609 return Extension(U,U1,U2,Dist,T);
610 }
611 if (U > myPolygon.Last().ParamOnBis()) {
612 return Extension(U,U1,U2,Dist,T);
613 }
614
615 //-------------------------------------------------------
0d969553 616 // Find start parameter by using <myPolygon>.
7fd59977 617 //-------------------------------------------------------
618 Standard_Integer IntervalIndex = myPolygon.Interval(U);
619 Standard_Real UMin = myPolygon.Value(IntervalIndex ).ParamOnBis();
620 Standard_Real UMax = myPolygon.Value(IntervalIndex + 1).ParamOnBis();
621 Standard_Real VMin = myPolygon.Value(IntervalIndex ).ParamOnC2();
622 Standard_Real VMax = myPolygon.Value(IntervalIndex + 1).ParamOnC2();
623 Standard_Real Alpha,VInit;
624
625 if (Abs(UMax - UMin) < gp::Resolution()) {
626 VInit = VMin;
627 }
628 else {
629 Alpha = (U - UMin)/(UMax - UMin);
630 VInit = VMin + Alpha*(VMax - VMin);
631 }
632
633 U1 = LinkBisCurve(U);
634 Standard_Real VTemp = Min(VMin,VMax);
635 VMax = Max(VMin,VMax); VMin = VTemp;
636 Standard_Boolean Valid = Standard_True;
637 //---------------------------------------------------------------
0d969553 638 // Calculate parameter U2 on curve C2 solution of H(u,v)=0
7fd59977 639 //---------------------------------------------------------------
640 gp_Pnt2d P1;
641 gp_Vec2d T1;
873c119f 642 Standard_Real EpsH = 1.E-9;
643 Standard_Real EpsH100 = 1.E-7;
7fd59977 644 curve1->D1 (U1,P1,T1);
645 gp_Vec2d N1(T1.Y(), - T1.X());
873c119f 646
7fd59977 647 if ((VMax - VMin) < Precision::PConfusion()) {
648 U2 = VInit;
649 }
650 else {
651 Bisector_FunctionH H (curve2,P1,sign1*sign2*T1);
652 Standard_Real FInit;
653 H.Value(VInit,FInit);
654 if (Abs(FInit) < EpsH) {
655 U2 = VInit;
656 }
859a47c3 657 else
658 {
659 math_BissecNewton aNewSolution(EpsH);
660 aNewSolution.Perform(H, VMin - EpsH100, VMax + EpsH100, 10);
661
662 if (aNewSolution.IsDone())
663 {
664 U2 = aNewSolution.Root();
7fd59977 665 }
859a47c3 666 else
667 {
873c119f 668 math_FunctionRoot SolRoot (H,VInit,EpsH,VMin - EpsH100,VMax + EpsH100);
859a47c3 669
670 if (SolRoot.IsDone())
873c119f 671 U2 = SolRoot.Root();
859a47c3 672 else
673 Valid = Standard_False;
7fd59977 674 }
675 }
676 }
677
678 gp_Pnt2d PBis = pointStart;
679 //----------------
680 // P(U) = F(U1,U2)
681 //----------------
682 if (Valid) {
683 gp_Pnt2d P2 = curve2->Value(U2);
684 gp_Vec2d P2P1(P1.X() - P2.X(),P1.Y() - P2.Y());
685 Standard_Real SquareP2P1 = P2P1.SquareMagnitude();
686 Standard_Real N1P2P1 = N1.Dot(P2P1);
687
688 if (P1.IsEqual(P2,Precision::Confusion())) {
689 PBis = P1 ;
690 Dist = 0.0;
691 }
692 else if (N1P2P1*sign1 < 0) {
693 Valid = Standard_False;
694 }
695 else {
873c119f 696 PBis = P1.Translated(- (0.5*SquareP2P1/N1P2P1)*N1);
697 Dist = P1.SquareDistance(PBis);
7fd59977 698 }
699 }
700
701 //----------------------------------------------------------------
0d969553
Y
702 // If the point is not valid
703 // calculate by intersection.
7fd59977 704 //----------------------------------------------------------------
705 if (!Valid) {
706 //--------------------------------------------------------------------
0d969553
Y
707 // Construction of the bisectrice point curve and of the straight line passing
708 // by P1 and carried by the normal. curve2 is limited by VMin and VMax.
7fd59977 709 //--------------------------------------------------------------------
710 Standard_Real DMin = Precision::Infinite();
711 gp_Pnt2d P;
712 Handle(Bisector_BisecPC) BisPC
713 = new Bisector_BisecPC(curve2, P1, sign2, VMin, VMax);
714 Handle(Geom2d_Line) NorLi = new Geom2d_Line (P1,N1);
715
716 Geom2dAdaptor_Curve ABisPC(BisPC);
717 Geom2dAdaptor_Curve ANorLi(NorLi);
718 //-------------------------------------------------------------------------
719 Geom2dInt_GInter Intersect(ABisPC,ANorLi,
873c119f 720 Precision::Confusion(),Precision::Confusion());
7fd59977 721 //-------------------------------------------------------------------------
722
723 if (Intersect.IsDone() && !Intersect.IsEmpty()) {
724 for (Standard_Integer i = 1; i <= Intersect.NbPoints(); i++) {
873c119f 725 if (Intersect.Point(i).ParamOnSecond()*sign1 < Precision::PConfusion()) {
726 P = Intersect.Point(i).Value();
727 if (P.SquareDistance(P1) < DMin) {
728 DMin = P.SquareDistance(P1);
729 PBis = P;
730 U2 = BisPC->LinkBisCurve(Intersect.Point(i).ParamOnFirst());
731 Dist = DMin;
732 }
733 }
7fd59977 734 }
735 }
736 }
737 return PBis;
738}
739
740//=============================================================================
741//function : ValueByInt
0d969553
Y
742//purpose : CALCULATE THE CURRENT POINT BY INTERSECTION.
743// -------------------------------------------
744// Calculate the current point, the distance from the current point
745// to two curves, the parameters on each curve of the projection of the
746// current point.
747// the current point with parameter U is the intersection of the
748// bissectrice point curve (P1,curve2) and of the straight line
749// passing through P1 of director vector N1.
750// P1 is the current point of parameter U on curve1 and N1 the
751// normal at this point.
7fd59977 752//=============================================================================
753gp_Pnt2d Bisector_BisecCC::ValueByInt (const Standard_Real U,
873c119f 754 Standard_Real& U1,
755 Standard_Real& U2,
756 Standard_Real& Dist) const
7fd59977 757{
758 //------------------------------------------------------------------
0d969553 759 // Return point, tangent, normal on C1 at parameter U.
7fd59977 760 //-------------------------------------------------------------------
761 U1 = LinkBisCurve(U);
762
763 gp_Pnt2d P1,P2,P,PSol;
764 gp_Vec2d Tan1,Tan2;
765 curve1->D1(U1,P1,Tan1);
766 gp_Vec2d N1( Tan1.Y(), - Tan1.X());
873c119f 767
7fd59977 768 //--------------------------------------------------------------------------
0d969553 769 // test confusion of P1 with extremity of curve2.
7fd59977 770 //--------------------------------------------------------------------------
771 if (P1.Distance(curve2->Value(curve2->FirstParameter())) < Precision::Confusion()) {
772 U2 = curve2->FirstParameter();
773 curve2->D1(U2,P2,Tan2);
774 if ( isConvex1 && isConvex2 ) {
775 Dist = 0.;
776 return P1;
777 }
778 if (! Tan1.IsParallel(Tan2,Precision::Angular())) {
779 Dist = 0.;
780 return P1;
781 }
782 }
783 if (P1.Distance(curve2->Value(curve2->LastParameter())) < Precision::Confusion()) {
784 U2 = curve2->LastParameter();
785 curve2->D1(U2,P2,Tan2);
786 if ( isConvex1 && isConvex2 ) {
787 Dist = 0.;
788 return P1;
789 }
790 if (! Tan1.IsParallel(Tan2,Precision::Angular())) {
791 Dist = 0.;
792 return P1;
793 }
794 }
795
796 Standard_Boolean YaSol = Standard_False;
797 Standard_Real DMin = Precision::Infinite();
798 Standard_Real USol;
799 Standard_Real EpsMax = 1.E-6;
800 Standard_Real EpsX;
801 Standard_Real EpsH = 1.E-8;
802 Standard_Real DistPP1;
803 Standard_Integer NbSamples =20;
804 Standard_Real UFirstOnC2 = curve2->FirstParameter();
805 Standard_Real ULastOnC2 = curve2->LastParameter();
873c119f 806
7fd59977 807 if (!myPolygon.IsEmpty()){
808 if (sign1 == sign2) { ULastOnC2 = myPolygon.Last().ParamOnC2();}
809 else { UFirstOnC2 = myPolygon.Last().ParamOnC2();}
810 }
811
812 if (Abs(ULastOnC2 - UFirstOnC2) < Precision::PConfusion()/100.) {
813 Dist = Precision::Infinite();
814 return P1;
815 }
816
817 DiscretPar(Abs(ULastOnC2 - UFirstOnC2),EpsH,EpsMax,2,20,EpsX,NbSamples);
818
819 Bisector_FunctionH H (curve2,P1,sign1*sign2*Tan1);
820 math_FunctionRoots SolRoot (H,
873c119f 821 UFirstOnC2,
822 ULastOnC2 ,
823 NbSamples,
824 EpsX,EpsH,EpsH);
7fd59977 825 if (SolRoot.IsDone()) {
826 for (Standard_Integer j = 1; j <= SolRoot.NbSolutions(); j++) {
827 USol = SolRoot.Value(j);
51740958 828 gp_Pnt2d P2Curve2 = curve2->Value(USol);
829 gp_Vec2d P2P1(P1.X() - P2Curve2.X(),P1.Y() - P2Curve2.Y());
7fd59977 830 Standard_Real SquareP2P1 = P2P1.SquareMagnitude();
831 Standard_Real N1P2P1 = N1.Dot(P2P1);
832
0d969553 833 // Test if the solution is at the proper side of the curves.
7fd59977 834 if (N1P2P1*sign1 > 0 ) {
873c119f 835 P = P1.Translated(- (0.5*SquareP2P1/N1P2P1)*N1);
836 DistPP1 = P1.SquareDistance(P);
837 if (DistPP1 < DMin) {
838 DMin = DistPP1;
839 PSol = P;
840 U2 = USol;
841 YaSol = Standard_True;
842 }
7fd59977 843 }
844 }
845 }
846
873c119f 847 /*
7fd59977 848 if (!YaSol) {
873c119f 849 //--------------------------------------------------------------------
850 // Construction de la bisectrice point courbe et de la droite passant
851 // par P1 et portee par la normale.
852 //--------------------------------------------------------------------
853 Handle(Bisector_BisecPC) BisPC
854 = new Bisector_BisecPC(curve2,P1,sign2,2*distMax);
855 //-------------------------------
856 // Test si la bissectrice existe.
857 //-------------------------------
858 if (BisPC->IsEmpty()) {
859 Dist = Precision::Infinite();
860 PSol = P1;
861 return PSol;
862 }
7fd59977 863
873c119f 864 Handle(Geom2d_Line) NorLi = new Geom2d_Line (P1,N1);
865 Geom2dAdaptor_Curve NorLiAd;
866 if (sign1 < 0.) {NorLiAd.Load(NorLi,0. ,distMax);}
867 else {NorLiAd.Load(NorLi,- distMax,0. );}
868
869 //-------------------------------------------------------------------------
870 Geom2dInt_GInter Intersect(BisPC,NorLiAd,
871 Precision::Confusion(),Precision::Confusion());
872 //-------------------------------------------------------------------------
873 if (Intersect.IsDone() && !Intersect.IsEmpty()) {
874 for (Standard_Integer i = 1; i <= Intersect.NbPoints(); i++) {
875 if (Intersect.Point(i).ParamOnSecond()*sign1< Precision::PConfusion()) {
876 P = Intersect.Point(i).Value();
877 DistPP1 = P.SquareDistance(P1);
878 if (DistPP1 < DMin) {
879 DMin = DistPP1;
880 PSol = P;
881 U2 = Intersect.Point(i).ParamOnFirst();
882 YaSol = Standard_True;
883 }
884 }
7fd59977 885 }
873c119f 886 }
887 }
888 */
7fd59977 889
890 if (YaSol) {
891 Dist = DMin;
892 //--------------------------------------------------------------
0d969553 893 // Point found => Test curve distance + Angular Test
7fd59977 894 //---------------------------------------------------------------
895 P2 = curve2->Value(U2);
896 gp_Vec2d PP1(P1.X() - PSol.X(),P1.Y() - PSol.Y());
897 gp_Vec2d PP2(P2.X() - PSol.X(),P2.Y() - PSol.Y());
873c119f 898
7fd59977 899 //-----------------------------------------------
0d969553 900 // Dist = product of norms = distance at the square.
7fd59977 901 //-----------------------------------------------
902 if (PP1.Dot(PP2) > (1. - Precision::Angular())*Dist) {
903 YaSol = Standard_False;
904 }
905 else {
906 if ( !isConvex1 ) {
873c119f 907 Standard_Real K1 = Curvature(curve1,U1,Precision::Confusion());
908 if (K1 != 0.) {
909 if (Dist > 1/(K1*K1)) YaSol = Standard_False;
910 }
7fd59977 911 }
912 if (YaSol) {
873c119f 913 if ( !isConvex2 ) {
914 Standard_Real K2 = Curvature(curve2,U2,Precision::Confusion());
915 if (K2 != 0.) {
916 if (Dist > 1/(K2*K2)) YaSol = Standard_False;
917 }
918 }
7fd59977 919 }
920 }
921 }
922 if (!YaSol) {
923 Dist = Precision::Infinite();
924 PSol = P1;
925 }
926 return PSol;
927}
928
929//=============================================================================
930//function : D0
931//purpose :
932//=============================================================================
933void Bisector_BisecCC::D0(const Standard_Real U,
873c119f 934 gp_Pnt2d& P) const
7fd59977 935{
936 Standard_Real U1,U2,Dist;
937
938 P = ValueAndDist(U,U1,U2,Dist);
939}
940
941//=============================================================================
942//function : D1
943//purpose :
944//=============================================================================
945void Bisector_BisecCC::D1(const Standard_Real U,
873c119f 946 gp_Pnt2d& P,
947 gp_Vec2d& V ) const
7fd59977 948{
949 V.SetCoord(0.,0.);
950 gp_Vec2d V2,V3;
951 Values(U,1,P,V,V2,V3);
952}
953
954//=============================================================================
955//function : D2
956//purpose :
957//=============================================================================
958void Bisector_BisecCC::D2(const Standard_Real U,
873c119f 959 gp_Pnt2d& P,
960 gp_Vec2d& V1,
961 gp_Vec2d& V2) const
7fd59977 962{
963 V1.SetCoord(0.,0.);
964 V2.SetCoord(0.,0.);
965 gp_Vec2d V3;
966 Values(U,2,P,V1,V2,V3);
967}
968
969//=============================================================================
970//function : D3
971//purpose :
972//=============================================================================
973void Bisector_BisecCC::D3(const Standard_Real U,
873c119f 974 gp_Pnt2d& P,
975 gp_Vec2d& V1,
976 gp_Vec2d& V2,
977 gp_Vec2d& V3) const
7fd59977 978{
979 V1.SetCoord(0.,0.);
980 V2.SetCoord(0.,0.);
981 V3.SetCoord(0.,0.);
982 Values(U,3,P,V1,V2,V3);
983}
984
985//=============================================================================
986//function : DN
987//purpose :
988//=============================================================================
989gp_Vec2d Bisector_BisecCC::DN(const Standard_Real U,
873c119f 990 const Standard_Integer N) const
7fd59977 991{
992 gp_Pnt2d P;
993 gp_Vec2d V1(0.,0.);
994 gp_Vec2d V2(0.,0.);
995 gp_Vec2d V3(0.,0.);
996 Values (U,N,P,V1,V2,V3);
997 switch (N) {
873c119f 998 case 1 : return V1;
999 case 2 : return V2;
1000 case 3 : return V3;
1001 default: {
1002 Standard_NotImplemented::Raise();
1003 }
7fd59977 1004 }
1005 return V1;
1006}
1007
1008//=============================================================================
1009//function : Values
0d969553 1010// purpose : the curve can be described by the following equations:
7fd59977 1011//
1012// B(u) = F(u,v0)
0d969553 1013// where v0 = Phi(u) is given by H (u,v) = 0.
7fd59977 1014//
0d969553 1015// with :
7fd59977 1016// ||P2(v0)P1(u)||**2
1017// F(u,v) = P1(u) - 1/2 *----------------* N(u)
1018// (N(u).P2(v0)P1(u))
1019//
1020// H(u,v) = (Tu.P1(u)P2(v))**2||Tv||**2 - (Tv.P1(u)P2(v))**2||Tu||**2
1021//
1022// => dB(u)/du = dF/du + dF/dv(- dH/du:dH/dv)
1023//
0d969553
Y
1024// Note : tangent to the bisectrice is bissectrice at
1025// tangents T1(u) and T2(v0)
7fd59977 1026//
1027//=============================================================================
1028void Bisector_BisecCC::Values (const Standard_Real U,
873c119f 1029 const Standard_Integer N,
1030 gp_Pnt2d& P,
1031 gp_Vec2d& V1,
1032 gp_Vec2d& V2,
1033 gp_Vec2d& V3) const
7fd59977 1034{
1035 V1 = gp_Vec2d(0.,0.);
1036 V2 = gp_Vec2d(0.,0.);
1037 V3 = gp_Vec2d(0.,0.);
1038 //-------------------------------------------------------------------------
0d969553
Y
1039 // Calculate the current point on the bisectrice and the parameters on each
1040 // curve.
7fd59977 1041 //-------------------------------------------------------------------------
1042 Standard_Real U0,V0,Dist;
873c119f 1043
7fd59977 1044 //-----------------------------------------------
0d969553 1045 // is the polygon reduced to a point or empty?
7fd59977 1046 //-----------------------------------------------
1047 if (myPolygon.Length() <= 1) {
1048 P = Extension(U,U0,V0,Dist,V1);
1049 }
1050 if (U < myPolygon.First().ParamOnBis()) {
1051 P = Extension(U,U0,V0,Dist,V1);
1052 return;
1053 }
1054 if (U > myPolygon.Last().ParamOnBis()) {
1055 P = Extension(U,U0,V0,Dist,V1);
1056 return;
1057 }
1058 P = ValueAndDist(U,U0,V0,Dist);
1059
1060 if (N == 0) return;
1061 //------------------------------------------------------------------
0d969553 1062 // Return point, tangent, normal to C1 by parameter U0.
7fd59977 1063 //-------------------------------------------------------------------
0d969553
Y
1064 gp_Pnt2d P1 ; // point on C1.
1065 gp_Vec2d Tu ; // tangent to C1 by U0.
1066 gp_Vec2d Tuu ; // second derivative to C1 by U0.
7fd59977 1067 curve1->D2(U0,P1,Tu,Tuu);
0d969553
Y
1068 gp_Vec2d Nor( - Tu .Y() , Tu .X()); // Normal by U0.
1069 gp_Vec2d Nu ( - Tuu.Y() , Tuu.X()); // derivative of the normal by U0.
7fd59977 1070
1071 //-------------------------------------------------------------------
0d969553 1072 // Return point, tangent, normale to C2 by parameter V0.
7fd59977 1073 //-------------------------------------------------------------------
0d969553
Y
1074 gp_Pnt2d P2 ; // point on C2.
1075 gp_Vec2d Tv ; // tangent to C2 by V.
1076 gp_Vec2d Tvv ; // second derivative to C2 by V.
7fd59977 1077 curve2->D2(V0,P2,Tv,Tvv);
1078
1079 gp_Vec2d PuPv(P2.X() - P1.X(), P2.Y() - P1.Y());
1080
1081 //-----------------------------
0d969553 1082 // Calculate dH/du and dH/dv.
7fd59977 1083 //-----------------------------
1084 Standard_Real TuTu,TvTv,TuTv;
1085 Standard_Real TuPuPv,TvPuPv ;
1086 Standard_Real TuuPuPv,TuTuu ;
1087 Standard_Real TvvPuPv,TvTvv ;
1088
1089 TuTu = Tu.Dot(Tu) ; TvTv = Tv.Dot(Tv) ; TuTv = Tu.Dot(Tv);
1090 TuPuPv = Tu.Dot(PuPv) ; TvPuPv = Tv.Dot(PuPv);
1091 TuuPuPv = Tuu.Dot(PuPv) ; TuTuu = Tu.Dot(Tuu) ;
1092 TvvPuPv = Tvv.Dot(PuPv) ; TvTvv = Tv.Dot(Tvv) ;
1093
1094 Standard_Real dHdu = 2*(TuPuPv*(TuuPuPv - TuTu)*TvTv +
873c119f 1095 TvPuPv*TuTv*TuTu -TuTuu*TvPuPv*TvPuPv);
7fd59977 1096 Standard_Real dHdv = 2*(TuPuPv*TuTv*TvTv + TvTvv*TuPuPv*TuPuPv -
873c119f 1097 TvPuPv*(TvvPuPv + TvTv)*TuTu);
7fd59977 1098
1099 //-----------------------------
0d969553 1100 // Calculate dF/du and dF/dv.
7fd59977 1101 //-----------------------------
1102 Standard_Real NorPuPv,NuPuPv,NorTv;
1103 Standard_Real A,B,dAdu,dAdv,dBdu,dBdv,BB;
873c119f 1104
7fd59977 1105 NorPuPv = Nor.Dot(PuPv);
1106 NuPuPv = Nu .Dot(PuPv);
1107 NorTv = Nor.Dot(Tv) ;
1108
1109 A = 0.5*PuPv.SquareMagnitude();
1110 B = - NorPuPv;
1111 BB = B*B;
1112 dAdu = - TuPuPv;
1113 dBdu = - NuPuPv ;
1114 dAdv = TvPuPv;
1115 dBdv = - NorTv;
873c119f 1116
7fd59977 1117 //---------------------------------------
1118 // F(u,v) = Pu - (A(u,v)/B(u,v))*Nor(u)
1119 //----------------------------------------
1120 if (BB < gp::Resolution()) {
1121 V1 = Tu.Normalized() + Tv.Normalized();
1122 V1 = 0.5*Tu.SquareMagnitude()*V1;
1123 }
1124 else {
1125 gp_Vec2d dFdu = Tu - (dAdu/B - dBdu*A/BB)*Nor - (A/B)*Nu;
1126 gp_Vec2d dFdv = ( - dAdv/B + dBdv*A/BB)*Nor ;
873c119f 1127
7fd59977 1128 if (Abs(dHdv) > gp::Resolution()) {
1129 V1 = dFdu + dFdv*( - dHdu / dHdv );
1130 }
1131 else {
1132 V1 = Tu;
1133 }
1134 }
1135 if (N == 1) return;
1136}
1137
1138//=============================================================================
1139//function : Extension
0d969553
Y
1140// purpose : Calculate the current point on the extensions
1141// by tangence of the curve.
7fd59977 1142//============================================================================
1143gp_Pnt2d Bisector_BisecCC::Extension (const Standard_Real U,
873c119f 1144 Standard_Real& U1,
1145 Standard_Real& U2,
1146 Standard_Real& Dist,
1147 gp_Vec2d& T ) const
7fd59977 1148{
1149 Bisector_PointOnBis PRef;
1150 gp_Pnt2d P,P1,P2,PBis;
1151 gp_Vec2d T1,Tang;
7fd59977 1152 Standard_Real dU = 0.;
7fd59977 1153 Standard_Boolean ExtensionTangent = Standard_False;
1154
1155 if (myPolygon.Length() == 0) {
1156 //---------------------------------------------
0d969553 1157 // Empty Polygon => segment (pointStart,pointEnd)
7fd59977 1158 //---------------------------------------------
1159 dU = U - startIntervals.First();
1160 P = pointStart;
1161 P1 = pointEnd;
1162 U1 = curve1->LastParameter();
1163 if (sign1 == sign2) { U2 = curve2->FirstParameter();}
1164 else { U2 = curve2->LastParameter() ;}
1165 Tang.SetCoord(P1.X() - P.X(),P1.Y() - P.Y());
1166 }
1167 else if (U < myPolygon.First().ParamOnBis()) {
1168 PRef = myPolygon.First();
1169 P = PRef.Point();
1170 dU = U - PRef.ParamOnBis();
1171 if (extensionStart) {
1172 //------------------------------------------------------------
0d969553 1173 // extension = segment (pointstart, first point of the polygon.)
7fd59977 1174 //------------------------------------------------------------
1175 P1 = pointStart;
1176 U1 = curve1->FirstParameter();
1177 if (sign1 == sign2) { U2 = curve2->LastParameter();}
1178 else { U2 = curve2->FirstParameter();}
1179 Tang.SetCoord(P.X() - P1.X(),P.Y() - P1.Y());
1180 }
1181 else {
1182 ExtensionTangent = Standard_True;
1183 }
1184 }
1185 else if (U > myPolygon.Last().ParamOnBis()) {
1186 PRef = myPolygon.Last();
1187 P = PRef.Point();
1188 dU = U - PRef.ParamOnBis();
1189 if (extensionEnd) {
1190 //------------------------------------------------------------
0d969553 1191 // extension = segment (last point of the polygon.pointEnd)
7fd59977 1192 //------------------------------------------------------------
1193 P1 = pointEnd;
1194 U1 = curve1->LastParameter();
1195 if (sign1 == sign2) { U2 = curve2->LastParameter();}
1196 else { U2 = curve2->FirstParameter();}
1197 Tang.SetCoord(P1.X() - P.X(),P1.Y() - P.Y());
1198 }
1199 else {
1200 ExtensionTangent = Standard_True;
1201 }
1202 }
1203
1204 if (ExtensionTangent) {
1205 //-----------------------------------------------------------
0d969553 1206 // If the la curve has no a extension, it is extended by tangency
7fd59977 1207 //------------------------------------------------------------
1208 U1 = PRef.ParamOnC1();
1209 U2 = PRef.ParamOnC2();
1210 P2 = curve2->Value(U2);
1211 curve1->D1(U1,P1,T1);
1212 Tang.SetCoord(2*P.X() - P1.X() - P2.X(), 2*P.Y() - P1.Y() - P2.Y());
1213 if (Tang.Magnitude() < Precision::Confusion()) {
1214 Tang = T1;
1215 }
1216 if (T1.Dot(Tang) < 0.) Tang = - Tang;
1217 }
873c119f 1218
7fd59977 1219 T = Tang.Normalized();
1220 PBis.SetCoord(P.X() + dU*T.X(),P.Y() + dU*T.Y());
1221 Dist = P1.Distance(PBis);
1222 return PBis;
1223}
1224
1225//=============================================================================
1226//function : PointByInt
1227// purpose :
1228//=============================================================================
1229static Standard_Boolean PointByInt(const Handle(Geom2d_Curve)& CA,
873c119f 1230 const Handle(Geom2d_Curve)& CB,
1231 const Standard_Real SignA,
1232 const Standard_Real SignB,
1233 const Standard_Real UOnA,
1234 Standard_Real& UOnB,
1235 Standard_Real& Dist)
7fd59977 1236{
1237 //------------------------------------------------------------------
0d969553 1238 // Return point,tangent, normal on CA with parameter UOnA.
7fd59977 1239 //-------------------------------------------------------------------
1240 gp_Pnt2d P1,P2,P,PSol;
1241 gp_Vec2d Tan1,Tan2;
1242 Standard_Boolean IsConvexA = Bisector::IsConvex(CA,SignA);
1243 Standard_Boolean IsConvexB = Bisector::IsConvex(CB,SignB);
1244
1245 CA->D1(UOnA,P1,Tan1);
1246 gp_Vec2d N1(Tan1.Y(), - Tan1.X());
873c119f 1247
7fd59977 1248 //--------------------------------------------------------------------------
0d969553 1249 // test of confusion of P1 with extremity of curve2.
7fd59977 1250 //--------------------------------------------------------------------------
1251 if (P1.Distance(CB->Value(CB->FirstParameter())) < Precision::Confusion()) {
1252 UOnB = CB->FirstParameter();
1253 CB->D1(UOnB,P2,Tan2);
1254 if ( IsConvexA && IsConvexB ) {
1255 Dist = 0.;
1256 return Standard_True;
1257 }
1258 if (! Tan1.IsParallel(Tan2,Precision::Angular())) {
1259 Dist = 0.;
1260 return Standard_False;
1261 }
1262 }
1263 if (P1.Distance(CB->Value(CB->LastParameter())) < Precision::Confusion()) {
1264 UOnB = CB->LastParameter();
1265 CB->D1(UOnB,P2,Tan2);
1266 if ( IsConvexA && IsConvexB ) {
1267 Dist = 0.;
1268 return Standard_True;
1269 }
1270 if (! Tan1.IsParallel(Tan2,Precision::Angular())) {
1271 Dist = 0.;
1272 return Standard_False;
1273 }
1274 }
1275
1276 Standard_Real DMin = Precision::Infinite();
1277 Standard_Real UPC;
1278 Standard_Boolean YaSol = Standard_False;
873c119f 1279 //--------------------------------------------------------------------
0d969553
Y
1280 // Construction of the bisectrice point curve and of the straight line passing
1281 // through P1 and carried by the normal.
7fd59977 1282 //--------------------------------------------------------------------
1283 Handle(Bisector_BisecPC) BisPC
1284 = new Bisector_BisecPC(CB,P1,SignB );
1285 //-------------------------------
0d969553 1286 // Test if the bissectrice exists.
7fd59977 1287 //-------------------------------
1288 if (BisPC->IsEmpty()) {
1289 Dist = Precision::Infinite();
1290 PSol = P1;
1291 return Standard_False;
1292 }
1293
1294 Handle(Geom2d_Line) NorLi = new Geom2d_Line (P1,N1);
1295
1296 Geom2dAdaptor_Curve ABisPC(BisPC);
1297 Geom2dAdaptor_Curve ANorLi(NorLi);
1298 //-------------------------------------------------------------------------
1299 Geom2dInt_GInter Intersect(ABisPC,ANorLi,
873c119f 1300 Precision::Confusion(),Precision::Confusion());
7fd59977 1301 //-------------------------------------------------------------------------
1302
1303 if (Intersect.IsDone() && !Intersect.IsEmpty()) {
1304 for (Standard_Integer i = 1; i <= Intersect.NbPoints(); i++) {
1305 if (Intersect.Point(i).ParamOnSecond()*SignA < Precision::PConfusion()) {
873c119f 1306 P = Intersect.Point(i).Value();
1307 if (P.SquareDistance(P1) < DMin) {
1308 DMin = P.SquareDistance(P1);
1309 PSol = P;
1310 UPC = Intersect.Point(i).ParamOnFirst();
1311 UOnB = BisPC->LinkBisCurve(UPC);
1312 Dist = DMin;
1313 YaSol = Standard_True;
1314 }
7fd59977 1315 }
1316 }
1317 }
1318 if (YaSol) {
1319 //--------------------------------------------------------------
0d969553 1320 // Point found => Test distance curvature + Angular test
7fd59977 1321 //---------------------------------------------------------------
1322 P2 = CB->Value(UOnB);
873c119f 1323 if(P1.SquareDistance(PSol) < 1.e-32)
1324 {
1325 YaSol = Standard_False;
1326 return YaSol;
1327 }
1328 if(P2.SquareDistance(PSol) < 1.e-32)
1329 {
1330 YaSol = Standard_False;
1331 return YaSol;
1332 }
1333
7fd59977 1334 gp_Dir2d PP1Unit(P1.X() - PSol.X(),P1.Y() - PSol.Y());
1335 gp_Dir2d PP2Unit(P2.X() - PSol.X(),P2.Y() - PSol.Y());
873c119f 1336
7fd59977 1337 if (PP1Unit*PP2Unit > 1. - Precision::Angular()) {
1338 YaSol = Standard_False;
1339 }
1340 else {
1341 Dist = sqrt(Dist);
1342 if ( !IsConvexA ) {
873c119f 1343 Standard_Real K1 = Curvature(CA,UOnA,Precision::Confusion());
1344 if (K1 != 0.) {
1345 if (Dist > Abs(1/K1)) YaSol = Standard_False;
1346 }
7fd59977 1347 }
1348 if (YaSol) {
873c119f 1349 if ( !IsConvexB ) {
1350 Standard_Real K2 = Curvature(CB,UOnB,Precision::Confusion());
1351 if (K2 != 0.) {
1352 if (Dist > Abs(1/K2)) YaSol = Standard_False;
1353 }
1354 }
7fd59977 1355 }
1356 }
1357 }
1358 return YaSol;
1359}
1360
1361//=============================================================================
1362//function : SupLastParameter
1363// purpose :
1364//=============================================================================
1365void Bisector_BisecCC::SupLastParameter()
1366{
1367 endIntervals.Append(curve1->LastParameter());
1368 // ----------------------------------------------------------------------
0d969553
Y
1369 // Calculate parameter on curve1 associated to one or the other of the extremities
1370 // of curve2 following the values of sign1 and sign2.
1371 // the bissectrice is limited by the obtained parameters.
7fd59977 1372 //------------------------------------------------------------------------
1373 Standard_Real UOnC1,UOnC2,Dist;
1374 if (sign1 == sign2) {
1375 UOnC2 = curve2->FirstParameter();
1376 }
1377 else {
1378 UOnC2 = curve2->LastParameter();
1379 }
1380 Standard_Boolean YaSol = PointByInt(curve2,curve1,sign2,sign1,UOnC2,UOnC1,Dist);
1381 if (YaSol) {
1382 if (UOnC1 > startIntervals.First() && UOnC1 < endIntervals.Last()) {
1383 endIntervals.SetValue(1,UOnC1);
1384 }
1385 }
1386}
1387
1388//=============================================================================
1389//function : Curve
1390// purpose :
1391//=============================================================================
1392Handle(Geom2d_Curve) Bisector_BisecCC::Curve(const Standard_Integer I) const
1393{
1394 if (I == 1) return curve1;
1395 else if (I == 2) return curve2;
1396 else Standard_OutOfRange::Raise();
1397 return curve1;
1398}
1399
1400//=============================================================================
1401//function : LinkBisCurve
1402//purpose :
1403//=============================================================================
1404Standard_Real Bisector_BisecCC::LinkBisCurve(const Standard_Real U) const
1405{
1406 return (U - shiftParameter);
1407}
1408
1409//=============================================================================
1410//function : LinkCurveBis
1411//purpose :
1412//=============================================================================
1413Standard_Real Bisector_BisecCC::LinkCurveBis(const Standard_Real U) const
1414{
1415 return (U + shiftParameter);
1416}
1417
1418//=============================================================================
1419//function : Indent
1420//purpose :
1421//=============================================================================
1422static void Indent(const Standard_Integer Offset) {
1423 if (Offset > 0) {
1424 for (Standard_Integer i = 0; i < Offset; i++) {cout << " ";}
1425 }
1426}
1427
1428//=============================================================================
1429//function : Polygon
1430// purpose :
1431//=============================================================================
1432const Bisector_PolyBis& Bisector_BisecCC::Polygon() const
1433{
1434 return myPolygon;
1435}
1436
1437//==========================================================================
1438//function : Parameter
1439//purpose :
1440//==========================================================================
1441Standard_Real Bisector_BisecCC::Parameter(const gp_Pnt2d& P) const
1442{
1443 Standard_Real UOnCurve;
1444
1445 if (P.IsEqual(Value(FirstParameter()),Precision::Confusion())) {
1446 UOnCurve = FirstParameter();
1447 }
1448 else if (P.IsEqual(Value(LastParameter()),Precision::Confusion())) {
1449 UOnCurve = LastParameter();
1450 }
91bb31f3 1451 else
1452 {
1453 ProjOnCurve(P, curve1, UOnCurve);
7fd59977 1454 }
91bb31f3 1455
7fd59977 1456 return UOnCurve;
1457}
1458
1459
1460//=============================================================================
1461//function : Dump
1462// purpose :
1463//=============================================================================
1464//void Bisector_BisecCC::Dump(const Standard_Integer Deep,
1465void Bisector_BisecCC::Dump(const Standard_Integer ,
873c119f 1466 const Standard_Integer Offset) const
7fd59977 1467{
1468 Indent (Offset);
1469 cout <<"Bisector_BisecCC :"<<endl;
1470 Indent (Offset);
873c119f 1471 // cout <<"Curve1 :"<<curve1<<endl;
1472 // cout <<"Curve2 :"<<curve2<<endl;
7fd59977 1473 cout <<"Sign1 :"<<sign1<<endl;
1474 cout <<"Sign2 :"<<sign2<<endl;
1475
1476 cout <<"Number Of Intervals :"<<startIntervals.Length()<<endl;
1477 for (Standard_Integer i = 1; i <= startIntervals.Length(); i++) {
1478 cout <<"Interval number :"<<i<<"Start :"<<startIntervals.Value(i)
873c119f 1479 <<" end :"<< endIntervals.Value(i)<<endl ;
7fd59977 1480 }
1481 cout <<"Index Current Interval :"<<currentInterval<<endl;
1482}
1483
1484//=============================================================================
1485//function : Curve
1486// purpose :
1487//=============================================================================
1488void Bisector_BisecCC::Curve(const Standard_Integer I,
873c119f 1489 const Handle(Geom2d_Curve)& C)
7fd59977 1490{
1491 if (I == 1) curve1 = C;
1492 else if (I == 2) curve2 = C;
1493 else Standard_OutOfRange::Raise();
1494}
1495
1496//=============================================================================
1497//function : Sign
1498// purpose :
1499//=============================================================================
1500void Bisector_BisecCC::Sign(const Standard_Integer I,
873c119f 1501 const Standard_Real S)
7fd59977 1502{
1503 if (I == 1) sign1 = S;
1504 else if (I == 2) sign2 = S;
1505 else Standard_OutOfRange::Raise();
1506}
1507
1508//=============================================================================
1509//function : Polygon
1510// purpose :
1511//=============================================================================
1512void Bisector_BisecCC::Polygon(const Bisector_PolyBis& P)
1513{
1514 myPolygon = P;
1515}
1516
1517//=============================================================================
1518//function : DistMax
1519// purpose :
1520//=============================================================================
1521void Bisector_BisecCC::DistMax(const Standard_Real D)
1522{
1523 distMax = D;
1524}
1525
1526//=============================================================================
1527//function : IsConvex
1528// purpose :
1529//=============================================================================
1530void Bisector_BisecCC::IsConvex(const Standard_Integer I,
873c119f 1531 const Standard_Boolean IsConvex)
7fd59977 1532{
1533 if (I == 1) isConvex1 = IsConvex;
1534 else if (I == 2) isConvex2 = IsConvex;
1535 else Standard_OutOfRange::Raise();
1536}
1537
1538//=============================================================================
1539//function : IsEmpty
1540// purpose :
1541//=============================================================================
1542void Bisector_BisecCC::IsEmpty ( const Standard_Boolean IsEmpty)
1543{
1544 isEmpty = IsEmpty;
1545}
1546
1547//=============================================================================
1548//function : ExtensionStart
1549// purpose :
1550//=============================================================================
1551void Bisector_BisecCC::ExtensionStart( const Standard_Boolean ExtensionStart)
1552{
1553 extensionStart = ExtensionStart;
1554}
1555
1556//=============================================================================
1557//function : ExtensionEnd
1558// purpose :
1559//=============================================================================
1560void Bisector_BisecCC::ExtensionEnd( const Standard_Boolean ExtensionEnd)
1561{
1562 extensionEnd = ExtensionEnd;
1563}
1564
1565//=============================================================================
1566//function : PointStart
1567// purpose :
1568//=============================================================================
1569void Bisector_BisecCC::PointStart( const gp_Pnt2d& Point)
1570{
1571 pointStart = Point;
1572}
1573
1574//=============================================================================
1575//function : PointEnd
1576// purpose :
1577//=============================================================================
1578void Bisector_BisecCC::PointEnd( const gp_Pnt2d& Point)
1579{
1580 pointEnd = Point;
1581}
1582
1583//=============================================================================
1584//function : StartIntervals
1585// purpose :
1586//=============================================================================
1587void Bisector_BisecCC::StartIntervals
1588 (const TColStd_SequenceOfReal& StartIntervals)
1589{
1590 startIntervals = StartIntervals;
1591}
1592
1593//=============================================================================
1594//function : EndIntervals
1595// purpose :
1596//=============================================================================
1597void Bisector_BisecCC::EndIntervals
1598 (const TColStd_SequenceOfReal& EndIntervals)
1599{
1600 endIntervals = EndIntervals;
1601}
1602
1603//=============================================================================
1604//function : FirstParameter
1605// purpose :
1606//=============================================================================
1607void Bisector_BisecCC::FirstParameter (const Standard_Real U)
1608{
1609 startIntervals.Append(U);
1610}
1611
1612//=============================================================================
1613//function : LastParameter
1614// purpose :
1615//=============================================================================
1616void Bisector_BisecCC::LastParameter (const Standard_Real U)
1617{
1618 endIntervals.Append(U);
1619}
1620
1621//=============================================================================
1622//function : SearchBound
1623// purpose :
1624//=============================================================================
1625Standard_Real Bisector_BisecCC::SearchBound (const Standard_Real U1,
873c119f 1626 const Standard_Real U2) const
7fd59977 1627{
1628 Standard_Real UMid,Dist1,Dist2,DistMid,U11,U22;
1629 Standard_Real UC1,UC2;
1630 gp_Pnt2d PBis,PBisPrec;
1631 Standard_Real TolPnt = Precision::Confusion();
1632 Standard_Real TolPar = Precision::PConfusion();
1633 U11 = U1; U22 = U2;
1634 PBisPrec = ValueByInt(U11,UC1,UC2,Dist1);
1635 PBis = ValueByInt(U22,UC1,UC2,Dist2);
873c119f 1636
7fd59977 1637 while ((U22 - U11) > TolPar ||
873c119f 1638 ((Dist1 < Precision::Infinite() &&
1639 Dist2 < Precision::Infinite() &&
1640 !PBis.IsEqual(PBisPrec,TolPnt)))) {
1641 PBisPrec = PBis;
1642 UMid = 0.5*( U22 + U11);
1643 PBis = ValueByInt(UMid,UC1,UC2,DistMid);
1644 if ((Dist1 < Precision::Infinite()) == (DistMid < Precision::Infinite())) {
1645 U11 = UMid;
1646 Dist1 = DistMid;
1647 }
1648 else {
1649 U22 = UMid;
1650 Dist2 = DistMid;
1651 }
7fd59977 1652 }
1653 PBis = ValueByInt(U11,UC1,UC2,Dist1);
1654 if (Dist1 < Precision::Infinite()) {
1655 UMid = U11;
1656 }
1657 else {
1658 UMid = U22;
1659 }
1660 return UMid;
1661}
1662
1663//=============================================================================
1664//function : ProjOnCurve
1665// purpose :
1666//=============================================================================
91bb31f3 1667static Standard_Boolean ProjOnCurve (const gp_Pnt2d& P,
873c119f 1668 const Handle(Geom2d_Curve)& C,
1669 Standard_Real& theParam)
7a06c690 1670{
91bb31f3 1671 //Standard_Real UOnCurve =0.;
1672 theParam = 0.0;
7fd59977 1673 gp_Pnt2d PF,PL;
1674 gp_Vec2d TF,TL;
1675
1676 C->D1(C->FirstParameter(),PF,TF);
1677 C->D1(C->LastParameter() ,PL,TL);
1678
91bb31f3 1679 if (P.IsEqual(PF ,Precision::Confusion()))
1680 {
1681 theParam = C->FirstParameter();
1682 return Standard_True;
7fd59977 1683 }
873c119f 1684
91bb31f3 1685 if (P.IsEqual(PL ,Precision::Confusion()))
1686 {
1687 theParam = C->LastParameter();
1688 return Standard_True;
7fd59977 1689 }
873c119f 1690
7fd59977 1691 gp_Vec2d PPF(PF.X() - P.X(), PF.Y() - P.Y());
1692 TF.Normalize();
873c119f 1693
91bb31f3 1694 if ( Abs (PPF.Dot(TF)) < Precision::Confusion())
1695 {
1696 theParam = C->FirstParameter();
1697 return Standard_True;
7fd59977 1698 }
1699 gp_Vec2d PPL (PL.X() - P.X(), PL.Y() - P.Y());
1700 TL.Normalize();
91bb31f3 1701 if ( Abs (PPL.Dot(TL)) < Precision::Confusion())
1702 {
1703 theParam = C->LastParameter();
1704 return Standard_True;
7fd59977 1705 }
1706 Geom2dAPI_ProjectPointOnCurve Proj(P,C,
873c119f 1707 C->FirstParameter(),
1708 C->LastParameter());
7fd59977 1709 if (Proj.NbPoints() > 0) {
91bb31f3 1710 theParam = Proj.LowerDistanceParameter();
7fd59977 1711 }
1712 else {
91bb31f3 1713 return Standard_False;
7fd59977 1714 }
91bb31f3 1715
1716 return Standard_True;
7fd59977 1717}
1718
1719//=============================================================================
1720//function : TestExtension
1721// purpose :
1722//=============================================================================
1723static Standard_Boolean TestExtension (const Handle(Geom2d_Curve)& C1,
873c119f 1724 const Handle(Geom2d_Curve)& C2,
1725 const Standard_Integer Start_End)
7fd59977 1726{
1727 gp_Pnt2d P1,P2;
1728 gp_Vec2d T1,T2;
1729 Standard_Boolean Test = Standard_False;
1730 if (Start_End == 1) {
1731 C1->D1(C1->FirstParameter(),P1,T1);
1732 }
1733 else {
1734 C1->D1(C1->LastParameter(),P1,T1);
1735 }
1736 C2->D1(C2->FirstParameter(),P2,T2);
1737 if (P1.IsEqual(P2,Precision::Confusion())) {
1738 T1.Normalize(); T2.Normalize();
1739 if (T1.Dot(T2) > 1.0 - Precision::Confusion()) {
1740 Test = Standard_True;
1741 }
1742 }
1743 else {
1744 C2->D1(C2->LastParameter(),P2,T2);
1745 if (P1.IsEqual(P2,Precision::Confusion())) {
1746 T2.Normalize();
1747 if (T1.Dot(T2) > 1.0 - Precision::Confusion()) {
873c119f 1748 Test = Standard_True;
7fd59977 1749 }
1750 }
1751 }
1752 return Test;
1753}
1754
1755//=============================================================================
1756//function : ComputePointEnd
1757// purpose :
1758//=============================================================================
1759void Bisector_BisecCC::ComputePointEnd ()
1760{
1761 Standard_Real U1,U2;
1762 Standard_Real KC,RC;
1763 U1 = curve1->FirstParameter();
1764 if (sign1 == sign2) {
1765 U2 = curve2->LastParameter();
1766 }
1767 else {
1768 U2 = curve2->FirstParameter();
1769 }
1770 Standard_Real K1 = Curvature(curve1,U1,Precision::Confusion());
1771 Standard_Real K2 = Curvature(curve2,U2,Precision::Confusion());
1772 if (!isConvex1 && !isConvex2) {
1773 if (K1 < K2) {KC = K1;} else {KC = K2;}
1774 }
1775 else if (!isConvex1) {KC = K1;}
1776 else {KC = K2;}
1777
1778 gp_Pnt2d PF;
1779 gp_Vec2d TF;
1780 curve1->D1(U1,PF,TF);
1781 TF.Normalize();
1782 if (KC != 0.) { RC = Abs(1/KC);}
1783 else { RC = Precision::Infinite();}
1784 pointEnd.SetCoord(PF.X() - sign1*RC*TF.Y(), PF.Y() + sign1*RC*TF.X());
1785
1786}
1787
1788//=============================================================================
1789//function : DiscretPar
1790// purpose :
1791//=============================================================================
1792static Standard_Boolean DiscretPar(const Standard_Real DU,
873c119f 1793 const Standard_Real EpsMin,
1794 const Standard_Real EpsMax,
1795 const Standard_Integer NbMin,
1796 const Standard_Integer NbMax,
1797 Standard_Real& Eps,
1798 Standard_Integer& Nb)
7fd59977 1799{
1800 if (DU <= NbMin*EpsMin) {
1801 Eps = DU/(NbMin + 1) ;
1802 Nb = NbMin;
1803 return Standard_False;
1804 }
1805
1806 Eps = Min (EpsMax,DU/NbMax);
1807
1808 if (Eps < EpsMin) {
1809 Eps = EpsMin;
1810 Nb = Standard_Integer(DU/EpsMin);
1811 }
1812 else { Nb = NbMax;}
1813
1814 return Standard_True;
1815}
1816
1817