projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023948: Wrong intersection between a surface of revolution and a plane.
[occt.git] / src / Law / Law_Interpolate.cxx
CommitLineData
b311480e 1// Created on: 1995-11-16
2// Created by: Laurent BOURESCHE
3// Copyright (c) 1995-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.
b311480e 16
7fd59977 17// Programme cree
18
19#include <Law_Interpolate.ixx>
20#include <Standard_ConstructionError.hxx>
21#include <PLib.hxx>
22#include <BSplCLib.hxx>
23#include <TColStd_Array1OfBoolean.hxx>
24#include <TColStd_Array1OfInteger.hxx>
25#include <gp_Pnt.hxx>
26
27//=======================================================================
28//function : CheckParameters
29//purpose :
30//=======================================================================
31
32static Standard_Boolean CheckParameters
33(const TColStd_Array1OfReal& Parameters)
34{
35 Standard_Integer ii;
36 Standard_Real distance;
37 Standard_Boolean result = Standard_True;
38 for (ii = Parameters.Lower() ; result && ii < Parameters.Upper() ; ii++) {
39 distance = Parameters.Value(ii+1) - Parameters.Value(ii);
40 result = (distance >= RealSmall());
41 }
42 return result;
43}
44
45//=======================================================================
46//function : BuildParameters
47//purpose :
48//=======================================================================
49static void BuildParameters(const Standard_Boolean PeriodicFlag,
50 const TColStd_Array1OfReal& PointsArray,
51 Handle(TColStd_HArray1OfReal)& ParametersPtr)
52{
53 Standard_Integer ii, index = 2;
54 Standard_Real distance;
55 Standard_Integer num_parameters = PointsArray.Length();
56 if (PeriodicFlag) {
57 num_parameters += 1;
58 }
59 ParametersPtr = new TColStd_HArray1OfReal(1, num_parameters);
60 ParametersPtr->SetValue(1,0.);
61 for (ii = PointsArray.Lower(); ii < PointsArray.Upper(); ii++) {
62 distance = Abs(PointsArray.Value(ii) - PointsArray.Value(ii+1));
63 ParametersPtr->SetValue(index, ParametersPtr->Value(ii) + distance);
64 index += 1 ;
65 }
66 if (PeriodicFlag) {
67 distance = Abs(PointsArray.Value(PointsArray.Upper()) -
68 PointsArray.Value(PointsArray.Lower()));
69 ParametersPtr->SetValue(index, ParametersPtr->Value(ii) + distance);
70 }
71}
72
73//=======================================================================
74//function : BuildPeriodicTangents
75//purpose :
76//=======================================================================
77
78static void BuildPeriodicTangent
79(const TColStd_Array1OfReal& PointsArray,
80 TColStd_Array1OfReal& TangentsArray,
81 TColStd_Array1OfBoolean& TangentFlags,
82 const TColStd_Array1OfReal& ParametersArray)
83{
84 Standard_Real point_array[3], parameter_array[3], eval_result[2];
85
86 if (PointsArray.Length() < 2) {
87 TangentFlags.SetValue(1,Standard_True);
88 TangentsArray.SetValue(1,0.);
89 }
90 else if (!TangentFlags.Value(1)){
91 //Pour les periodiques on evalue la tangente du point de fermeture
92 //par une interpolation de degre 2 entre le dernier point, le point
93 //de fermeture et le deuxieme point.
94 Standard_Integer degree = 2;
95 Standard_Real period = (ParametersArray.Value(ParametersArray.Upper()) -
96 ParametersArray.Value(ParametersArray.Lower()));
97 point_array [0] = PointsArray.Value(PointsArray.Upper());
98 point_array [1] = PointsArray.Value(PointsArray.Lower());
99 point_array [2] = PointsArray.Value(PointsArray.Lower()+1);
100 parameter_array[0] = ParametersArray.Value(ParametersArray.Upper() - 1) - period;
101 parameter_array[1] = ParametersArray.Value(ParametersArray.Lower());
102 parameter_array[2] = ParametersArray.Value(ParametersArray.Lower() + 1);
103 TangentFlags.SetValue(1,Standard_True);
104 PLib::EvalLagrange(parameter_array[1],
105 1,
106 degree,
107 1,
108 point_array[0],
109 parameter_array[0],
110 eval_result[0]);
111 TangentsArray.SetValue(1,eval_result[1]);
112 }
113}
114
115//=======================================================================
116//function : BuildTangents
117//purpose :
118//=======================================================================
119
120static void BuildTangents(const TColStd_Array1OfReal& PointsArray,
121 TColStd_Array1OfReal& TangentsArray,
122 TColStd_Array1OfBoolean& TangentFlags,
123 const TColStd_Array1OfReal& ParametersArray)
124{
125 Standard_Integer degree = 3;//,ii;
126 Standard_Real *point_array, *parameter_array, eval_result[2];
127
128 if ( PointsArray.Length() < 3) {
129 Standard_ConstructionError::Raise();
130 }
131 if (PointsArray.Length() == 3) {
132 degree = 2;
133 }
134 if (!TangentFlags.Value(1)) {
135 point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower());
136 parameter_array = (Standard_Real *) &ParametersArray.Value(1);
137 TangentFlags.SetValue(1,Standard_True);
138 PLib::EvalLagrange(ParametersArray.Value(1),
139 1,
140 degree,
141 1,
142 point_array[0],
143 parameter_array[0],
144 eval_result[0]);
145 TangentsArray.SetValue(1,eval_result[1]);
146 }
147 if (! TangentFlags.Value(TangentFlags.Upper())) {
148 point_array =
149 (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree);
150 TangentFlags.SetValue(TangentFlags.Upper(),Standard_True);
151 Standard_Integer iup = ParametersArray.Upper() - degree;
152 parameter_array = (Standard_Real *) &ParametersArray.Value(iup);
153 PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()),
154 1,
155 degree,
156 1,
157 point_array[0],
158 parameter_array[0],
159 eval_result[0]);
160 TangentsArray.SetValue(TangentsArray.Upper(),eval_result[1]);
161 }
162}
163
164//=======================================================================
165//function : Law_Interpolate
166//purpose :
167//=======================================================================
168
169Law_Interpolate::Law_Interpolate
170(const Handle(TColStd_HArray1OfReal)& PointsPtr,
171 const Standard_Boolean PeriodicFlag,
172 const Standard_Real Tolerance) :
173 myTolerance(Tolerance),
174 myPoints(PointsPtr),
175 myIsDone(Standard_False),
176 myPeriodic(PeriodicFlag),
177 myTangentRequest(Standard_False)
178
179{
180//Standard_Integer ii;
181 myTangents = new TColStd_HArray1OfReal (myPoints->Lower(),
182 myPoints->Upper());
183 myTangentFlags = new TColStd_HArray1OfBoolean(myPoints->Lower(),
184 myPoints->Upper());
185
186 BuildParameters(PeriodicFlag,
187 PointsPtr->Array1(),
188 myParameters);
189 myTangentFlags->Init(Standard_False);
190}
191
192//=======================================================================
193//function : Law_Interpolate
194//purpose :
195//=======================================================================
196
197Law_Interpolate::Law_Interpolate
198(const Handle(TColStd_HArray1OfReal)& PointsPtr,
199 const Handle(TColStd_HArray1OfReal)& ParametersPtr,
200 const Standard_Boolean PeriodicFlag,
201 const Standard_Real Tolerance) :
202 myTolerance(Tolerance),
203 myPoints(PointsPtr),
204 myIsDone(Standard_False),
205 myParameters(ParametersPtr),
206 myPeriodic(PeriodicFlag),
207 myTangentRequest(Standard_False)
208{
209//Standard_Integer ii;
210 if (PeriodicFlag) {
211 if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) {
212 Standard_ConstructionError::Raise();
213 }
214 }
215 myTangents = new TColStd_HArray1OfReal(myPoints->Lower(),
216 myPoints->Upper());
217 myTangentFlags = new TColStd_HArray1OfBoolean(myPoints->Lower(),
218 myPoints->Upper());
219 Standard_Boolean result = CheckParameters(ParametersPtr->Array1());
220 if (!result) {
221 Standard_ConstructionError::Raise();
222 }
223 myTangentFlags->Init(Standard_False);
224}
225
226//=======================================================================
227//function : Load
228//purpose :
229//=======================================================================
230
231void Law_Interpolate::Load
232(const TColStd_Array1OfReal& Tangents,
233 const Handle(TColStd_HArray1OfBoolean)& TangentFlagsPtr)
234{
235//Standard_Boolean result;
236 Standard_Integer ii;
237 myTangentRequest = Standard_True;
238 myTangentFlags = TangentFlagsPtr;
239 if (Tangents.Length() != myPoints->Length() ||
240 TangentFlagsPtr->Length() != myPoints->Length()) {
241 Standard_ConstructionError::Raise();
242 }
243 myTangents = new TColStd_HArray1OfReal(Tangents.Lower(),Tangents.Upper());
244 for (ii = Tangents.Lower() ; ii <= Tangents.Upper() ; ii++ ) {
245 myTangents->SetValue(ii,Tangents.Value(ii));
246 }
247}
248
249//=======================================================================
250//function : Load
251//purpose :
252//=======================================================================
253
254void Law_Interpolate::Load(const Standard_Real InitialTangent,
255 const Standard_Real FinalTangent)
256{
257//Standard_Boolean result;
258 myTangentRequest = Standard_True;
259 myTangentFlags->SetValue(1,Standard_True);
260 myTangents->SetValue(1,InitialTangent);
261 myTangentFlags->SetValue(myPoints->Length(),Standard_True);
262 myTangents->SetValue(myPoints->Length(),FinalTangent) ;
263}
264
265//=======================================================================
266//function : Perform
267//purpose :
268//=======================================================================
269
270void Law_Interpolate::Perform()
271{
272 if (myPeriodic) {
273 PerformPeriodic();
274 }
275 else {
276 PerformNonPeriodic();
277 }
278}
279
280//=======================================================================
281//function : PerformPeriodic
282//purpose :
283//=======================================================================
284
285void Law_Interpolate::PerformPeriodic()
286{
287 Standard_Integer degree,
288 ii,
289//jj,
290 index,
291 index1,
292//index2,
293 mult_index,
294 half_order,
295 inversion_problem,
296 num_points,
297 num_distinct_knots,
298 num_poles;
299
300 Standard_Real period;
301
302//gp_Pnt a_point;
303
304 num_points = myPoints->Length();
305 period = myParameters->Value(myParameters->Upper()) -
306 myParameters->Value(myParameters->Lower()) ;
307 num_poles = num_points + 1 ;
308 num_distinct_knots = num_points + 1;
309 half_order = 2;
310 degree = 3;
311 num_poles += 2;
312 if (myTangentRequest) {
313 for (ii = myTangentFlags->Lower() + 1;
314 ii <= myTangentFlags->Upper(); ii++) {
315 if (myTangentFlags->Value(ii)) {
316 num_poles += 1;
317 }
318 }
319 }
320 TColStd_Array1OfReal parameters(1,num_poles);
321 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1);
322 TColStd_Array1OfInteger mults(1,num_distinct_knots);
323 TColStd_Array1OfInteger contact_order_array(1, num_poles);
324 TColStd_Array1OfReal poles(1,num_poles);
325
326 for (ii = 1 ; ii <= half_order ; ii++) {
327 flatknots.SetValue(ii,myParameters->Value(myParameters->Upper() -1) -
328 period);
329 flatknots.SetValue(ii + half_order,myParameters->
330 Value(myParameters->Lower()));
331 flatknots.SetValue(num_poles + ii,
332 myParameters->Value(myParameters->Upper()));
333 flatknots.SetValue(num_poles + half_order + ii,
334 myParameters->Value(half_order) + period);
335 }
336 for (ii = 1 ; ii <= num_poles ; ii++) {
337 contact_order_array.SetValue(ii,0) ;
338 }
339 for (ii = 2; ii < num_distinct_knots; ii++) {
340 mults.SetValue(ii,1);
341 }
342 mults.SetValue(1,half_order);
343 mults.SetValue(num_distinct_knots ,half_order);
344
345 BuildPeriodicTangent(myPoints->Array1(),
346 myTangents->ChangeArray1(),
347 myTangentFlags->ChangeArray1(),
348 myParameters->Array1());
349
350 contact_order_array.SetValue(2,1);
351 parameters.SetValue(1,myParameters->Value(1));
352 parameters.SetValue(2,myParameters->Value(1));
353 poles.SetValue(1,myPoints->Value(1));
354 poles.SetValue(2,myTangents->Value(1));
355 mult_index = 2;
356 index = 3;
357 index1 = degree + 2;
358 if (myTangentRequest) {
359 for (ii = myTangentFlags->Lower() + 1;
360 ii <= myTangentFlags->Upper(); ii++) {
361 parameters.SetValue(index,myParameters->Value(ii));
362 flatknots.SetValue(index1,myParameters->Value(ii));
363 poles.SetValue(index,myPoints->Value(ii));
364 index += 1;
365 index1 += 1;
366 if (myTangentFlags->Value(ii)) {
367 mults.SetValue(mult_index,mults.Value(mult_index) + 1);
368 contact_order_array(index) = 1;
369
370 parameters.SetValue(index,
371 myParameters->Value(ii));
372 flatknots.SetValue(index1,myParameters->Value(ii));
373 poles.SetValue(index,myTangents->Value(ii));
374 index += 1;
375 index1 += 1;
376 }
377 mult_index += 1;
378 }
379 }
380 else {
381 index = degree + 1;
382 index1 = 2 ;
383 for(ii = myParameters->Lower(); ii <= myParameters->Upper(); ii++) {
384 parameters.SetValue(index1, myParameters->Value(ii));
385 flatknots.SetValue(index, myParameters->Value(ii));
386 index += 1;
387 index1 += 1;
388 }
389 index = 3;
390 for (ii = myPoints->Lower() + 1; ii <= myPoints->Upper(); ii++) {
391 //
392 // copy all the given points since the last one will be initialized
393 // below by the first point in the array myPoints
394 //
395 poles.SetValue(index, myPoints->Value(ii));
396 index += 1;
397 }
398 }
399 contact_order_array.SetValue(num_poles - 1, 1);
400 parameters.SetValue(num_poles-1,
401 myParameters->Value(myParameters->Upper()));
402 //
403 // for the periodic curve ONLY the tangent of the first point
404 // will be used since the curve should close itself at the first
405 // point See BuildPeriodicTangent
406 //
407 poles.SetValue(num_poles-1,myTangents->Value(1));
408 parameters.SetValue(num_poles,
409 myParameters->Value(myParameters->Upper()));
410 poles.SetValue(num_poles,
411 myPoints->Value(1));
412
413 BSplCLib::Interpolate(degree,
414 flatknots,
415 parameters,
416 contact_order_array,
417 1,
418 poles(1),
419 inversion_problem);
420 if (!inversion_problem) {
421 TColStd_Array1OfReal newpoles(poles.Value(1),
422 1,
423 num_poles - 2) ;
424 myCurve = new Law_BSpline(newpoles,
425 myParameters->Array1(),
426 mults,
427 degree,
428 myPeriodic);
429 myIsDone = Standard_True;
430 }
431}
432
433
434//=======================================================================
435//function : PerformNonPeriodic
436//purpose :
437//=======================================================================
438
439void Law_Interpolate::PerformNonPeriodic()
440{
441 Standard_Integer degree,
442 ii,
443//jj,
444 index,
445 index1,
446 index2,
447 index3,
448 mult_index,
449 inversion_problem,
450 num_points,
451 num_distinct_knots,
452 num_poles;
453
454 num_points =
455 num_distinct_knots =
456 num_poles = myPoints->Length();
457 if (num_poles == 2 && !myTangentRequest) {
458 degree = 1;
459 }
460 else if (num_poles == 3 && !myTangentRequest) {
461 degree = 2;
462 num_distinct_knots = 2;
463 }
464 else {
465 degree = 3;
466 num_poles += 2;
467 if (myTangentRequest) {
468 for (ii = myTangentFlags->Lower() + 1;
469 ii < myTangentFlags->Upper(); ii++) {
470 if (myTangentFlags->Value(ii)) {
471 num_poles += 1;
472 }
473 }
474 }
475 }
476 TColStd_Array1OfReal parameters(1,num_poles) ;
477 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1) ;
478 TColStd_Array1OfInteger mults(1,num_distinct_knots) ;
479 TColStd_Array1OfReal knots(1,num_distinct_knots) ;
480 TColStd_Array1OfInteger contact_order_array(1, num_poles) ;
481 TColStd_Array1OfReal poles(1,num_poles) ;
482
483 for (ii = 1 ; ii <= degree + 1 ; ii++) {
484 flatknots.SetValue(ii,myParameters->Value(1));
485 flatknots.SetValue(ii + num_poles,
486 myParameters->Value(num_points));
487 }
488 for (ii = 1 ; ii <= num_poles ; ii++) {
489 contact_order_array.SetValue(ii,0);
490 }
491 for (ii = 2 ; ii < num_distinct_knots ; ii++) {
492 mults.SetValue(ii,1);
493 }
494 mults.SetValue(1,degree + 1);
495 mults.SetValue(num_distinct_knots ,degree + 1);
496
497 switch (degree) {
498 case 1:
499 for (ii = 1 ; ii <= num_poles ; ii++) {
500 poles.SetValue(ii ,myPoints->Value(ii));
501 }
502 myCurve =
503 new Law_BSpline(poles,
504 myParameters->Array1(),
505 mults,
506 degree) ;
507 myIsDone = Standard_True ;
508 break ;
509 case 2:
510 knots.SetValue(1,myParameters->Value(1)) ;
511 knots.SetValue(2,myParameters->Value(num_poles)) ;
512 for (ii = 1 ; ii <= num_poles ; ii++) {
513 poles.SetValue(ii,myPoints->Value(ii)) ;
514
515 }
516 BSplCLib::Interpolate(degree,
517 flatknots,
518 myParameters->Array1(),
519 contact_order_array,
520 1,
521 poles(1),
522 inversion_problem) ;
523 if (!inversion_problem) {
524 myCurve = new Law_BSpline(poles,
525 knots,
526 mults,
527 degree);
528 myIsDone = Standard_True;
529 }
530 break;
531 case 3:
532 //
533 // check if the boundary conditions are set
534 //
535 if (num_points >= 3) {
536 //
537 // cannot build the tangents with degree 3 with only 2 points
538 // if those where not given in advance
539 //
540 BuildTangents(myPoints->Array1(),
541 myTangents->ChangeArray1(),
542 myTangentFlags->ChangeArray1(),
543 myParameters->Array1()) ;
544 }
545 contact_order_array.SetValue(2,1);
546 parameters.SetValue(1,myParameters->Value(1));
547 parameters.SetValue(2,myParameters->Value(1));
548 poles.SetValue(1,myPoints->Value(1));
549 poles.SetValue(2,myTangents->Value(1));
550 mult_index = 2 ;
551 index = 3 ;
552 index1 = 2 ;
553 index2 = myPoints->Lower() + 1 ;
554 index3 = degree + 2 ;
555 if (myTangentRequest) {
556 for (ii = myParameters->Lower() + 1 ;
557 ii < myParameters->Upper() ; ii++) {
558 parameters.SetValue(index,myParameters->Value(ii)) ;
559 poles.SetValue(index,myPoints->Value(index2)) ;
560 flatknots.SetValue(index3,myParameters->Value(ii)) ;
561 index += 1 ;
562 index3 += 1 ;
563 if (myTangentFlags->Value(index1)) {
564 //
565 // set the multiplicities, the order of the contact, the
566 // the flatknots,
567 //
568 mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
569 contact_order_array(index) = 1 ;
570 flatknots.SetValue(index3, myParameters->Value(ii)) ;
571 parameters.SetValue(index,
572 myParameters->Value(ii));
573 poles.SetValue(index,myTangents->Value(ii));
574 index += 1 ;
575 index3 += 1 ;
576 }
577 mult_index += 1 ;
578 index1 += 1 ;
579 index2 += 1 ;
580 }
581 }
582 else {
583 index1 = 2 ;
584 for(ii = myParameters->Lower(); ii <= myParameters->Upper(); ii++) {
585 parameters.SetValue(index1, myParameters->Value(ii));
586 index1 += 1;
587 }
588 index = 3 ;
589 for (ii = myPoints->Lower() + 1; ii <= myPoints->Upper() - 1; ii++){
590 poles.SetValue(index, myPoints->Value(ii));
591 index += 1 ;
592 }
593 index = degree + 1;
594 for(ii = myParameters->Lower(); ii <= myParameters->Upper(); ii++){
595 flatknots.SetValue(index, myParameters->Value(ii));
596 index += 1 ;
597 }
598 }
599 poles.SetValue(num_poles-1, myTangents->Value(num_points));
600
601 contact_order_array.SetValue(num_poles - 1, 1);
602 parameters.SetValue(num_poles,
603 myParameters->Value(myParameters->Upper())) ;
604 parameters.SetValue(num_poles -1,
605 myParameters->Value(myParameters->Upper())) ;
606
607 poles.SetValue(num_poles,
608 myPoints->Value(num_points)) ;
609
610 BSplCLib::Interpolate(degree,
611 flatknots,
612 parameters,
613 contact_order_array,
614 1,
615 poles(1),
616 inversion_problem) ;
617 if (!inversion_problem) {
618 myCurve = new Law_BSpline(poles,
619 myParameters->Array1(),
620 mults,
621 degree) ;
622 myIsDone = Standard_True;
623 }
624 break ;
625 }
626}
627
628//=======================================================================
857ffd5e 629//function : Handle(Geom_BSplineCurve)&
7fd59977 630//purpose :
631//=======================================================================
632
633const Handle(Law_BSpline)& Law_Interpolate::Curve() const
634{
635 if ( !myIsDone)
636 StdFail_NotDone::Raise(" ");
637 return myCurve;
638}
639
640//=======================================================================
641//function : IsDone
642//purpose :
643//=======================================================================
644
645Standard_Boolean Law_Interpolate::IsDone() const
646{
647 return myIsDone;
648}
Open CASCADE Technology repository