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