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