1 // Created on: 1993-07-07
2 // Created by: Jean Claude VAUTHIER
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 // Version:
18 //pmn 24/09/96 Ajout du prolongement de courbe.
19 // jct 15/04/97 Ajout du prolongement de surface.
20 // jct 24/04/97 simplification ou suppression de calculs
21 // inutiles dans ExtendSurfByLength
22 // correction de Tbord et Continuity=0 accepte
23 // correction du calcul de lambda et appel a
24 // TangExtendToConstraint avec lambmin au lieu de 1.
25 // correction du passage Sr rat --> BSp nD
26 // xab 26/06/97 treatement partiel anulation des derivees
27 // partiels du denonimateur des Surfaces BSplines Rationnelles
28 // dans le cas de valeurs proportionnelles des denominateurs
29 // en umin umax et/ou vmin vmax.
30 // pmn 4/07/97 Gestion de la continuite dans BuildCurve3d (PRO9097)
31 // xab 10/07/97 on revient en arriere sur l'ajout du 26/06/97
32 // pmn 26/09/97 Ajout des parametres d'approx dans BuildCurve3d
33 // xab 29/09/97 on reintegre l'ajout du 26/06/97
34 // pmn 31/10/97 Ajoute AdjustExtremity
35 // jct 26/11/98 blindage dans ExtendSurf qd NTgte = 0 (CTS21288)
36 // jct 19/01/99 traitement de la periodicite dans ExtendSurf
37 // Design:
38 // Warning: None
39 // References: None
40 // Language: C++2.0
41 // Purpose:
42 // Declarations:
43
44 #include <GeomLib.hxx>
45
46 #include <Adaptor2d_Curve2d.hxx>
47 #include <Adaptor3d_Curve.hxx>
48 #include <Adaptor3d_CurveOnSurface.hxx>
49 #include <Adaptor3d_Surface.hxx>
50 #include <AdvApprox_ApproxAFunction.hxx>
51 #include <AdvApprox_PrefAndRec.hxx>
52 #include <Approx_CurveOnSurface.hxx>
53 #include <BSplCLib.hxx>
54 #include <BSplSLib.hxx>
55 #include <CSLib.hxx>
56 #include <CSLib_NormalStatus.hxx>
57 #include <ElCLib.hxx>
58 #include <Geom2d_BezierCurve.hxx>
59 #include <Geom2d_BSplineCurve.hxx>
60 #include <Geom2d_Circle.hxx>
61 #include <Geom2d_Curve.hxx>
62 #include <Geom2d_Ellipse.hxx>
63 #include <Geom2d_Hyperbola.hxx>
64 #include <Geom2d_Line.hxx>
65 #include <Geom2d_OffsetCurve.hxx>
66 #include <Geom2d_Parabola.hxx>
67 #include <Geom2d_TrimmedCurve.hxx>
68 #include <Geom2dAdaptor_Curve.hxx>
69 #include <Geom2dConvert.hxx>
70 #include <Geom_BezierCurve.hxx>
71 #include <Geom_BezierSurface.hxx>
72 #include <Geom_BoundedCurve.hxx>
73 #include <Geom_BoundedSurface.hxx>
74 #include <Geom_BSplineCurve.hxx>
75 #include <Geom_BSplineSurface.hxx>
76 #include <Geom_Circle.hxx>
77 #include <Geom_Curve.hxx>
78 #include <Geom_Ellipse.hxx>
79 #include <Geom_Hyperbola.hxx>
80 #include <Geom_Line.hxx>
81 #include <Geom_OffsetCurve.hxx>
82 #include <Geom_Parabola.hxx>
83 #include <Geom_Plane.hxx>
84 #include <Geom_RectangularTrimmedSurface.hxx>
85 #include <Geom_Surface.hxx>
86 #include <Geom_TrimmedCurve.hxx>
87 #include <GeomAdaptor_Surface.hxx>
88 #include <GeomConvert.hxx>
89 #include <GeomConvert_ApproxSurface.hxx>
90 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
91 #include <GeomLib_DenominatorMultiplier.hxx>
92 #include <GeomLib_DenominatorMultiplierPtr.hxx>
93 #include <GeomLib_LogSample.hxx>
94 #include <GeomLib_MakeCurvefromApprox.hxx>
95 #include <GeomLib_PolyFunc.hxx>
96 #include <gp_Ax2.hxx>
97 #include <gp_Circ.hxx>
98 #include <gp_Circ2d.hxx>
99 #include <gp_Dir.hxx>
100 #include <gp_Elips.hxx>
101 #include <gp_Elips2d.hxx>
102 #include <gp_GTrsf2d.hxx>
103 #include <gp_Hypr.hxx>
104 #include <gp_Hypr2d.hxx>
105 #include <gp_Lin.hxx>
106 #include <gp_Lin2d.hxx>
107 #include <gp_Parab.hxx>
108 #include <gp_Parab2d.hxx>
109 #include <gp_Pnt.hxx>
110 #include <gp_Pnt2d.hxx>
111 #include <gp_Trsf2d.hxx>
112 #include <gp_TrsfForm.hxx>
113 #include <gp_Vec.hxx>
114 #include <Hermit.hxx>
115 #include <math.hxx>
116 #include <math_FunctionAllRoots.hxx>
117 #include <math_FunctionSample.hxx>
118 #include <math_Jacobi.hxx>
119 #include <math_Matrix.hxx>
120 #include <math_Vector.hxx>
121 #include <PLib.hxx>
122 #include <Precision.hxx>
123 #include <Standard_ConstructionError.hxx>
124 #include <Standard_NotImplemented.hxx>
125 #include <TColgp_Array1OfPnt.hxx>
126 #include <TColgp_Array1OfPnt2d.hxx>
127 #include <TColgp_Array1OfVec.hxx>
128 #include <TColgp_Array1OfXYZ.hxx>
129 #include <TColgp_Array2OfPnt.hxx>
130 #include <TColgp_HArray2OfPnt.hxx>
131 #include <TColStd_Array1OfInteger.hxx>
132 #include <TColStd_Array1OfReal.hxx>
133 #include <TColStd_Array2OfReal.hxx>
134 #include <TColStd_HArray1OfReal.hxx>
135 #include <TColStd_HArray2OfReal.hxx>
136 //
137 static Standard_Boolean CompareWeightPoles(const TColgp_Array1OfPnt& thePoles1,
138 const TColStd_Array1OfReal* const theW1,
139 const TColgp_Array1OfPnt& thePoles2,
140 const TColStd_Array1OfReal* const theW2,
141 const Standard_Real theTol);
142
143 //=======================================================================
144 //function : ComputeLambda
145 //purpose : Calcul le facteur lambda qui minimise la variation de vittesse
146 // sur une interpolation d'hermite d'ordre (i,0)
147 //=======================================================================
ComputeLambda(const math_Matrix & Constraint,const math_Matrix & Hermit,const Standard_Real Length,Standard_Real & Lambda)148 static void ComputeLambda(const math_Matrix& Constraint,
149 const math_Matrix& Hermit,
150 const Standard_Real Length,
151 Standard_Real& Lambda )
152 {
153 Standard_Integer size = Hermit.RowNumber();
154 Standard_Integer Continuity = size-2;
155 Standard_Integer ii, jj, ip, pp;
156
157 //Minimization
158 math_Matrix HDer(1, size-1, 1, size);
159 for (jj=1; jj<=size; jj++) {
160 for (ii=1; ii<size;ii++) {
161 HDer(ii, jj) = ii*Hermit(jj, ii+1);
162 }
163 }
164
165 math_Vector V(1, size);
166 math_Vector Vec1(1, Constraint.RowNumber());
167 math_Vector Vec2(1, Constraint.RowNumber());
168 math_Vector Vec3(1, Constraint.RowNumber());
169 math_Vector Vec4(1, Constraint.RowNumber());
170
171 Standard_Real * polynome = &HDer(1,1);
172 Standard_Real * valhder = &V(1);
173 Vec2 = Constraint.Col(2);
174 Vec2 /= Length;
175 Standard_Real t, squared1 = Vec2.Norm2(), GW;
176 // math_Matrix Vec(1, Constraint.RowNumber(), 1, size-1);
177 // gp_Vec Vfirst(p0.XYZ()), Vlast(Point.XYZ());
178 // TColgp_Array1OfVec Der(2, 4);
179 // Der(2) = d1; Der(3) = d2; Der(4) = d3;
180
181 Standard_Integer GOrdre = 4 + 4*Continuity,
182 DDim=Continuity*(Continuity+2);
183 math_Vector GaussP(1, GOrdre), GaussW(1, GOrdre),
184 pol2(1, 2*Continuity+1),
185 pol4(1, 4*Continuity+1);
186 math::GaussPoints(GOrdre, GaussP);
187 math::GaussWeights (GOrdre, GaussW);
188 pol4.Init(0.);
189
190 for (ip=1; ip<=GOrdre; ip++) {
191 t = (GaussP(ip)+1.)/2;
192 GW = GaussW(ip);
193 PLib::NoDerivativeEvalPolynomial(t , Continuity, Continuity+2, DDim,
194 polynome[0], valhder[0]);
195 V /= Length; //Normalisation
196
197 // i
198 // C'(t) = SUM Vi*Lambda
199 Vec1 = Constraint.Col(1);
200 Vec1 *= V(1);
201 Vec1 += V(size)*Constraint.Col(size);
202 Vec2 = Constraint.Col(2);
203 Vec2 *= V(2);
204 if (Continuity > 1) {
205 Vec3 = Constraint.Col(3);
206 Vec3 *= V(3);
207 if (Continuity > 2) {
208 Vec4 = Constraint.Col(4);
209 Vec4 *= V(4);
210 }
211 }
212
213
214 // 2 2
215 // C'(t) - C'(0)
216
217 pol2(1) = Vec1.Norm2();
218 pol2(2) = 2*(Vec1.Multiplied(Vec2));
219 pol2(3) = Vec2.Norm2() - squared1;
220 if (Continuity>1) {
221 pol2(3) += 2*(Vec1.Multiplied(Vec3));
222 pol2(4) = 2*(Vec2.Multiplied(Vec3));
223 pol2(5) = Vec3.Norm2();
224 if (Continuity>2) {
225 pol2(4)+= 2*(Vec1.Multiplied(Vec4));
226 pol2(5)+= 2*(Vec2.Multiplied(Vec4));
227 pol2(6) = 2*(Vec3.Multiplied(Vec4));
228 pol2(7) = Vec4.Norm2();
229 }
230 }
231
232 // 2 2 2
233 // Integrale de ( C'(t) - C'(0) )
234 for (ii=1; ii<=pol2.Length(); ii++) {
235 pp = ii;
236 for(jj=1; jj<ii; jj++, pp++) {
237 pol4(pp) += 2*GW*pol2(ii)*pol2(jj);
238 }
239 pol4(2*ii-1) += GW*Pow(pol2(ii), 2);
240 }
241 }
242
243 Standard_Real EMin, E;
244 PLib::NoDerivativeEvalPolynomial(Lambda , pol4.Length()-1, 1,
245 pol4.Length()-1,
246 pol4(1), EMin);
247
248 if (EMin > Precision::Confusion()) {
249 // Recheche des extrema de la fonction
250 GeomLib_PolyFunc FF(pol4);
251 GeomLib_LogSample S(Lambda/1000, 50*Lambda, 100);
252 math_FunctionAllRoots Solve(FF, S, Precision::Confusion(),
253 Precision::Confusion()*(Length+1),
254 1.e-15);
255 if (Solve.IsDone()) {
256 for (ii=1; ii<=Solve.NbPoints(); ii++) {
257 t = Solve.GetPoint(ii);
258 PLib::NoDerivativeEvalPolynomial(t , pol4.Length()-1, 1,
259 pol4.Length()-1,
260 pol4(1), E);
261 if (E < EMin) {
262 Lambda = t;
263 EMin = E;
264 }
265 }
266 }
267 }
268 }
269
270 #include <Extrema_LocateExtPC.hxx>
271 #include <Geom2d_Curve.hxx>
272 //=======================================================================
273 //function : RemovePointsFromArray
274 //purpose :
275 //=======================================================================
276
RemovePointsFromArray(const Standard_Integer NumPoints,const TColStd_Array1OfReal & InParameters,Handle (TColStd_HArray1OfReal)& OutParameters)277 void GeomLib::RemovePointsFromArray(const Standard_Integer NumPoints,
278 const TColStd_Array1OfReal& InParameters,
279 Handle(TColStd_HArray1OfReal)& OutParameters)
280 {
281 Standard_Integer ii,
282 jj,
283 add_one_point,
284 loc_num_points,
285 num_points,
286 index ;
287 Standard_Real delta,
288 current_parameter ;
289
290 loc_num_points = Max(0,NumPoints-2) ;
291 delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
292 delta /= (Standard_Real) (loc_num_points + 1) ;
293 num_points = 1 ;
294 current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
295 ii = InParameters.Lower() + 1 ;
296 for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
297 add_one_point = 0 ;
298 while ( ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
299 ii += 1 ;
300 add_one_point = 1 ;
301 }
302 num_points += add_one_point ;
303 current_parameter += delta ;
304 }
305 if (NumPoints <= 2) {
306 num_points = 2 ;
307 }
308 index = 2 ;
309 current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
310 OutParameters =
311 new TColStd_HArray1OfReal(1,num_points) ;
312 OutParameters->ChangeArray1()(1) = InParameters(InParameters.Lower()) ;
313 ii = InParameters.Lower() + 1 ;
314 for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
315 add_one_point = 0 ;
316 while (ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
317 ii += 1 ;
318 add_one_point = 1 ;
319 }
320 if (add_one_point && index <= num_points) {
321 OutParameters->ChangeArray1()(index) = InParameters(ii-1) ;
322 index += 1 ;
323 }
324 current_parameter += delta ;
325 }
326 OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
327 }
328 //=======================================================================
329 //function : DensifyArray1OfReal
330 //purpose :
331 //=======================================================================
332
DensifyArray1OfReal(const Standard_Integer MinNumPoints,const TColStd_Array1OfReal & InParameters,Handle (TColStd_HArray1OfReal)& OutParameters)333 void GeomLib::DensifyArray1OfReal(const Standard_Integer MinNumPoints,
334 const TColStd_Array1OfReal& InParameters,
335 Handle(TColStd_HArray1OfReal)& OutParameters)
336 {
337 Standard_Integer ii,
338 in_order,
339 num_points,
340 num_parameters_to_add,
341 index ;
342 Standard_Real delta,
343 current_parameter ;
344
345 in_order = 1 ;
346 if (MinNumPoints > InParameters.Length()) {
347
348 //
349 // checks the parameters are in increasing order
350 //
351 for (ii = InParameters.Lower() ; ii < InParameters.Upper() ; ii++) {
352 if (InParameters(ii) > InParameters(ii+1)) {
353 in_order = 0 ;
354 break ;
355 }
356 }
357 if (in_order) {
358 num_parameters_to_add = MinNumPoints - InParameters.Length() ;
359 delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
360 delta /= (Standard_Real) (num_parameters_to_add + 1) ;
361 num_points = MinNumPoints ;
362 OutParameters =
363 new TColStd_HArray1OfReal(1,num_points) ;
364 index = 1 ;
365 current_parameter = InParameters(InParameters.Lower()) ;
366 OutParameters->ChangeArray1()(index) = current_parameter ;
367 index += 1 ;
368 current_parameter += delta ;
369 for (ii = InParameters.Lower() + 1 ; index <= num_points && ii <= InParameters.Upper() ; ii++) {
370 while (current_parameter < InParameters(ii) && index <= num_points) {
371 OutParameters->ChangeArray1()(index) = current_parameter ;
372 index += 1 ;
373 current_parameter += delta ;
374 }
375 if (index <= num_points) {
376 OutParameters->ChangeArray1()(index) = InParameters(ii) ;
377 }
378 index += 1 ;
379 }
380 //
381 // beware of roundoff !
382 //
383 OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
384 }
385 else {
386 index = 1 ;
387 num_points = InParameters.Length() ;
388 OutParameters =
389 new TColStd_HArray1OfReal(1,num_points) ;
390 for (ii = InParameters.Lower() ; ii <= InParameters.Upper() ; ii++) {
391 OutParameters->ChangeArray1()(index) = InParameters(ii) ;
392 index += 1 ;
393 }
394 }
395 }
396 else {
397 index = 1 ;
398 num_points = InParameters.Length() ;
399 OutParameters =
400 new TColStd_HArray1OfReal(1,num_points) ;
401 for (ii = InParameters.Lower() ; ii <= InParameters.Upper() ; ii++) {
402 OutParameters->ChangeArray1()(index) = InParameters(ii) ;
403 index += 1 ;
404 }
405 }
406 }
407
408 //=======================================================================
409 //function : FuseIntervals
410 //purpose :
411 //=======================================================================
FuseIntervals(const TColStd_Array1OfReal & I1,const TColStd_Array1OfReal & I2,TColStd_SequenceOfReal & Seq,const Standard_Real Epspar,const Standard_Boolean IsAdjustToFirstInterval)412 void GeomLib::FuseIntervals(const TColStd_Array1OfReal& I1,
413 const TColStd_Array1OfReal& I2,
414 TColStd_SequenceOfReal& Seq,
415 const Standard_Real Epspar,
416 const Standard_Boolean IsAdjustToFirstInterval)
417 {
418 Standard_Integer ind1=1, ind2=1;
419 Standard_Real v1, v2;
420 // Initialisations : les IND1 et IND2 pointent sur le 1er element
421 // de chacune des 2 tables a traiter.INDS pointe sur le dernier
422 // element cree de TABSOR
423
424
425 //--- On remplit TABSOR en parcourant TABLE1 et TABLE2 simultanement ---
426 //------------------ en eliminant les occurrences multiples ------------
427
428 while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
429 v1 = I1(ind1);
430 v2 = I2(ind2);
431 if (Abs(v1-v2)<= Epspar) {
432 // Ici les elements de I1 et I2 conviennent .
433 if (IsAdjustToFirstInterval)
434 {
435 Seq.Append(v1);
436 }
437 else
438 {
439 Seq.Append((v1 + v2) / 2);
440 }
441 ind1++;
442 ind2++;
443 }
444 else if (v1 < v2) {
445 // Ici l' element de I1 convient.
446 Seq.Append(v1);
447 ind1++;
448 }
449 else {
450 // Ici l' element de TABLE2 convient.
451 Seq.Append(v2);
452 ind2++;
453 }
454 }
455
456 if (ind1>I1.Upper()) {
457 //----- Ici I1 est epuise, on complete avec la fin de TABLE2 -------
458
459 for (; ind2<=I2.Upper(); ind2++) {
460 Seq.Append(I2(ind2));
461 }
462 }
463
464 if (ind2>I2.Upper()) {
465 //----- Ici I2 est epuise, on complete avec la fin de I1 -------
466 for (; ind1<=I1.Upper(); ind1++) {
467 Seq.Append(I1(ind1));
468 }
469 }
470 }
471
472
473 //=======================================================================
474 //function : EvalMaxParametricDistance
475 //purpose :
476 //=======================================================================
477
EvalMaxParametricDistance(const Adaptor3d_Curve & ACurve,const Adaptor3d_Curve & AReferenceCurve,const Standard_Real,const TColStd_Array1OfReal & Parameters,Standard_Real & MaxDistance)478 void GeomLib::EvalMaxParametricDistance(const Adaptor3d_Curve& ACurve,
479 const Adaptor3d_Curve& AReferenceCurve,
480 // const Standard_Real Tolerance,
481 const Standard_Real ,
482 const TColStd_Array1OfReal& Parameters,
483 Standard_Real& MaxDistance)
484 {
485 Standard_Integer ii ;
486
487 Standard_Real max_squared = 0.0e0,
488 // tolerance_squared,
489 local_distance_squared ;
490
491 // tolerance_squared = Tolerance * Tolerance ;
492 gp_Pnt Point1 ;
493 gp_Pnt Point2 ;
494 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
495 ACurve.D0(Parameters(ii),
496 Point1) ;
497 AReferenceCurve.D0(Parameters(ii),
498 Point2) ;
499 local_distance_squared =
500 Point1.SquareDistance (Point2) ;
501 max_squared = Max(max_squared,local_distance_squared) ;
502 }
503 if (max_squared > 0.0e0) {
504 MaxDistance = sqrt(max_squared) ;
505 }
506 else {
507 MaxDistance = 0.0e0 ;
508 }
509
510 }
511 //=======================================================================
512 //function : EvalMaxDistanceAlongParameter
513 //purpose :
514 //=======================================================================
515
EvalMaxDistanceAlongParameter(const Adaptor3d_Curve & ACurve,const Adaptor3d_Curve & AReferenceCurve,const Standard_Real Tolerance,const TColStd_Array1OfReal & Parameters,Standard_Real & MaxDistance)516 void GeomLib::EvalMaxDistanceAlongParameter(const Adaptor3d_Curve& ACurve,
517 const Adaptor3d_Curve& AReferenceCurve,
518 const Standard_Real Tolerance,
519 const TColStd_Array1OfReal& Parameters,
520 Standard_Real& MaxDistance)
521 {
522 Standard_Integer ii ;
523 Standard_Real max_squared = 0.0e0,
524 tolerance_squared = Tolerance * Tolerance,
525 other_parameter,
526 para_tolerance,
527 local_distance_squared ;
528 gp_Pnt Point1 ;
529 gp_Pnt Point2 ;
530
531
532
533 para_tolerance =
534 AReferenceCurve.Resolution(Tolerance) ;
535 other_parameter = Parameters(Parameters.Lower()) ;
536 ACurve.D0(other_parameter,
537 Point1) ;
538 Extrema_LocateExtPC a_projector(Point1,
539 AReferenceCurve,
540 other_parameter,
541 para_tolerance) ;
542 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
543 ACurve.D0(Parameters(ii),
544 Point1) ;
545 AReferenceCurve.D0(Parameters(ii),
546 Point2) ;
547 local_distance_squared =
548 Point1.SquareDistance (Point2) ;
549
550 local_distance_squared =
551 Point1.SquareDistance (Point2) ;
552
553
554 if (local_distance_squared > tolerance_squared) {
555
556
557 a_projector.Perform(Point1,
558 other_parameter) ;
559 if (a_projector.IsDone()) {
560 other_parameter =
561 a_projector.Point().Parameter() ;
562 AReferenceCurve.D0(other_parameter,
563 Point2) ;
564 local_distance_squared =
565 Point1.SquareDistance (Point2) ;
566 }
567 else {
568 local_distance_squared = 0.0e0 ;
569 other_parameter = Parameters(ii) ;
570 }
571 }
572 else {
573 other_parameter = Parameters(ii) ;
574 }
575
576
577 max_squared = Max(max_squared,local_distance_squared) ;
578 }
579 if (max_squared > tolerance_squared) {
580 MaxDistance = sqrt(max_squared) ;
581 }
582 else {
583 MaxDistance = Tolerance ;
584 }
585 }
586
587
588
589 // Aliases:
590
591 // Global data definitions:
592
593 // Methods :
594
595
596 //=======================================================================
597 //function : To3d
598 //purpose :
599 //=======================================================================
600
Handle(Geom_Curve)601 Handle(Geom_Curve) GeomLib::To3d (const gp_Ax2& Position,
602 const Handle(Geom2d_Curve)& Curve2d ) {
603 Handle(Geom_Curve) Curve3d;
604 Handle(Standard_Type) KindOfCurve = Curve2d->DynamicType();
605
606 if (KindOfCurve == STANDARD_TYPE (Geom2d_TrimmedCurve)) {
607 Handle(Geom2d_TrimmedCurve) Ct =
608 Handle(Geom2d_TrimmedCurve)::DownCast(Curve2d);
609 Standard_Real U1 = Ct->FirstParameter ();
610 Standard_Real U2 = Ct->LastParameter ();
611 Handle(Geom2d_Curve) CBasis2d = Ct->BasisCurve();
612 Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
613 Curve3d = new Geom_TrimmedCurve (CC, U1, U2);
614 }
615 else if (KindOfCurve == STANDARD_TYPE (Geom2d_OffsetCurve)) {
616 Handle(Geom2d_OffsetCurve) Co =
617 Handle(Geom2d_OffsetCurve)::DownCast(Curve2d);
618 Standard_Real Offset = Co->Offset();
619 Handle(Geom2d_Curve) CBasis2d = Co->BasisCurve();
620 Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
621 Curve3d = new Geom_OffsetCurve (CC, Offset, Position.Direction());
622 }
623 else if (KindOfCurve == STANDARD_TYPE (Geom2d_BezierCurve)) {
624 Handle(Geom2d_BezierCurve) CBez2d =
625 Handle(Geom2d_BezierCurve)::DownCast (Curve2d);
626 Standard_Integer Nbpoles = CBez2d->NbPoles ();
627 TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
628 CBez2d->Poles (Poles2d);
629 TColgp_Array1OfPnt Poles3d (1, Nbpoles);
630 for (Standard_Integer i = 1; i <= Nbpoles; i++) {
631 Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
632 }
633 Handle(Geom_BezierCurve) CBez3d;
634 if (CBez2d->IsRational()) {
635 TColStd_Array1OfReal TheWeights (1, Nbpoles);
636 CBez2d->Weights (TheWeights);
637 CBez3d = new Geom_BezierCurve (Poles3d, TheWeights);
638 }
639 else {
640 CBez3d = new Geom_BezierCurve (Poles3d);
641 }
642 Curve3d = CBez3d;
643 }
644 else if (KindOfCurve == STANDARD_TYPE (Geom2d_BSplineCurve)) {
645 Handle(Geom2d_BSplineCurve) CBSpl2d =
646 Handle(Geom2d_BSplineCurve)::DownCast (Curve2d);
647 Standard_Integer Nbpoles = CBSpl2d->NbPoles ();
648 Standard_Integer Nbknots = CBSpl2d->NbKnots ();
649 Standard_Integer TheDegree = CBSpl2d->Degree ();
650 Standard_Boolean IsPeriodic = CBSpl2d->IsPeriodic();
651 TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
652 CBSpl2d->Poles (Poles2d);
653 TColgp_Array1OfPnt Poles3d (1, Nbpoles);
654 for (Standard_Integer i = 1; i <= Nbpoles; i++) {
655 Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
656 }
657 TColStd_Array1OfReal TheKnots (1, Nbknots);
658 TColStd_Array1OfInteger TheMults (1, Nbknots);
659 CBSpl2d->Knots (TheKnots);
660 CBSpl2d->Multiplicities (TheMults);
661 Handle(Geom_BSplineCurve) CBSpl3d;
662 if (CBSpl2d->IsRational()) {
663 TColStd_Array1OfReal TheWeights (1, Nbpoles);
664 CBSpl2d->Weights (TheWeights);
665 CBSpl3d = new Geom_BSplineCurve (Poles3d, TheWeights, TheKnots, TheMults, TheDegree, IsPeriodic);
666 }
667 else {
668 CBSpl3d = new Geom_BSplineCurve (Poles3d, TheKnots, TheMults, TheDegree, IsPeriodic);
669 }
670 Curve3d = CBSpl3d;
671 }
672 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Line)) {
673 Handle(Geom2d_Line) Line2d = Handle(Geom2d_Line)::DownCast (Curve2d);
674 gp_Lin2d L2d = Line2d->Lin2d();
675 gp_Lin L3d = ElCLib::To3d (Position, L2d);
676 Handle(Geom_Line) GeomL3d = new Geom_Line (L3d);
677 Curve3d = GeomL3d;
678 }
679 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Circle)) {
680 Handle(Geom2d_Circle) Circle2d =
681 Handle(Geom2d_Circle)::DownCast (Curve2d);
682 gp_Circ2d C2d = Circle2d->Circ2d();
683 gp_Circ C3d = ElCLib::To3d (Position, C2d);
684 Handle(Geom_Circle) GeomC3d = new Geom_Circle (C3d);
685 Curve3d = GeomC3d;
686 }
687 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Ellipse)) {
688 Handle(Geom2d_Ellipse) Ellipse2d =
689 Handle(Geom2d_Ellipse)::DownCast (Curve2d);
690 gp_Elips2d E2d = Ellipse2d->Elips2d ();
691 gp_Elips E3d = ElCLib::To3d (Position, E2d);
692 Handle(Geom_Ellipse) GeomE3d = new Geom_Ellipse (E3d);
693 Curve3d = GeomE3d;
694 }
695 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Parabola)) {
696 Handle(Geom2d_Parabola) Parabola2d =
697 Handle(Geom2d_Parabola)::DownCast (Curve2d);
698 gp_Parab2d Prb2d = Parabola2d->Parab2d ();
699 gp_Parab Prb3d = ElCLib::To3d (Position, Prb2d);
700 Handle(Geom_Parabola) GeomPrb3d = new Geom_Parabola (Prb3d);
701 Curve3d = GeomPrb3d;
702 }
703 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Hyperbola)) {
704 Handle(Geom2d_Hyperbola) Hyperbola2d =
705 Handle(Geom2d_Hyperbola)::DownCast (Curve2d);
706 gp_Hypr2d H2d = Hyperbola2d->Hypr2d ();
707 gp_Hypr H3d = ElCLib::To3d (Position, H2d);
708 Handle(Geom_Hyperbola) GeomH3d = new Geom_Hyperbola (H3d);
709 Curve3d = GeomH3d;
710 }
711 else {
712 throw Standard_NotImplemented();
713 }
714
715 return Curve3d;
716 }
717
718
719
720 //=======================================================================
721 //function : GTransform
722 //purpose :
723 //=======================================================================
724
Handle(Geom2d_Curve)725 Handle(Geom2d_Curve) GeomLib::GTransform(const Handle(Geom2d_Curve)& Curve,
726 const gp_GTrsf2d& GTrsf)
727 {
728 gp_TrsfForm Form = GTrsf.Form();
729
730 if ( Form != gp_Other) {
731
732 // Alors, la GTrsf est en fait une Trsf.
733 // La geometrie des courbes sera alors inchangee.
734
735 Handle(Geom2d_Curve) C =
736 Handle(Geom2d_Curve)::DownCast(Curve->Transformed(GTrsf.Trsf2d()));
737 return C;
738 }
739 else {
740
741 // Alors, la GTrsf est une other Transformation.
742 // La geometrie des courbes est alors changee, et les conics devront
743 // etre converties en BSplines.
744
745 Handle(Standard_Type) TheType = Curve->DynamicType();
746
747 if ( TheType == STANDARD_TYPE(Geom2d_TrimmedCurve)) {
748
749 // On va recurer sur la BasisCurve
750
751 Handle(Geom2d_TrimmedCurve) C =
752 Handle(Geom2d_TrimmedCurve)::DownCast(Curve->Copy());
753
754 Handle(Standard_Type) TheBasisType = (C->BasisCurve())->DynamicType();
755
756 if (TheBasisType == STANDARD_TYPE(Geom2d_BSplineCurve) ||
757 TheBasisType == STANDARD_TYPE(Geom2d_BezierCurve) ) {
758
759 // Dans ces cas le parametrage est conserve sur la courbe transformee
760 // on peut donc la trimmer avec les parametres de la courbe de base.
761
762 Standard_Real U1 = C->FirstParameter();
763 Standard_Real U2 = C->LastParameter();
764
765 Handle(Geom2d_TrimmedCurve) result =
766 new Geom2d_TrimmedCurve(GTransform(C->BasisCurve(), GTrsf), U1,U2);
767 return result;
768 }
769 else if ( TheBasisType == STANDARD_TYPE(Geom2d_Line)) {
770
771 // Dans ce cas, le parametrage n`est plus conserve.
772 // Il faut recalculer les parametres de Trimming sur la courbe
773 // resultante. ( Calcul par projection ( ElCLib) des points debut
774 // et fin transformes)
775
776 Handle(Geom2d_Line) L =
777 Handle(Geom2d_Line)::DownCast(GTransform(C->BasisCurve(), GTrsf));
778 gp_Lin2d Lin = L->Lin2d();
779
780 gp_Pnt2d P1 = C->StartPoint();
781 gp_Pnt2d P2 = C->EndPoint();
782 P1.SetXY(GTrsf.Transformed(P1.XY()));
783 P2.SetXY(GTrsf.Transformed(P2.XY()));
784 Standard_Real U1 = ElCLib::Parameter(Lin,P1);
785 Standard_Real U2 = ElCLib::Parameter(Lin,P2);
786
787 Handle(Geom2d_TrimmedCurve) result =
788 new Geom2d_TrimmedCurve(L,U1,U2);
789 return result;
790 }
791 else if (TheBasisType == STANDARD_TYPE(Geom2d_Circle) ||
792 TheBasisType == STANDARD_TYPE(Geom2d_Ellipse) ||
793 TheBasisType == STANDARD_TYPE(Geom2d_Parabola) ||
794 TheBasisType == STANDARD_TYPE(Geom2d_Hyperbola) ) {
795
796 // Dans ces cas, la geometrie de la courbe n`est pas conservee
797 // on la convertir en BSpline avant de lui appliquer la Trsf.
798
799 Handle(Geom2d_BSplineCurve) BS =
800 Geom2dConvert::CurveToBSplineCurve(C);
801 return GTransform(BS,GTrsf);
802 }
803 else {
804
805 // La transformee d`une OffsetCurve vaut ????? Sais pas faire !!
806
807 Handle(Geom2d_Curve) dummy;
808 return dummy;
809 }
810 }
811 else if ( TheType == STANDARD_TYPE(Geom2d_Line)) {
812
813 Handle(Geom2d_Line) L =
814 Handle(Geom2d_Line)::DownCast(Curve->Copy());
815 gp_Lin2d Lin = L->Lin2d();
816 gp_Pnt2d P = Lin.Location();
817 gp_Pnt2d PP = L->Value(10.); // pourquoi pas !!
818 P.SetXY(GTrsf.Transformed(P.XY()));
819 PP.SetXY(GTrsf.Transformed(PP.XY()));
820 L->SetLocation(P);
821 gp_Vec2d V(P,PP);
822 L->SetDirection(gp_Dir2d(V));
823 return L;
824 }
825 else if ( TheType == STANDARD_TYPE(Geom2d_BezierCurve)) {
826
827 // Les GTrsf etant des operation lineaires, la transformee d`une courbe
828 // a poles est la courbe dont les poles sont la transformee des poles
829 // de la courbe de base.
830
831 Handle(Geom2d_BezierCurve) C =
832 Handle(Geom2d_BezierCurve)::DownCast(Curve->Copy());
833 Standard_Integer NbPoles = C->NbPoles();
834 TColgp_Array1OfPnt2d Poles(1,NbPoles);
835 C->Poles(Poles);
836 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
837 Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
838 C->SetPole(i,Poles(i));
839 }
840 return C;
841 }
842 else if ( TheType == STANDARD_TYPE(Geom2d_BSplineCurve)) {
843
844 // Voir commentaire pour les Bezier.
845
846 Handle(Geom2d_BSplineCurve) C =
847 Handle(Geom2d_BSplineCurve)::DownCast(Curve->Copy());
848 Standard_Integer NbPoles = C->NbPoles();
849 TColgp_Array1OfPnt2d Poles(1,NbPoles);
850 C->Poles(Poles);
851 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
852 Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
853 C->SetPole(i,Poles(i));
854 }
855 return C;
856 }
857 else if ( TheType == STANDARD_TYPE(Geom2d_Circle) ||
858 TheType == STANDARD_TYPE(Geom2d_Ellipse) ) {
859
860 // Dans ces cas, la geometrie de la courbe n`est pas conservee
861 // on la convertir en BSpline avant de lui appliquer la Trsf.
862
863 Handle(Geom2d_BSplineCurve) C =
864 Geom2dConvert::CurveToBSplineCurve(Curve);
865 return GTransform(C, GTrsf);
866 }
867 else if ( TheType == STANDARD_TYPE(Geom2d_Parabola) ||
868 TheType == STANDARD_TYPE(Geom2d_Hyperbola) ||
869 TheType == STANDARD_TYPE(Geom2d_OffsetCurve) ) {
870
871 // On ne sait pas faire : return a null Handle;
872
873 Handle(Geom2d_Curve) dummy;
874 return dummy;
875 }
876 }
877
878 Handle(Geom2d_Curve) WNT__; // portage Windows.
879 return WNT__;
880 }
881
882
883 //=======================================================================
884 //function : SameRange
885 //purpose :
886 //=======================================================================
SameRange(const Standard_Real Tolerance,const Handle (Geom2d_Curve)& CurvePtr,const Standard_Real FirstOnCurve,const Standard_Real LastOnCurve,const Standard_Real RequestedFirst,const Standard_Real RequestedLast,Handle (Geom2d_Curve)& NewCurvePtr)887 void GeomLib::SameRange(const Standard_Real Tolerance,
888 const Handle(Geom2d_Curve)& CurvePtr,
889 const Standard_Real FirstOnCurve,
890 const Standard_Real LastOnCurve,
891 const Standard_Real RequestedFirst,
892 const Standard_Real RequestedLast,
893 Handle(Geom2d_Curve)& NewCurvePtr)
894 {
895 if(CurvePtr.IsNull()) throw Standard_Failure();
896 if (Abs(LastOnCurve - RequestedLast) <= Tolerance &&
897 Abs(FirstOnCurve - RequestedFirst) <= Tolerance)
898 {
899 NewCurvePtr = CurvePtr;
900 return;
901 }
902
903 // the parametrisation length must at least be the same.
904 if (Abs(LastOnCurve - FirstOnCurve - RequestedLast + RequestedFirst)
905 <= Tolerance)
906 {
907 if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Line)))
908 {
909 Handle(Geom2d_Line) Line =
910 Handle(Geom2d_Line)::DownCast(CurvePtr->Copy());
911 Standard_Real dU = FirstOnCurve - RequestedFirst;
912 gp_Dir2d D = Line->Direction() ;
913 Line->Translate(dU * gp_Vec2d(D));
914 NewCurvePtr = Line;
915 }
916 else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Circle)))
917 {
918 gp_Trsf2d Trsf;
919 NewCurvePtr = Handle(Geom2d_Curve)::DownCast(CurvePtr->Copy());
920 Handle(Geom2d_Circle) Circ =
921 Handle(Geom2d_Circle)::DownCast(NewCurvePtr);
922 gp_Pnt2d P = Circ->Location();
923 Standard_Real dU;
924 if (Circ->Circ2d().IsDirect()) {
925 dU = FirstOnCurve - RequestedFirst;
926 }
927 else {
928 dU = RequestedFirst - FirstOnCurve;
929 }
930 Trsf.SetRotation(P,dU);
931 NewCurvePtr->Transform(Trsf) ;
932 }
933 else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
934 {
935 Handle(Geom2d_TrimmedCurve) TC =
936 Handle(Geom2d_TrimmedCurve)::DownCast(CurvePtr);
937 GeomLib::SameRange(Tolerance,
938 TC->BasisCurve(),
939 FirstOnCurve , LastOnCurve,
940 RequestedFirst, RequestedLast,
941 NewCurvePtr);
942 NewCurvePtr = new Geom2d_TrimmedCurve( NewCurvePtr, RequestedFirst, RequestedLast );
943 }
944 //
945 // attention a des problemes de limitation : utiliser le MEME test que dans
946 // Geom2d_TrimmedCurve::SetTrim car sinon comme on risque de relimite sur
947 // RequestedFirst et RequestedLast on aura un probleme
948 //
949 //
950 else if (Abs(LastOnCurve - FirstOnCurve) > Precision::PConfusion() ||
951 Abs(RequestedLast + RequestedFirst) > Precision::PConfusion())
952 {
953
954 Handle(Geom2d_TrimmedCurve) TC =
955 new Geom2d_TrimmedCurve(CurvePtr,FirstOnCurve,LastOnCurve);
956
957 Handle(Geom2d_BSplineCurve) BS =
958 Geom2dConvert::CurveToBSplineCurve(TC);
959 TColStd_Array1OfReal Knots(1,BS->NbKnots());
960 BS->Knots(Knots);
961
962 BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
963
964 BS->SetKnots(Knots);
965 NewCurvePtr = BS;
966 }
967 }
968 else
969 { // On segmente le resultat
970 Handle(Geom2d_TrimmedCurve) TC;
971 Handle(Geom2d_Curve) aCCheck = CurvePtr;
972
973 if(aCCheck->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
974 {
975 aCCheck = Handle(Geom2d_TrimmedCurve)::DownCast(aCCheck)->BasisCurve();
976 }
977
978 if(aCCheck->IsPeriodic())
979 {
980 if(Abs(LastOnCurve - FirstOnCurve) > Precision::PConfusion())
981 {
982 TC = new Geom2d_TrimmedCurve( CurvePtr, FirstOnCurve, LastOnCurve );
983 }
984 else
985 {
986 TC = new Geom2d_TrimmedCurve( CurvePtr, CurvePtr->FirstParameter(), CurvePtr->LastParameter() );
987 }
988 }
989 else
990 {
991 const Standard_Real Udeb = Max(CurvePtr->FirstParameter(), FirstOnCurve);
992 const Standard_Real Ufin = Min(CurvePtr->LastParameter(), LastOnCurve);
993 if(Abs(Ufin - Udeb) > Precision::PConfusion())
994 {
995 TC = new Geom2d_TrimmedCurve( CurvePtr, Udeb, Ufin );
996 }
997 else
998 {
999 TC = new Geom2d_TrimmedCurve( CurvePtr, CurvePtr->FirstParameter(), CurvePtr->LastParameter());
1000 }
1001 }
1002
1003 //
1004 Handle(Geom2d_BSplineCurve) BS =
1005 Geom2dConvert::CurveToBSplineCurve(TC);
1006 TColStd_Array1OfReal Knots(1,BS->NbKnots());
1007 BS->Knots(Knots);
1008
1009 BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
1010
1011 BS->SetKnots(Knots);
1012 NewCurvePtr = BS;
1013 }
1014 }
1015
1016 //=======================================================================
1017 //class : GeomLib_CurveOnSurfaceEvaluator
1018 //purpose: The evaluator for the Curve 3D building
1019 //=======================================================================
1020
1021 class GeomLib_CurveOnSurfaceEvaluator : public AdvApprox_EvaluatorFunction
1022 {
1023 public:
GeomLib_CurveOnSurfaceEvaluator(Adaptor3d_CurveOnSurface & theCurveOnSurface,Standard_Real theFirst,Standard_Real theLast)1024 GeomLib_CurveOnSurfaceEvaluator (Adaptor3d_CurveOnSurface& theCurveOnSurface,
1025 Standard_Real theFirst, Standard_Real theLast)
1026 : CurveOnSurface(theCurveOnSurface), FirstParam(theFirst), LastParam(theLast) {}
1027
1028 virtual void Evaluate (Standard_Integer *Dimension,
1029 Standard_Real StartEnd[2],
1030 Standard_Real *Parameter,
1031 Standard_Integer *DerivativeRequest,
1032 Standard_Real *Result, // [Dimension]
1033 Standard_Integer *ErrorCode);
1034
1035 private:
1036 Adaptor3d_CurveOnSurface& CurveOnSurface;
1037 Standard_Real FirstParam;
1038 Standard_Real LastParam;
1039
1040 Handle(Adaptor3d_Curve) TrimCurve;
1041 };
1042
Evaluate(Standard_Integer *,Standard_Real DebutFin[2],Standard_Real * Parameter,Standard_Integer * DerivativeRequest,Standard_Real * Result,Standard_Integer * ReturnCode)1043 void GeomLib_CurveOnSurfaceEvaluator::Evaluate (Standard_Integer *,/*Dimension*/
1044 Standard_Real DebutFin[2],
1045 Standard_Real *Parameter,
1046 Standard_Integer *DerivativeRequest,
1047 Standard_Real *Result,// [Dimension]
1048 Standard_Integer *ReturnCode)
1049 {
1050 gp_Pnt Point;
1051
1052 //Gestion des positionnements gauche / droite
1053 if ((DebutFin[0] != FirstParam) || (DebutFin[1] != LastParam))
1054 {
1055 TrimCurve = CurveOnSurface.Trim(DebutFin[0], DebutFin[1], Precision::PConfusion());
1056 FirstParam = DebutFin[0];
1057 LastParam = DebutFin[1];
1058 }
1059
1060 //Positionemment
1061 if (*DerivativeRequest == 0)
1062 {
1063 TrimCurve->D0((*Parameter), Point) ;
1064
1065 for (Standard_Integer ii = 0 ; ii < 3 ; ii++)
1066 Result[ii] = Point.Coord(ii + 1);
1067 }
1068 if (*DerivativeRequest == 1)
1069 {
1070 gp_Vec Vector;
1071 TrimCurve->D1((*Parameter), Point, Vector);
1072 for (Standard_Integer ii = 0 ; ii < 3 ; ii++)
1073 Result[ii] = Vector.Coord(ii + 1) ;
1074 }
1075 if (*DerivativeRequest == 2)
1076 {
1077 gp_Vec Vector, VecBis;
1078 TrimCurve->D2((*Parameter), Point, VecBis, Vector);
1079 for (Standard_Integer ii = 0 ; ii < 3 ; ii++)
1080 Result[ii] = Vector.Coord(ii + 1) ;
1081 }
1082 ReturnCode[0] = 0;
1083 }
1084
1085 //=======================================================================
1086 //function : BuildCurve3d
1087 //purpose :
1088 //=======================================================================
1089
BuildCurve3d(const Standard_Real Tolerance,Adaptor3d_CurveOnSurface & Curve,const Standard_Real FirstParameter,const Standard_Real LastParameter,Handle (Geom_Curve)& NewCurvePtr,Standard_Real & MaxDeviation,Standard_Real & AverageDeviation,const GeomAbs_Shape Continuity,const Standard_Integer MaxDegree,const Standard_Integer MaxSegment)1090 void GeomLib::BuildCurve3d(const Standard_Real Tolerance,
1091 Adaptor3d_CurveOnSurface& Curve,
1092 const Standard_Real FirstParameter,
1093 const Standard_Real LastParameter,
1094 Handle(Geom_Curve)& NewCurvePtr,
1095 Standard_Real& MaxDeviation,
1096 Standard_Real& AverageDeviation,
1097 const GeomAbs_Shape Continuity,
1098 const Standard_Integer MaxDegree,
1099 const Standard_Integer MaxSegment)
1100
1101 {
1102
1103
1104 MaxDeviation = 0.0e0 ;
1105 AverageDeviation = 0.0e0 ;
1106 Handle(GeomAdaptor_Surface) geom_adaptor_surface_ptr (Handle(GeomAdaptor_Surface)::DownCast(Curve.GetSurface()) );
1107 Handle(Geom2dAdaptor_Curve) geom_adaptor_curve_ptr (Handle(Geom2dAdaptor_Curve)::DownCast(Curve.GetCurve()) );
1108
1109 if (! geom_adaptor_curve_ptr.IsNull() &&
1110 ! geom_adaptor_surface_ptr.IsNull()) {
1111 Handle(Geom_Plane) P ;
1112 const GeomAdaptor_Surface& geom_surface = *geom_adaptor_surface_ptr;
1113
1114 Handle(Geom_RectangularTrimmedSurface) RT = Handle(Geom_RectangularTrimmedSurface)::DownCast(geom_surface.Surface());
1115 if ( RT.IsNull()) {
1116 P = Handle(Geom_Plane)::DownCast(geom_surface.Surface());
1117 }
1118 else {
1119 P = Handle(Geom_Plane)::DownCast(RT->BasisSurface());
1120 }
1121
1122
1123 if (! P.IsNull()) {
1124 // compute the 3d curve
1125 gp_Ax2 axes = P->Position().Ax2();
1126 const Geom2dAdaptor_Curve& geom2d_curve = *geom_adaptor_curve_ptr;
1127 NewCurvePtr =
1128 GeomLib::To3d(axes,
1129 geom2d_curve.Curve());
1130 return;
1131
1132 }
1133
1134 Handle(Adaptor2d_Curve2d) TrimmedC2D = geom_adaptor_curve_ptr->Trim (FirstParameter, LastParameter, Precision::PConfusion());
1135
1136 Standard_Boolean isU, isForward;
1137 Standard_Real aParam;
1138 if (isIsoLine(TrimmedC2D, isU, aParam, isForward))
1139 {
1140 NewCurvePtr = buildC3dOnIsoLine (TrimmedC2D, geom_adaptor_surface_ptr, FirstParameter, LastParameter, Tolerance, isU, aParam, isForward);
1141 if (!NewCurvePtr.IsNull())
1142 {
1143 return;
1144 }
1145 }
1146 }
1147
1148 //
1149 // Entree
1150 //
1151 Handle(TColStd_HArray1OfReal) Tolerance1DPtr,Tolerance2DPtr;
1152 Handle(TColStd_HArray1OfReal) Tolerance3DPtr =
1153 new TColStd_HArray1OfReal(1,1) ;
1154 Tolerance3DPtr->SetValue(1,Tolerance);
1155
1156 // Recherche des discontinuitees
1157 Standard_Integer NbIntervalC2 = Curve.NbIntervals(GeomAbs_C2);
1158 TColStd_Array1OfReal Param_de_decoupeC2 (1, NbIntervalC2+1);
1159 Curve.Intervals(Param_de_decoupeC2, GeomAbs_C2);
1160
1161 Standard_Integer NbIntervalC3 = Curve.NbIntervals(GeomAbs_C3);
1162 TColStd_Array1OfReal Param_de_decoupeC3 (1, NbIntervalC3+1);
1163 Curve.Intervals(Param_de_decoupeC3, GeomAbs_C3);
1164
1165 // Note extension of the parameteric range
1166 // Pour forcer le Trim au premier appel de l'evaluateur
1167 GeomLib_CurveOnSurfaceEvaluator ev (Curve, FirstParameter - 1., LastParameter + 1.);
1168
1169 // Approximation avec decoupe preferentiel
1170 AdvApprox_PrefAndRec Preferentiel(Param_de_decoupeC2,
1171 Param_de_decoupeC3);
1172 AdvApprox_ApproxAFunction anApproximator(0,
1173 0,
1174 1,
1175 Tolerance1DPtr,
1176 Tolerance2DPtr,
1177 Tolerance3DPtr,
1178 FirstParameter,
1179 LastParameter,
1180 Continuity,
1181 MaxDegree,
1182 MaxSegment,
1183 ev,
1184 // CurveOnSurfaceEvaluator,
1185 Preferentiel) ;
1186
1187 if (anApproximator.HasResult()) {
1188 GeomLib_MakeCurvefromApprox
1189 aCurveBuilder(anApproximator) ;
1190
1191 Handle(Geom_BSplineCurve) aCurvePtr =
1192 aCurveBuilder.Curve(1) ;
1193 // On rend les resultats de l'approx
1194 MaxDeviation = anApproximator.MaxError(3,1) ;
1195 AverageDeviation = anApproximator.AverageError(3,1) ;
1196 NewCurvePtr = aCurvePtr ;
1197 }
1198 }
1199
1200 //=======================================================================
1201 //function : AdjustExtremity
1202 //purpose :
1203 //=======================================================================
1204
AdjustExtremity(Handle (Geom_BoundedCurve)& Curve,const gp_Pnt & P1,const gp_Pnt & P2,const gp_Vec & T1,const gp_Vec & T2)1205 void GeomLib::AdjustExtremity(Handle(Geom_BoundedCurve)& Curve,
1206 const gp_Pnt& P1,
1207 const gp_Pnt& P2,
1208 const gp_Vec& T1,
1209 const gp_Vec& T2)
1210 {
1211 // il faut Convertir l'entree (en preservant si possible le parametrage)
1212 Handle(Geom_BSplineCurve) aIn, aDef;
1213 aIn = GeomConvert::CurveToBSplineCurve(Curve, Convert_QuasiAngular);
1214
1215 Standard_Integer ii, jj;
1216 gp_Pnt P;
1217 gp_Vec V, Vtan, DV;
1218 TColgp_Array1OfPnt PolesDef(1,4), Coeffs(1,4);
1219 TColStd_Array1OfReal FK(1, 8);
1220 TColStd_Array1OfReal Ti(1, 4);
1221 TColStd_Array1OfInteger Contact(1, 4);
1222
1223 Ti(1) = Ti(2) = aIn->FirstParameter();
1224 Ti(3) = Ti(4) = aIn->LastParameter();
1225 Contact(1) = Contact(3) = 0;
1226 Contact(2) = Contact(4) = 1;
1227 for (ii=1; ii<=4; ii++) {
1228 FK(ii) = aIn->FirstParameter();
1229 FK(ii) = aIn->LastParameter();
1230 }
1231
1232 // Calculs des contraintes de deformations
1233 aIn->D1(Ti(1), P, V);
1234 PolesDef(1).ChangeCoord() = P1.XYZ()-P.XYZ();
1235 Vtan = T1;
1236 Vtan.Normalize();
1237 DV = Vtan * (Vtan * V) - V;
1238 PolesDef(2).ChangeCoord() = (Ti(4)-Ti(1))*DV.XYZ();
1239
1240 aIn->D1(Ti(4), P, V);
1241 PolesDef(3).ChangeCoord() = P2.XYZ()-P.XYZ();
1242 Vtan = T2;
1243 Vtan.Normalize();
1244 DV = Vtan * (Vtan * V) - V;
1245 PolesDef(4).ChangeCoord() = (Ti(4)-Ti(1))* DV.XYZ();
1246
1247 // Interpolation des contraintes
1248 math_Matrix Mat(1, 4, 1, 4);
1249 if (!PLib::HermiteCoefficients(0., 1., 1, 1, Mat))
1250 throw Standard_ConstructionError();
1251
1252 for (jj=1; jj<=4; jj++) {
1253 gp_XYZ aux(0.,0.,0.);
1254 for (ii=1; ii<=4; ii++) {
1255 aux.SetLinearForm(Mat(ii,jj), PolesDef(ii).XYZ(), aux);
1256 }
1257 Coeffs(jj).SetXYZ(aux);
1258 }
1259
1260 PLib::CoefficientsPoles(Coeffs, PLib::NoWeights(),
1261 PolesDef, PLib::NoWeights());
1262
1263 // Ajout de la deformation
1264 TColStd_Array1OfReal K(1, 2);
1265 TColStd_Array1OfInteger M(1, 2);
1266 K(1) = Ti(1);
1267 K(2) = Ti(4);
1268 M.Init(4);
1269
1270 aDef = new (Geom_BSplineCurve) (PolesDef, K, M, 3);
1271 if (aIn->Degree() < 3) aIn->IncreaseDegree(3);
1272 else aDef->IncreaseDegree(aIn->Degree());
1273
1274 for (ii=2; ii<aIn->NbKnots(); ii++) {
1275 aDef->InsertKnot(aIn->Knot(ii), aIn->Multiplicity(ii));
1276 }
1277
1278 if (aDef->NbPoles() != aIn->NbPoles())
1279 throw Standard_ConstructionError("Inconsistent poles's number");
1280
1281 for (ii=1; ii<=aDef->NbPoles(); ii++) {
1282 P = aIn->Pole(ii);
1283 P.ChangeCoord() += aDef->Pole(ii).XYZ();
1284 aIn->SetPole(ii, P);
1285 }
1286 Curve = aIn;
1287 }
1288 //=======================================================================
1289 //function : ExtendCurveToPoint
1290 //purpose :
1291 //=======================================================================
1292
ExtendCurveToPoint(Handle (Geom_BoundedCurve)& Curve,const gp_Pnt & Point,const Standard_Integer Continuity,const Standard_Boolean After)1293 void GeomLib::ExtendCurveToPoint(Handle(Geom_BoundedCurve)& Curve,
1294 const gp_Pnt& Point,
1295 const Standard_Integer Continuity,
1296 const Standard_Boolean After)
1297 {
1298 if(Continuity < 1 || Continuity > 3) return;
1299 Standard_Integer size = Continuity + 2;
1300 Standard_Real Ubord, Tol=1.e-6;
1301 math_Matrix MatCoefs(1,size, 1,size);
1302 Standard_Real Lambda, L1;
1303 Standard_Integer ii, jj;
1304 gp_Vec d1, d2, d3;
1305 gp_Pnt p0;
1306 // il faut Convertir l'entree (en preservant si possible le parametrage)
1307 GeomConvert_CompCurveToBSplineCurve Concat(Curve, Convert_QuasiAngular);
1308
1309 // Les contraintes de constructions
1310 TColgp_Array1OfXYZ Cont(1,size);
1311 if (After) {
1312 Ubord = Curve->LastParameter();
1313
1314 }
1315 else {
1316 Ubord = Curve->FirstParameter();
1317 }
1318 PLib::HermiteCoefficients(0, 1, // Les Bornes
1319 Continuity, 0, // Les Ordres de contraintes
1320 MatCoefs);
1321
1322 Curve->D3(Ubord, p0, d1, d2, d3);
1323 if (!After) { // Inversion du parametrage
1324 d1 *= -1;
1325 d3 *= -1;
1326 }
1327
1328 L1 = p0.Distance(Point);
1329 if (L1 > Tol) {
1330 // Lambda est le ratio qu'il faut appliquer a la derive de la courbe
1331 // pour obtenir la derive du prolongement (fixe arbitrairement a la
1332 // longueur du segment bout de la courbe - point cible.
1333 // On essai d'avoir sur le prolongement la vitesse moyenne que l'on
1334 // a sur la courbe.
1335 gp_Vec daux;
1336 gp_Pnt pp;
1337 Standard_Real f= Curve->FirstParameter(), t, dt, norm;
1338 dt = (Curve->LastParameter()-f)/9;
1339 norm = d1.Magnitude();
1340 for (ii=1, t=f+dt; ii<=8; ii++, t+=dt) {
1341 Curve->D1(t, pp, daux);
1342 norm += daux.Magnitude();
1343 }
1344 norm /= 9;
1345 dt = d1.Magnitude() / norm;
1346 if ((dt<1.5) && (dt>0.75)) { // Le bord est dans la moyenne on le garde
1347 Lambda = ((Standard_Real)1) / Max (d1.Magnitude() / L1, Tol);
1348 }
1349 else {
1350 Lambda = ((Standard_Real)1) / Max (norm / L1, Tol);
1351 }
1352 }
1353 else {
1354 return; // Pas d'extension
1355 }
1356
1357 // Optimisation du Lambda
1358 math_Matrix Cons(1, 3, 1, size);
1359 Cons(1,1) = p0.X(); Cons(2,1) = p0.Y(); Cons(3,1) = p0.Z();
1360 Cons(1,2) = d1.X(); Cons(2,2) = d1.Y(); Cons(3,2) = d1.Z();
1361 Cons(1,size) = Point.X(); Cons(2,size) = Point.Y(); Cons(3,size) = Point.Z();
1362 if (Continuity >= 2) {
1363 Cons(1,3) = d2.X(); Cons(2,3) = d2.Y(); Cons(3,3) = d2.Z();
1364 }
1365 if (Continuity >= 3) {
1366 Cons(1,4) = d3.X(); Cons(2,4) = d3.Y(); Cons(3,4) = d3.Z();
1367 }
1368 ComputeLambda(Cons, MatCoefs, L1, Lambda);
1369
1370 // Construction dans la Base Polynomiale
1371 Cont(1) = p0.XYZ();
1372 Cont(2) = d1.XYZ() * Lambda;
1373 if(Continuity >= 2) Cont(3) = d2.XYZ() * Pow(Lambda,2);
1374 if(Continuity >= 3) Cont(4) = d3.XYZ() * Pow(Lambda,3);
1375 Cont(size) = Point.XYZ();
1376
1377
1378 TColgp_Array1OfPnt ExtrapPoles(1, size);
1379 TColgp_Array1OfPnt ExtraCoeffs(1, size);
1380
1381 gp_Pnt PNull(0.,0.,0.);
1382 ExtraCoeffs.Init(PNull);
1383 for (ii=1; ii<=size; ii++) {
1384 for (jj=1; jj<=size; jj++) {
1385 ExtraCoeffs(jj).ChangeCoord() += MatCoefs(ii,jj)*Cont(ii);
1386 }
1387 }
1388
1389 // Convertion Dans la Base de Bernstein
1390 PLib::CoefficientsPoles(ExtraCoeffs, PLib::NoWeights(),
1391 ExtrapPoles, PLib::NoWeights());
1392
1393 Handle(Geom_BezierCurve) Bezier = new (Geom_BezierCurve) (ExtrapPoles);
1394
1395 Standard_Real dist = ExtrapPoles(1).Distance(p0);
1396 Standard_Boolean Ok;
1397 Tol += dist;
1398
1399 // Concatenation
1400 Ok = Concat.Add(Bezier, Tol, After);
1401 if (!Ok) throw Standard_ConstructionError("ExtendCurveToPoint");
1402
1403 Curve = Concat.BSplineCurve();
1404 }
1405
1406
1407 //=======================================================================
1408 //function : ExtendKPart
1409 //purpose : Extension par longueur des surfaces cannonique
1410 //=======================================================================
1411 static Standard_Boolean
ExtendKPart(Handle (Geom_RectangularTrimmedSurface)& Surface,const Standard_Real Length,const Standard_Boolean InU,const Standard_Boolean After)1412 ExtendKPart(Handle(Geom_RectangularTrimmedSurface)& Surface,
1413 const Standard_Real Length,
1414 const Standard_Boolean InU,
1415 const Standard_Boolean After)
1416 {
1417
1418 if (Surface.IsNull()) return Standard_False;
1419
1420 Standard_Boolean Ok=Standard_True;
1421 Standard_Real Uf, Ul, Vf, Vl;
1422 Handle(Geom_Surface) Support = Surface->BasisSurface();
1423 GeomAbs_SurfaceType Type;
1424
1425 Surface->Bounds(Uf, Ul, Vf, Vl);
1426 GeomAdaptor_Surface AS(Surface);
1427 Type = AS.GetType();
1428
1429 if (InU) {
1430 switch(Type) {
1431 case GeomAbs_Plane :
1432 {
1433 if (After) Ul+=Length;
1434 else Uf-=Length;
1435 Surface = new (Geom_RectangularTrimmedSurface)
1436 (Support, Uf, Ul, Vf, Vl);
1437 break;
1438 }
1439
1440 default:
1441 Ok = Standard_False;
1442 }
1443 }
1444 else {
1445 switch(Type) {
1446 case GeomAbs_Plane :
1447 case GeomAbs_Cylinder :
1448 case GeomAbs_SurfaceOfExtrusion :
1449 {
1450 if (After) Vl+=Length;
1451 else Vf-=Length;
1452 Surface = new (Geom_RectangularTrimmedSurface)
1453 (Support, Uf, Ul, Vf, Vl);
1454 break;
1455 }
1456 default:
1457 Ok = Standard_False;
1458 }
1459 }
1460
1461 return Ok;
1462 }
1463
1464 //=======================================================================
1465 //function : ExtendSurfByLength
1466 //purpose :
1467 //=======================================================================
ExtendSurfByLength(Handle (Geom_BoundedSurface)& Surface,const Standard_Real Length,const Standard_Integer Continuity,const Standard_Boolean InU,const Standard_Boolean After)1468 void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
1469 const Standard_Real Length,
1470 const Standard_Integer Continuity,
1471 const Standard_Boolean InU,
1472 const Standard_Boolean After)
1473 {
1474 if(Continuity < 0 || Continuity > 3) return;
1475 Standard_Integer Cont = Continuity;
1476
1477 // Kpart ?
1478 Handle(Geom_RectangularTrimmedSurface) TS =
1479 Handle(Geom_RectangularTrimmedSurface)::DownCast (Surface);
1480 if (ExtendKPart(TS,Length, InU, After) ) {
1481 Surface = TS;
1482 return;
1483 }
1484
1485 // format BSplineSurface avec un degre suffisant pour la continuite voulue
1486 Handle(Geom_BSplineSurface) BS =
1487 Handle(Geom_BSplineSurface)::DownCast (Surface);
1488 if (BS.IsNull()) {
1489 //BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1490 Standard_Real Tol = Precision::Confusion(); //1.e-4;
1491 GeomAbs_Shape UCont = GeomAbs_C1, VCont = GeomAbs_C1;
1492 Standard_Integer degU = 14, degV = 14;
1493 Standard_Integer nmax = 16;
1494 Standard_Integer thePrec = 1;
1495 const Handle(Geom_Surface)& aSurf = Surface; // to resolve ambiguity
1496 GeomConvert_ApproxSurface theApprox(aSurf,Tol,UCont,VCont,degU,degV,nmax,thePrec);
1497 if (theApprox.HasResult())
1498 BS = theApprox.Surface();
1499 else
1500 BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1501 }
1502 if (InU&&(BS->UDegree()<Continuity+1))
1503 BS->IncreaseDegree(Continuity+1,BS->VDegree());
1504 if (!InU&&(BS->VDegree()<Continuity+1))
1505 BS->IncreaseDegree(BS->UDegree(),Continuity+1);
1506
1507 // si BS etait periodique dans le sens de l'extension, elle ne le sera plus
1508 if ( (InU&&(BS->IsUPeriodic())) || (!InU&&(BS->IsVPeriodic())) ) {
1509 Standard_Real U0,U1,V0,V1;
1510 BS->Bounds(U0,U1,V0,V1);
1511 BS->Segment(U0,U1,V0,V1);
1512 }
1513
1514
1515 // IFV Fix OCC bug 0022694 - wrong result extrapolating rational surfaces
1516 // Standard_Boolean rational = ( InU && BS->IsURational() )
1517 // || ( !InU && BS->IsVRational() ) ;
1518 Standard_Boolean rational = (BS->IsURational() || BS->IsVRational());
1519 Standard_Boolean NullWeight;
1520 Standard_Real EpsW = 10*Precision::PConfusion();
1521 Standard_Integer gap = 3;
1522 if ( rational ) gap++;
1523
1524
1525
1526 Standard_Integer Cdeg = 0, Cdim = 0, NbP = 0, Ksize = 0, Psize = 1;
1527 Standard_Integer ii, jj, ipole, Kount;
1528 Standard_Real Tbord, lambmin=Length;
1529 Standard_Real * Padr = NULL;
1530 Standard_Boolean Ok;
1531 Handle(TColStd_HArray1OfReal) FKnots, Point, lambda, Tgte, Poles;
1532
1533
1534
1535
1536 for (Kount=0, Ok=Standard_False; Kount<=2 && !Ok; Kount++) {
1537 // transformation de la surface en une BSpline non rationnelle a une variable
1538 // de degre UDegree ou VDegree et de dimension 3 ou 4 x NbVpoles ou NbUpoles
1539 // le nombre de poles egal a NbUpoles ou NbVpoles
1540 // ATTENTION : dans le cas rationnel, un point de coordonnees (x,y,z)
1541 // et de poids w devient un point de coordonnees (wx, wy, wz, w )
1542
1543
1544 if (InU) {
1545 Cdeg = BS->UDegree();
1546 NbP = BS->NbUPoles();
1547 Cdim = BS->NbVPoles() * gap;
1548 }
1549 else {
1550 Cdeg = BS->VDegree();
1551 NbP = BS->NbVPoles();
1552 Cdim = BS->NbUPoles() * gap;
1553 }
1554
1555 // les noeuds plats
1556 Ksize = NbP + Cdeg + 1;
1557 FKnots = new (TColStd_HArray1OfReal) (1,Ksize);
1558 if (InU)
1559 BS->UKnotSequence(FKnots->ChangeArray1());
1560 else
1561 BS->VKnotSequence(FKnots->ChangeArray1());
1562
1563 // le parametre du noeud de raccord
1564 if (After)
1565 Tbord = FKnots->Value(FKnots->Upper()-Cdeg);
1566 else
1567 Tbord = FKnots->Value(FKnots->Lower()+Cdeg);
1568
1569 // les poles
1570 Psize = Cdim * NbP;
1571 Poles = new (TColStd_HArray1OfReal) (1,Psize);
1572
1573 if (InU) {
1574 for (ii=1,ipole=1; ii<=NbP; ii++) {
1575 for (jj=1;jj<=BS->NbVPoles();jj++) {
1576 Poles->SetValue(ipole, BS->Pole(ii,jj).X());
1577 Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1578 Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1579 if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1580 ipole+=gap;
1581 }
1582 }
1583 }
1584 else {
1585 for (jj=1,ipole=1; jj<=NbP; jj++) {
1586 for (ii=1;ii<=BS->NbUPoles();ii++) {
1587 Poles->SetValue(ipole, BS->Pole(ii,jj).X());
1588 Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1589 Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1590 if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1591 ipole+=gap;
1592 }
1593 }
1594 }
1595 Padr = (Standard_Real *) &Poles->ChangeValue(1);
1596
1597 // calcul du point de raccord et de la tangente
1598 Point = new (TColStd_HArray1OfReal)(1,Cdim);
1599 Tgte = new (TColStd_HArray1OfReal)(1,Cdim);
1600 lambda = new (TColStd_HArray1OfReal)(1,Cdim);
1601
1602 Standard_Boolean periodic_flag = Standard_False ;
1603 Standard_Integer extrap_mode[2], derivative_request = Max(Continuity,1);
1604 extrap_mode[0] = extrap_mode[1] = Cdeg;
1605 TColStd_Array1OfReal Result(1, Cdim * (derivative_request+1)) ;
1606
1607 TColStd_Array1OfReal& tgte = Tgte->ChangeArray1();
1608 TColStd_Array1OfReal& point = Point->ChangeArray1();
1609 TColStd_Array1OfReal& lamb = lambda->ChangeArray1();
1610
1611 Standard_Real * Radr = (Standard_Real *) &Result(1) ;
1612
1613 BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
1614 Cdeg,FKnots->Array1(),Cdim,*Padr,*Radr);
1615 Ok = Standard_True;
1616 for (ii=1;ii<=Cdim;ii++) {
1617 point(ii) = Result(ii);
1618 tgte(ii) = Result(ii+Cdim);
1619 }
1620
1621 // calcul de la contrainte a atteindre
1622
1623 gp_Vec CurT, OldT;
1624
1625 Standard_Real NTgte, val, Tgtol = 1.e-12, OldN = 0.0;
1626 if (rational) {
1627 for (ii=gap;ii<=Cdim;ii+=gap) {
1628 tgte(ii) = 0.;
1629 }
1630 for (ii=gap;ii<=Cdim;ii+=gap) {
1631 CurT.SetCoord(tgte(ii-3),tgte(ii-2), tgte(ii-1));
1632 NTgte=CurT.Magnitude();
1633 if (NTgte>Tgtol) {
1634 val = Length/NTgte;
1635 // Attentions aux Cas ou le segment donne par les poles
1636 // est oppose au sens de la derive
1637 // Exemple: Certaine portions de tore.
1638 if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1639 Ok = Standard_False;
1640 }
1641
1642 lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = val;
1643 lamb(ii) = 0.;
1644 lambmin = Min(lambmin, val);
1645 }
1646 else {
1647 lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = 0.;
1648 lamb(ii) = 0.;
1649 }
1650 OldT = CurT;
1651 OldN = NTgte;
1652 }
1653 }
1654 else {
1655 for (ii=gap;ii<=Cdim;ii+=gap) {
1656 CurT.SetCoord(tgte(ii-2),tgte(ii-1), tgte(ii));
1657 NTgte=CurT.Magnitude();
1658 if (NTgte>Tgtol) {
1659 val = Length/NTgte;
1660 // Attentions aux Cas ou le segment donne par les poles
1661 // est oppose au sens de la derive
1662 // Exemple: Certaine portion de tore.
1663 if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1664 Ok = Standard_False;
1665 }
1666 lamb(ii) = lamb(ii-1) = lamb(ii-2) = val;
1667 lambmin = Min(lambmin, val);
1668 }
1669 else {
1670 lamb(ii) =lamb(ii-1) = lamb(ii-2) = 0.;
1671 }
1672 OldT = CurT;
1673 OldN = NTgte;
1674 }
1675 }
1676 if (!Ok && Kount<2) {
1677 // On augmente le degre de l'iso bord afin de rapprocher les poles de la surface
1678 // Et on ressaye
1679 if (InU) BS->IncreaseDegree(BS->UDegree(), BS->VDegree()+2);
1680 else BS->IncreaseDegree(BS->UDegree()+2, BS->VDegree());
1681 }
1682 }
1683
1684
1685 TColStd_Array1OfReal ConstraintPoint(1,Cdim);
1686 if (After) {
1687 for (ii=1;ii<=Cdim;ii++) {
1688 ConstraintPoint(ii) = Point->Value(ii) + lambda->Value(ii)*Tgte->Value(ii);
1689 }
1690 }
1691 else {
1692 for (ii=1;ii<=Cdim;ii++) {
1693 ConstraintPoint(ii) = Point->Value(ii) - lambda->Value(ii)*Tgte->Value(ii);
1694 }
1695 }
1696
1697 // cas particulier du rationnel
1698 if (rational) {
1699 for (ipole=1;ipole<=Psize;ipole+=gap) {
1700 Poles->ChangeValue(ipole) *= Poles->Value(ipole+3);
1701 Poles->ChangeValue(ipole+1) *= Poles->Value(ipole+3);
1702 Poles->ChangeValue(ipole+2) *= Poles->Value(ipole+3);
1703 }
1704 for (ii=1;ii<=Cdim;ii+=gap) {
1705 ConstraintPoint(ii) *= ConstraintPoint(ii+3);
1706 ConstraintPoint(ii+1) *= ConstraintPoint(ii+3);
1707 ConstraintPoint(ii+2) *= ConstraintPoint(ii+3);
1708 }
1709 }
1710
1711 // tableaux necessaires pour l'extension
1712 Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots = 0;
1713 TColStd_Array1OfReal FK(1, Ksize2) ;
1714 Standard_Real * FKRadr = &FK(1);
1715
1716 Standard_Integer Psize2 = Psize+Cdeg*Cdim;
1717 TColStd_Array1OfReal PRes(1, Psize2) ;
1718 Standard_Real * PRadr = &PRes(1);
1719 Standard_Real ww;
1720 Standard_Boolean ExtOk = Standard_False;
1721 Handle(TColgp_HArray2OfPnt) NewPoles;
1722 Handle(TColStd_HArray2OfReal) NewWeights;
1723
1724
1725 for (Kount=1; Kount<=5 && !ExtOk; Kount++) {
1726 // extension
1727 BSplCLib::TangExtendToConstraint(FKnots->Array1(),
1728 lambmin,NbP,*Padr,
1729 Cdim,Cdeg,
1730 ConstraintPoint, Cont, After,
1731 NbPoles, NbKnots,*FKRadr, *PRadr);
1732
1733 // recopie des poles du resultat sous forme de points 3D et de poids
1734 Standard_Integer NU, NV, indice ;
1735 if (InU) {
1736 NU = NbPoles;
1737 NV = BS->NbVPoles();
1738 }
1739 else {
1740 NU = BS->NbUPoles();
1741 NV = NbPoles;
1742 }
1743
1744 NewPoles = new (TColgp_HArray2OfPnt)(1,NU,1,NV);
1745 TColgp_Array2OfPnt& NewP = NewPoles->ChangeArray2();
1746 NewWeights = new (TColStd_HArray2OfReal) (1,NU,1,NV);
1747 TColStd_Array2OfReal& NewW = NewWeights->ChangeArray2();
1748
1749 if (!rational) NewW.Init(1.);
1750 NullWeight= Standard_False;
1751
1752 if (InU) {
1753 for (ii=1; ii<=NU && !NullWeight; ii++) {
1754 for (jj=1; jj<=NV && !NullWeight; jj++) {
1755 indice = 1+(ii-1)*Cdim+(jj-1)*gap;
1756 NewP(ii,jj).SetCoord(1,PRes(indice));
1757 NewP(ii,jj).SetCoord(2,PRes(indice+1));
1758 NewP(ii,jj).SetCoord(3,PRes(indice+2));
1759 if (rational) {
1760 ww = PRes(indice+3);
1761 if (Abs(ww - 1.0) < EpsW)
1762 ww = 1.0;
1763 if (ww < EpsW) {
1764 NullWeight = Standard_True;
1765 }
1766 else {
1767 NewW(ii,jj) = ww;
1768 NewP(ii,jj).ChangeCoord() /= ww;
1769 }
1770 }
1771 }
1772 }
1773 }
1774 else {
1775 for (jj=1; jj<=NV && !NullWeight; jj++) {
1776 for (ii=1; ii<=NU && !NullWeight; ii++) {
1777 indice = 1+(ii-1)*gap+(jj-1)*Cdim;
1778 NewP(ii,jj).SetCoord(1,PRes(indice));
1779 NewP(ii,jj).SetCoord(2,PRes(indice+1));
1780 NewP(ii,jj).SetCoord(3,PRes(indice+2));
1781 if (rational) {
1782 ww = PRes(indice+3);
1783 if (Abs(ww - 1.0) < EpsW)
1784 ww = 1.0;
1785 if (ww < EpsW) {
1786 NullWeight = Standard_True;
1787 }
1788 else {
1789 NewW(ii,jj) = ww;
1790 NewP(ii,jj).ChangeCoord() /= ww;
1791 }
1792 }
1793 }
1794 }
1795 }
1796
1797 if (NullWeight) {
1798 #ifdef OCCT_DEBUG
1799 std::cout << "Echec de l'Extension rationnelle" << std::endl;
1800 #endif
1801 lambmin /= 3.;
1802 NullWeight = Standard_False;
1803 }
1804 else {
1805 ExtOk = Standard_True;
1806 }
1807 }
1808
1809
1810 // recopie des noeuds plats sous forme de noeuds avec leurs multiplicites
1811 // calcul des degres du resultat
1812 Standard_Integer Usize = BS->NbUKnots(), Vsize = BS->NbVKnots(), UDeg, VDeg;
1813 if (InU)
1814 Usize++;
1815 else
1816 Vsize++;
1817 TColStd_Array1OfReal UKnots(1,Usize);
1818 TColStd_Array1OfReal VKnots(1,Vsize);
1819 TColStd_Array1OfInteger UMults(1,Usize);
1820 TColStd_Array1OfInteger VMults(1,Vsize);
1821 TColStd_Array1OfReal FKRes(1, NbKnots);
1822
1823 for (ii=1; ii<=NbKnots; ii++)
1824 FKRes(ii) = FK(ii);
1825
1826 if (InU) {
1827 BSplCLib::Knots(FKRes, UKnots, UMults);
1828 UDeg = Cdeg;
1829 UMults(Usize) = UDeg+1; // Petite verrue utile quand la continuite
1830 // n'est pas ok.
1831 BS->VKnots(VKnots);
1832 BS->VMultiplicities(VMults);
1833 VDeg = BS->VDegree();
1834 }
1835 else {
1836 BSplCLib::Knots(FKRes, VKnots, VMults);
1837 VDeg = Cdeg;
1838 VMults(Vsize) = VDeg+1;
1839 BS->UKnots(UKnots);
1840 BS->UMultiplicities(UMults);
1841 UDeg = BS->UDegree();
1842 }
1843
1844 // construction de la surface BSpline resultat
1845 Handle(Geom_BSplineSurface) Res =
1846 new (Geom_BSplineSurface) (NewPoles->Array2(),
1847 NewWeights->Array2(),
1848 UKnots,VKnots,
1849 UMults,VMults,
1850 UDeg,VDeg,
1851 BS->IsUPeriodic(),
1852 BS->IsVPeriodic());
1853 Surface = Res;
1854 }
1855
1856 //=======================================================================
1857 //function : Inertia
1858 //purpose :
1859 //=======================================================================
Inertia(const TColgp_Array1OfPnt & Points,gp_Pnt & Bary,gp_Dir & XDir,gp_Dir & YDir,Standard_Real & Xgap,Standard_Real & Ygap,Standard_Real & Zgap)1860 void GeomLib::Inertia(const TColgp_Array1OfPnt& Points,
1861 gp_Pnt& Bary,
1862 gp_Dir& XDir,
1863 gp_Dir& YDir,
1864 Standard_Real& Xgap,
1865 Standard_Real& Ygap,
1866 Standard_Real& Zgap)
1867 {
1868 gp_XYZ GB(0., 0., 0.), Diff;
1869 // gp_Vec A,B,C,D;
1870
1871 Standard_Integer i,nb=Points.Length();
1872 GB.SetCoord(0.,0.,0.);
1873 for (i=1; i<=nb; i++)
1874 GB += Points(i).XYZ();
1875
1876 GB /= nb;
1877
1878 math_Matrix M (1, 3, 1, 3);
1879 M.Init(0.);
1880 for (i=1; i<=nb; i++) {
1881 Diff.SetLinearForm(-1, Points(i).XYZ(), GB);
1882 M(1,1) += Diff.X() * Diff.X();
1883 M(2,2) += Diff.Y() * Diff.Y();
1884 M(3,3) += Diff.Z() * Diff.Z();
1885 M(1,2) += Diff.X() * Diff.Y();
1886 M(1,3) += Diff.X() * Diff.Z();
1887 M(2,3) += Diff.Y() * Diff.Z();
1888 }
1889
1890 M(2,1)=M(1,2) ;
1891 M(3,1)=M(1,3) ;
1892 M(3,2)=M(2,3) ;
1893
1894 M /= nb;
1895
1896 math_Jacobi J(M);
1897 if (!J.IsDone()) {
1898 #ifdef OCCT_DEBUG
1899 std::cout << "Erreur dans Jacobbi" << std::endl;
1900 M.Dump(std::cout);
1901 #endif
1902 }
1903
1904 Standard_Real n1,n2,n3;
1905
1906 n1=J.Value(1);
1907 n2=J.Value(2);
1908 n3=J.Value(3);
1909
1910 Standard_Real r1 = Min(Min(n1,n2),n3), r2;
1911 Standard_Integer m1, m2, m3;
1912 if (r1==n1) {
1913 m1 = 1;
1914 r2 = Min(n2,n3);
1915 if (r2==n2) {
1916 m2 = 2;
1917 m3 = 3;
1918 }
1919 else {
1920 m2 = 3;
1921 m3 = 2;
1922 }
1923 }
1924 else {
1925 if (r1==n2) {
1926 m1 = 2 ;
1927 r2 = Min(n1,n3);
1928 if (r2==n1) {
1929 m2 = 1;
1930 m3 = 3;
1931 }
1932 else {
1933 m2 = 3;
1934 m3 = 1;
1935 }
1936 }
1937 else {
1938 m1 = 3 ;
1939 r2 = Min(n1,n2);
1940 if (r2==n1) {
1941 m2 = 1;
1942 m3 = 2;
1943 }
1944 else {
1945 m2 = 2;
1946 m3 = 1;
1947 }
1948 }
1949 }
1950
1951 math_Vector V2(1,3),V3(1,3);
1952 J.Vector(m2,V2);
1953 J.Vector(m3,V3);
1954
1955 Bary.SetXYZ(GB);
1956 XDir.SetCoord(V3(1),V3(2),V3(3));
1957 YDir.SetCoord(V2(1),V2(2),V2(3));
1958
1959 Zgap = sqrt(Abs(J.Value(m1)));
1960 Ygap = sqrt(Abs(J.Value(m2)));
1961 Xgap = sqrt(Abs(J.Value(m3)));
1962 }
1963 //=======================================================================
1964 //function : AxeOfInertia
1965 //purpose :
1966 //=======================================================================
AxeOfInertia(const TColgp_Array1OfPnt & Points,gp_Ax2 & Axe,Standard_Boolean & IsSingular,const Standard_Real Tol)1967 void GeomLib::AxeOfInertia(const TColgp_Array1OfPnt& Points,
1968 gp_Ax2& Axe,
1969 Standard_Boolean& IsSingular,
1970 const Standard_Real Tol)
1971 {
1972 gp_Pnt Bary;
1973 gp_Dir OX,OY,OZ;
1974 Standard_Real gx, gy, gz;
1975
1976 GeomLib::Inertia(Points, Bary, OX, OY, gx, gy, gz);
1977
1978 if (gy*Points.Length()<=Tol) {
1979 gp_Ax2 axe (Bary, OX);
1980 OY = axe.XDirection();
1981 IsSingular = Standard_True;
1982 }
1983 else {
1984 IsSingular = Standard_False;
1985 }
1986
1987 OZ = OX^OY;
1988 gp_Ax2 TheAxe(Bary, OZ, OX);
1989 Axe = TheAxe;
1990 }
1991
1992 //=======================================================================
1993 //function : CanBeTreated
1994 //purpose : indicates if the surface can be treated(if the conditions are
1995 // filled) and need to be treated(if the surface hasn't been yet
1996 // treated or if the surface is rationnal and non periodic)
1997 //=======================================================================
1998
CanBeTreated(Handle (Geom_BSplineSurface)& BSurf)1999 static Standard_Boolean CanBeTreated(Handle(Geom_BSplineSurface)& BSurf)
2000
2001 {Standard_Integer i;
2002 Standard_Real lambda; //proportionnality coefficient
2003 Standard_Boolean AlreadyTreated=Standard_True;
2004
2005 if (!BSurf->IsURational()||(BSurf->IsUPeriodic()))
2006 return Standard_False;
2007 else {
2008 lambda=(BSurf->Weight(1,1)/BSurf->Weight(BSurf->NbUPoles(),1));
2009 for (i=1;i<=BSurf->NbVPoles();i++) //test of the proportionnality of the denominator on the boundaries
2010 if ((BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))<(1-Precision::Confusion()))||
2011 (BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))>(1+Precision::Confusion())))
2012 return Standard_False;
2013 i=1;
2014 while ((AlreadyTreated) && (i<=BSurf->NbVPoles())){ //tests if the surface has already been treated
2015 if (((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))<(1-Precision::Confusion()))||
2016 ((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))>(1+Precision::Confusion()))||
2017 ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))<(1-Precision::Confusion()))||
2018 ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))>(1+Precision::Confusion())))
2019 AlreadyTreated=Standard_False;
2020 i++;
2021 }
2022 if (AlreadyTreated)
2023 return Standard_False;
2024 }
2025 return Standard_True;
2026 }
2027
2028 //=======================================================================
2029 //class : law_evaluator
2030 //purpose : useful to estimate the value of a function of 2 variables
2031 //=======================================================================
2032
2033 class law_evaluator : public BSplSLib_EvaluatorFunction
2034 {
2035
2036 public:
2037
law_evaluator(const GeomLib_DenominatorMultiplierPtr theDenominatorPtr)2038 law_evaluator (const GeomLib_DenominatorMultiplierPtr theDenominatorPtr)
2039 : myDenominator (theDenominatorPtr) {}
2040
Evaluate(const Standard_Integer theDerivativeRequest,const Standard_Real theUParameter,const Standard_Real theVParameter,Standard_Real & theResult,Standard_Integer & theErrorCode) const2041 virtual void Evaluate (const Standard_Integer theDerivativeRequest,
2042 const Standard_Real theUParameter,
2043 const Standard_Real theVParameter,
2044 Standard_Real& theResult,
2045 Standard_Integer& theErrorCode) const
2046 {
2047 if ((myDenominator != NULL) && (theDerivativeRequest == 0))
2048 {
2049 theResult = myDenominator->Value (theUParameter, theVParameter);
2050 theErrorCode = 0;
2051 }
2052 else
2053 {
2054 theErrorCode = 1;
2055 }
2056 }
2057
2058 private:
2059
2060 GeomLib_DenominatorMultiplierPtr myDenominator;
2061
2062 };
2063
2064 //=======================================================================
2065 //function : CheckIfKnotExists
2066 //purpose : true if the knot already exists in the knot sequence
2067 //=======================================================================
2068
CheckIfKnotExists(const TColStd_Array1OfReal & surface_knots,const Standard_Real knot)2069 static Standard_Boolean CheckIfKnotExists(const TColStd_Array1OfReal& surface_knots,
2070 const Standard_Real knot)
2071
2072 {Standard_Integer i;
2073 for (i=1;i<=surface_knots.Length();i++)
2074 if ((surface_knots(i)-Precision::Confusion()<=knot)&&(surface_knots(i)+Precision::Confusion()>=knot))
2075 return Standard_True;
2076 return Standard_False;
2077 }
2078
2079 //=======================================================================
2080 //function : AddAKnot
2081 //purpose : add a knot and its multiplicity to the knot sequence. This knot
2082 // will be C2 and the degree is increased of deltasurface_degree
2083 //=======================================================================
2084
AddAKnot(const TColStd_Array1OfReal & knots,const TColStd_Array1OfInteger & mults,const Standard_Real knotinserted,const Standard_Integer deltasurface_degree,const Standard_Integer finalsurfacedegree,Handle (TColStd_HArray1OfReal)& newknots,Handle (TColStd_HArray1OfInteger)& newmults)2085 static void AddAKnot(const TColStd_Array1OfReal& knots,
2086 const TColStd_Array1OfInteger& mults,
2087 const Standard_Real knotinserted,
2088 const Standard_Integer deltasurface_degree,
2089 const Standard_Integer finalsurfacedegree,
2090 Handle(TColStd_HArray1OfReal) & newknots,
2091 Handle(TColStd_HArray1OfInteger) & newmults)
2092
2093 {Standard_Integer i;
2094
2095 newknots=new TColStd_HArray1OfReal(1,knots.Length()+1);
2096 newmults=new TColStd_HArray1OfInteger(1,knots.Length()+1);
2097 i=1;
2098 while (knots(i)<knotinserted){
2099 newknots->SetValue(i,knots(i));
2100 newmults->SetValue(i,mults(i)+deltasurface_degree);
2101 i++;
2102 }
2103 newknots->SetValue(i,knotinserted); //insertion of the new knot
2104 newmults->SetValue(i,finalsurfacedegree-2);
2105 i++;
2106 while (i<=newknots->Length()){
2107 newknots->SetValue(i,knots(i-1));
2108 newmults->SetValue(i,mults(i-1)+deltasurface_degree);
2109 i++;
2110 }
2111 }
2112
2113 //=======================================================================
2114 //function : Sort
2115 //purpose : give the new flat knots(u or v) of the surface
2116 //=======================================================================
2117
BuildFlatKnot(const TColStd_Array1OfReal & surface_knots,const TColStd_Array1OfInteger & surface_mults,const Standard_Integer deltasurface_degree,const Standard_Integer finalsurface_degree,const Standard_Real knotmin,const Standard_Real knotmax,Handle (TColStd_HArray1OfReal)& ResultKnots,Handle (TColStd_HArray1OfInteger)& ResultMults)2118 static void BuildFlatKnot(const TColStd_Array1OfReal& surface_knots,
2119 const TColStd_Array1OfInteger& surface_mults,
2120 const Standard_Integer deltasurface_degree,
2121 const Standard_Integer finalsurface_degree,
2122 const Standard_Real knotmin,
2123 const Standard_Real knotmax,
2124 Handle(TColStd_HArray1OfReal)& ResultKnots,
2125 Handle(TColStd_HArray1OfInteger)& ResultMults)
2126
2127 {
2128 Standard_Integer i;
2129
2130 if (CheckIfKnotExists(surface_knots,knotmin) &&
2131 CheckIfKnotExists(surface_knots,knotmax)){
2132 ResultKnots=new TColStd_HArray1OfReal(1,surface_knots.Length());
2133 ResultMults=new TColStd_HArray1OfInteger(1,surface_knots.Length());
2134 for (i=1;i<=surface_knots.Length();i++){
2135 ResultKnots->SetValue(i,surface_knots(i));
2136 ResultMults->SetValue(i,surface_mults(i)+deltasurface_degree);
2137 }
2138 }
2139 else{
2140 if ((CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax)))
2141 AddAKnot(surface_knots,surface_mults,knotmax,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2142 else{
2143 if ((!CheckIfKnotExists(surface_knots,knotmin))&&(CheckIfKnotExists(surface_knots,knotmax)))
2144 AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2145 else{
2146 if ((!CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax))&&
2147 (knotmin==knotmax)){
2148 AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2149 }
2150 else{
2151 Handle(TColStd_HArray1OfReal) IntermedKnots;
2152 Handle(TColStd_HArray1OfInteger) IntermedMults;
2153 AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,IntermedKnots,IntermedMults);
2154 AddAKnot(IntermedKnots->ChangeArray1(),IntermedMults->ChangeArray1(),knotmax,0,finalsurface_degree,ResultKnots,ResultMults);
2155 }
2156 }
2157 }
2158 }
2159 }
2160
2161 //=======================================================================
2162 //function : FunctionMultiply
2163 //purpose : multiply the surface BSurf by a(u,v) (law_evaluator) on its
2164 // numerator and denominator
2165 //=======================================================================
2166
FunctionMultiply(Handle (Geom_BSplineSurface)& BSurf,const Standard_Real knotmin,const Standard_Real knotmax)2167 static void FunctionMultiply(Handle(Geom_BSplineSurface)& BSurf,
2168 const Standard_Real knotmin,
2169 const Standard_Real knotmax)
2170
2171 {TColStd_Array1OfReal surface_u_knots(1,BSurf->NbUKnots()) ;
2172 TColStd_Array1OfInteger surface_u_mults(1,BSurf->NbUKnots()) ;
2173 TColStd_Array1OfReal surface_v_knots(1,BSurf->NbVKnots()) ;
2174 TColStd_Array1OfInteger surface_v_mults(1,BSurf->NbVKnots()) ;
2175 TColgp_Array2OfPnt surface_poles(1,BSurf->NbUPoles(),
2176 1,BSurf->NbVPoles()) ;
2177 TColStd_Array2OfReal surface_weights(1,BSurf->NbUPoles(),
2178 1,BSurf->NbVPoles()) ;
2179 Standard_Integer i,j,k,status,new_num_u_poles,new_num_v_poles,length=0;
2180 Handle(TColStd_HArray1OfReal) newuknots,newvknots;
2181 Handle(TColStd_HArray1OfInteger) newumults,newvmults;
2182
2183 BSurf->UKnots(surface_u_knots) ;
2184 BSurf->UMultiplicities(surface_u_mults) ;
2185 BSurf->VKnots(surface_v_knots) ;
2186 BSurf->VMultiplicities(surface_v_mults) ;
2187 BSurf->Poles(surface_poles) ;
2188 BSurf->Weights(surface_weights) ;
2189
2190 TColStd_Array1OfReal Knots(1,2);
2191 TColStd_Array1OfInteger Mults(1,2);
2192 Handle(TColStd_HArray1OfReal) NewKnots;
2193 Handle(TColStd_HArray1OfInteger) NewMults;
2194
2195 Knots(1)=0;
2196 Knots(2)=1;
2197 Mults(1)=4;
2198 Mults(2)=4;
2199 BuildFlatKnot(Knots,Mults,0,3,knotmin,knotmax,NewKnots,NewMults);
2200
2201 for (i=1;i<=NewMults->Length();i++)
2202 length+=NewMults->Value(i);
2203 TColStd_Array1OfReal FlatKnots(1,length);
2204 BSplCLib::KnotSequence(NewKnots->ChangeArray1(),NewMults->ChangeArray1(),FlatKnots);
2205
2206 GeomLib_DenominatorMultiplier aDenominator (BSurf, FlatKnots);
2207
2208 BuildFlatKnot(surface_u_knots,
2209 surface_u_mults,
2210 3,
2211 BSurf->UDegree()+3,
2212 knotmin,
2213 knotmax,
2214 newuknots,
2215 newumults);
2216 BuildFlatKnot(surface_v_knots,
2217 surface_v_mults,
2218 BSurf->VDegree(),
2219 2*(BSurf->VDegree()),
2220 1.0,
2221 0.0,
2222 newvknots,
2223 newvmults);
2224 length=0;
2225 for (i=1;i<=newumults->Length();i++)
2226 length+=newumults->Value(i);
2227 new_num_u_poles=(length-BSurf->UDegree()-3-1);
2228 TColStd_Array1OfReal newuflatknots(1,length);
2229 length=0;
2230 for (i=1;i<=newvmults->Length();i++)
2231 length+=newvmults->Value(i);
2232 new_num_v_poles=(length-2*BSurf->VDegree()-1);
2233 TColStd_Array1OfReal newvflatknots(1,length);
2234
2235 TColgp_Array2OfPnt NewNumerator(1,new_num_u_poles,1,new_num_v_poles);
2236 TColStd_Array2OfReal NewDenominator(1,new_num_u_poles,1,new_num_v_poles);
2237
2238 BSplCLib::KnotSequence(newuknots->ChangeArray1(),newumults->ChangeArray1(),newuflatknots);
2239 BSplCLib::KnotSequence(newvknots->ChangeArray1(),newvmults->ChangeArray1(),newvflatknots);
2240 //POP pour WNT
2241 law_evaluator ev (&aDenominator);
2242 // BSplSLib::FunctionMultiply(law_evaluator, //multiplication
2243 BSplSLib::FunctionMultiply(ev, //multiplication
2244 BSurf->UDegree(),
2245 BSurf->VDegree(),
2246 surface_u_knots,
2247 surface_v_knots,
2248 &surface_u_mults,
2249 &surface_v_mults,
2250 surface_poles,
2251 &surface_weights,
2252 newuflatknots,
2253 newvflatknots,
2254 BSurf->UDegree()+3,
2255 2*(BSurf->VDegree()),
2256 NewNumerator,
2257 NewDenominator,
2258 status);
2259 if (status!=0)
2260 throw Standard_ConstructionError("GeomLib Multiplication Error") ;
2261 for (i = 1 ; i <= new_num_u_poles ; i++) {
2262 for (j = 1 ; j <= new_num_v_poles ; j++) {
2263 for (k = 1 ; k <= 3 ; k++) {
2264 NewNumerator(i,j).SetCoord(k,NewNumerator(i,j).Coord(k)/NewDenominator(i,j)) ;
2265 }
2266 }
2267 }
2268 BSurf= new Geom_BSplineSurface(NewNumerator,
2269 NewDenominator,
2270 newuknots->ChangeArray1(),
2271 newvknots->ChangeArray1(),
2272 newumults->ChangeArray1(),
2273 newvmults->ChangeArray1(),
2274 BSurf->UDegree()+3,
2275 2*(BSurf->VDegree()) );
2276 }
2277
2278 //=======================================================================
2279 //function : CancelDenominatorDerivative1D
2280 //purpose : cancel the denominator derivative in one direction
2281 //=======================================================================
2282
CancelDenominatorDerivative1D(Handle (Geom_BSplineSurface)& BSurf)2283 static void CancelDenominatorDerivative1D(Handle(Geom_BSplineSurface) & BSurf)
2284
2285 {Standard_Integer i,j;
2286 Standard_Real uknotmin=1.0,uknotmax=0.0,
2287 x,y,
2288 startu_value,
2289 endu_value;
2290 TColStd_Array1OfReal BSurf_u_knots(1,BSurf->NbUKnots()) ;
2291
2292 startu_value=BSurf->UKnot(1);
2293 endu_value=BSurf->UKnot(BSurf->NbUKnots());
2294 BSurf->UKnots(BSurf_u_knots) ;
2295 BSplCLib::Reparametrize(0.0,1.0,BSurf_u_knots);
2296 BSurf->SetUKnots(BSurf_u_knots); //reparametrisation of the surface
2297 Handle(Geom_BSplineCurve) BCurve;
2298 TColStd_Array1OfReal BCurveWeights(1,BSurf->NbUPoles());
2299 TColgp_Array1OfPnt BCurvePoles(1,BSurf->NbUPoles());
2300 TColStd_Array1OfReal BCurveKnots(1,BSurf->NbUKnots());
2301 TColStd_Array1OfInteger BCurveMults(1,BSurf->NbUKnots());
2302
2303 if (CanBeTreated(BSurf)){
2304 for (i=1;i<=BSurf->NbVPoles();i++){ //loop on each pole function
2305 x=1.0;y=0.0;
2306 for (j=1;j<=BSurf->NbUPoles();j++){
2307 BCurveWeights(j)=BSurf->Weight(j,i);
2308 BCurvePoles(j)=BSurf->Pole(j,i);
2309 }
2310 BSurf->UKnots(BCurveKnots);
2311 BSurf->UMultiplicities(BCurveMults);
2312 BCurve = new Geom_BSplineCurve(BCurvePoles, //building of a pole function
2313 BCurveWeights,
2314 BCurveKnots,
2315 BCurveMults,
2316 BSurf->UDegree());
2317 Hermit::Solutionbis(BCurve,x,y,Precision::Confusion(),Precision::Confusion());
2318 if (x<uknotmin)
2319 uknotmin=x; //uknotmin,uknotmax:extremal knots
2320 if ((x!=1.0)&&(x>uknotmax))
2321 uknotmax=x;
2322 if ((y!=0.0)&&(y<uknotmin))
2323 uknotmin=y;
2324 if (y>uknotmax)
2325 uknotmax=y;
2326 }
2327
2328 FunctionMultiply(BSurf,uknotmin,uknotmax); //multiplication
2329
2330 BSurf->UKnots(BSurf_u_knots) ;
2331 BSplCLib::Reparametrize(startu_value,endu_value,BSurf_u_knots);
2332 BSurf->SetUKnots(BSurf_u_knots);
2333 }
2334 }
2335
2336 //=======================================================================
2337 //function : CancelDenominatorDerivative
2338 //purpose :
2339 //=======================================================================
2340
CancelDenominatorDerivative(Handle (Geom_BSplineSurface)& BSurf,const Standard_Boolean udirection,const Standard_Boolean vdirection)2341 void GeomLib::CancelDenominatorDerivative(Handle(Geom_BSplineSurface) & BSurf,
2342 const Standard_Boolean udirection,
2343 const Standard_Boolean vdirection)
2344
2345 {if (udirection && !vdirection)
2346 CancelDenominatorDerivative1D(BSurf);
2347 else{
2348 if (!udirection && vdirection) {
2349 BSurf->ExchangeUV();
2350 CancelDenominatorDerivative1D(BSurf);
2351 BSurf->ExchangeUV();
2352 }
2353 else{
2354 if (udirection && vdirection){ //optimize the treatment
2355 if (BSurf->UDegree()<=BSurf->VDegree()){
2356 CancelDenominatorDerivative1D(BSurf);
2357 BSurf->ExchangeUV();
2358 CancelDenominatorDerivative1D(BSurf);
2359 BSurf->ExchangeUV();
2360 }
2361 else{
2362 BSurf->ExchangeUV();
2363 CancelDenominatorDerivative1D(BSurf);
2364 BSurf->ExchangeUV();
2365 CancelDenominatorDerivative1D(BSurf);
2366 }
2367 }
2368 }
2369 }
2370 }
2371
2372 //=======================================================================
2373 //function : NormEstim
2374 //purpose :
2375 //=======================================================================
2376
NormEstim(const Handle (Geom_Surface)& S,const gp_Pnt2d & UV,const Standard_Real Tol,gp_Dir & N)2377 Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S,
2378 const gp_Pnt2d& UV,
2379 const Standard_Real Tol, gp_Dir& N)
2380 {
2381 gp_Vec DU, DV;
2382 gp_Pnt DummyPnt;
2383 Standard_Real aTol2 = Square(Tol);
2384
2385 S->D1(UV.X(), UV.Y(), DummyPnt, DU, DV);
2386
2387 Standard_Real MDU = DU.SquareMagnitude(), MDV = DV.SquareMagnitude();
2388
2389 if(MDU >= aTol2 && MDV >= aTol2) {
2390 gp_Vec Norm = DU^DV;
2391 Standard_Real Magn = Norm.SquareMagnitude();
2392 if(Magn < aTol2) return 3;
2393
2394 //Magn = sqrt(Magn);
2395 N.SetXYZ(Norm.XYZ());
2396
2397 return 0;
2398 }
2399 else {
2400 gp_Vec D2U, D2V, D2UV;
2401 Standard_Boolean isDone;
2402 CSLib_NormalStatus aStatus;
2403 gp_Dir aNormal;
2404
2405 S->D2(UV.X(), UV.Y(), DummyPnt, DU, DV, D2U, D2V, D2UV);
2406 CSLib::Normal(DU, DV, D2U, D2V, D2UV, Tol, isDone, aStatus, aNormal);
2407
2408 if (isDone) {
2409 Standard_Real Umin, Umax, Vmin, Vmax;
2410 Standard_Real step = 1.0e-5;
2411 Standard_Real eps = 1.0e-16;
2412 Standard_Real sign = -1.0;
2413
2414 S->Bounds(Umin, Umax, Vmin, Vmax);
2415
2416 // check for cone apex singularity point
2417 if ((UV.Y() > Vmin + step) && (UV.Y() < Vmax - step))
2418 {
2419 gp_Dir aNormal1, aNormal2;
2420 Standard_Real aConeSingularityAngleEps = 1.0e-4;
2421 S->D1(UV.X(), UV.Y() - sign * step, DummyPnt, DU, DV);
2422 if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
2423 aNormal1 = DU^DV;
2424 S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
2425 if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
2426 aNormal2 = DU^DV;
2427 if (aNormal1.IsOpposite(aNormal2, aConeSingularityAngleEps))
2428 return 2;
2429 }
2430 }
2431 }
2432
2433 // Along V
2434 if(MDU < aTol2 && MDV >= aTol2) {
2435 if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
2436 sign = 1.0;
2437 S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
2438 gp_Vec Norm = DU^DV;
2439 if (Norm.SquareMagnitude() < eps) {
2440 Standard_Real sign1 = -1.0;
2441 if ((Umax - UV.X()) > (UV.X() - Umin))
2442 sign1 = 1.0;
2443 S->D1(UV.X() + sign1 * step, UV.Y() + sign * step, DummyPnt, DU, DV);
2444 Norm = DU^DV;
2445 }
2446 if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
2447 aNormal.Reverse();
2448 }
2449
2450 // Along U
2451 if(MDV < aTol2 && MDU >= aTol2) {
2452 if ((Umax - UV.X()) > (UV.X() - Umin))
2453 sign = 1.0;
2454 S->D1(UV.X() + sign * step, UV.Y(), DummyPnt, DU, DV);
2455 gp_Vec Norm = DU^DV;
2456 if (Norm.SquareMagnitude() < eps) {
2457 Standard_Real sign1 = -1.0;
2458 if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
2459 sign1 = 1.0;
2460 S->D1(UV.X() + sign * step, UV.Y() + sign1 * step, DummyPnt, DU, DV);
2461 Norm = DU^DV;
2462 }
2463 if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
2464 aNormal.Reverse();
2465 }
2466
2467 // quasysingular
2468 if ((aStatus == CSLib_D1NuIsNull) || (aStatus == CSLib_D1NvIsNull) ||
2469 (aStatus == CSLib_D1NuIsParallelD1Nv)) {
2470 N.SetXYZ(aNormal.XYZ());
2471 return 1;
2472 }
2473 // conical
2474 if (aStatus == CSLib_InfinityOfSolutions)
2475 return 2;
2476 }
2477 // computation is impossible
2478 else {
2479 // conical
2480 if (aStatus == CSLib_D1NIsNull) {
2481 return 2;
2482 }
2483 return 3;
2484 }
2485 }
2486 return 3;
2487 }
2488
2489 //=======================================================================
2490 //function : IsClosed
2491 //purpose :
2492 //=======================================================================
IsClosed(const Handle (Geom_Surface)& S,const Standard_Real Tol,Standard_Boolean & isUClosed,Standard_Boolean & isVClosed)2493 void GeomLib::IsClosed (const Handle(Geom_Surface)& S,
2494 const Standard_Real Tol,
2495 Standard_Boolean& isUClosed, Standard_Boolean& isVClosed)
2496 {
2497 isUClosed = Standard_False;
2498 isVClosed = Standard_False;
2499 //
2500 GeomAdaptor_Surface aGAS(S);
2501 GeomAbs_SurfaceType aSType = aGAS.GetType();
2502 //
2503 Standard_Real u1, u2, v1, v2;
2504 u1 = aGAS.FirstUParameter();
2505 u2 = aGAS.LastUParameter();
2506 v1 = aGAS.FirstVParameter();
2507 v2 = aGAS.LastVParameter();
2508 //
2509 Standard_Real Tol2 = Tol * Tol;
2510 switch (aSType)
2511 {
2512 case GeomAbs_Plane:
2513 {
2514 return;
2515 }
2516 case GeomAbs_SurfaceOfExtrusion:
2517 {
2518 if (Precision::IsInfinite(u1) || Precision::IsInfinite(u2)) {
2519 // not closed
2520 return;
2521 }
2522 }
2523 Standard_FALLTHROUGH
2524 case GeomAbs_Cylinder:
2525 {
2526 if(Precision::IsInfinite(v1))
2527 v1 = 0.;
2528 gp_Pnt p1 = aGAS.Value(u1, v1);
2529 gp_Pnt p2 = aGAS.Value(u2, v1);
2530 isUClosed = p1.SquareDistance(p2) <= Tol2;
2531 return;
2532 }
2533 case GeomAbs_Cone:
2534 {
2535 //find v with maximal distance from axis
2536 if(!(Precision::IsInfinite(v1) || Precision::IsInfinite(v2)))
2537 {
2538 gp_Cone aCone = aGAS.Cone();
2539 gp_Pnt anApex = aCone.Apex();
2540 gp_Pnt P1 = aGAS.Value(u1, v1);
2541 gp_Pnt P2 = aGAS.Value(u1, v2);
2542 if(P2.SquareDistance(anApex) > P1.SquareDistance(anApex))
2543 {
2544 v1 = v2;
2545 }
2546 }
2547 else
2548 {
2549 v1 = 0.;
2550 }
2551 gp_Pnt p1 = aGAS.Value(u1, v1);
2552 gp_Pnt p2 = aGAS.Value(u2, v1);
2553 isUClosed = p1.SquareDistance(p2) <= Tol2;
2554 return;
2555 }
2556 case GeomAbs_Sphere:
2557 {
2558 //find v with maximal distance from axis
2559 if(v1*v2 <= 0.)
2560 {
2561 v1 = 0.;
2562 }
2563 else
2564 {
2565 if(v1 < 0.)
2566 {
2567 v1 = v2;
2568 }
2569 }
2570 gp_Pnt p1 = aGAS.Value(u1, v1);
2571 gp_Pnt p2 = aGAS.Value(u2, v1);
2572 isUClosed = p1.SquareDistance(p2) <= Tol2;
2573 return;
2574 }
2575 case GeomAbs_Torus:
2576 {
2577 Standard_Real ures = aGAS.UResolution(Tol);
2578 Standard_Real vres = aGAS.VResolution(Tol);
2579 //
2580 isUClosed = (u2 - u1) >= aGAS.UPeriod() - ures;
2581 isVClosed = (v2 - v1) >= aGAS.VPeriod() - vres;
2582 return;
2583 }
2584 case GeomAbs_BSplineSurface:
2585 {
2586 Handle(Geom_BSplineSurface) aBSpl = aGAS.BSpline();
2587 isUClosed = GeomLib::IsBSplUClosed(aBSpl, u1, u2, Tol);
2588 isVClosed = GeomLib::IsBSplVClosed(aBSpl, v1, v2, Tol);
2589 return;
2590 }
2591 case GeomAbs_BezierSurface:
2592 {
2593 Handle(Geom_BezierSurface) aBz = aGAS.Bezier();
2594 isUClosed = GeomLib::IsBzUClosed(aBz, u1, u2, Tol);
2595 isVClosed = GeomLib::IsBzVClosed(aBz, v1, v2, Tol);
2596 return;
2597 }
2598 case GeomAbs_SurfaceOfRevolution:
2599 case GeomAbs_OffsetSurface:
2600 case GeomAbs_OtherSurface:
2601 {
2602 Standard_Integer nbp = 23;
2603 if(Precision::IsInfinite(v1))
2604 {
2605 v1 = Sign(1., v1);
2606 }
2607 if(Precision::IsInfinite(v2))
2608 {
2609 v2 = Sign(1., v2);
2610 }
2611 //
2612 if(aSType == GeomAbs_OffsetSurface ||
2613 aSType == GeomAbs_OtherSurface)
2614 {
2615 if(Precision::IsInfinite(u1))
2616 {
2617 u1 = Sign(1., u1);
2618 }
2619 if(Precision::IsInfinite(u2))
2620 {
2621 u2 = Sign(1., u2);
2622 }
2623 }
2624 isUClosed = Standard_True;
2625 Standard_Real dt = (v2 - v1) / (nbp - 1);
2626 Standard_Real res = Max(aGAS.UResolution(Tol), Precision::PConfusion());
2627 if(dt <= res)
2628 {
2629 nbp = RealToInt((v2 - v1) /(2.*res)) + 1;
2630 nbp = Max(nbp, 2);
2631 dt = (v2 - v1) / (nbp - 1);
2632 }
2633 Standard_Real t;
2634 Standard_Integer i;
2635 for(i = 0; i < nbp; ++i)
2636 {
2637 t = (i == nbp-1 ? v2 : v1 + i * dt);
2638 gp_Pnt p1 = aGAS.Value(u1, t);
2639 gp_Pnt p2 = aGAS.Value(u2, t);
2640 if(p1.SquareDistance(p2) > Tol2)
2641 {
2642 isUClosed = Standard_False;
2643 break;
2644 }
2645 }
2646 //
2647 nbp = 23;
2648 isVClosed = Standard_True;
2649 dt = (u2 - u1) / (nbp - 1);
2650 res = Max(aGAS.VResolution(Tol), Precision::PConfusion());
2651 if(dt <= res)
2652 {
2653 nbp = RealToInt((u2 - u1) /(2.*res)) + 1;
2654 nbp = Max(nbp, 2);
2655 dt = (u2 - u1) / (nbp - 1);
2656 }
2657 for(i = 0; i < nbp; ++i)
2658 {
2659 t = (i == nbp-1 ? u2 : u1 + i * dt);
2660 gp_Pnt p1 = aGAS.Value(t, v1);
2661 gp_Pnt p2 = aGAS.Value(t, v2);
2662 if(p1.SquareDistance(p2) > Tol2)
2663 {
2664 isVClosed = Standard_False;
2665 break;
2666 }
2667 }
2668 return;
2669 }
2670 default:
2671 {
2672 return;
2673 }
2674 }
2675 }
2676
2677 //=======================================================================
2678 //function : IsBSplUClosed
2679 //purpose :
2680 //=======================================================================
IsBSplUClosed(const Handle (Geom_BSplineSurface)& S,const Standard_Real U1,const Standard_Real U2,const Standard_Real Tol)2681 Standard_Boolean GeomLib::IsBSplUClosed (const Handle(Geom_BSplineSurface)& S,
2682 const Standard_Real U1,
2683 const Standard_Real U2,
2684 const Standard_Real Tol)
2685 {
2686 Handle(Geom_Curve) aCUF = S->UIso( U1 );
2687 Handle(Geom_Curve) aCUL = S->UIso( U2 );
2688 if(aCUF.IsNull() || aCUL.IsNull())
2689 return Standard_False;
2690 Standard_Real Tol2 = 2.*Tol;
2691 Handle(Geom_BSplineCurve) aBsF = Handle(Geom_BSplineCurve)::DownCast(aCUF);
2692 Handle(Geom_BSplineCurve) aBsL = Handle(Geom_BSplineCurve)::DownCast(aCUL);
2693 const TColgp_Array1OfPnt& aPF = aBsF->Poles();
2694 const TColgp_Array1OfPnt& aPL = aBsL->Poles();
2695 const TColStd_Array1OfReal* WF = aBsF->Weights();
2696 const TColStd_Array1OfReal* WL = aBsL->Weights();
2697 return CompareWeightPoles(aPF, WF, aPL, WL, Tol2);
2698 }
2699
2700 //=======================================================================
2701 //function : IsBSplVClosed
2702 //purpose :
2703 //=======================================================================
IsBSplVClosed(const Handle (Geom_BSplineSurface)& S,const Standard_Real V1,const Standard_Real V2,const Standard_Real Tol)2704 Standard_Boolean GeomLib::IsBSplVClosed (const Handle(Geom_BSplineSurface)& S,
2705 const Standard_Real V1,
2706 const Standard_Real V2,
2707 const Standard_Real Tol)
2708 {
2709 Handle(Geom_Curve) aCVF = S->VIso( V1 );
2710 Handle(Geom_Curve) aCVL = S->VIso( V2 );
2711 if(aCVF.IsNull() || aCVL.IsNull())
2712 return Standard_False;
2713 Standard_Real Tol2 = 2.*Tol;
2714 Handle(Geom_BSplineCurve) aBsF = Handle(Geom_BSplineCurve)::DownCast(aCVF);
2715 Handle(Geom_BSplineCurve) aBsL = Handle(Geom_BSplineCurve)::DownCast(aCVL);
2716 const TColgp_Array1OfPnt& aPF = aBsF->Poles();
2717 const TColgp_Array1OfPnt& aPL = aBsL->Poles();
2718 const TColStd_Array1OfReal* WF = aBsF->Weights();
2719 const TColStd_Array1OfReal* WL = aBsL->Weights();
2720 return CompareWeightPoles(aPF, WF, aPL, WL, Tol2);
2721 }
2722 //=======================================================================
2723 //function : IsBzUClosed
2724 //purpose :
2725 //=======================================================================
IsBzUClosed(const Handle (Geom_BezierSurface)& S,const Standard_Real U1,const Standard_Real U2,const Standard_Real Tol)2726 Standard_Boolean GeomLib::IsBzUClosed (const Handle(Geom_BezierSurface)& S,
2727 const Standard_Real U1,
2728 const Standard_Real U2,
2729 const Standard_Real Tol)
2730 {
2731 Handle(Geom_Curve) aCUF = S->UIso( U1 );
2732 Handle(Geom_Curve) aCUL = S->UIso( U2 );
2733 if(aCUF.IsNull() || aCUL.IsNull())
2734 return Standard_False;
2735 Standard_Real Tol2 = 2.*Tol;
2736 Handle(Geom_BezierCurve) aBzF = Handle(Geom_BezierCurve)::DownCast(aCUF);
2737 Handle(Geom_BezierCurve) aBzL = Handle(Geom_BezierCurve)::DownCast(aCUL);
2738 const TColgp_Array1OfPnt& aPF = aBzF->Poles();
2739 const TColgp_Array1OfPnt& aPL = aBzL->Poles();
2740 //
2741 return CompareWeightPoles(aPF, 0, aPL, 0, Tol2);
2742 }
2743
2744 //=======================================================================
2745 //function : IsBzVClosed
2746 //purpose :
2747 //=======================================================================
IsBzVClosed(const Handle (Geom_BezierSurface)& S,const Standard_Real V1,const Standard_Real V2,const Standard_Real Tol)2748 Standard_Boolean GeomLib::IsBzVClosed (const Handle(Geom_BezierSurface)& S,
2749 const Standard_Real V1,
2750 const Standard_Real V2,
2751 const Standard_Real Tol)
2752 {
2753 Handle(Geom_Curve) aCVF = S->VIso( V1 );
2754 Handle(Geom_Curve) aCVL = S->VIso( V2 );
2755 if(aCVF.IsNull() || aCVL.IsNull())
2756 return Standard_False;
2757 Standard_Real Tol2 = 2.*Tol;
2758 Handle(Geom_BezierCurve) aBzF = Handle(Geom_BezierCurve)::DownCast(aCVF);
2759 Handle(Geom_BezierCurve) aBzL = Handle(Geom_BezierCurve)::DownCast(aCVL);
2760 const TColgp_Array1OfPnt& aPF = aBzF->Poles();
2761 const TColgp_Array1OfPnt& aPL = aBzL->Poles();
2762 //
2763 return CompareWeightPoles(aPF, 0, aPL, 0, Tol2);
2764 }
2765
2766 //=======================================================================
2767 //function : CompareWeightPoles
2768 //purpose : Checks if thePoles1(i)*theW1(i) is equal to thePoles2(i)*theW2(i)
2769 // with tolerance theTol.
2770 // It is necessary for non-rational B-splines and Bezier curves
2771 // to set theW1 and theW2 addresses to zero.
2772 //=======================================================================
CompareWeightPoles(const TColgp_Array1OfPnt & thePoles1,const TColStd_Array1OfReal * const theW1,const TColgp_Array1OfPnt & thePoles2,const TColStd_Array1OfReal * const theW2,const Standard_Real theTol)2773 static Standard_Boolean CompareWeightPoles(const TColgp_Array1OfPnt& thePoles1,
2774 const TColStd_Array1OfReal* const theW1,
2775 const TColgp_Array1OfPnt& thePoles2,
2776 const TColStd_Array1OfReal* const theW2,
2777 const Standard_Real theTol)
2778 {
2779 if(thePoles1.Length() != thePoles2.Length())
2780 {
2781 return Standard_False;
2782 }
2783 //
2784 Standard_Integer i = 1;
2785 for( i = 1 ; i <= thePoles1.Length(); i++ )
2786 {
2787 const Standard_Real aW1 = (theW1 == 0) ? 1.0 : theW1->Value(i);
2788 const Standard_Real aW2 = (theW2 == 0) ? 1.0 : theW2->Value(i);
2789
2790 gp_XYZ aPole1 = thePoles1.Value(i).XYZ() * aW1;
2791 gp_XYZ aPole2 = thePoles2.Value(i).XYZ() * aW2;
2792 if(!aPole1.IsEqual(aPole2, theTol))
2793 return Standard_False;
2794 }
2795 //
2796 return Standard_True;
2797 }
2798
2799 //=============================================================================
2800 //function : isIsoLine
2801 //purpose :
2802 //=============================================================================
isIsoLine(const Handle (Adaptor2d_Curve2d)theC2D,Standard_Boolean & theIsU,Standard_Real & theParam,Standard_Boolean & theIsForward)2803 Standard_Boolean GeomLib::isIsoLine (const Handle(Adaptor2d_Curve2d) theC2D,
2804 Standard_Boolean& theIsU,
2805 Standard_Real& theParam,
2806 Standard_Boolean& theIsForward)
2807 {
2808 // These variables are used to check line state (vertical or horizontal).
2809 Standard_Boolean isAppropriateType = Standard_False;
2810 gp_Pnt2d aLoc2d;
2811 gp_Dir2d aDir2d;
2812
2813 // Test type.
2814 const GeomAbs_CurveType aType = theC2D->GetType();
2815 if (aType == GeomAbs_Line)
2816 {
2817 gp_Lin2d aLin2d = theC2D->Line();
2818 aLoc2d = aLin2d.Location();
2819 aDir2d = aLin2d.Direction();
2820 isAppropriateType = Standard_True;
2821 }
2822 else if (aType == GeomAbs_BSplineCurve)
2823 {
2824 Handle(Geom2d_BSplineCurve) aBSpline2d = theC2D->BSpline();
2825 if (aBSpline2d->Degree() != 1 || aBSpline2d->NbPoles() != 2)
2826 return Standard_False; // Not a line or uneven parameterization.
2827
2828 aLoc2d = aBSpline2d->Pole(1);
2829
2830 // Vector should be non-degenerated.
2831 gp_Vec2d aVec2d(aBSpline2d->Pole(1), aBSpline2d->Pole(2));
2832 if (aVec2d.SquareMagnitude() < Precision::Confusion())
2833 return Standard_False; // Degenerated spline.
2834 aDir2d = aVec2d;
2835
2836 isAppropriateType = Standard_True;
2837 }
2838 else if (aType == GeomAbs_BezierCurve)
2839 {
2840 Handle(Geom2d_BezierCurve) aBezier2d = theC2D->Bezier();
2841 if (aBezier2d->Degree() != 1 || aBezier2d->NbPoles() != 2)
2842 return Standard_False; // Not a line or uneven parameterization.
2843
2844 aLoc2d = aBezier2d->Pole(1);
2845
2846 // Vector should be non-degenerated.
2847 gp_Vec2d aVec2d(aBezier2d->Pole(1), aBezier2d->Pole(2));
2848 if (aVec2d.SquareMagnitude() < Precision::Confusion())
2849 return Standard_False; // Degenerated spline.
2850 aDir2d = aVec2d;
2851
2852 isAppropriateType = Standard_True;
2853 }
2854
2855 if (!isAppropriateType)
2856 return Standard_False;
2857
2858 // Check line to be vertical or horizontal.
2859 if (aDir2d.IsParallel(gp::DX2d(), Precision::Angular()))
2860 {
2861 // Horizontal line. V = const.
2862 theIsU = Standard_False;
2863 theParam = aLoc2d.Y();
2864 theIsForward = aDir2d.Dot(gp::DX2d()) > 0.0;
2865 return Standard_True;
2866 }
2867 else if (aDir2d.IsParallel(gp::DY2d(), Precision::Angular()))
2868 {
2869 // Vertical line. U = const.
2870 theIsU = Standard_True;
2871 theParam = aLoc2d.X();
2872 theIsForward = aDir2d.Dot(gp::DY2d()) > 0.0;
2873 return Standard_True;
2874 }
2875
2876 return Standard_False;
2877 }
2878
2879 //=============================================================================
2880 //function : buildC3dOnIsoLine
2881 //purpose :
2882 //=============================================================================
Handle(Geom_Curve)2883 Handle(Geom_Curve) GeomLib::buildC3dOnIsoLine (const Handle(Adaptor2d_Curve2d) theC2D,
2884 const Handle(Adaptor3d_Surface) theSurf,
2885 const Standard_Real theFirst,
2886 const Standard_Real theLast,
2887 const Standard_Real theTolerance,
2888 const Standard_Boolean theIsU,
2889 const Standard_Real theParam,
2890 const Standard_Boolean theIsForward)
2891 {
2892 // Convert adapter to the appropriate type.
2893 Handle(GeomAdaptor_Surface) aGeomAdapter = Handle(GeomAdaptor_Surface)::DownCast(theSurf);
2894 if (aGeomAdapter.IsNull())
2895 return Handle(Geom_Curve)();
2896
2897 if (theSurf->GetType() == GeomAbs_Sphere)
2898 return Handle(Geom_Curve)();
2899
2900 // Extract isoline
2901 Handle(Geom_Surface) aSurf = aGeomAdapter->Surface();
2902 Handle(Geom_Curve) aC3d;
2903
2904 gp_Pnt2d aF2d = theC2D->Value(theC2D->FirstParameter());
2905 gp_Pnt2d aL2d = theC2D->Value(theC2D->LastParameter());
2906
2907 Standard_Boolean isToTrim = Standard_True;
2908 Standard_Real U1, U2, V1, V2;
2909 aSurf->Bounds(U1, U2, V1, V2);
2910
2911 if (theIsU)
2912 {
2913 Standard_Real aV1Param = Min(aF2d.Y(), aL2d.Y());
2914 Standard_Real aV2Param = Max(aF2d.Y(), aL2d.Y());
2915 if (aV2Param < V1 - theTolerance || aV1Param > V2 + theTolerance)
2916 {
2917 return Handle(Geom_Curve)();
2918 }
2919 else if (Precision::IsInfinite(V1) || Precision::IsInfinite(V2))
2920 {
2921 if (Abs(aV2Param - aV1Param) < Precision::PConfusion())
2922 {
2923 return Handle(Geom_Curve)();
2924 }
2925 aSurf = new Geom_RectangularTrimmedSurface(aSurf, U1, U2, aV1Param, aV2Param);
2926 isToTrim = Standard_False;
2927 }
2928 else
2929 {
2930 aV1Param = Max(aV1Param, V1);
2931 aV2Param = Min(aV2Param, V2);
2932 if (Abs(aV2Param - aV1Param) < Precision::PConfusion())
2933 {
2934 return Handle(Geom_Curve)();
2935 }
2936 }
2937 aC3d = aSurf->UIso(theParam);
2938 if (isToTrim)
2939 aC3d = new Geom_TrimmedCurve(aC3d, aV1Param, aV2Param);
2940 }
2941 else
2942 {
2943 Standard_Real aU1Param = Min(aF2d.X(), aL2d.X());
2944 Standard_Real aU2Param = Max(aF2d.X(), aL2d.X());
2945 if (aU2Param < U1 - theTolerance || aU1Param > U2 + theTolerance)
2946 {
2947 return Handle(Geom_Curve)();
2948 }
2949 else if (Precision::IsInfinite(U1) || Precision::IsInfinite(U2))
2950 {
2951 if (Abs(aU2Param - aU1Param) < Precision::PConfusion())
2952 {
2953 return Handle(Geom_Curve)();
2954 }
2955 aSurf = new Geom_RectangularTrimmedSurface(aSurf, aU1Param, aU2Param, V1, V2);
2956 isToTrim = Standard_False;
2957 }
2958 else
2959 {
2960 aU1Param = Max(aU1Param, U1);
2961 aU2Param = Min(aU2Param, U2);
2962 if (Abs(aU2Param - aU1Param) < Precision::PConfusion())
2963 {
2964 return Handle(Geom_Curve)();
2965 }
2966 }
2967 aC3d = aSurf->VIso(theParam);
2968 if (isToTrim)
2969 aC3d = new Geom_TrimmedCurve(aC3d, aU1Param, aU2Param);
2970 }
2971
2972 // Convert arbitrary curve type to the b-spline.
2973 Handle(Geom_BSplineCurve) aCurve3d = GeomConvert::CurveToBSplineCurve(aC3d, Convert_QuasiAngular);
2974 if (!theIsForward)
2975 aCurve3d->Reverse();
2976
2977 // Rebuild parameterization for the 3d curve to have the same parameterization with
2978 // a two-dimensional curve.
2979 TColStd_Array1OfReal aKnots = aCurve3d->Knots();
2980 BSplCLib::Reparametrize(theC2D->FirstParameter(), theC2D->LastParameter(), aKnots);
2981 aCurve3d->SetKnots(aKnots);
2982
2983 // Evaluate error.
2984 Standard_Real anError3d = 0.0;
2985
2986 const Standard_Real aParF = theFirst;
2987 const Standard_Real aParL = theLast;
2988 const Standard_Integer aNbPnt = 23;
2989 for (Standard_Integer anIdx = 0; anIdx <= aNbPnt; ++anIdx)
2990 {
2991 const Standard_Real aPar = aParF + ((aParL - aParF) * anIdx) / aNbPnt;
2992
2993 const gp_Pnt2d aPnt2d = theC2D->Value(aPar);
2994
2995 const gp_Pnt aPntC3D = aCurve3d->Value(aPar);
2996 const gp_Pnt aPntC2D = theSurf->Value(aPnt2d.X(), aPnt2d.Y());
2997
2998 const Standard_Real aSqDeviation = aPntC3D.SquareDistance(aPntC2D);
2999 anError3d = Max (aSqDeviation, anError3d);
3000 }
3001
3002 anError3d = Sqrt(anError3d);
3003
3004 // Target tolerance is not obtained. This situation happens for isolines on the sphere.
3005 // OCCT is unable to convert it keeping original parameterization, while the geometric
3006 // form of the result is entirely identical. In that case, it is better to utilize
3007 // a general-purpose approach.
3008 if (anError3d > theTolerance)
3009 return Handle(Geom_Curve)();
3010
3011 return aCurve3d;
3012 }
3013