7fd59977 |
1 | // File: BSplCLib.cxx |
2 | // Created: Fri Aug 9 16:32:46 1991 |
3 | // Author: JCV |
4 | |
5 | // |
6 | // Modified RLE 9 Sep 1993 |
7 | // |
8 | // pmn : modified 28-01-97 : fixed a mistake in LocateParameter (PRO6973) |
9 | // pmn : modified 4-11-96 : fixed a mistake in BuildKnots (PRO6124) |
10 | // pmn : modified 28-Jun-96 : fixed a mistake in AntiBoorScheme |
11 | // xab : modified 15-Jun-95 : fixed a mistake in IsRational |
12 | // xab : modified 15-Mar-95 : removed Epsilon comparison in IsRational |
13 | // added RationalDerivatives. |
14 | // xab : 30-Mar-95 : fixed coupling with lti in RationalDerivatives |
15 | // xab : 15-Mar-96 : fixed a typo in Eval with extrapolation |
16 | // jct : 15-Apr-97 : added TangExtendToConstraint |
17 | // jct : 24-Apr-97 : correction on computation of Tbord and NewFlatKnots |
18 | // in TangExtendToConstraint; Continuity can be equal to 0 |
19 | // |
20 | |
21 | #define No_Standard_RangeError |
22 | #define No_Standard_OutOfRange |
23 | |
24 | #include <BSplCLib.ixx> |
25 | #include <PLib.hxx> |
26 | #include <Precision.hxx> |
27 | #include <Standard_NotImplemented.hxx> |
28 | |
29 | typedef gp_Pnt Pnt; |
30 | typedef gp_Vec Vec; |
31 | typedef TColgp_Array1OfPnt Array1OfPnt; |
32 | typedef TColStd_Array1OfReal Array1OfReal; |
33 | typedef TColStd_Array1OfInteger Array1OfInteger; |
34 | |
35 | //======================================================================= |
36 | //class : BSplCLib_LocalMatrix |
37 | //purpose: Auxiliary class optimizing creation of matrix buffer for |
38 | // evaluation of bspline (using stack allocation for main matrix) |
39 | //======================================================================= |
40 | |
41 | class BSplCLib_LocalMatrix : public math_Matrix |
42 | { |
43 | public: |
44 | BSplCLib_LocalMatrix (Standard_Integer DerivativeRequest, Standard_Integer Order) |
45 | : math_Matrix (myBuffer, 1, DerivativeRequest + 1, 1, Order) |
46 | { |
47 | if ( DerivativeRequest > BSplCLib::MaxDegree() || |
48 | Order > BSplCLib::MaxDegree()+1 || |
49 | BSplCLib::MaxDegree() > 25 ) |
50 | Standard_OutOfRange::Raise ("BSplCLib: bspline degree is greater than maximum supported"); |
51 | } |
52 | |
53 | private: |
54 | // local buffer, to be sufficient for addressing by index [Degree+1][Degree+1] |
55 | // (see math_Matrix implementation) |
56 | Standard_Real myBuffer[27*27]; |
57 | }; |
58 | |
59 | //======================================================================= |
60 | //class : BSplCLib_LocalArray |
61 | //purpose: Auxiliary class optimizing creation of array buffer for |
62 | // evaluation of bspline (using stack allocation for small arrays) |
63 | //======================================================================= |
64 | |
65 | #define LOCARRAY_BUFFER 1024 |
66 | class BSplCLib_LocalArray |
67 | { |
68 | public: |
69 | BSplCLib_LocalArray (Standard_Integer Size) |
70 | : myPtr(myBuffer) |
71 | { |
72 | if ( Size > LOCARRAY_BUFFER ) |
73 | myPtr = (Standard_Real*)Standard::Allocate (Size * sizeof(Standard_Real)); |
74 | } |
75 | |
76 | ~BSplCLib_LocalArray () |
77 | { |
78 | if ( myPtr != myBuffer ) |
79 | Standard::Free (*(Standard_Address*)&myPtr); |
80 | } |
81 | |
82 | Standard_Real& operator [] (int i) { return myPtr[i]; } |
83 | |
84 | private: |
85 | Standard_Real myBuffer[LOCARRAY_BUFFER]; |
86 | Standard_Real* myPtr; |
87 | }; |
88 | |
89 | //======================================================================= |
90 | //function : Hunt |
91 | //purpose : |
92 | //======================================================================= |
93 | |
94 | void BSplCLib::Hunt (const Array1OfReal& XX, |
95 | const Standard_Real X, |
96 | Standard_Integer& Ilc) |
97 | { |
98 | // replaced by simple dichotomy (RLE) |
99 | Ilc = XX.Lower(); |
100 | const Standard_Real *px = &XX(Ilc); |
101 | px -= Ilc; |
102 | |
103 | if (X < px[Ilc]) { |
104 | Ilc--; |
105 | return; |
106 | } |
107 | Standard_Integer Ihi = XX.Upper(); |
108 | if (X > px[Ihi]) { |
109 | Ilc = Ihi + 1; |
110 | return; |
111 | } |
112 | Standard_Integer Im; |
113 | |
114 | while (Ihi - Ilc != 1) { |
115 | Im = (Ihi + Ilc) >> 1; |
116 | if (X > px[Im]) Ilc = Im; |
117 | else Ihi = Im; |
118 | } |
119 | } |
120 | |
121 | //======================================================================= |
122 | //function : FirstUKnotIndex |
123 | //purpose : |
124 | //======================================================================= |
125 | |
126 | Standard_Integer BSplCLib::FirstUKnotIndex (const Standard_Integer Degree, |
127 | const TColStd_Array1OfInteger& Mults) |
128 | { |
129 | Standard_Integer Index = Mults.Lower(); |
130 | Standard_Integer SigmaMult = Mults(Index); |
131 | |
132 | while (SigmaMult <= Degree) { |
133 | Index++; |
134 | SigmaMult += Mults (Index); |
135 | } |
136 | return Index; |
137 | } |
138 | |
139 | //======================================================================= |
140 | //function : LastUKnotIndex |
141 | //purpose : |
142 | //======================================================================= |
143 | |
144 | Standard_Integer BSplCLib::LastUKnotIndex (const Standard_Integer Degree, |
145 | const Array1OfInteger& Mults) |
146 | { |
147 | Standard_Integer Index = Mults.Upper(); |
148 | Standard_Integer SigmaMult = Mults(Index); |
149 | |
150 | while (SigmaMult <= Degree) { |
151 | Index--; |
152 | SigmaMult += Mults.Value (Index); |
153 | } |
154 | return Index; |
155 | } |
156 | |
157 | //======================================================================= |
158 | //function : FlatIndex |
159 | //purpose : |
160 | //======================================================================= |
161 | |
162 | Standard_Integer BSplCLib::FlatIndex |
163 | (const Standard_Integer Degree, |
164 | const Standard_Integer Index, |
165 | const TColStd_Array1OfInteger& Mults, |
166 | const Standard_Boolean Periodic) |
167 | { |
168 | Standard_Integer i, index = Index; |
169 | const Standard_Integer MLower = Mults.Lower(); |
170 | const Standard_Integer *pmu = &Mults(MLower); |
171 | pmu -= MLower; |
172 | |
173 | for (i = MLower + 1; i <= Index; i++) |
174 | index += pmu[i] - 1; |
175 | if ( Periodic) |
176 | index += Degree; |
177 | else |
178 | index += pmu[MLower] - 1; |
179 | return index; |
180 | } |
181 | |
182 | //======================================================================= |
183 | //function : LocateParameter |
184 | //purpose : Traitement des noeuds avec multiplicites |
185 | //pmn 28-01-97 -> calcule eventuel de la periode. |
186 | //======================================================================= |
187 | |
188 | void BSplCLib::LocateParameter |
189 | (const Standard_Integer , //Degree, |
190 | const Array1OfReal& Knots, |
191 | const Array1OfInteger& , //Mults, |
192 | const Standard_Real U, |
193 | const Standard_Boolean IsPeriodic, |
194 | const Standard_Integer FromK1, |
195 | const Standard_Integer ToK2, |
196 | Standard_Integer& KnotIndex, |
197 | Standard_Real& NewU) |
198 | { |
199 | Standard_Real uf = 0, ul=1; |
200 | if (IsPeriodic) { |
201 | uf = Knots(Knots.Lower()); |
202 | ul = Knots(Knots.Upper()); |
203 | } |
204 | BSplCLib::LocateParameter(Knots,U,IsPeriodic,FromK1,ToK2, |
205 | KnotIndex,NewU, uf, ul); |
206 | } |
207 | |
208 | //======================================================================= |
209 | //function : LocateParameter |
210 | //purpose : Pour des noeuds plats |
211 | // pmn 28-01-97 -> On a bel est bien besoin du degree pour calculer |
212 | // la periode eventuelle |
213 | //======================================================================= |
214 | |
215 | void BSplCLib::LocateParameter |
216 | (const Standard_Integer Degree, |
217 | const Array1OfReal& Knots, |
218 | const Standard_Real U, |
219 | const Standard_Boolean IsPeriodic, |
220 | const Standard_Integer FromK1, |
221 | const Standard_Integer ToK2, |
222 | Standard_Integer& KnotIndex, |
223 | Standard_Real& NewU) |
224 | { |
225 | if (IsPeriodic) |
226 | BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2, |
227 | KnotIndex, NewU, |
228 | Knots(Knots.Lower() + Degree), |
229 | Knots(Knots.Upper() - Degree)); |
230 | else |
231 | BSplCLib::LocateParameter(Knots, U, IsPeriodic, FromK1, ToK2, |
232 | KnotIndex, NewU, |
233 | 0., |
234 | 1.); |
235 | } |
236 | |
237 | //======================================================================= |
238 | //function : LocateParameter |
239 | //purpose : Claculs effectifs |
240 | // pmn 28-01-97 : Ajoute les bornes de la periode en argument d'entree, car il est |
241 | // car il est imposible de les inventer a ce niveaux. |
242 | //======================================================================= |
243 | |
244 | void BSplCLib::LocateParameter |
245 | (const TColStd_Array1OfReal& Knots, |
246 | const Standard_Real U, |
247 | const Standard_Boolean IsPeriodic, |
248 | const Standard_Integer FromK1, |
249 | const Standard_Integer ToK2, |
250 | Standard_Integer& KnotIndex, |
251 | Standard_Real& NewU, |
252 | const Standard_Real UFirst, |
253 | const Standard_Real ULast) |
254 | { |
255 | Standard_Integer First,Last; |
256 | if (FromK1 < ToK2) { |
257 | First = FromK1; |
258 | Last = ToK2; |
259 | } |
260 | else { |
261 | First = ToK2; |
262 | Last = FromK1; |
263 | } |
264 | Standard_Integer Last1 = Last - 1; |
265 | NewU = U; |
266 | if (IsPeriodic) { |
267 | Standard_Real Period = ULast - UFirst; |
268 | |
269 | while (NewU > ULast ) |
270 | NewU -= Period; |
271 | |
272 | while (NewU < UFirst) |
273 | NewU += Period; |
274 | } |
275 | |
276 | BSplCLib::Hunt (Knots, NewU, KnotIndex); |
277 | |
278 | Standard_Real Eps = Epsilon(U); |
279 | Standard_Real val; |
280 | if (Eps < 0) Eps = - Eps; |
281 | Standard_Integer KLower = Knots.Lower(); |
282 | const Standard_Real *knots = &Knots(KLower); |
283 | knots -= KLower; |
284 | if ( KnotIndex < Knots.Upper()) { |
285 | val = NewU - knots[KnotIndex + 1]; |
286 | if (val < 0) val = - val; |
287 | // <= pour etre coherant avec les Segment ou Eps correspond a un bit d'erreur. |
288 | if (val <= Eps) KnotIndex++; |
289 | } |
290 | if (KnotIndex < First) KnotIndex = First; |
291 | if (KnotIndex > Last1) KnotIndex = Last1; |
292 | |
293 | if (KnotIndex != Last1) { |
294 | Standard_Real K1 = knots[KnotIndex]; |
295 | Standard_Real K2 = knots[KnotIndex + 1]; |
296 | val = K2 - K1; |
297 | if (val < 0) val = - val; |
298 | |
299 | while (val <= Eps) { |
300 | KnotIndex++; |
301 | K1 = K2; |
302 | K2 = knots[KnotIndex + 1]; |
303 | val = K2 - K1; |
304 | if (val < 0) val = - val; |
305 | } |
306 | } |
307 | } |
308 | |
309 | //======================================================================= |
310 | //function : LocateParameter |
311 | //purpose : the index is recomputed only if out of range |
312 | //pmn 28-01-97 -> calcule eventuel de la periode. |
313 | //======================================================================= |
314 | |
315 | void BSplCLib::LocateParameter |
316 | (const Standard_Integer Degree, |
317 | const TColStd_Array1OfReal& Knots, |
318 | const TColStd_Array1OfInteger& Mults, |
319 | const Standard_Real U, |
320 | const Standard_Boolean Periodic, |
321 | Standard_Integer& KnotIndex, |
322 | Standard_Real& NewU) |
323 | { |
324 | Standard_Integer first,last; |
325 | if (&Mults) { |
326 | if (Periodic) { |
327 | first = Knots.Lower(); |
328 | last = Knots.Upper(); |
329 | } |
330 | else { |
331 | first = FirstUKnotIndex(Degree,Mults); |
332 | last = LastUKnotIndex (Degree,Mults); |
333 | } |
334 | } |
335 | else { |
336 | first = Knots.Lower() + Degree; |
337 | last = Knots.Upper() - Degree; |
338 | } |
339 | if ( KnotIndex < first || KnotIndex > last) |
340 | BSplCLib::LocateParameter(Knots, U, Periodic, first, last, |
341 | KnotIndex, NewU, Knots(first), Knots(last)); |
342 | else |
343 | NewU = U; |
344 | } |
345 | |
346 | //======================================================================= |
347 | //function : MaxKnotMult |
348 | //purpose : |
349 | //======================================================================= |
350 | |
351 | Standard_Integer BSplCLib::MaxKnotMult |
352 | (const Array1OfInteger& Mults, |
353 | const Standard_Integer FromK1, |
354 | const Standard_Integer ToK2) |
355 | { |
356 | Standard_Integer MLower = Mults.Lower(); |
357 | const Standard_Integer *pmu = &Mults(MLower); |
358 | pmu -= MLower; |
359 | Standard_Integer MaxMult = pmu[FromK1]; |
360 | |
361 | for (Standard_Integer i = FromK1; i <= ToK2; i++) { |
362 | if (MaxMult < pmu[i]) MaxMult = pmu[i]; |
363 | } |
364 | return MaxMult; |
365 | } |
366 | |
367 | //======================================================================= |
368 | //function : MinKnotMult |
369 | //purpose : |
370 | //======================================================================= |
371 | |
372 | Standard_Integer BSplCLib::MinKnotMult |
373 | (const Array1OfInteger& Mults, |
374 | const Standard_Integer FromK1, |
375 | const Standard_Integer ToK2) |
376 | { |
377 | Standard_Integer MLower = Mults.Lower(); |
378 | const Standard_Integer *pmu = &Mults(MLower); |
379 | pmu -= MLower; |
380 | Standard_Integer MinMult = pmu[FromK1]; |
381 | |
382 | for (Standard_Integer i = FromK1; i <= ToK2; i++) { |
383 | if (MinMult > pmu[i]) MinMult = pmu[i]; |
384 | } |
385 | return MinMult; |
386 | } |
387 | |
388 | //======================================================================= |
389 | //function : NbPoles |
390 | //purpose : |
391 | //======================================================================= |
392 | |
393 | Standard_Integer BSplCLib::NbPoles(const Standard_Integer Degree, |
394 | const Standard_Boolean Periodic, |
395 | const TColStd_Array1OfInteger& Mults) |
396 | { |
397 | Standard_Integer i,sigma = 0; |
398 | Standard_Integer f = Mults.Lower(); |
399 | Standard_Integer l = Mults.Upper(); |
400 | const Standard_Integer * pmu = &Mults(f); |
401 | pmu -= f; |
402 | Standard_Integer Mf = pmu[f]; |
403 | Standard_Integer Ml = pmu[l]; |
404 | if (Mf <= 0) return 0; |
405 | if (Ml <= 0) return 0; |
406 | if (Periodic) { |
407 | if (Mf > Degree) return 0; |
408 | if (Ml > Degree) return 0; |
409 | if (Mf != Ml ) return 0; |
410 | sigma = Mf; |
411 | } |
412 | else { |
413 | Standard_Integer Deg1 = Degree + 1; |
414 | if (Mf > Deg1) return 0; |
415 | if (Ml > Deg1) return 0; |
416 | sigma = Mf + Ml - Deg1; |
417 | } |
418 | |
419 | for (i = f + 1; i < l; i++) { |
420 | if (pmu[i] <= 0 ) return 0; |
421 | if (pmu[i] > Degree) return 0; |
422 | sigma += pmu[i]; |
423 | } |
424 | return sigma; |
425 | } |
426 | |
427 | //======================================================================= |
428 | //function : KnotSequenceLength |
429 | //purpose : |
430 | //======================================================================= |
431 | |
432 | Standard_Integer BSplCLib::KnotSequenceLength |
433 | (const TColStd_Array1OfInteger& Mults, |
434 | const Standard_Integer Degree, |
435 | const Standard_Boolean Periodic) |
436 | { |
437 | Standard_Integer i,l = 0; |
438 | Standard_Integer MLower = Mults.Lower(); |
439 | Standard_Integer MUpper = Mults.Upper(); |
440 | const Standard_Integer * pmu = &Mults(MLower); |
441 | pmu -= MLower; |
442 | |
443 | for (i = MLower; i <= MUpper; i++) |
444 | l += pmu[i]; |
445 | if (Periodic) l += 2 * (Degree + 1 - pmu[MLower]); |
446 | return l; |
447 | } |
448 | |
449 | //======================================================================= |
450 | //function : KnotSequence |
451 | //purpose : |
452 | //======================================================================= |
453 | |
454 | void BSplCLib::KnotSequence |
455 | (const TColStd_Array1OfReal& Knots, |
456 | const TColStd_Array1OfInteger& Mults, |
457 | TColStd_Array1OfReal& KnotSeq) |
458 | { |
459 | BSplCLib::KnotSequence(Knots,Mults,0,Standard_False,KnotSeq); |
460 | } |
461 | |
462 | //======================================================================= |
463 | //function : KnotSequence |
464 | //purpose : |
465 | //======================================================================= |
466 | |
467 | void BSplCLib::KnotSequence |
468 | (const TColStd_Array1OfReal& Knots, |
469 | const TColStd_Array1OfInteger& Mults, |
470 | const Standard_Integer Degree, |
471 | const Standard_Boolean Periodic, |
472 | TColStd_Array1OfReal& KnotSeq) |
473 | { |
474 | Standard_Real K; |
475 | Standard_Integer Mult; |
476 | Standard_Integer MLower = Mults.Lower(); |
477 | const Standard_Integer * pmu = &Mults(MLower); |
478 | pmu -= MLower; |
479 | Standard_Integer KLower = Knots.Lower(); |
480 | Standard_Integer KUpper = Knots.Upper(); |
481 | const Standard_Real * pkn = &Knots(KLower); |
482 | pkn -= KLower; |
483 | Standard_Integer M1 = Degree + 1 - pmu[MLower]; // for periodic |
484 | Standard_Integer i,j,index = Periodic ? M1 + 1 : 1; |
485 | |
486 | for (i = KLower; i <= KUpper; i++) { |
487 | Mult = pmu[i]; |
488 | K = pkn[i]; |
489 | |
490 | for (j = 1; j <= Mult; j++) { |
491 | KnotSeq (index) = K; |
492 | index++; |
493 | } |
494 | } |
495 | if (Periodic) { |
496 | Standard_Real period = pkn[KUpper] - pkn[KLower]; |
497 | Standard_Integer m; |
498 | m = 1; |
499 | j = KUpper - 1; |
500 | |
501 | for (i = M1; i >= 1; i--) { |
502 | KnotSeq(i) = pkn[j] - period; |
503 | m++; |
504 | if (m > pmu[j]) { |
505 | j--; |
506 | m = 1; |
507 | } |
508 | } |
509 | m = 1; |
510 | j = KLower + 1; |
511 | |
512 | for (i = index; i <= KnotSeq.Upper(); i++) { |
513 | KnotSeq(i) = pkn[j] + period; |
514 | m++; |
515 | if (m > pmu[j]) { |
516 | j++; |
517 | m = 1; |
518 | } |
519 | } |
520 | } |
521 | } |
522 | |
523 | //======================================================================= |
524 | //function : KnotsLength |
525 | //purpose : |
526 | //======================================================================= |
527 | Standard_Integer BSplCLib::KnotsLength(const TColStd_Array1OfReal& SeqKnots, |
528 | // const Standard_Boolean Periodic) |
529 | const Standard_Boolean ) |
530 | { |
531 | Standard_Integer sizeMult = 1; |
532 | Standard_Real val = SeqKnots(1); |
533 | for (Standard_Integer jj=2; |
534 | jj<=SeqKnots.Length();jj++) |
535 | { |
536 | // test on strict equality on nodes |
537 | if (SeqKnots(jj)!=val) |
538 | { |
539 | val = SeqKnots(jj); |
540 | sizeMult++; |
541 | } |
542 | } |
543 | return sizeMult; |
544 | } |
545 | |
546 | //======================================================================= |
547 | //function : Knots |
548 | //purpose : |
549 | //======================================================================= |
550 | void BSplCLib::Knots(const TColStd_Array1OfReal& SeqKnots, |
551 | TColStd_Array1OfReal &knots, |
552 | TColStd_Array1OfInteger &mult, |
553 | // const Standard_Boolean Periodic) |
554 | const Standard_Boolean ) |
555 | { |
556 | Standard_Real val = SeqKnots(1); |
557 | Standard_Integer kk=1; |
558 | knots(kk) = val; |
559 | mult(kk) = 1; |
560 | |
561 | for (Standard_Integer jj=2;jj<=SeqKnots.Length();jj++) |
562 | { |
563 | // test on strict equality on nodes |
564 | if (SeqKnots(jj)!=val) |
565 | { |
566 | val = SeqKnots(jj); |
567 | kk++; |
568 | knots(kk) = val; |
569 | mult(kk) = 1; |
570 | } |
571 | else |
572 | { |
573 | mult(kk)++; |
574 | } |
575 | } |
576 | } |
577 | |
578 | //======================================================================= |
579 | //function : KnotForm |
580 | //purpose : |
581 | //======================================================================= |
582 | |
583 | BSplCLib_KnotDistribution BSplCLib::KnotForm |
584 | (const Array1OfReal& Knots, |
585 | const Standard_Integer FromK1, |
586 | const Standard_Integer ToK2) |
587 | { |
588 | Standard_Real DU0,DU1,Ui,Uj,Eps0,val; |
589 | BSplCLib_KnotDistribution KForm = BSplCLib_Uniform; |
590 | |
591 | Standard_Integer KLower = Knots.Lower(); |
592 | const Standard_Real * pkn = &Knots(KLower); |
593 | pkn -= KLower; |
594 | Ui = pkn[FromK1]; |
595 | if (Ui < 0) Ui = - Ui; |
596 | Uj = pkn[FromK1 + 1]; |
597 | if (Uj < 0) Uj = - Uj; |
598 | DU0 = Uj - Ui; |
599 | if (DU0 < 0) DU0 = - DU0; |
600 | Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0); |
601 | Standard_Integer i = FromK1 + 1; |
602 | |
603 | while (KForm != BSplCLib_NonUniform && i < ToK2) { |
604 | Ui = pkn[i]; |
605 | if (Ui < 0) Ui = - Ui; |
606 | i++; |
607 | Uj = pkn[i]; |
608 | if (Uj < 0) Uj = - Uj; |
609 | DU1 = Uj - Ui; |
610 | if (DU1 < 0) DU1 = - DU1; |
611 | val = DU1 - DU0; |
612 | if (val < 0) val = -val; |
613 | if (val > Eps0) KForm = BSplCLib_NonUniform; |
614 | DU0 = DU1; |
615 | Eps0 = Epsilon (Ui) + Epsilon (Uj) + Epsilon (DU0); |
616 | } |
617 | return KForm; |
618 | } |
619 | |
620 | //======================================================================= |
621 | //function : MultForm |
622 | //purpose : |
623 | //======================================================================= |
624 | |
625 | BSplCLib_MultDistribution BSplCLib::MultForm |
626 | (const Array1OfInteger& Mults, |
627 | const Standard_Integer FromK1, |
628 | const Standard_Integer ToK2) |
629 | { |
630 | Standard_Integer First,Last; |
631 | if (FromK1 < ToK2) { |
632 | First = FromK1; |
633 | Last = ToK2; |
634 | } |
635 | else { |
636 | First = ToK2; |
637 | Last = FromK1; |
638 | } |
639 | Standard_Integer MLower = Mults.Lower(); |
640 | const Standard_Integer *pmu = &Mults(MLower); |
641 | pmu -= MLower; |
642 | Standard_Integer FirstMult = pmu[First]; |
643 | BSplCLib_MultDistribution MForm = BSplCLib_Constant; |
644 | Standard_Integer i = First + 1; |
645 | Standard_Integer Mult = pmu[i]; |
646 | |
647 | // while (MForm != BSplCLib_NonUniform && i <= Last) { ???????????JR???????? |
648 | while (MForm != BSplCLib_NonConstant && i <= Last) { |
649 | if (i == First + 1) { |
650 | if (Mult != FirstMult) MForm = BSplCLib_QuasiConstant; |
651 | } |
652 | else if (i == Last) { |
653 | if (MForm == BSplCLib_QuasiConstant) { |
654 | if (FirstMult != pmu[i]) MForm = BSplCLib_NonConstant; |
655 | } |
656 | else { |
657 | if (Mult != pmu[i]) MForm = BSplCLib_NonConstant; |
658 | } |
659 | } |
660 | else { |
661 | if (Mult != pmu[i]) MForm = BSplCLib_NonConstant; |
662 | Mult = pmu[i]; |
663 | } |
664 | i++; |
665 | } |
666 | return MForm; |
667 | } |
668 | |
669 | //======================================================================= |
670 | //function : Reparametrize |
671 | //purpose : |
672 | //======================================================================= |
673 | |
674 | void BSplCLib::Reparametrize |
675 | (const Standard_Real U1, |
676 | const Standard_Real U2, |
677 | Array1OfReal& Knots) |
678 | { |
679 | Standard_Integer Lower = Knots.Lower(); |
680 | Standard_Integer Upper = Knots.Upper(); |
681 | Standard_Real UFirst = Min (U1, U2); |
682 | Standard_Real ULast = Max (U1, U2); |
683 | Standard_Real NewLength = ULast - UFirst; |
684 | BSplCLib_KnotDistribution KSet = BSplCLib::KnotForm (Knots, Lower, Upper); |
685 | if (KSet == BSplCLib_Uniform) { |
686 | Standard_Real DU = NewLength / (Upper - Lower); |
687 | Knots (Lower) = UFirst; |
688 | |
689 | for (Standard_Integer i = Lower + 1; i <= Upper; i++) { |
690 | Knots (i) = Knots (i-1) + DU; |
691 | } |
692 | } |
693 | else { |
694 | Standard_Real K2; |
695 | Standard_Real Ratio; |
696 | Standard_Real K1 = Knots (Lower); |
697 | Standard_Real Length = Knots (Upper) - Knots (Lower); |
698 | Knots (Lower) = UFirst; |
699 | |
700 | for (Standard_Integer i = Lower + 1; i <= Upper; i++) { |
701 | K2 = Knots (i); |
702 | Ratio = (K2 - K1) / Length; |
703 | Knots (i) = Knots (i-1) + (NewLength * Ratio); |
704 | |
705 | //for CheckCurveData |
706 | Standard_Real Eps = Epsilon( Abs(Knots(i-1)) ); |
707 | if (Knots(i) - Knots(i-1) <= Eps) |
708 | Knots(i) += 1.1*Eps; |
709 | |
710 | K1 = K2; |
711 | } |
712 | } |
713 | } |
714 | |
715 | //======================================================================= |
716 | //function : Reverse |
717 | //purpose : |
718 | //======================================================================= |
719 | |
720 | void BSplCLib::Reverse(TColStd_Array1OfReal& Knots) |
721 | { |
722 | Standard_Integer first = Knots.Lower(); |
723 | Standard_Integer last = Knots.Upper(); |
724 | Standard_Real kfirst = Knots(first); |
725 | Standard_Real klast = Knots(last); |
726 | Standard_Real tfirst = kfirst; |
727 | Standard_Real tlast = klast; |
728 | first++; |
729 | last--; |
730 | |
731 | while (first <= last) { |
732 | tfirst += klast - Knots(last); |
733 | tlast -= Knots(first) - kfirst; |
734 | kfirst = Knots(first); |
735 | klast = Knots(last); |
736 | Knots(first) = tfirst; |
737 | Knots(last) = tlast; |
738 | first++; |
739 | last--; |
740 | } |
741 | } |
742 | |
743 | //======================================================================= |
744 | //function : Reverse |
745 | //purpose : |
746 | //======================================================================= |
747 | |
748 | void BSplCLib::Reverse(TColStd_Array1OfInteger& Mults) |
749 | { |
750 | Standard_Integer first = Mults.Lower(); |
751 | Standard_Integer last = Mults.Upper(); |
752 | Standard_Integer temp; |
753 | |
754 | while (first < last) { |
755 | temp = Mults(first); |
756 | Mults(first) = Mults(last); |
757 | Mults(last) = temp; |
758 | first++; |
759 | last--; |
760 | } |
761 | } |
762 | |
763 | //======================================================================= |
764 | //function : Reverse |
765 | //purpose : |
766 | //======================================================================= |
767 | |
768 | void BSplCLib::Reverse(TColStd_Array1OfReal& Weights, |
769 | const Standard_Integer L) |
770 | { |
771 | Standard_Integer i, l = L; |
772 | l = Weights.Lower()+(l-Weights.Lower())%(Weights.Upper()-Weights.Lower()+1); |
773 | |
774 | TColStd_Array1OfReal temp(0,Weights.Length()-1); |
775 | |
776 | for (i = Weights.Lower(); i <= l; i++) |
777 | temp(l-i) = Weights(i); |
778 | |
779 | for (i = l+1; i <= Weights.Upper(); i++) |
780 | temp(l-Weights.Lower()+Weights.Upper()-i+1) = Weights(i); |
781 | |
782 | for (i = Weights.Lower(); i <= Weights.Upper(); i++) |
783 | Weights(i) = temp(i-Weights.Lower()); |
784 | } |
785 | |
786 | //======================================================================= |
787 | //function : IsRational |
788 | //purpose : |
789 | //======================================================================= |
790 | |
791 | Standard_Boolean BSplCLib::IsRational(const TColStd_Array1OfReal& Weights, |
792 | const Standard_Integer I1, |
793 | const Standard_Integer I2, |
794 | // const Standard_Real Epsi) |
795 | const Standard_Real ) |
796 | { |
797 | Standard_Integer i, f = Weights.Lower(), l = Weights.Length(); |
798 | Standard_Integer I3 = I2 - f; |
799 | const Standard_Real * WG = &Weights(f); |
800 | WG -= f; |
801 | |
802 | for (i = I1 - f; i < I3; i++) { |
803 | if (WG[f + (i % l)] != WG[f + ((i + 1) % l)]) return Standard_True; |
804 | } |
805 | return Standard_False ; |
806 | } |
807 | |
808 | //======================================================================= |
809 | //function : Eval |
810 | //purpose : evaluate point and derivatives |
811 | //======================================================================= |
812 | |
813 | void BSplCLib::Eval(const Standard_Real U, |
814 | const Standard_Integer Degree, |
815 | Standard_Real& Knots, |
816 | const Standard_Integer Dimension, |
817 | Standard_Real& Poles) |
818 | { |
819 | Standard_Integer step,i,Dms,Dm1,Dpi,Sti; |
820 | Standard_Real X, Y, *poles, *knots = &Knots; |
821 | Dm1 = Dms = Degree; |
822 | Dm1--; |
823 | Dms++; |
824 | switch (Dimension) { |
825 | |
826 | case 1 : { |
827 | |
828 | for (step = - 1; step < Dm1; step++) { |
829 | Dms--; |
830 | poles = &Poles; |
831 | Dpi = Dm1; |
832 | Sti = step; |
833 | |
834 | for (i = 0; i < Dms; i++) { |
835 | Dpi++; |
836 | Sti++; |
837 | X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]); |
838 | Y = 1 - X; |
839 | poles[0] *= X; poles[0] += Y * poles[1]; |
840 | poles += 1; |
841 | } |
842 | } |
843 | break; |
844 | } |
845 | case 2 : { |
846 | |
847 | for (step = - 1; step < Dm1; step++) { |
848 | Dms--; |
849 | poles = &Poles; |
850 | Dpi = Dm1; |
851 | Sti = step; |
852 | |
853 | for (i = 0; i < Dms; i++) { |
854 | Dpi++; |
855 | Sti++; |
856 | X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]); |
857 | Y = 1 - X; |
858 | poles[0] *= X; poles[0] += Y * poles[2]; |
859 | poles[1] *= X; poles[1] += Y * poles[3]; |
860 | poles += 2; |
861 | } |
862 | } |
863 | break; |
864 | } |
865 | case 3 : { |
866 | |
867 | for (step = - 1; step < Dm1; step++) { |
868 | Dms--; |
869 | poles = &Poles; |
870 | Dpi = Dm1; |
871 | Sti = step; |
872 | |
873 | for (i = 0; i < Dms; i++) { |
874 | Dpi++; |
875 | Sti++; |
876 | X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]); |
877 | Y = 1 - X; |
878 | poles[0] *= X; poles[0] += Y * poles[3]; |
879 | poles[1] *= X; poles[1] += Y * poles[4]; |
880 | poles[2] *= X; poles[2] += Y * poles[5]; |
881 | poles += 3; |
882 | } |
883 | } |
884 | break; |
885 | } |
886 | case 4 : { |
887 | |
888 | for (step = - 1; step < Dm1; step++) { |
889 | Dms--; |
890 | poles = &Poles; |
891 | Dpi = Dm1; |
892 | Sti = step; |
893 | |
894 | for (i = 0; i < Dms; i++) { |
895 | Dpi++; |
896 | Sti++; |
897 | X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]); |
898 | Y = 1 - X; |
899 | poles[0] *= X; poles[0] += Y * poles[4]; |
900 | poles[1] *= X; poles[1] += Y * poles[5]; |
901 | poles[2] *= X; poles[2] += Y * poles[6]; |
902 | poles[3] *= X; poles[3] += Y * poles[7]; |
903 | poles += 4; |
904 | } |
905 | } |
906 | break; |
907 | } |
908 | default : { |
909 | Standard_Integer k; |
910 | |
911 | for (step = - 1; step < Dm1; step++) { |
912 | Dms--; |
913 | poles = &Poles; |
914 | Dpi = Dm1; |
915 | Sti = step; |
916 | |
917 | for (i = 0; i < Dms; i++) { |
918 | Dpi++; |
919 | Sti++; |
920 | X = (knots[Dpi] - U) / (knots[Dpi] - knots[Sti]); |
921 | Y = 1 - X; |
922 | |
923 | for (k = 0; k < Dimension; k++) { |
924 | poles[k] *= X; |
925 | poles[k] += Y * poles[k + Dimension]; |
926 | } |
927 | poles += Dimension; |
928 | } |
929 | } |
930 | } |
931 | } |
932 | } |
933 | |
934 | //======================================================================= |
935 | //function : BoorScheme |
936 | //purpose : |
937 | //======================================================================= |
938 | |
939 | void BSplCLib::BoorScheme(const Standard_Real U, |
940 | const Standard_Integer Degree, |
941 | Standard_Real& Knots, |
942 | const Standard_Integer Dimension, |
943 | Standard_Real& Poles, |
944 | const Standard_Integer Depth, |
945 | const Standard_Integer Length) |
946 | { |
947 | // |
948 | // Compute the values |
949 | // |
950 | // P(i,j) (U). |
951 | // |
952 | // for i = 0 to Depth, |
953 | // j = 0 to Length - i |
954 | // |
955 | // The Boor scheme is : |
956 | // |
957 | // P(0,j) = Pole(j) |
958 | // P(i,j) = x * P(i-1,j) + (1-x) * P(i-1,j+1) |
959 | // |
960 | // where x = (knot(i+j+Degree) - U) / (knot(i+j+Degree) - knot(i+j)) |
961 | // |
962 | // |
963 | // The values are stored in the array Poles |
964 | // They are alternatively written if the odd and even positions. |
965 | // |
966 | // The successives contents of the array are |
967 | // ***** means unitialised, l = Degree + Length |
968 | // |
969 | // P(0,0) ****** P(0,1) ...... P(0,l-1) ******** P(0,l) |
970 | // P(0,0) P(1,0) P(0,1) ...... P(0,l-1) P(1,l-1) P(0,l) |
971 | // P(0,0) P(1,0) P(2,0) ...... P(2,l-1) P(1,l-1) P(0,l) |
972 | // |
973 | |
974 | Standard_Integer i,k,step; |
975 | Standard_Real *knots = &Knots; |
976 | Standard_Real *pole, *firstpole = &Poles - 2 * Dimension; |
977 | // the steps of recursion |
978 | |
979 | for (step = 0; step < Depth; step++) { |
980 | firstpole += Dimension; |
981 | pole = firstpole; |
982 | // compute the new row of poles |
983 | |
984 | for (i = step; i < Length; i++) { |
985 | pole += 2 * Dimension; |
986 | // coefficient |
987 | Standard_Real X = (knots[i+Degree-step] - U) |
988 | / (knots[i+Degree-step] - knots[i]); |
989 | Standard_Real Y = 1. - X; |
990 | // Value |
991 | // P(i,j) = X * P(i-1,j) + (1-X) * P(i-1,j+1) |
992 | |
993 | for (k = 0; k < Dimension; k++) |
994 | pole[k] = X * pole[k - Dimension] + Y * pole[k + Dimension]; |
995 | } |
996 | } |
997 | } |
998 | |
999 | //======================================================================= |
1000 | //function : AntiBoorScheme |
1001 | //purpose : |
1002 | //======================================================================= |
1003 | |
1004 | Standard_Boolean BSplCLib::AntiBoorScheme(const Standard_Real U, |
1005 | const Standard_Integer Degree, |
1006 | Standard_Real& Knots, |
1007 | const Standard_Integer Dimension, |
1008 | Standard_Real& Poles, |
1009 | const Standard_Integer Depth, |
1010 | const Standard_Integer Length, |
1011 | const Standard_Real Tolerance) |
1012 | { |
1013 | // do the Boor scheme reverted. |
1014 | |
1015 | Standard_Integer i,k,step, half_length; |
1016 | Standard_Real *knots = &Knots; |
1017 | Standard_Real z,X,Y,*pole, *firstpole = &Poles + (Depth-1) * Dimension; |
1018 | |
1019 | // Test the special case length = 1 |
1020 | // only verification of the central point |
1021 | |
1022 | if (Length == 1) { |
1023 | X = (knots[Degree] - U) / (knots[Degree] - knots[0]); |
1024 | Y = 1. - X; |
1025 | |
1026 | for (k = 0; k < Dimension; k++) { |
1027 | z = X * firstpole[k] + Y * firstpole[k+2*Dimension]; |
1028 | if (Abs(z - firstpole[k+Dimension]) > Tolerance) |
1029 | return Standard_False; |
1030 | } |
1031 | return Standard_True; |
1032 | } |
1033 | |
1034 | // General case |
1035 | // the steps of recursion |
1036 | |
1037 | for (step = Depth-1; step >= 0; step--) { |
1038 | firstpole -= Dimension; |
1039 | pole = firstpole; |
1040 | |
1041 | // first step from left to right |
1042 | |
1043 | for (i = step; i < Length-1; i++) { |
1044 | pole += 2 * Dimension; |
1045 | |
1046 | X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]); |
1047 | Y = 1. - X; |
1048 | |
1049 | for (k = 0; k < Dimension; k++) |
1050 | pole[k+Dimension] = (pole[k] - X*pole[k-Dimension]) / Y; |
1051 | |
1052 | } |
1053 | |
1054 | // second step from right to left |
1055 | pole += 4* Dimension; |
1056 | half_length = (Length - 1 + step) / 2 ; |
1057 | // |
1058 | // only do half of the way from right to left |
1059 | // otherwise it start degenerating because of |
1060 | // overflows |
1061 | // |
1062 | |
1063 | for (i = Length-1; i > half_length ; i--) { |
1064 | pole -= 2 * Dimension; |
1065 | |
1066 | // coefficient |
1067 | X = (knots[i+Degree-step] - U) / (knots[i+Degree-step] - knots[i]); |
1068 | Y = 1. - X; |
1069 | |
1070 | for (k = 0; k < Dimension; k++) { |
1071 | z = (pole[k] - Y * pole[k+Dimension]) / X; |
1072 | if (Abs(z-pole[k-Dimension]) > Tolerance) |
1073 | return Standard_False; |
1074 | pole[k-Dimension] += z; |
1075 | pole[k-Dimension] /= 2.; |
1076 | } |
1077 | } |
1078 | } |
1079 | return Standard_True; |
1080 | } |
1081 | |
1082 | //======================================================================= |
1083 | //function : Derivative |
1084 | //purpose : |
1085 | //======================================================================= |
1086 | |
1087 | void BSplCLib::Derivative(const Standard_Integer Degree, |
1088 | Standard_Real& Knots, |
1089 | const Standard_Integer Dimension, |
1090 | const Standard_Integer Length, |
1091 | const Standard_Integer Order, |
1092 | Standard_Real& Poles) |
1093 | { |
1094 | Standard_Integer i,k,step,span = Degree; |
1095 | Standard_Real *knot = &Knots; |
1096 | |
1097 | for (step = 1; step <= Order; step++) { |
1098 | Standard_Real* pole = &Poles; |
1099 | |
1100 | for (i = step; i < Length; i++) { |
1101 | Standard_Real coef = - span / (knot[i+span] - knot[i]); |
1102 | |
1103 | for (k = 0; k < Dimension; k++) { |
1104 | pole[k] -= pole[k+Dimension]; |
1105 | pole[k] *= coef; |
1106 | } |
1107 | pole += Dimension; |
1108 | } |
1109 | span--; |
1110 | } |
1111 | } |
1112 | |
1113 | //======================================================================= |
1114 | //function : Bohm |
1115 | //purpose : |
1116 | //======================================================================= |
1117 | |
1118 | void BSplCLib::Bohm(const Standard_Real U, |
1119 | const Standard_Integer Degree, |
1120 | const Standard_Integer N, |
1121 | Standard_Real& Knots, |
1122 | const Standard_Integer Dimension, |
1123 | Standard_Real& Poles) |
1124 | { |
1125 | // First phase independant of U, compute the poles of the derivatives |
1126 | Standard_Integer i,j,iDim,min,Dmi,DDmi,jDmi,Degm1; |
1127 | Standard_Real *knot = &Knots, *pole, coef, *tbis, *psav, *psDD, *psDDmDim; |
1128 | psav = &Poles; |
1129 | if (N < Degree) min = N; |
1130 | else min = Degree; |
1131 | Degm1 = Degree - 1; |
1132 | DDmi = (Degree << 1) + 1; |
1133 | switch (Dimension) { |
1134 | case 1 : { |
1135 | psDD = psav + Degree; |
1136 | psDDmDim = psDD - 1; |
1137 | |
1138 | for (i = 0; i < Degree; i++) { |
1139 | DDmi--; |
1140 | pole = psDD; |
1141 | tbis = psDDmDim; |
1142 | jDmi = DDmi; |
1143 | |
1144 | for (j = Degm1; j >= i; j--) { |
1145 | jDmi--; |
1146 | *pole -= *tbis; *pole /= (knot[jDmi] - knot[j]); |
1147 | pole--; |
1148 | tbis--; |
1149 | } |
1150 | } |
1151 | // Second phase, dependant of U |
1152 | iDim = - 1; |
1153 | |
1154 | for (i = 0; i < Degree; i++) { |
1155 | iDim += 1; |
1156 | pole = psav + iDim; |
1157 | tbis = pole + 1; |
1158 | coef = U - knot[i]; |
1159 | |
1160 | for (j = i; j >= 0; j--) { |
1161 | *pole += coef * (*tbis); |
1162 | pole--; |
1163 | tbis--; |
1164 | } |
1165 | } |
1166 | // multiply by the degrees |
1167 | coef = Degree; |
1168 | Dmi = Degree; |
1169 | pole = psav + 1; |
1170 | |
1171 | for (i = 1; i <= min; i++) { |
1172 | *pole *= coef; pole++; |
1173 | Dmi--; |
1174 | coef *= Dmi; |
1175 | } |
1176 | break; |
1177 | } |
1178 | case 2 : { |
1179 | psDD = psav + (Degree << 1); |
1180 | psDDmDim = psDD - 2; |
1181 | |
1182 | for (i = 0; i < Degree; i++) { |
1183 | DDmi--; |
1184 | pole = psDD; |
1185 | tbis = psDDmDim; |
1186 | jDmi = DDmi; |
1187 | |
1188 | for (j = Degm1; j >= i; j--) { |
1189 | jDmi--; |
1190 | coef = 1. / (knot[jDmi] - knot[j]); |
1191 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1192 | *pole -= *tbis; *pole *= coef; |
1193 | pole -= 3; |
1194 | tbis -= 3; |
1195 | } |
1196 | } |
1197 | // Second phase, dependant of U |
1198 | iDim = - 2; |
1199 | |
1200 | for (i = 0; i < Degree; i++) { |
1201 | iDim += 2; |
1202 | pole = psav + iDim; |
1203 | tbis = pole + 2; |
1204 | coef = U - knot[i]; |
1205 | |
1206 | for (j = i; j >= 0; j--) { |
1207 | *pole += coef * (*tbis); pole++; tbis++; |
1208 | *pole += coef * (*tbis); |
1209 | pole -= 3; |
1210 | tbis -= 3; |
1211 | } |
1212 | } |
1213 | // multiply by the degrees |
1214 | coef = Degree; |
1215 | Dmi = Degree; |
1216 | pole = psav + 2; |
1217 | |
1218 | for (i = 1; i <= min; i++) { |
1219 | *pole *= coef; pole++; |
1220 | *pole *= coef; pole++; |
1221 | Dmi--; |
1222 | coef *= Dmi; |
1223 | } |
1224 | break; |
1225 | } |
1226 | case 3 : { |
1227 | psDD = psav + (Degree << 1) + Degree; |
1228 | psDDmDim = psDD - 3; |
1229 | |
1230 | for (i = 0; i < Degree; i++) { |
1231 | DDmi--; |
1232 | pole = psDD; |
1233 | tbis = psDDmDim; |
1234 | jDmi = DDmi; |
1235 | |
1236 | for (j = Degm1; j >= i; j--) { |
1237 | jDmi--; |
1238 | coef = 1. / (knot[jDmi] - knot[j]); |
1239 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1240 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1241 | *pole -= *tbis; *pole *= coef; |
1242 | pole -= 5; |
1243 | tbis -= 5; |
1244 | } |
1245 | } |
1246 | // Second phase, dependant of U |
1247 | iDim = - 3; |
1248 | |
1249 | for (i = 0; i < Degree; i++) { |
1250 | iDim += 3; |
1251 | pole = psav + iDim; |
1252 | tbis = pole + 3; |
1253 | coef = U - knot[i]; |
1254 | |
1255 | for (j = i; j >= 0; j--) { |
1256 | *pole += coef * (*tbis); pole++; tbis++; |
1257 | *pole += coef * (*tbis); pole++; tbis++; |
1258 | *pole += coef * (*tbis); |
1259 | pole -= 5; |
1260 | tbis -= 5; |
1261 | } |
1262 | } |
1263 | // multiply by the degrees |
1264 | coef = Degree; |
1265 | Dmi = Degree; |
1266 | pole = psav + 3; |
1267 | |
1268 | for (i = 1; i <= min; i++) { |
1269 | *pole *= coef; pole++; |
1270 | *pole *= coef; pole++; |
1271 | *pole *= coef; pole++; |
1272 | Dmi--; |
1273 | coef *= Dmi; |
1274 | } |
1275 | break; |
1276 | } |
1277 | case 4 : { |
1278 | psDD = psav + (Degree << 2); |
1279 | psDDmDim = psDD - 4; |
1280 | |
1281 | for (i = 0; i < Degree; i++) { |
1282 | DDmi--; |
1283 | pole = psDD; |
1284 | tbis = psDDmDim; |
1285 | jDmi = DDmi; |
1286 | |
1287 | for (j = Degm1; j >= i; j--) { |
1288 | jDmi--; |
1289 | coef = 1. / (knot[jDmi] - knot[j]); |
1290 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1291 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1292 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1293 | *pole -= *tbis; *pole *= coef; |
1294 | pole -= 7; |
1295 | tbis -= 7; |
1296 | } |
1297 | } |
1298 | // Second phase, dependant of U |
1299 | iDim = - 4; |
1300 | |
1301 | for (i = 0; i < Degree; i++) { |
1302 | iDim += 4; |
1303 | pole = psav + iDim; |
1304 | tbis = pole + 4; |
1305 | coef = U - knot[i]; |
1306 | |
1307 | for (j = i; j >= 0; j--) { |
1308 | *pole += coef * (*tbis); pole++; tbis++; |
1309 | *pole += coef * (*tbis); pole++; tbis++; |
1310 | *pole += coef * (*tbis); pole++; tbis++; |
1311 | *pole += coef * (*tbis); |
1312 | pole -= 7; |
1313 | tbis -= 7; |
1314 | } |
1315 | } |
1316 | // multiply by the degrees |
1317 | coef = Degree; |
1318 | Dmi = Degree; |
1319 | pole = psav + 4; |
1320 | |
1321 | for (i = 1; i <= min; i++) { |
1322 | *pole *= coef; pole++; |
1323 | *pole *= coef; pole++; |
1324 | *pole *= coef; pole++; |
1325 | *pole *= coef; pole++; |
1326 | Dmi--; |
1327 | coef *= Dmi; |
1328 | } |
1329 | break; |
1330 | } |
1331 | default : { |
1332 | Standard_Integer k; |
1333 | Standard_Integer Dim2 = Dimension << 1; |
1334 | psDD = psav + Degree * Dimension; |
1335 | psDDmDim = psDD - Dimension; |
1336 | |
1337 | for (i = 0; i < Degree; i++) { |
1338 | DDmi--; |
1339 | pole = psDD; |
1340 | tbis = psDDmDim; |
1341 | jDmi = DDmi; |
1342 | |
1343 | for (j = Degm1; j >= i; j--) { |
1344 | jDmi--; |
1345 | coef = 1. / (knot[jDmi] - knot[j]); |
1346 | |
1347 | for (k = 0; k < Dimension; k++) { |
1348 | *pole -= *tbis; *pole *= coef; pole++; tbis++; |
1349 | } |
1350 | pole -= Dim2; |
1351 | tbis -= Dim2; |
1352 | } |
1353 | } |
1354 | // Second phase, dependant of U |
1355 | iDim = - Dimension; |
1356 | |
1357 | for (i = 0; i < Degree; i++) { |
1358 | iDim += Dimension; |
1359 | pole = psav + iDim; |
1360 | tbis = pole + Dimension; |
1361 | coef = U - knot[i]; |
1362 | |
1363 | for (j = i; j >= 0; j--) { |
1364 | |
1365 | for (k = 0; k < Dimension; k++) { |
1366 | *pole += coef * (*tbis); pole++; tbis++; |
1367 | } |
1368 | pole -= Dim2; |
1369 | tbis -= Dim2; |
1370 | } |
1371 | } |
1372 | // multiply by the degrees |
1373 | coef = Degree; |
1374 | Dmi = Degree; |
1375 | pole = psav + Dimension; |
1376 | |
1377 | for (i = 1; i <= min; i++) { |
1378 | |
1379 | for (k = 0; k < Dimension; k++) { |
1380 | *pole *= coef; pole++; |
1381 | } |
1382 | Dmi--; |
1383 | coef *= Dmi; |
1384 | } |
1385 | } |
1386 | } |
1387 | } |
1388 | |
1389 | //======================================================================= |
1390 | //function : BuildKnots |
1391 | //purpose : |
1392 | //======================================================================= |
1393 | |
1394 | void BSplCLib::BuildKnots(const Standard_Integer Degree, |
1395 | const Standard_Integer Index, |
1396 | const Standard_Boolean Periodic, |
1397 | const TColStd_Array1OfReal& Knots, |
1398 | const TColStd_Array1OfInteger& Mults, |
1399 | Standard_Real& LK) |
1400 | { |
1401 | Standard_Integer KLower = Knots.Lower(); |
1402 | const Standard_Real * pkn = &Knots(KLower); |
1403 | pkn -= KLower; |
1404 | Standard_Real *knot = &LK; |
1405 | if (&Mults == NULL) { |
1406 | switch (Degree) { |
1407 | case 1 : { |
1408 | Standard_Integer j = Index ; |
1409 | knot[0] = pkn[j]; j++; |
1410 | knot[1] = pkn[j]; |
1411 | break; |
1412 | } |
1413 | case 2 : { |
1414 | Standard_Integer j = Index - 1; |
1415 | knot[0] = pkn[j]; j++; |
1416 | knot[1] = pkn[j]; j++; |
1417 | knot[2] = pkn[j]; j++; |
1418 | knot[3] = pkn[j]; |
1419 | break; |
1420 | } |
1421 | case 3 : { |
1422 | Standard_Integer j = Index - 2; |
1423 | knot[0] = pkn[j]; j++; |
1424 | knot[1] = pkn[j]; j++; |
1425 | knot[2] = pkn[j]; j++; |
1426 | knot[3] = pkn[j]; j++; |
1427 | knot[4] = pkn[j]; j++; |
1428 | knot[5] = pkn[j]; |
1429 | break; |
1430 | } |
1431 | case 4 : { |
1432 | Standard_Integer j = Index - 3; |
1433 | knot[0] = pkn[j]; j++; |
1434 | knot[1] = pkn[j]; j++; |
1435 | knot[2] = pkn[j]; j++; |
1436 | knot[3] = pkn[j]; j++; |
1437 | knot[4] = pkn[j]; j++; |
1438 | knot[5] = pkn[j]; j++; |
1439 | knot[6] = pkn[j]; j++; |
1440 | knot[7] = pkn[j]; |
1441 | break; |
1442 | } |
1443 | case 5 : { |
1444 | Standard_Integer j = Index - 4; |
1445 | knot[0] = pkn[j]; j++; |
1446 | knot[1] = pkn[j]; j++; |
1447 | knot[2] = pkn[j]; j++; |
1448 | knot[3] = pkn[j]; j++; |
1449 | knot[4] = pkn[j]; j++; |
1450 | knot[5] = pkn[j]; j++; |
1451 | knot[6] = pkn[j]; j++; |
1452 | knot[7] = pkn[j]; j++; |
1453 | knot[8] = pkn[j]; j++; |
1454 | knot[9] = pkn[j]; |
1455 | break; |
1456 | } |
1457 | case 6 : { |
1458 | Standard_Integer j = Index - 5; |
1459 | knot[ 0] = pkn[j]; j++; |
1460 | knot[ 1] = pkn[j]; j++; |
1461 | knot[ 2] = pkn[j]; j++; |
1462 | knot[ 3] = pkn[j]; j++; |
1463 | knot[ 4] = pkn[j]; j++; |
1464 | knot[ 5] = pkn[j]; j++; |
1465 | knot[ 6] = pkn[j]; j++; |
1466 | knot[ 7] = pkn[j]; j++; |
1467 | knot[ 8] = pkn[j]; j++; |
1468 | knot[ 9] = pkn[j]; j++; |
1469 | knot[10] = pkn[j]; j++; |
1470 | knot[11] = pkn[j]; |
1471 | break; |
1472 | } |
1473 | default : { |
1474 | Standard_Integer i,j; |
1475 | Standard_Integer Deg2 = Degree << 1; |
1476 | j = Index - Degree; |
1477 | |
1478 | for (i = 0; i < Deg2; i++) { |
1479 | j++; |
1480 | knot[i] = pkn[j]; |
1481 | } |
1482 | } |
1483 | } |
1484 | } |
1485 | else { |
1486 | Standard_Integer i; |
1487 | Standard_Integer Deg1 = Degree - 1; |
1488 | Standard_Integer KUpper = Knots.Upper(); |
1489 | Standard_Integer MLower = Mults.Lower(); |
1490 | Standard_Integer MUpper = Mults.Upper(); |
1491 | const Standard_Integer * pmu = &Mults(MLower); |
1492 | pmu -= MLower; |
1493 | Standard_Real dknot = 0; |
1494 | Standard_Integer ilow = Index , mlow = 0; |
1495 | Standard_Integer iupp = Index + 1, mupp = 0; |
1496 | Standard_Real loffset = 0., uoffset = 0.; |
1497 | Standard_Boolean getlow = Standard_True, getupp = Standard_True; |
1498 | if (Periodic) { |
1499 | dknot = pkn[KUpper] - pkn[KLower]; |
1500 | if (iupp > MUpper) { |
1501 | iupp = MLower + 1; |
1502 | uoffset = dknot; |
1503 | } |
1504 | } |
1505 | // Find the knots around Index |
1506 | |
1507 | for (i = 0; i < Degree; i++) { |
1508 | if (getlow) { |
1509 | mlow++; |
1510 | if (mlow > pmu[ilow]) { |
1511 | mlow = 1; |
1512 | ilow--; |
1513 | getlow = (ilow >= MLower); |
1514 | if (Periodic && !getlow) { |
1515 | ilow = MUpper - 1; |
1516 | loffset = dknot; |
1517 | getlow = Standard_True; |
1518 | } |
1519 | } |
1520 | if (getlow) |
1521 | knot[Deg1 - i] = pkn[ilow] - loffset; |
1522 | } |
1523 | if (getupp) { |
1524 | mupp++; |
1525 | if (mupp > pmu[iupp]) { |
1526 | mupp = 1; |
1527 | iupp++; |
1528 | getupp = (iupp <= MUpper); |
1529 | if (Periodic && !getupp) { |
1530 | iupp = MLower + 1; |
1531 | uoffset = dknot; |
1532 | getupp = Standard_True; |
1533 | } |
1534 | } |
1535 | if (getupp) |
1536 | knot[Degree + i] = pkn[iupp] + uoffset; |
1537 | } |
1538 | } |
1539 | } |
1540 | } |
1541 | |
1542 | //======================================================================= |
1543 | //function : PoleIndex |
1544 | //purpose : |
1545 | //======================================================================= |
1546 | |
1547 | Standard_Integer BSplCLib::PoleIndex(const Standard_Integer Degree, |
1548 | const Standard_Integer Index, |
1549 | const Standard_Boolean Periodic, |
1550 | const TColStd_Array1OfInteger& Mults) |
1551 | { |
1552 | Standard_Integer i, pindex = 0; |
1553 | |
1554 | for (i = Mults.Lower(); i <= Index; i++) |
1555 | pindex += Mults(i); |
1556 | if (Periodic) |
1557 | pindex -= Mults(Mults.Lower()); |
1558 | else |
1559 | pindex -= Degree + 1; |
1560 | |
1561 | return pindex; |
1562 | } |
1563 | |
1564 | //======================================================================= |
1565 | //function : BuildBoor |
1566 | //purpose : builds the local array for boor |
1567 | //======================================================================= |
1568 | |
1569 | void BSplCLib::BuildBoor(const Standard_Integer Index, |
1570 | const Standard_Integer Length, |
1571 | const Standard_Integer Dimension, |
1572 | const TColStd_Array1OfReal& Poles, |
1573 | Standard_Real& LP) |
1574 | { |
1575 | Standard_Real *poles = &LP; |
1576 | Standard_Integer i,k, ip = Poles.Lower() + Index * Dimension; |
1577 | |
1578 | for (i = 0; i < Length+1; i++) { |
1579 | |
1580 | for (k = 0; k < Dimension; k++) { |
1581 | poles[k] = Poles(ip); |
1582 | ip++; |
1583 | if (ip > Poles.Upper()) ip = Poles.Lower(); |
1584 | } |
1585 | poles += 2 * Dimension; |
1586 | } |
1587 | } |
1588 | |
1589 | //======================================================================= |
1590 | //function : BoorIndex |
1591 | //purpose : |
1592 | //======================================================================= |
1593 | |
1594 | Standard_Integer BSplCLib::BoorIndex(const Standard_Integer Index, |
1595 | const Standard_Integer Length, |
1596 | const Standard_Integer Depth) |
1597 | { |
1598 | if (Index <= Depth) return Index; |
1599 | if (Index <= Length) return 2 * Index - Depth; |
1600 | return Length + Index - Depth; |
1601 | } |
1602 | |
1603 | //======================================================================= |
1604 | //function : GetPole |
1605 | //purpose : |
1606 | //======================================================================= |
1607 | |
1608 | void BSplCLib::GetPole(const Standard_Integer Index, |
1609 | const Standard_Integer Length, |
1610 | const Standard_Integer Depth, |
1611 | const Standard_Integer Dimension, |
1612 | Standard_Real& LP, |
1613 | Standard_Integer& Position, |
1614 | TColStd_Array1OfReal& Pole) |
1615 | { |
1616 | Standard_Integer k; |
1617 | Standard_Real* pole = &LP + BoorIndex(Index,Length,Depth) * Dimension; |
1618 | |
1619 | for (k = 0; k < Dimension; k++) { |
1620 | Pole(Position) = pole[k]; |
1621 | Position++; |
1622 | } |
1623 | if (Position > Pole.Upper()) Position = Pole.Lower(); |
1624 | } |
1625 | |
1626 | //======================================================================= |
1627 | //function : PrepareInsertKnots |
1628 | //purpose : |
1629 | //======================================================================= |
1630 | |
1631 | Standard_Boolean BSplCLib::PrepareInsertKnots |
1632 | (const Standard_Integer Degree, |
1633 | const Standard_Boolean Periodic, |
1634 | const TColStd_Array1OfReal& Knots, |
1635 | const TColStd_Array1OfInteger& Mults, |
1636 | const TColStd_Array1OfReal& AddKnots, |
1637 | const TColStd_Array1OfInteger& AddMults, |
1638 | Standard_Integer& NbPoles, |
1639 | Standard_Integer& NbKnots, |
1640 | const Standard_Real Tolerance, |
1641 | const Standard_Boolean Add) |
1642 | { |
1643 | Standard_Boolean addflat = &AddMults == NULL; |
1644 | |
1645 | Standard_Integer first,last; |
1646 | if (Periodic) { |
1647 | first = Knots.Lower(); |
1648 | last = Knots.Upper(); |
1649 | } |
1650 | else { |
1651 | first = FirstUKnotIndex(Degree,Mults); |
1652 | last = LastUKnotIndex(Degree,Mults); |
1653 | } |
1654 | Standard_Real adeltaK1 = Knots(first)-AddKnots(AddKnots.Lower()); |
1655 | Standard_Real adeltaK2 = AddKnots(AddKnots.Upper())-Knots(last); |
1656 | if (adeltaK1 > Tolerance) return Standard_False; |
1657 | if (adeltaK2 > Tolerance) return Standard_False; |
1658 | |
1659 | Standard_Integer sigma = 0, mult, amult, lastmult = 0; |
1660 | NbKnots = 0; |
1661 | Standard_Integer k = Knots.Lower() - 1; |
1662 | Standard_Integer ak = AddKnots.Lower(); |
1663 | |
1664 | if(Periodic && AddKnots.Length() > 1) |
1665 | { |
1666 | //gka for case when segments was produced on full period only one knot |
1667 | //was added in the end of curve |
1668 | if(fabs(adeltaK1) <= Precision::PConfusion() && |
1669 | fabs(adeltaK2) <= Precision::PConfusion()) |
1670 | ak++; |
1671 | } |
1672 | |
1673 | Standard_Real au,oldau = AddKnots(ak),Eps; |
1674 | |
1675 | while (ak <= AddKnots.Upper()) { |
1676 | au = AddKnots(ak); |
1677 | if (au < oldau) return Standard_False; |
1678 | oldau = au; |
1679 | |
1680 | Eps = Max(Tolerance,Epsilon(au)); |
1681 | |
1682 | while ((k < Knots.Upper()) && (Knots(k+1) - au <= Eps)) { |
1683 | k++; |
1684 | NbKnots++; |
1685 | sigma += Mults(k); |
1686 | } |
1687 | |
1688 | if (addflat) amult = 1; |
1689 | else amult = Max(0,AddMults(ak)); |
1690 | |
1691 | while ((ak < AddKnots.Upper()) && |
1692 | (Abs(au - AddKnots(ak+1)) <= Eps)) { |
1693 | ak++; |
1694 | if (Add) { |
1695 | if (addflat) amult++; |
1696 | else amult += Max(0,AddMults(ak)); |
1697 | } |
1698 | } |
1699 | |
1700 | |
1701 | if (Abs(au - Knots(k)) <= Eps) { |
1702 | // identic to existing knot |
1703 | mult = Mults(k); |
1704 | lastmult = mult;//gka |
1705 | if (Add) { |
1706 | if (mult + amult > Degree) |
1707 | amult = Max(0,Degree - mult); |
1708 | sigma += amult; |
1709 | //lastmult = mult + amult; |
1710 | |
1711 | } |
1712 | else if (amult > mult) { |
1713 | if (amult > Degree) amult = Degree; |
1714 | sigma += amult - mult; |
1715 | //lastmult = amult;//gka modified |
1716 | } |
1717 | /* |
1718 | // on periodic curves if this is the last knot |
1719 | // the multiplicity is added twice to account for the first knot |
1720 | if (k == Knots.Upper() && Periodic) { |
1721 | if (Add) |
1722 | sigma += amult; |
1723 | else |
1724 | sigma += amult - mult; |
1725 | } |
1726 | */ |
1727 | } |
1728 | else { |
1729 | // not identic to existing knot |
1730 | if (amult > 0) { |
1731 | if (amult > Degree) amult = Degree; |
1732 | NbKnots++; |
1733 | //lastmult = amult; |
1734 | sigma += amult; |
1735 | } |
1736 | } |
1737 | |
1738 | ak++; |
1739 | } |
1740 | |
1741 | // count the last knots |
1742 | if (lastmult == 0)// || k < Knots.Upper()) |
1743 | lastmult = Mults(Knots.Upper()); |
1744 | |
1745 | while (k < Knots.Upper()) { |
1746 | k++; |
1747 | NbKnots++; |
1748 | sigma += Mults(k); |
1749 | } |
1750 | |
1751 | if (Periodic) { |
1752 | NbPoles = sigma - lastmult; |
1753 | } |
1754 | else { |
1755 | NbPoles = sigma - Degree - 1; |
1756 | } |
1757 | |
1758 | return Standard_True; |
1759 | } |
1760 | |
1761 | //======================================================================= |
1762 | //function : Copy |
1763 | //purpose : copy reals from an array to an other |
1764 | // |
1765 | // NbValues are copied from OldPoles(OldFirst) |
1766 | // to NewPoles(NewFirst) |
1767 | // |
1768 | // Periodicity is handled. |
1769 | // OldFirst and NewFirst are updated |
1770 | // to the position after the last copied pole. |
1771 | // |
1772 | //======================================================================= |
1773 | |
1774 | static void Copy(const Standard_Integer NbPoles, |
1775 | Standard_Integer& OldFirst, |
1776 | const TColStd_Array1OfReal& OldPoles, |
1777 | Standard_Integer& NewFirst, |
1778 | TColStd_Array1OfReal& NewPoles) |
1779 | { |
1780 | // reset the index in the range for periodicity |
1781 | |
1782 | OldFirst = OldPoles.Lower() + |
1783 | (OldFirst - OldPoles.Lower()) % (OldPoles.Upper() - OldPoles.Lower() + 1); |
1784 | |
1785 | NewFirst = NewPoles.Lower() + |
1786 | (NewFirst - NewPoles.Lower()) % (NewPoles.Upper() - NewPoles.Lower() + 1); |
1787 | |
1788 | // copy |
1789 | Standard_Integer i; |
1790 | |
1791 | for (i = 1; i <= NbPoles; i++) { |
1792 | NewPoles(NewFirst) = OldPoles(OldFirst); |
1793 | OldFirst++; |
1794 | if (OldFirst > OldPoles.Upper()) OldFirst = OldPoles.Lower(); |
1795 | NewFirst++; |
1796 | if (NewFirst > NewPoles.Upper()) NewFirst = NewPoles.Lower(); |
1797 | } |
1798 | } |
1799 | |
1800 | //======================================================================= |
1801 | //function : InsertKnots |
1802 | //purpose : insert an array of knots and multiplicities |
1803 | //======================================================================= |
1804 | |
1805 | void BSplCLib::InsertKnots |
1806 | (const Standard_Integer Degree, |
1807 | const Standard_Boolean Periodic, |
1808 | const Standard_Integer Dimension, |
1809 | const TColStd_Array1OfReal& Poles, |
1810 | const TColStd_Array1OfReal& Knots, |
1811 | const TColStd_Array1OfInteger& Mults, |
1812 | const TColStd_Array1OfReal& AddKnots, |
1813 | const TColStd_Array1OfInteger& AddMults, |
1814 | TColStd_Array1OfReal& NewPoles, |
1815 | TColStd_Array1OfReal& NewKnots, |
1816 | TColStd_Array1OfInteger& NewMults, |
1817 | const Standard_Real Tolerance, |
1818 | const Standard_Boolean Add) |
1819 | { |
1820 | Standard_Boolean addflat = &AddMults == NULL; |
1821 | |
1822 | Standard_Integer i,k,mult,firstmult; |
1823 | Standard_Integer index,kn,curnk,curk; |
1824 | Standard_Integer p,np, curp, curnp, length, depth; |
1825 | Standard_Real u; |
1826 | Standard_Integer need; |
1827 | Standard_Real Eps; |
1828 | |
1829 | // ------------------- |
1830 | // create local arrays |
1831 | // ------------------- |
1832 | |
1833 | Standard_Real *knots = new Standard_Real[2*Degree]; |
1834 | Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension]; |
1835 | |
1836 | //---------------------------- |
1837 | // loop on the knots to insert |
1838 | //---------------------------- |
1839 | |
1840 | curk = Knots.Lower()-1; // current position in Knots |
1841 | curnk = NewKnots.Lower()-1; // current position in NewKnots |
1842 | curp = Poles.Lower(); // current position in Poles |
1843 | curnp = NewPoles.Lower(); // current position in NewPoles |
1844 | |
1845 | // NewKnots, NewMults, NewPoles contains the current state of the curve |
1846 | |
1847 | // index is the first pole of the current curve for insertion schema |
1848 | |
1849 | if (Periodic) index = -Mults(Mults.Lower()); |
1850 | else index = -Degree-1; |
1851 | |
1852 | // on Periodic curves the first knot and the last knot are inserted later |
1853 | // (they are the same knot) |
1854 | firstmult = 0; // multiplicity of the first-last knot for periodic |
1855 | |
1856 | |
1857 | // kn current knot to insert in AddKnots |
1858 | |
1859 | for (kn = AddKnots.Lower(); kn <= AddKnots.Upper(); kn++) { |
1860 | |
1861 | u = AddKnots(kn); |
1862 | Eps = Max(Tolerance,Epsilon(u)); |
1863 | |
1864 | //----------------------------------- |
1865 | // find the position in the old knots |
1866 | // and copy to the new knots |
1867 | //----------------------------------- |
1868 | |
1869 | while (curk < Knots.Upper() && Knots(curk+1) - u <= Eps) { |
1870 | curk++; curnk++; |
1871 | NewKnots(curnk) = Knots(curk); |
1872 | index += NewMults(curnk) = Mults(curk); |
1873 | } |
1874 | |
1875 | //----------------------------------- |
1876 | // Slice the knots and the mults |
1877 | // to the current size of the new curve |
1878 | //----------------------------------- |
1879 | |
1880 | i = curnk + Knots.Upper() - curk; |
1881 | TColStd_Array1OfReal nknots(NewKnots(NewKnots.Lower()),NewKnots.Lower(),i); |
1882 | TColStd_Array1OfInteger nmults(NewMults(NewMults.Lower()),NewMults.Lower(),i); |
1883 | |
1884 | //----------------------------------- |
1885 | // copy enough knots |
1886 | // to compute the insertion schema |
1887 | //----------------------------------- |
1888 | |
1889 | k = curk; |
1890 | i = curnk; |
1891 | mult = 0; |
1892 | |
1893 | while (mult < Degree && k < Knots.Upper()) { |
1894 | k++; i++; |
1895 | nknots(i) = Knots(k); |
1896 | mult += nmults(i) = Mults(k); |
1897 | } |
1898 | |
1899 | // copy knots at the end for periodic curve |
1900 | if (Periodic) { |
1901 | mult = 0; |
1902 | k = Knots.Upper(); |
1903 | i = nknots.Upper(); |
1904 | |
1905 | while (mult < Degree && i > curnk) { |
1906 | nknots(i) = Knots(k); |
1907 | mult += nmults(i) = Mults(k); |
1908 | k--; |
1909 | i--; |
1910 | } |
1911 | nmults(nmults.Upper()) = nmults(nmults.Lower()); |
1912 | } |
1913 | |
1914 | |
1915 | |
1916 | //------------------------------------ |
1917 | // do the boor scheme on the new curve |
1918 | // to insert the new knot |
1919 | //------------------------------------ |
1920 | |
1921 | Standard_Boolean sameknot = (Abs(u-NewKnots(curnk)) <= Eps); |
1922 | |
1923 | if (sameknot) length = Max(0,Degree - NewMults(curnk)); |
1924 | else length = Degree; |
1925 | |
1926 | if (addflat) depth = 1; |
1927 | else depth = Min(Degree,AddMults(kn)); |
1928 | |
1929 | if (sameknot) { |
1930 | if (Add) { |
1931 | if ((NewMults(curnk) + depth) > Degree) |
1932 | depth = Degree - NewMults(curnk); |
1933 | } |
1934 | else { |
1935 | depth = Max(0,depth-NewMults(curnk)); |
1936 | } |
1937 | |
1938 | if (Periodic) { |
1939 | // on periodic curve the first and last knot are delayed to the end |
1940 | if (curk == Knots.Lower() || (curk == Knots.Upper())) { |
1941 | firstmult += depth; |
1942 | depth = 0; |
1943 | } |
1944 | } |
1945 | } |
1946 | if (depth <= 0) continue; |
1947 | |
1948 | BuildKnots(Degree,curnk,Periodic,nknots,nmults,*knots); |
1949 | |
1950 | // copy the poles |
1951 | |
1952 | need = NewPoles.Lower()+(index+length+1)*Dimension - curnp; |
1953 | need = Min(need,Poles.Upper() - curp + 1); |
1954 | |
1955 | p = curp; |
1956 | np = curnp; |
1957 | Copy(need,p,Poles,np,NewPoles); |
1958 | curp += need; |
1959 | curnp += need; |
1960 | |
1961 | // slice the poles to the current number of poles in case of periodic |
1962 | TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1); |
1963 | |
1964 | BuildBoor(index,length,Dimension,npoles,*poles); |
1965 | BoorScheme(u,Degree,*knots,Dimension,*poles,depth,length); |
1966 | |
1967 | //------------------- |
1968 | // copy the new poles |
1969 | //------------------- |
1970 | |
1971 | curnp += depth * Dimension; // number of poles is increased by depth |
1972 | TColStd_Array1OfReal ThePoles(NewPoles(NewPoles.Lower()),NewPoles.Lower(),curnp-1); |
1973 | np = NewKnots.Lower()+(index+1)*Dimension; |
1974 | |
1975 | for (i = 1; i <= length + depth; i++) |
1976 | GetPole(i,length,depth,Dimension,*poles,np,ThePoles); |
1977 | |
1978 | //------------------- |
1979 | // insert the knot |
1980 | //------------------- |
1981 | |
1982 | index += depth; |
1983 | if (sameknot) { |
1984 | NewMults(curnk) += depth; |
1985 | } |
1986 | else { |
1987 | curnk++; |
1988 | NewKnots(curnk) = u; |
1989 | NewMults(curnk) = depth; |
1990 | } |
1991 | } |
1992 | |
1993 | //------------------------------ |
1994 | // copy the last poles and knots |
1995 | //------------------------------ |
1996 | |
1997 | Copy(Poles.Upper() - curp + 1,curp,Poles,curnp,NewPoles); |
1998 | |
1999 | while (curk < Knots.Upper()) { |
2000 | curk++; curnk++; |
2001 | NewKnots(curnk) = Knots(curk); |
2002 | NewMults(curnk) = Mults(curk); |
2003 | } |
2004 | |
2005 | //------------------------------ |
2006 | // process the first-last knot |
2007 | // on periodic curves |
2008 | //------------------------------ |
2009 | |
2010 | if (firstmult > 0) { |
2011 | curnk = NewKnots.Lower(); |
2012 | if (NewMults(curnk) + firstmult > Degree) { |
2013 | firstmult = Degree - NewMults(curnk); |
2014 | } |
2015 | if (firstmult > 0) { |
2016 | |
2017 | length = Degree - NewMults(curnk); |
2018 | depth = firstmult; |
2019 | |
2020 | BuildKnots(Degree,curnk,Periodic,NewKnots,NewMults,*knots); |
2021 | TColStd_Array1OfReal npoles(NewPoles(NewPoles.Lower()), |
2022 | NewPoles.Lower(), |
2023 | NewPoles.Upper()-depth*Dimension); |
2024 | BuildBoor(0,length,Dimension,npoles,*poles); |
2025 | BoorScheme(NewKnots(curnk),Degree,*knots,Dimension,*poles,depth,length); |
2026 | |
2027 | //--------------------------- |
2028 | // copy the new poles |
2029 | // but rotate them with depth |
2030 | //--------------------------- |
2031 | |
2032 | np = NewPoles.Lower(); |
2033 | |
2034 | for (i = depth; i < length + depth; i++) |
2035 | GetPole(i,length,depth,Dimension,*poles,np,NewPoles); |
2036 | |
2037 | np = NewPoles.Upper() - depth*Dimension + 1; |
2038 | |
2039 | for (i = 0; i < depth; i++) |
2040 | GetPole(i,length,depth,Dimension,*poles,np,NewPoles); |
2041 | |
2042 | NewMults(NewMults.Lower()) += depth; |
2043 | NewMults(NewMults.Upper()) += depth; |
2044 | } |
2045 | } |
2046 | // free local arrays |
2047 | delete [] knots; |
2048 | delete [] poles; |
2049 | } |
2050 | |
2051 | //======================================================================= |
2052 | //function : RemoveKnot |
2053 | //purpose : |
2054 | //======================================================================= |
2055 | |
2056 | Standard_Boolean BSplCLib::RemoveKnot |
2057 | (const Standard_Integer Index, |
2058 | const Standard_Integer Mult, |
2059 | const Standard_Integer Degree, |
2060 | const Standard_Boolean Periodic, |
2061 | const Standard_Integer Dimension, |
2062 | const TColStd_Array1OfReal& Poles, |
2063 | const TColStd_Array1OfReal& Knots, |
2064 | const TColStd_Array1OfInteger& Mults, |
2065 | TColStd_Array1OfReal& NewPoles, |
2066 | TColStd_Array1OfReal& NewKnots, |
2067 | TColStd_Array1OfInteger& NewMults, |
2068 | const Standard_Real Tolerance) |
2069 | { |
2070 | Standard_Integer index,i,j,k,p,np; |
2071 | |
2072 | Standard_Integer TheIndex = Index; |
2073 | |
2074 | // protection |
2075 | Standard_Integer first,last; |
2076 | if (Periodic) { |
2077 | first = Knots.Lower(); |
2078 | last = Knots.Upper(); |
2079 | } |
2080 | else { |
2081 | first = BSplCLib::FirstUKnotIndex(Degree,Mults) + 1; |
2082 | last = BSplCLib::LastUKnotIndex(Degree,Mults) - 1; |
2083 | } |
2084 | if (Index < first) return Standard_False; |
2085 | if (Index > last) return Standard_False; |
2086 | |
2087 | if ( Periodic && (Index == first)) TheIndex = last; |
2088 | |
2089 | Standard_Integer depth = Mults(TheIndex) - Mult; |
2090 | Standard_Integer length = Degree - Mult; |
2091 | |
2092 | // ------------------- |
2093 | // create local arrays |
2094 | // ------------------- |
2095 | |
2096 | Standard_Real *knots = new Standard_Real[4*Degree]; |
2097 | Standard_Real *poles = new Standard_Real[(2*Degree+1)*Dimension]; |
2098 | |
2099 | |
2100 | // ------------------------------------ |
2101 | // build the knots for anti Boor Scheme |
2102 | // ------------------------------------ |
2103 | |
2104 | // the new sequence of knots |
2105 | // is obtained from the knots at Index-1 and Index |
2106 | |
2107 | BSplCLib::BuildKnots(Degree,TheIndex-1,Periodic,Knots,Mults,*knots); |
2108 | index = PoleIndex(Degree,TheIndex-1,Periodic,Mults); |
2109 | BSplCLib::BuildKnots(Degree,TheIndex,Periodic,Knots,Mults,knots[2*Degree]); |
2110 | |
2111 | index += Mult; |
2112 | |
2113 | for (i = 0; i < Degree - Mult; i++) |
2114 | knots[i] = knots[i+Mult]; |
2115 | |
2116 | for (i = Degree-Mult; i < 2*Degree; i++) |
2117 | knots[i] = knots[2*Degree+i]; |
2118 | |
2119 | |
2120 | // ------------------------------------ |
2121 | // build the poles for anti Boor Scheme |
2122 | // ------------------------------------ |
2123 | |
2124 | p = Poles.Lower()+index * Dimension; |
2125 | |
2126 | for (i = 0; i <= length + depth; i++) { |
2127 | j = Dimension * BoorIndex(i,length,depth); |
2128 | |
2129 | for (k = 0; k < Dimension; k++) { |
2130 | poles[j+k] = Poles(p+k); |
2131 | } |
2132 | p += Dimension; |
2133 | if (p > Poles.Upper()) p = Poles.Lower(); |
2134 | } |
2135 | |
2136 | |
2137 | // ---------------- |
2138 | // Anti Boor Scheme |
2139 | // ---------------- |
2140 | |
2141 | Standard_Boolean result = AntiBoorScheme(Knots(TheIndex),Degree,*knots, |
2142 | Dimension,*poles, |
2143 | depth,length,Tolerance); |
2144 | |
2145 | // ---------------- |
2146 | // copy the results |
2147 | // ---------------- |
2148 | |
2149 | if (result) { |
2150 | |
2151 | // poles |
2152 | |
2153 | p = Poles.Lower(); |
2154 | np = NewPoles.Lower(); |
2155 | |
2156 | // unmodified poles before |
2157 | Copy((index+1)*Dimension,p,Poles,np,NewPoles); |
2158 | |
2159 | |
2160 | // modified |
2161 | |
2162 | for (i = 1; i <= length; i++) |
2163 | BSplCLib::GetPole(i,length,0,Dimension,*poles,np,NewPoles); |
2164 | p += (length + depth) * Dimension ; |
2165 | |
2166 | // unmodified poles after |
2167 | if (p != Poles.Lower()) { |
2168 | i = Poles.Upper() - p + 1; |
2169 | Copy(i,p,Poles,np,NewPoles); |
2170 | } |
2171 | |
2172 | // knots and mults |
2173 | |
2174 | if (Mult > 0) { |
2175 | NewKnots = Knots; |
2176 | NewMults = Mults; |
2177 | NewMults(TheIndex) = Mult; |
2178 | if (Periodic) { |
2179 | if (TheIndex == first) NewMults(last) = Mult; |
2180 | if (TheIndex == last) NewMults(first) = Mult; |
2181 | } |
2182 | } |
2183 | else { |
2184 | if (!Periodic || (TheIndex != first && TheIndex != last)) { |
2185 | |
2186 | for (i = Knots.Lower(); i < TheIndex; i++) { |
2187 | NewKnots(i) = Knots(i); |
2188 | NewMults(i) = Mults(i); |
2189 | } |
2190 | |
2191 | for (i = TheIndex+1; i <= Knots.Upper(); i++) { |
2192 | NewKnots(i-1) = Knots(i); |
2193 | NewMults(i-1) = Mults(i); |
2194 | } |
2195 | } |
2196 | else { |
2197 | // The interesting case of a Periodic curve |
2198 | // where the first and last knot is removed. |
2199 | |
2200 | for (i = first; i < last-1; i++) { |
2201 | NewKnots(i) = Knots(i+1); |
2202 | NewMults(i) = Mults(i+1); |
2203 | } |
2204 | NewKnots(last-1) = NewKnots(first) + Knots(last) - Knots(first); |
2205 | NewMults(last-1) = NewMults(first); |
2206 | } |
2207 | } |
2208 | } |
2209 | |
2210 | |
2211 | // free local arrays |
2212 | delete [] knots; |
2213 | delete [] poles; |
2214 | |
2215 | return result; |
2216 | } |
2217 | |
2218 | //======================================================================= |
2219 | //function : IncreaseDegreeCountKnots |
2220 | //purpose : |
2221 | //======================================================================= |
2222 | |
2223 | Standard_Integer BSplCLib::IncreaseDegreeCountKnots |
2224 | (const Standard_Integer Degree, |
2225 | const Standard_Integer NewDegree, |
2226 | const Standard_Boolean Periodic, |
2227 | const TColStd_Array1OfInteger& Mults) |
2228 | { |
2229 | if (Periodic) return Mults.Length(); |
2230 | Standard_Integer f = FirstUKnotIndex(Degree,Mults); |
2231 | Standard_Integer l = LastUKnotIndex(Degree,Mults); |
2232 | Standard_Integer m,i,removed = 0, step = NewDegree - Degree; |
2233 | |
2234 | i = Mults.Lower(); |
2235 | m = Degree + (f - i + 1) * step + 1; |
2236 | |
2237 | while (m > NewDegree+1) { |
2238 | removed++; |
2239 | m -= Mults(i) + step; |
2240 | i++; |
2241 | } |
2242 | if (m < NewDegree+1) removed--; |
2243 | |
2244 | i = Mults.Upper(); |
2245 | m = Degree + (i - l + 1) * step + 1; |
2246 | |
2247 | while (m > NewDegree+1) { |
2248 | removed++; |
2249 | m -= Mults(i) + step; |
2250 | i--; |
2251 | } |
2252 | if (m < NewDegree+1) removed--; |
2253 | |
2254 | return Mults.Length() - removed; |
2255 | } |
2256 | |
2257 | //======================================================================= |
2258 | //function : IncreaseDegree |
2259 | //purpose : |
2260 | //======================================================================= |
2261 | |
2262 | void BSplCLib::IncreaseDegree |
2263 | (const Standard_Integer Degree, |
2264 | const Standard_Integer NewDegree, |
2265 | const Standard_Boolean Periodic, |
2266 | const Standard_Integer Dimension, |
2267 | const TColStd_Array1OfReal& Poles, |
2268 | const TColStd_Array1OfReal& Knots, |
2269 | const TColStd_Array1OfInteger& Mults, |
2270 | TColStd_Array1OfReal& NewPoles, |
2271 | TColStd_Array1OfReal& NewKnots, |
2272 | TColStd_Array1OfInteger& NewMults) |
2273 | { |
2274 | // Degree elevation of a BSpline Curve |
2275 | |
2276 | // This algorithms loops on degree incrementation from Degree to NewDegree. |
2277 | // The variable curDeg is the current degree to increment. |
2278 | |
2279 | // Before degree incrementations a "working curve" is created. |
2280 | // It has the same knot, poles and multiplicities. |
2281 | |
2282 | // If the curve is periodic knots are added on the working curve before |
2283 | // and after the existing knots to make it a non-periodic curves. |
2284 | // The poles are also copied. |
2285 | |
2286 | // The first and last multiplicity of the working curve are set to Degree+1, |
2287 | // null poles are added if necessary. |
2288 | |
2289 | // Then the degree is incremented on the working curve. |
2290 | // The knots are unchanged but all multiplicities will be incremented. |
2291 | |
2292 | // Each degree incrementation is achieved by averaging curDeg+1 curves. |
2293 | |
2294 | // See : Degree elevation of B-spline curves |
2295 | // Hartmut PRAUTZSCH |
2296 | // CAGD 1 (1984) |
2297 | |
2298 | |
2299 | //------------------------- |
2300 | // create the working curve |
2301 | //------------------------- |
2302 | |
2303 | Standard_Integer i,k,f,l,m,pf,pl,firstknot; |
2304 | |
2305 | pf = 0; // number of null poles added at beginning |
2306 | pl = 0; // number of null poles added at end |
2307 | |
2308 | Standard_Integer nbwknots = Knots.Length(); |
2309 | f = FirstUKnotIndex(Degree,Mults); |
2310 | l = LastUKnotIndex (Degree,Mults); |
2311 | |
2312 | if (Periodic) { |
2313 | // Periodic curves are transformed in non-periodic curves |
2314 | |
2315 | nbwknots += f - Mults.Lower(); |
2316 | |
2317 | pf = -Degree - 1; |
2318 | |
2319 | for (i = Mults.Lower(); i <= f; i++) |
2320 | pf += Mults(i); |
2321 | |
2322 | nbwknots += Mults.Upper() - l; |
2323 | |
2324 | pl = -Degree - 1; |
2325 | |
2326 | for (i = l; i <= Mults.Upper(); i++) |
2327 | pl += Mults(i); |
2328 | } |
2329 | |
2330 | // copy the knots and multiplicities |
2331 | TColStd_Array1OfReal wknots(1,nbwknots); |
2332 | TColStd_Array1OfInteger wmults(1,nbwknots); |
2333 | if (!Periodic) { |
2334 | wknots = Knots; |
2335 | wmults = Mults; |
2336 | } |
2337 | else { |
2338 | // copy the knots for a periodic curve |
2339 | Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower()); |
2340 | i = 0; |
2341 | |
2342 | for (k = l; k < Knots.Upper(); k++) { |
2343 | i++; |
2344 | wknots(i) = Knots(k) - period; |
2345 | wmults(i) = Mults(k); |
2346 | } |
2347 | |
2348 | for (k = Knots.Lower(); k <= Knots.Upper(); k++) { |
2349 | i++; |
2350 | wknots(i) = Knots(k); |
2351 | wmults(i) = Mults(k); |
2352 | } |
2353 | |
2354 | for (k = Knots.Lower()+1; k <= f; k++) { |
2355 | i++; |
2356 | wknots(i) = Knots(k) + period; |
2357 | wmults(i) = Mults(k); |
2358 | } |
2359 | } |
2360 | |
2361 | // set the first and last mults to Degree+1 |
2362 | // and add null poles |
2363 | |
2364 | pf += Degree + 1 - wmults(1); |
2365 | wmults(1) = Degree + 1; |
2366 | pl += Degree + 1 - wmults(nbwknots); |
2367 | wmults(nbwknots) = Degree + 1; |
2368 | |
2369 | //--------------------------- |
2370 | // poles of the working curve |
2371 | //--------------------------- |
2372 | |
2373 | Standard_Integer nbwpoles = 0; |
2374 | |
2375 | for (i = 1; i <= nbwknots; i++) nbwpoles += wmults(i); |
2376 | nbwpoles -= Degree + 1; |
2377 | |
2378 | // we provide space for degree elevation |
2379 | TColStd_Array1OfReal |
2380 | wpoles(1,(nbwpoles + (nbwknots-1) * (NewDegree - Degree)) * Dimension); |
2381 | |
2382 | for (i = 1; i <= pf * Dimension; i++) |
2383 | wpoles(i) = 0; |
2384 | |
2385 | k = Poles.Lower(); |
2386 | |
2387 | for (i = pf * Dimension + 1; i <= (nbwpoles - pl) * Dimension; i++) { |
2388 | wpoles(i) = Poles(k); |
2389 | k++; |
2390 | if (k > Poles.Upper()) k = Poles.Lower(); |
2391 | } |
2392 | |
2393 | for (i = (nbwpoles-pl)*Dimension+1; i <= nbwpoles*Dimension; i++) |
2394 | wpoles(i) = 0; |
2395 | |
2396 | |
2397 | //----------------------------------------------------------- |
2398 | // Declare the temporary arrays used in degree incrementation |
2399 | //----------------------------------------------------------- |
2400 | |
2401 | Standard_Integer nbwp = nbwpoles + (nbwknots-1) * (NewDegree - Degree); |
2402 | // Arrays for storing the temporary curves |
2403 | TColStd_Array1OfReal tempc1(1,nbwp * Dimension); |
2404 | TColStd_Array1OfReal tempc2(1,nbwp * Dimension); |
2405 | |
2406 | // Array for storing the knots to insert |
2407 | TColStd_Array1OfReal iknots(1,nbwknots); |
2408 | |
2409 | // Arrays for receiving the knots after insertion |
2410 | TColStd_Array1OfReal nknots(1,nbwknots); |
2411 | |
2412 | |
2413 | |
2414 | //------------------------------ |
2415 | // Loop on degree incrementation |
2416 | //------------------------------ |
2417 | |
2418 | Standard_Integer step,curDeg; |
2419 | Standard_Integer nbp = nbwpoles; |
2420 | nbwp = nbp; |
2421 | |
2422 | for (curDeg = Degree; curDeg < NewDegree; curDeg++) { |
2423 | |
2424 | nbp = nbwp; // current number of poles |
2425 | nbwp = nbp + nbwknots - 1; // new number of poles |
2426 | |
2427 | // For the averaging |
2428 | TColStd_Array1OfReal nwpoles(1,nbwp * Dimension); |
2429 | nwpoles.Init(0.0e0) ; |
2430 | |
2431 | |
2432 | for (step = 0; step <= curDeg; step++) { |
2433 | |
2434 | // Compute the bspline of rank step. |
2435 | |
2436 | // if not the first time, decrement the multiplicities back |
2437 | if (step != 0) { |
2438 | for (i = 1; i <= nbwknots; i++) |
2439 | wmults(i)--; |
2440 | } |
2441 | |
2442 | // Poles are the current poles |
2443 | // but the poles congruent to step are duplicated. |
2444 | |
2445 | Standard_Integer offset = 0; |
2446 | |
2447 | for (i = 0; i < nbp; i++) { |
2448 | offset++; |
2449 | |
2450 | for (k = 0; k < Dimension; k++) { |
2451 | tempc1((offset-1)*Dimension+k+1) = |
2452 | wpoles(NewPoles.Lower()+i*Dimension+k); |
2453 | } |
2454 | if (i % (curDeg+1) == step) { |
2455 | offset++; |
2456 | |
2457 | for (k = 0; k < Dimension; k++) { |
2458 | tempc1((offset-1)*Dimension+k+1) = |
2459 | wpoles(NewPoles.Lower()+i*Dimension+k); |
2460 | } |
2461 | } |
2462 | } |
2463 | |
2464 | // Knots multiplicities are increased |
2465 | // For each knot where the sum of multiplicities is congruent to step |
2466 | |
2467 | Standard_Integer stepmult = step+1; |
2468 | Standard_Integer nbknots = 0; |
2469 | Standard_Integer smult = 0; |
2470 | |
2471 | for (k = 1; k <= nbwknots; k++) { |
2472 | smult += wmults(k); |
2473 | if (smult >= stepmult) { |
2474 | // this knot is increased |
2475 | stepmult += curDeg+1; |
2476 | wmults(k)++; |
2477 | } |
2478 | else { |
2479 | // this knot is inserted |
2480 | nbknots++; |
2481 | iknots(nbknots) = wknots(k); |
2482 | } |
2483 | } |
2484 | |
2485 | // the curve is obtained by inserting the knots |
2486 | // to raise the multiplicities |
2487 | |
2488 | // we build "slices" of the arrays to set the correct size |
2489 | if (nbknots > 0) { |
2490 | TColStd_Array1OfReal aknots(iknots(1),1,nbknots); |
2491 | TColStd_Array1OfReal curve (tempc1(1),1,offset * Dimension); |
2492 | TColStd_Array1OfReal ncurve(tempc2(1),1,nbwp * Dimension); |
2493 | // InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults, |
2494 | // aknots,NoMults(),ncurve,nknots,wmults,Epsilon(1.)); |
2495 | |
2496 | InsertKnots(curDeg+1,Standard_False,Dimension,curve,wknots,wmults, |
2497 | aknots,NoMults(),ncurve,nknots,wmults,0.0); |
2498 | |
2499 | // add to the average |
2500 | |
2501 | for (i = 1; i <= nbwp * Dimension; i++) |
2502 | nwpoles(i) += ncurve(i); |
2503 | } |
2504 | else { |
2505 | // add to the average |
2506 | |
2507 | for (i = 1; i <= nbwp * Dimension; i++) |
2508 | nwpoles(i) += tempc1(i); |
2509 | } |
2510 | } |
2511 | |
2512 | // The result is the average |
2513 | |
2514 | for (i = 1; i <= nbwp * Dimension; i++) { |
2515 | wpoles(i) = nwpoles(i) / (curDeg+1); |
2516 | } |
2517 | } |
2518 | |
2519 | //----------------- |
2520 | // Copy the results |
2521 | //----------------- |
2522 | |
2523 | // index in new knots of the first knot of the curve |
2524 | if (Periodic) |
2525 | firstknot = Mults.Upper() - l + 1; |
2526 | else |
2527 | firstknot = f; |
2528 | |
2529 | // the new curve starts at index firstknot |
2530 | // so we must remove knots until the sum of multiplicities |
2531 | // from the first to the start is NewDegree+1 |
2532 | |
2533 | // m is the current sum of multiplicities |
2534 | m = 0; |
2535 | |
2536 | for (k = 1; k <= firstknot; k++) |
2537 | m += wmults(k); |
2538 | |
2539 | // compute the new first knot (k), pf will be the index of the first pole |
2540 | k = 1; |
2541 | pf = 0; |
2542 | |
2543 | while (m > NewDegree+1) { |
2544 | k++; |
2545 | m -= wmults(k); |
2546 | pf += wmults(k); |
2547 | } |
2548 | if (m < NewDegree+1) { |
2549 | k--; |
2550 | wmults(k) += m - NewDegree - 1; |
2551 | pf += m - NewDegree - 1; |
2552 | } |
2553 | |
2554 | // on a periodic curve the knots start with firstknot |
2555 | if (Periodic) |
2556 | k = firstknot; |
2557 | |
2558 | // copy knots |
2559 | |
2560 | for (i = NewKnots.Lower(); i <= NewKnots.Upper(); i++) { |
2561 | NewKnots(i) = wknots(k); |
2562 | NewMults(i) = wmults(k); |
2563 | k++; |
2564 | } |
2565 | |
2566 | // copy poles |
2567 | pf *= Dimension; |
2568 | |
2569 | for (i = NewPoles.Lower(); i <= NewPoles.Upper(); i++) { |
2570 | pf++; |
2571 | NewPoles(i) = wpoles(pf); |
2572 | } |
2573 | } |
2574 | |
2575 | //======================================================================= |
2576 | //function : PrepareUnperiodize |
2577 | //purpose : |
2578 | //======================================================================= |
2579 | |
2580 | void BSplCLib::PrepareUnperiodize |
2581 | (const Standard_Integer Degree, |
2582 | const TColStd_Array1OfInteger& Mults, |
2583 | Standard_Integer& NbKnots, |
2584 | Standard_Integer& NbPoles) |
2585 | { |
2586 | Standard_Integer i; |
2587 | // initialize NbKnots and NbPoles |
2588 | NbKnots = Mults.Length(); |
2589 | NbPoles = - Degree - 1; |
2590 | |
2591 | for (i = Mults.Lower(); i <= Mults.Upper(); i++) |
2592 | NbPoles += Mults(i); |
2593 | |
2594 | Standard_Integer sigma, k; |
2595 | // Add knots at the beginning of the curve to raise Multiplicities |
2596 | // to Degre + 1; |
2597 | sigma = Mults(Mults.Lower()); |
2598 | k = Mults.Upper() - 1; |
2599 | |
2600 | while ( sigma < Degree + 1) { |
2601 | sigma += Mults(k); |
2602 | NbPoles += Mults(k); |
2603 | k--; |
2604 | NbKnots++; |
2605 | } |
2606 | // We must add exactly until Degree + 1 -> |
2607 | // Supress the excedent. |
2608 | if ( sigma > Degree + 1) |
2609 | NbPoles -= sigma - Degree - 1; |
2610 | |
2611 | // Add knots at the end of the curve to raise Multiplicities |
2612 | // to Degre + 1; |
2613 | sigma = Mults(Mults.Upper()); |
2614 | k = Mults.Lower() + 1; |
2615 | |
2616 | while ( sigma < Degree + 1) { |
2617 | sigma += Mults(k); |
2618 | NbPoles += Mults(k); |
2619 | k++; |
2620 | NbKnots++; |
2621 | } |
2622 | // We must add exactly until Degree + 1 -> |
2623 | // Supress the excedent. |
2624 | if ( sigma > Degree + 1) |
2625 | NbPoles -= sigma - Degree - 1; |
2626 | } |
2627 | |
2628 | //======================================================================= |
2629 | //function : Unperiodize |
2630 | //purpose : |
2631 | //======================================================================= |
2632 | |
2633 | void BSplCLib::Unperiodize |
2634 | (const Standard_Integer Degree, |
2635 | const Standard_Integer , // Dimension, |
2636 | const TColStd_Array1OfInteger& Mults, |
2637 | const TColStd_Array1OfReal& Knots, |
2638 | const TColStd_Array1OfReal& Poles, |
2639 | TColStd_Array1OfInteger& NewMults, |
2640 | TColStd_Array1OfReal& NewKnots, |
2641 | TColStd_Array1OfReal& NewPoles) |
2642 | { |
2643 | Standard_Integer sigma, k, index = 0; |
2644 | // evaluation of index : number of knots to insert before knot(1) to |
2645 | // raise sum of multiplicities to <Degree + 1> |
2646 | sigma = Mults(Mults.Lower()); |
2647 | k = Mults.Upper() - 1; |
2648 | |
2649 | while ( sigma < Degree + 1) { |
2650 | sigma += Mults(k); |
2651 | k--; |
2652 | index++; |
2653 | } |
2654 | |
2655 | Standard_Real period = Knots(Knots.Upper()) - Knots(Knots.Lower()); |
2656 | |
2657 | // set the 'interior' knots; |
2658 | |
2659 | for ( k = 1; k <= Knots.Length(); k++) { |
2660 | NewKnots ( k + index ) = Knots( k); |
2661 | NewMults ( k + index ) = Mults( k); |
2662 | } |
2663 | |
2664 | // set the 'starting' knots; |
2665 | |
2666 | for ( k = 1; k <= index; k++) { |
2667 | NewKnots( k) = NewKnots( k + Knots.Length() - 1) - period; |
2668 | NewMults( k) = NewMults( k + Knots.Length() - 1); |
2669 | } |
2670 | NewMults( 1) -= sigma - Degree -1; |
2671 | |
2672 | // set the 'ending' knots; |
2673 | sigma = NewMults( index + Knots.Length() ); |
2674 | |
2675 | for ( k = Knots.Length() + index + 1; k <= NewKnots.Length(); k++) { |
2676 | NewKnots( k) = NewKnots( k - Knots.Length() + 1) + period; |
2677 | NewMults( k) = NewMults( k - Knots.Length() + 1); |
2678 | sigma += NewMults( k - Knots.Length() + 1); |
2679 | } |
2680 | NewMults(NewMults.Length()) -= sigma - Degree - 1; |
2681 | |
2682 | for ( k = 1 ; k <= NewPoles.Length(); k++) { |
2683 | NewPoles(k ) = Poles( (k-1) % Poles.Length() + 1); |
2684 | } |
2685 | } |
2686 | |
2687 | //======================================================================= |
2688 | //function : PrepareTrimming |
2689 | //purpose : |
2690 | //======================================================================= |
2691 | |
2692 | void BSplCLib::PrepareTrimming(const Standard_Integer Degree, |
2693 | const Standard_Boolean Periodic, |
2694 | const TColStd_Array1OfReal& Knots, |
2695 | const TColStd_Array1OfInteger& Mults, |
2696 | const Standard_Real U1, |
2697 | const Standard_Real U2, |
2698 | Standard_Integer& NbKnots, |
2699 | Standard_Integer& NbPoles) |
2700 | { |
2701 | Standard_Integer i; |
2702 | Standard_Real NewU1, NewU2; |
2703 | Standard_Integer index1 = 0, index2 = 0; |
2704 | |
2705 | // Eval index1, index2 : position of U1 and U2 in the Array Knots |
2706 | // such as Knots(index1-1) <= U1 < Knots(index1) |
2707 | // Knots(index2-1) <= U2 < Knots(index2) |
2708 | LocateParameter( Degree, Knots, Mults, U1, Periodic, |
2709 | Knots.Lower(), Knots.Upper(), index1, NewU1); |
2710 | LocateParameter( Degree, Knots, Mults, U2, Periodic, |
2711 | Knots.Lower(), Knots.Upper(), index2, NewU2); |
2712 | index1++; |
2713 | if ( Abs(Knots(index2) - U2) <= Epsilon( U1)) |
2714 | index2--; |
2715 | |
2716 | // eval NbKnots: |
2717 | NbKnots = index2 - index1 + 3; |
2718 | |
2719 | // eval NbPoles: |
2720 | NbPoles = Degree + 1; |
2721 | |
2722 | for ( i = index1; i <= index2; i++) |
2723 | NbPoles += Mults(i); |
2724 | } |
2725 | |
2726 | //======================================================================= |
2727 | //function : Trimming |
2728 | //purpose : |
2729 | //======================================================================= |
2730 | |
2731 | void BSplCLib::Trimming(const Standard_Integer Degree, |
2732 | const Standard_Boolean Periodic, |
2733 | const Standard_Integer Dimension, |
2734 | const TColStd_Array1OfReal& Knots, |
2735 | const TColStd_Array1OfInteger& Mults, |
2736 | const TColStd_Array1OfReal& Poles, |
2737 | const Standard_Real U1, |
2738 | const Standard_Real U2, |
2739 | TColStd_Array1OfReal& NewKnots, |
2740 | TColStd_Array1OfInteger& NewMults, |
2741 | TColStd_Array1OfReal& NewPoles) |
2742 | { |
2743 | Standard_Integer i, nbpoles, nbknots; |
2744 | Standard_Real kk[2]; |
2745 | Standard_Integer mm[2]; |
2746 | TColStd_Array1OfReal K( kk[0], 1, 2 ); |
2747 | TColStd_Array1OfInteger M( mm[0], 1, 2 ); |
2748 | |
2749 | K(1) = U1; K(2) = U2; |
2750 | mm[0] = mm[1] = Degree; |
2751 | if (!PrepareInsertKnots( Degree, Periodic, Knots, Mults, K, M, |
2752 | nbpoles, nbknots, Epsilon( U1), 0)) |
2753 | Standard_OutOfRange::Raise(); |
2754 | |
2755 | TColStd_Array1OfReal TempPoles(1, nbpoles*Dimension); |
2756 | TColStd_Array1OfReal TempKnots(1, nbknots); |
2757 | TColStd_Array1OfInteger TempMults(1, nbknots); |
2758 | |
2759 | // |
2760 | // do not allow the multiplicities to Add : they must be less than Degree |
2761 | // |
2762 | InsertKnots(Degree, Periodic, Dimension, Poles, Knots, Mults, |
2763 | K, M, TempPoles, TempKnots, TempMults, Epsilon(U1), |
2764 | Standard_False); |
2765 | |
2766 | // find in TempPoles the index of the pole corresponding to U1 |
2767 | Standard_Integer Kindex = 0, Pindex; |
2768 | Standard_Real NewU1; |
2769 | LocateParameter( Degree, TempKnots, TempMults, U1, Periodic, |
2770 | TempKnots.Lower(), TempKnots.Upper(), Kindex, NewU1); |
2771 | Pindex = PoleIndex ( Degree, Kindex, Periodic, TempMults); |
2772 | Pindex *= Dimension; |
2773 | |
2774 | for ( i = 1; i <= NewPoles.Length(); i++) NewPoles(i) = TempPoles(Pindex + i); |
2775 | |
2776 | for ( i = 1; i <= NewKnots.Length(); i++) { |
2777 | NewKnots(i) = TempKnots( Kindex+i-1); |
2778 | NewMults(i) = TempMults( Kindex+i-1); |
2779 | } |
2780 | NewMults(1) = Min(Degree, NewMults(1)) + 1 ; |
2781 | NewMults(NewMults.Length())= Min(Degree, NewMults(NewMults.Length())) + 1 ; |
2782 | } |
2783 | |
2784 | //======================================================================= |
2785 | //function : Solves a LU factored Matrix |
2786 | //purpose : |
2787 | //======================================================================= |
2788 | |
2789 | Standard_Integer |
2790 | BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, |
2791 | const Standard_Integer UpperBandWidth, |
2792 | const Standard_Integer LowerBandWidth, |
2793 | const Standard_Integer ArrayDimension, |
2794 | Standard_Real& Array) |
2795 | { |
2796 | Standard_Integer ii, |
2797 | jj, |
2798 | kk, |
2799 | MinIndex, |
2800 | MaxIndex, |
2801 | ReturnCode = 0 ; |
2802 | |
2803 | Standard_Real *PolesArray = &Array ; |
2804 | Standard_Real Inverse ; |
2805 | |
2806 | |
2807 | if (Matrix.LowerCol() != 1 || |
2808 | Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) { |
2809 | ReturnCode = 1 ; |
2810 | goto FINISH ; |
2811 | } |
2812 | |
2813 | for (ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow() ; ii++) { |
2814 | MinIndex = (ii - LowerBandWidth >= Matrix.LowerRow() ? |
2815 | ii - LowerBandWidth : Matrix.LowerRow()) ; |
2816 | |
2817 | for ( jj = MinIndex ; jj < ii ; jj++) { |
2818 | |
2819 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
2820 | PolesArray[(ii-1) * ArrayDimension + kk] += |
2821 | PolesArray[(jj-1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1) ; |
2822 | } |
2823 | } |
2824 | } |
2825 | |
2826 | for (ii = Matrix.UpperRow() ; ii >= Matrix.LowerRow() ; ii--) { |
2827 | MaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ? |
2828 | ii + UpperBandWidth : Matrix.UpperRow()) ; |
2829 | |
2830 | for (jj = MaxIndex ; jj > ii ; jj--) { |
2831 | |
2832 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
2833 | PolesArray[(ii-1) * ArrayDimension + kk] -= |
2834 | PolesArray[(jj - 1) * ArrayDimension + kk] * |
2835 | Matrix(ii, jj - ii + LowerBandWidth + 1) ; |
2836 | } |
2837 | } |
2838 | |
2839 | //fixing a bug PRO18577 to avoid divizion by zero |
2840 | |
2841 | Standard_Real divizor = Matrix(ii,LowerBandWidth + 1) ; |
2842 | Standard_Real Toler = 1.0e-16; |
2843 | if ( Abs(divizor) > Toler ) |
2844 | Inverse = 1.0e0 / divizor ; |
2845 | else { |
2846 | Inverse = 1.0e0; |
2847 | // cout << " BSplCLib::SolveBandedSystem() : zero determinant " << endl; |
2848 | ReturnCode = 1; |
2849 | goto FINISH; |
2850 | } |
2851 | |
2852 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
2853 | PolesArray[(ii-1) * ArrayDimension + kk] *= Inverse ; |
2854 | } |
2855 | } |
2856 | FINISH : |
2857 | return (ReturnCode) ; |
2858 | } |
2859 | |
2860 | //======================================================================= |
2861 | //function : Solves a LU factored Matrix |
2862 | //purpose : |
2863 | //======================================================================= |
2864 | |
2865 | Standard_Integer |
2866 | BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, |
2867 | const Standard_Integer UpperBandWidth, |
2868 | const Standard_Integer LowerBandWidth, |
2869 | const Standard_Boolean HomogeneousFlag, |
2870 | const Standard_Integer ArrayDimension, |
2871 | Standard_Real& Poles, |
2872 | Standard_Real& Weights) |
2873 | { |
2874 | Standard_Integer ii, |
2875 | kk, |
2876 | ErrorCode = 0, |
2877 | ReturnCode = 0 ; |
2878 | |
2879 | Standard_Real Inverse, |
2880 | *PolesArray = &Poles, |
2881 | *WeightsArray = &Weights ; |
2882 | |
2883 | if (Matrix.LowerCol() != 1 || |
2884 | Matrix.UpperCol() != UpperBandWidth + LowerBandWidth + 1) { |
2885 | ReturnCode = 1 ; |
2886 | goto FINISH ; |
2887 | } |
2888 | if (HomogeneousFlag == Standard_False) { |
2889 | |
2890 | for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1; ii++) { |
2891 | |
2892 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
2893 | PolesArray[ii * ArrayDimension + kk] *= |
2894 | WeightsArray[ii] ; |
2895 | } |
2896 | } |
2897 | } |
2898 | ErrorCode = |
2899 | BSplCLib::SolveBandedSystem(Matrix, |
2900 | UpperBandWidth, |
2901 | LowerBandWidth, |
2902 | ArrayDimension, |
2903 | Poles) ; |
2904 | if (ErrorCode != 0) { |
2905 | ReturnCode = 2 ; |
2906 | goto FINISH ; |
2907 | } |
2908 | ErrorCode = |
2909 | BSplCLib::SolveBandedSystem(Matrix, |
2910 | UpperBandWidth, |
2911 | LowerBandWidth, |
2912 | 1, |
2913 | Weights) ; |
2914 | if (ErrorCode != 0) { |
2915 | ReturnCode = 3 ; |
2916 | goto FINISH ; |
2917 | } |
2918 | if (HomogeneousFlag == Standard_False) { |
2919 | |
2920 | for (ii = 0 ; ii < Matrix.UpperRow() - Matrix.LowerRow() + 1 ; ii++) { |
2921 | Inverse = 1.0e0 / WeightsArray[ii] ; |
2922 | |
2923 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
2924 | PolesArray[ii * ArrayDimension + kk] *= Inverse ; |
2925 | } |
2926 | } |
2927 | } |
2928 | FINISH : return (ReturnCode) ; |
2929 | } |
2930 | |
2931 | //======================================================================= |
2932 | //function : BuildSchoenbergPoints |
2933 | //purpose : |
2934 | //======================================================================= |
2935 | |
2936 | void BSplCLib::BuildSchoenbergPoints(const Standard_Integer Degree, |
2937 | const TColStd_Array1OfReal& FlatKnots, |
2938 | TColStd_Array1OfReal& Parameters) |
2939 | { |
2940 | Standard_Integer ii, |
2941 | jj ; |
2942 | Standard_Real Inverse ; |
2943 | Inverse = 1.0e0 / (Standard_Real)Degree ; |
2944 | |
2945 | for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) { |
2946 | Parameters(ii) = 0.0e0 ; |
2947 | |
2948 | for (jj = 1 ; jj <= Degree ; jj++) { |
2949 | Parameters(ii) += FlatKnots(jj + ii) ; |
2950 | } |
2951 | Parameters(ii) *= Inverse ; |
2952 | } |
2953 | } |
2954 | |
2955 | //======================================================================= |
2956 | //function : Interpolate |
2957 | //purpose : |
2958 | //======================================================================= |
2959 | |
2960 | void BSplCLib::Interpolate(const Standard_Integer Degree, |
2961 | const TColStd_Array1OfReal& FlatKnots, |
2962 | const TColStd_Array1OfReal& Parameters, |
2963 | const TColStd_Array1OfInteger& ContactOrderArray, |
2964 | const Standard_Integer ArrayDimension, |
2965 | Standard_Real& Poles, |
2966 | Standard_Integer& InversionProblem) |
2967 | { |
2968 | Standard_Integer ErrorCode, |
2969 | UpperBandWidth, |
2970 | LowerBandWidth ; |
2971 | // Standard_Real *PolesArray = &Poles ; |
2972 | math_Matrix InterpolationMatrix(1, Parameters.Length(), |
2973 | 1, 2 * Degree + 1) ; |
2974 | ErrorCode = |
2975 | BSplCLib::BuildBSpMatrix(Parameters, |
2976 | ContactOrderArray, |
2977 | FlatKnots, |
2978 | Degree, |
2979 | InterpolationMatrix, |
2980 | UpperBandWidth, |
2981 | LowerBandWidth) ; |
2982 | Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ; |
2983 | |
2984 | ErrorCode = |
2985 | BSplCLib::FactorBandedMatrix(InterpolationMatrix, |
2986 | UpperBandWidth, |
2987 | LowerBandWidth, |
2988 | InversionProblem) ; |
2989 | Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ; |
2990 | |
2991 | ErrorCode = |
2992 | BSplCLib::SolveBandedSystem(InterpolationMatrix, |
2993 | UpperBandWidth, |
2994 | LowerBandWidth, |
2995 | ArrayDimension, |
2996 | Poles) ; |
2997 | |
2998 | Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ; |
2999 | } |
3000 | |
3001 | //======================================================================= |
3002 | //function : Interpolate |
3003 | //purpose : |
3004 | //======================================================================= |
3005 | |
3006 | void BSplCLib::Interpolate(const Standard_Integer Degree, |
3007 | const TColStd_Array1OfReal& FlatKnots, |
3008 | const TColStd_Array1OfReal& Parameters, |
3009 | const TColStd_Array1OfInteger& ContactOrderArray, |
3010 | const Standard_Integer ArrayDimension, |
3011 | Standard_Real& Poles, |
3012 | Standard_Real& Weights, |
3013 | Standard_Integer& InversionProblem) |
3014 | { |
3015 | Standard_Integer ErrorCode, |
3016 | UpperBandWidth, |
3017 | LowerBandWidth ; |
3018 | |
3019 | math_Matrix InterpolationMatrix(1, Parameters.Length(), |
3020 | 1, 2 * Degree + 1) ; |
3021 | ErrorCode = |
3022 | BSplCLib::BuildBSpMatrix(Parameters, |
3023 | ContactOrderArray, |
3024 | FlatKnots, |
3025 | Degree, |
3026 | InterpolationMatrix, |
3027 | UpperBandWidth, |
3028 | LowerBandWidth) ; |
3029 | Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ; |
3030 | |
3031 | ErrorCode = |
3032 | BSplCLib::FactorBandedMatrix(InterpolationMatrix, |
3033 | UpperBandWidth, |
3034 | LowerBandWidth, |
3035 | InversionProblem) ; |
3036 | Standard_OutOfRange_Raise_if (ErrorCode != 0, "BSplCLib::Interpolate") ; |
3037 | |
3038 | ErrorCode = |
3039 | BSplCLib::SolveBandedSystem(InterpolationMatrix, |
3040 | UpperBandWidth, |
3041 | LowerBandWidth, |
3042 | Standard_False, |
3043 | ArrayDimension, |
3044 | Poles, |
3045 | Weights) ; |
3046 | |
3047 | Standard_OutOfRange_Raise_if (ErrorCode != 0,"BSplCLib::Interpolate") ; |
3048 | } |
3049 | |
3050 | //======================================================================= |
3051 | //function : Evaluates a Bspline function : uses the ExtrapMode |
3052 | //purpose : the function is extrapolated using the Taylor expansion |
3053 | // of degree ExtrapMode[0] to the left and the Taylor |
3054 | // expansion of degree ExtrapMode[1] to the right |
3055 | // this evaluates the numerator by multiplying by the weights |
3056 | // and evaluating it but does not call RationalDerivatives after |
3057 | //======================================================================= |
3058 | |
3059 | void BSplCLib::Eval |
3060 | (const Standard_Real Parameter, |
3061 | const Standard_Boolean PeriodicFlag, |
3062 | const Standard_Integer DerivativeRequest, |
3063 | Standard_Integer& ExtrapMode, |
3064 | const Standard_Integer Degree, |
3065 | const TColStd_Array1OfReal& FlatKnots, |
3066 | const Standard_Integer ArrayDimension, |
3067 | Standard_Real& Poles, |
3068 | Standard_Real& Weights, |
3069 | Standard_Real& PolesResults, |
3070 | Standard_Real& WeightsResults) |
3071 | { |
3072 | Standard_Integer ii, |
3073 | jj, |
3074 | kk=0, |
3075 | Index, |
3076 | Index1, |
3077 | Index2, |
3078 | *ExtrapModeArray, |
3079 | Modulus, |
3080 | NewRequest, |
3081 | ExtrapolatingFlag[2], |
3082 | ErrorCode, |
3083 | ReturnCode, |
3084 | Order = Degree + 1, |
3085 | FirstNonZeroBsplineIndex, |
3086 | LocalRequest = DerivativeRequest ; |
3087 | Standard_Real *PResultArray, |
3088 | *WResultArray, |
3089 | *PolesArray, |
3090 | *WeightsArray, |
3091 | LocalParameter, |
3092 | Period, |
3093 | Inverse, |
3094 | Delta ; |
3095 | PolesArray = &Poles ; |
3096 | WeightsArray = &Weights ; |
3097 | ExtrapModeArray = &ExtrapMode ; |
3098 | PResultArray = &PolesResults ; |
3099 | WResultArray = &WeightsResults ; |
3100 | LocalParameter = Parameter ; |
3101 | ExtrapolatingFlag[0] = |
3102 | ExtrapolatingFlag[1] = 0 ; |
3103 | // |
3104 | // check if we are extrapolating to a degree which is smaller than |
3105 | // the degree of the Bspline |
3106 | // |
3107 | if (PeriodicFlag) { |
3108 | Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ; |
3109 | |
3110 | while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) { |
3111 | LocalParameter -= Period ; |
3112 | } |
3113 | |
3114 | while (LocalParameter < FlatKnots(2)) { |
3115 | LocalParameter += Period ; |
3116 | } |
3117 | } |
3118 | if (Parameter < FlatKnots(2) && |
3119 | LocalRequest < ExtrapModeArray[0] && |
3120 | ExtrapModeArray[0] < Degree) { |
3121 | LocalRequest = ExtrapModeArray[0] ; |
3122 | LocalParameter = FlatKnots(2) ; |
3123 | ExtrapolatingFlag[0] = 1 ; |
3124 | } |
3125 | if (Parameter > FlatKnots(FlatKnots.Upper()-1) && |
3126 | LocalRequest < ExtrapModeArray[1] && |
3127 | ExtrapModeArray[1] < Degree) { |
3128 | LocalRequest = ExtrapModeArray[1] ; |
3129 | LocalParameter = FlatKnots(FlatKnots.Upper()-1) ; |
3130 | ExtrapolatingFlag[1] = 1 ; |
3131 | } |
3132 | Delta = Parameter - LocalParameter ; |
3133 | if (LocalRequest >= Order) { |
3134 | LocalRequest = Degree ; |
3135 | } |
3136 | if (PeriodicFlag) { |
3137 | Modulus = FlatKnots.Length() - Degree -1 ; |
3138 | } |
3139 | else { |
3140 | Modulus = FlatKnots.Length() - Degree ; |
3141 | } |
3142 | |
3143 | BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order); |
3144 | ErrorCode = |
3145 | BSplCLib::EvalBsplineBasis(1, |
3146 | LocalRequest, |
3147 | Order, |
3148 | FlatKnots, |
3149 | LocalParameter, |
3150 | FirstNonZeroBsplineIndex, |
3151 | BsplineBasis) ; |
3152 | if (ErrorCode != 0) { |
3153 | ReturnCode = 1 ; |
3154 | goto FINISH ; |
3155 | } |
3156 | if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) { |
3157 | Index = 0 ; |
3158 | Index2 = 0 ; |
3159 | |
3160 | for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) { |
3161 | Index1 = FirstNonZeroBsplineIndex ; |
3162 | |
3163 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3164 | PResultArray[Index + kk] = 0.0e0 ; |
3165 | } |
3166 | WResultArray[Index] = 0.0e0 ; |
3167 | |
3168 | for (jj = 1 ; jj <= Order ; jj++) { |
3169 | |
3170 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3171 | PResultArray[Index + kk] += |
3172 | PolesArray[(Index1-1) * ArrayDimension + kk] |
3173 | * WeightsArray[Index1-1] * BsplineBasis(ii,jj) ; |
3174 | } |
3175 | WResultArray[Index2] += WeightsArray[Index1-1] * BsplineBasis(ii,jj) ; |
3176 | |
3177 | Index1 = Index1 % Modulus ; |
3178 | Index1 += 1 ; |
3179 | } |
3180 | Index += ArrayDimension ; |
3181 | Index2 += 1 ; |
3182 | } |
3183 | } |
3184 | else { |
3185 | // |
3186 | // store Taylor expansion in LocalRealArray |
3187 | // |
3188 | NewRequest = DerivativeRequest ; |
3189 | if (NewRequest > Degree) { |
3190 | NewRequest = Degree ; |
3191 | } |
3192 | BSplCLib_LocalArray LocalRealArray((LocalRequest + 1)*ArrayDimension); |
3193 | Index = 0 ; |
3194 | Inverse = 1.0e0 ; |
3195 | |
3196 | for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) { |
3197 | Index1 = FirstNonZeroBsplineIndex ; |
3198 | |
3199 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3200 | LocalRealArray[Index + kk] = 0.0e0 ; |
3201 | } |
3202 | |
3203 | for (jj = 1 ; jj <= Order ; jj++) { |
3204 | |
3205 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3206 | LocalRealArray[Index + kk] += |
3207 | PolesArray[(Index1-1)*ArrayDimension + kk] * |
3208 | WeightsArray[Index1-1] * BsplineBasis(ii,jj) ; |
3209 | } |
3210 | Index1 = Index1 % Modulus ; |
3211 | Index1 += 1 ; |
3212 | } |
3213 | |
3214 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3215 | LocalRealArray[Index + kk] *= Inverse ; |
3216 | } |
3217 | Index += ArrayDimension ; |
3218 | Inverse /= (Standard_Real) ii ; |
3219 | } |
3220 | PLib::EvalPolynomial(Delta, |
3221 | NewRequest, |
3222 | Degree, |
3223 | ArrayDimension, |
3224 | LocalRealArray[0], |
3225 | PolesResults) ; |
3226 | Index = 0 ; |
3227 | Inverse = 1.0e0 ; |
3228 | |
3229 | for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) { |
3230 | Index1 = FirstNonZeroBsplineIndex ; |
3231 | LocalRealArray[Index] = 0.0e0 ; |
3232 | |
3233 | for (jj = 1 ; jj <= Order ; jj++) { |
3234 | LocalRealArray[Index] += |
3235 | WeightsArray[Index1-1] * BsplineBasis(ii,jj) ; |
3236 | Index1 = Index1 % Modulus ; |
3237 | Index1 += 1 ; |
3238 | } |
3239 | LocalRealArray[Index + kk] *= Inverse ; |
3240 | Index += 1 ; |
3241 | Inverse /= (Standard_Real) ii ; |
3242 | } |
3243 | PLib::EvalPolynomial(Delta, |
3244 | NewRequest, |
3245 | Degree, |
3246 | 1, |
3247 | LocalRealArray[0], |
3248 | WeightsResults) ; |
3249 | } |
3250 | FINISH : ; |
3251 | } |
3252 | |
3253 | //======================================================================= |
3254 | //function : Evaluates a Bspline function : uses the ExtrapMode |
3255 | //purpose : the function is extrapolated using the Taylor expansion |
3256 | // of degree ExtrapMode[0] to the left and the Taylor |
3257 | // expansion of degree ExtrapMode[1] to the right |
3258 | // WARNING : the array Results is supposed to have at least |
3259 | // (DerivativeRequest + 1) * ArrayDimension slots and the |
3260 | // |
3261 | //======================================================================= |
3262 | |
3263 | void BSplCLib::Eval |
3264 | (const Standard_Real Parameter, |
3265 | const Standard_Boolean PeriodicFlag, |
3266 | const Standard_Integer DerivativeRequest, |
3267 | Standard_Integer& ExtrapMode, |
3268 | const Standard_Integer Degree, |
3269 | const TColStd_Array1OfReal& FlatKnots, |
3270 | const Standard_Integer ArrayDimension, |
3271 | Standard_Real& Poles, |
3272 | Standard_Real& Results) |
3273 | { |
3274 | Standard_Integer ii, |
3275 | jj, |
3276 | kk, |
3277 | Index, |
3278 | Index1, |
3279 | *ExtrapModeArray, |
3280 | Modulus, |
3281 | NewRequest, |
3282 | ExtrapolatingFlag[2], |
3283 | ErrorCode, |
3284 | ReturnCode, |
3285 | Order = Degree + 1, |
3286 | FirstNonZeroBsplineIndex, |
3287 | LocalRequest = DerivativeRequest ; |
3288 | |
3289 | Standard_Real *ResultArray, |
3290 | *PolesArray, |
3291 | LocalParameter, |
3292 | Period, |
3293 | Inverse, |
3294 | Delta ; |
3295 | |
3296 | PolesArray = &Poles ; |
3297 | ExtrapModeArray = &ExtrapMode ; |
3298 | ResultArray = &Results ; |
3299 | LocalParameter = Parameter ; |
3300 | ExtrapolatingFlag[0] = |
3301 | ExtrapolatingFlag[1] = 0 ; |
3302 | // |
3303 | // check if we are extrapolating to a degree which is smaller than |
3304 | // the degree of the Bspline |
3305 | // |
3306 | if (PeriodicFlag) { |
3307 | Period = FlatKnots(FlatKnots.Upper() - 1) - FlatKnots(2) ; |
3308 | |
3309 | while (LocalParameter > FlatKnots(FlatKnots.Upper() - 1)) { |
3310 | LocalParameter -= Period ; |
3311 | } |
3312 | |
3313 | while (LocalParameter < FlatKnots(2)) { |
3314 | LocalParameter += Period ; |
3315 | } |
3316 | } |
3317 | if (Parameter < FlatKnots(2) && |
3318 | LocalRequest < ExtrapModeArray[0] && |
3319 | ExtrapModeArray[0] < Degree) { |
3320 | LocalRequest = ExtrapModeArray[0] ; |
3321 | LocalParameter = FlatKnots(2) ; |
3322 | ExtrapolatingFlag[0] = 1 ; |
3323 | } |
3324 | if (Parameter > FlatKnots(FlatKnots.Upper()-1) && |
3325 | LocalRequest < ExtrapModeArray[1] && |
3326 | ExtrapModeArray[1] < Degree) { |
3327 | LocalRequest = ExtrapModeArray[1] ; |
3328 | LocalParameter = FlatKnots(FlatKnots.Upper()-1) ; |
3329 | ExtrapolatingFlag[1] = 1 ; |
3330 | } |
3331 | Delta = Parameter - LocalParameter ; |
3332 | if (LocalRequest >= Order) { |
3333 | LocalRequest = Degree ; |
3334 | } |
3335 | |
3336 | if (PeriodicFlag) { |
3337 | Modulus = FlatKnots.Length() - Degree -1 ; |
3338 | } |
3339 | else { |
3340 | Modulus = FlatKnots.Length() - Degree ; |
3341 | } |
3342 | |
3343 | BSplCLib_LocalMatrix BsplineBasis (LocalRequest, Order); |
3344 | |
3345 | ErrorCode = |
3346 | BSplCLib::EvalBsplineBasis(1, |
3347 | LocalRequest, |
3348 | Order, |
3349 | FlatKnots, |
3350 | LocalParameter, |
3351 | FirstNonZeroBsplineIndex, |
3352 | BsplineBasis); |
3353 | if (ErrorCode != 0) { |
3354 | ReturnCode = 1 ; |
3355 | goto FINISH ; |
3356 | } |
3357 | if (ExtrapolatingFlag[0] == 0 && ExtrapolatingFlag[1] == 0) { |
3358 | Index = 0 ; |
3359 | |
3360 | for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) { |
3361 | Index1 = FirstNonZeroBsplineIndex ; |
3362 | |
3363 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3364 | ResultArray[Index + kk] = 0.0e0 ; |
3365 | } |
3366 | |
3367 | for (jj = 1 ; jj <= Order ; jj++) { |
3368 | |
3369 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3370 | ResultArray[Index + kk] += |
3371 | PolesArray[(Index1-1) * ArrayDimension + kk] * BsplineBasis(ii,jj) ; |
3372 | } |
3373 | Index1 = Index1 % Modulus ; |
3374 | Index1 += 1 ; |
3375 | } |
3376 | Index += ArrayDimension ; |
3377 | } |
3378 | } |
3379 | else { |
3380 | // |
3381 | // store Taylor expansion in LocalRealArray |
3382 | // |
3383 | NewRequest = DerivativeRequest ; |
3384 | if (NewRequest > Degree) { |
3385 | NewRequest = Degree ; |
3386 | } |
3387 | BSplCLib_LocalArray LocalRealArray((LocalRequest + 1)*ArrayDimension); |
3388 | |
3389 | Index = 0 ; |
3390 | Inverse = 1.0e0 ; |
3391 | |
3392 | for (ii = 1 ; ii <= LocalRequest + 1 ; ii++) { |
3393 | Index1 = FirstNonZeroBsplineIndex ; |
3394 | |
3395 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3396 | LocalRealArray[Index + kk] = 0.0e0 ; |
3397 | } |
3398 | |
3399 | for (jj = 1 ; jj <= Order ; jj++) { |
3400 | |
3401 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3402 | LocalRealArray[Index + kk] += |
3403 | PolesArray[(Index1-1)*ArrayDimension + kk] * BsplineBasis(ii,jj) ; |
3404 | } |
3405 | Index1 = Index1 % Modulus ; |
3406 | Index1 += 1 ; |
3407 | } |
3408 | |
3409 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
3410 | LocalRealArray[Index + kk] *= Inverse ; |
3411 | } |
3412 | Index += ArrayDimension ; |
3413 | Inverse /= (Standard_Real) ii ; |
3414 | } |
3415 | PLib::EvalPolynomial(Delta, |
3416 | NewRequest, |
3417 | Degree, |
3418 | ArrayDimension, |
3419 | LocalRealArray[0], |
3420 | Results) ; |
3421 | } |
3422 | FINISH : ; |
3423 | } |
3424 | |
3425 | //======================================================================= |
3426 | //function : TangExtendToConstraint |
3427 | //purpose : Extends a Bspline function using the tangency map |
3428 | // WARNING : |
3429 | // |
3430 | // |
3431 | //======================================================================= |
3432 | |
3433 | void BSplCLib::TangExtendToConstraint |
3434 | (const TColStd_Array1OfReal& FlatKnots, |
3435 | const Standard_Real C1Coefficient, |
3436 | const Standard_Integer NumPoles, |
3437 | Standard_Real& Poles, |
3438 | const Standard_Integer CDimension, |
3439 | const Standard_Integer CDegree, |
3440 | const TColStd_Array1OfReal& ConstraintPoint, |
3441 | const Standard_Integer Continuity, |
3442 | const Standard_Boolean After, |
3443 | Standard_Integer& NbPolesResult, |
3444 | Standard_Integer& NbKnotsResult, |
3445 | Standard_Real& KnotsResult, |
3446 | Standard_Real& PolesResult) |
3447 | { |
3448 | #if DEB |
3449 | if (CDegree<Continuity+1) { |
3450 | cout<<"The BSpline degree must be greater than the order of continuity"<<endl; |
3451 | } |
3452 | #endif |
3453 | Standard_Real * Padr = &Poles ; |
3454 | Standard_Real * KRadr = &KnotsResult ; |
3455 | Standard_Real * PRadr = &PolesResult ; |
3456 | |
3457 | //////////////////////////////////////////////////////////////////////// |
3458 | // |
3459 | // 1. calcul du prolongement nD |
3460 | // |
3461 | //////////////////////////////////////////////////////////////////////// |
3462 | |
3463 | // matrice d'Hermite |
3464 | Standard_Integer Csize = Continuity + 2; |
3465 | math_Matrix MatCoefs(1,Csize, 1,Csize); |
3466 | if (After) { |
3467 | PLib::HermiteCoefficients(0, 1, // Les Bornes |
3468 | Continuity, 0, // Les Ordres de contraintes |
3469 | MatCoefs); |
3470 | } |
3471 | else { |
3472 | PLib::HermiteCoefficients(0, 1, // Les Bornes |
3473 | 0, Continuity, // Les Ordres de contraintes |
3474 | MatCoefs); |
3475 | } |
3476 | |
3477 | |
3478 | // positionnement au noeud de raccord |
3479 | Standard_Real Tbord ; |
3480 | if (After) { |
3481 | Tbord = FlatKnots(FlatKnots.Upper()-CDegree); |
3482 | } |
3483 | else { |
3484 | Tbord = FlatKnots(FlatKnots.Lower()+CDegree); |
3485 | } |
3486 | Standard_Boolean periodic_flag = Standard_False ; |
3487 | Standard_Integer ipos, extrap_mode[2], derivative_request = Max(Continuity,1); |
3488 | extrap_mode[0] = extrap_mode[1] = CDegree; |
3489 | TColStd_Array1OfReal EvalBS(1, CDimension * (derivative_request+1)) ; |
3490 | Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ; |
3491 | BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0], |
3492 | CDegree,FlatKnots,CDimension,Poles,*Eadr); |
3493 | |
3494 | // norme de la tangente au noeud de raccord |
3495 | math_Vector Tgte(1,CDimension); |
3496 | |
3497 | for (ipos=1;ipos<=CDimension;ipos++) { |
3498 | Tgte(ipos) = EvalBS(ipos+CDimension); |
3499 | } |
3500 | Standard_Real L1=Tgte.Norm(); |
3501 | |
3502 | |
3503 | // matrice de contraintes |
3504 | math_Matrix Contraintes(1,Csize,1,CDimension); |
3505 | if (After) { |
3506 | |
3507 | for (ipos=1;ipos<=CDimension;ipos++) { |
3508 | Contraintes(1,ipos) = EvalBS(ipos); |
3509 | Contraintes(2,ipos) = C1Coefficient * EvalBS(ipos+CDimension); |
3510 | if(Continuity >= 2) Contraintes(3,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2); |
3511 | if(Continuity >= 3) Contraintes(4,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3); |
3512 | Contraintes(Continuity+2,ipos) = ConstraintPoint(ipos); |
3513 | } |
3514 | } |
3515 | else { |
3516 | |
3517 | for (ipos=1;ipos<=CDimension;ipos++) { |
3518 | Contraintes(1,ipos) = ConstraintPoint(ipos); |
3519 | Contraintes(2,ipos) = EvalBS(ipos); |
3520 | if(Continuity >= 1) Contraintes(3,ipos) = C1Coefficient * EvalBS(ipos+CDimension); |
3521 | if(Continuity >= 2) Contraintes(4,ipos) = EvalBS(ipos+2*CDimension) * Pow(C1Coefficient,2); |
3522 | if(Continuity >= 3) Contraintes(5,ipos) = EvalBS(ipos+3*CDimension) * Pow(C1Coefficient,3); |
3523 | } |
3524 | } |
3525 | |
3526 | // calcul des coefficients du prolongement |
3527 | Standard_Integer ii, jj, kk; |
3528 | TColStd_Array1OfReal ExtraCoeffs(1,Csize*CDimension); |
3529 | ExtraCoeffs.Init(0.); |
3530 | |
3531 | for (ii=1; ii<=Csize; ii++) { |
3532 | |
3533 | for (jj=1; jj<=Csize; jj++) { |
3534 | |
3535 | for (kk=1; kk<=CDimension; kk++) { |
3536 | ExtraCoeffs(kk+(jj-1)*CDimension) += MatCoefs(ii,jj)*Contraintes(ii,kk); |
3537 | } |
3538 | } |
3539 | } |
3540 | |
3541 | // calcul des poles du prolongement |
3542 | TColStd_Array1OfReal ExtrapPoles(1,Csize*CDimension); |
3543 | Standard_Real * EPadr = &ExtrapPoles(1) ; |
3544 | PLib::CoefficientsPoles(CDimension, |
3545 | ExtraCoeffs, PLib::NoWeights(), |
3546 | ExtrapPoles, PLib::NoWeights()); |
3547 | |
3548 | // calcul des noeuds du prolongement avec leurs multiplicites |
3549 | TColStd_Array1OfReal ExtrapNoeuds(1,2); |
3550 | ExtrapNoeuds(1) = 0.; |
3551 | ExtrapNoeuds(2) = 1.; |
3552 | TColStd_Array1OfInteger ExtrapMults(1,2); |
3553 | ExtrapMults(1) = Csize; |
3554 | ExtrapMults(2) = Csize; |
3555 | |
3556 | // noeuds plats du prolongement |
3557 | TColStd_Array1OfReal FK2(1, Csize*2); |
3558 | BSplCLib::KnotSequence(ExtrapNoeuds,ExtrapMults,FK2); |
3559 | |
3560 | // norme de la tangente au point de raccord |
3561 | if (After) { |
3562 | BSplCLib::Eval(0.,periodic_flag,1,extrap_mode[0], |
3563 | Csize-1,FK2,CDimension,*EPadr,*Eadr); |
3564 | } |
3565 | else { |
3566 | BSplCLib::Eval(1.,periodic_flag,1,extrap_mode[0], |
3567 | Csize-1,FK2,CDimension,*EPadr,*Eadr); |
3568 | } |
3569 | |
3570 | for (ipos=1;ipos<=CDimension;ipos++) { |
3571 | Tgte(ipos) = EvalBS(ipos+CDimension); |
3572 | } |
3573 | Standard_Real L2 = Tgte.Norm(); |
3574 | |
3575 | // harmonisation des degres |
3576 | TColStd_Array1OfReal NewP2(1, (CDegree+1)*CDimension); |
3577 | TColStd_Array1OfReal NewK2(1, 2); |
3578 | TColStd_Array1OfInteger NewM2(1, 2); |
3579 | if (Csize-1<CDegree) { |
3580 | BSplCLib::IncreaseDegree(Csize-1,CDegree,Standard_False,CDimension, |
3581 | ExtrapPoles,ExtrapNoeuds,ExtrapMults, |
3582 | NewP2,NewK2,NewM2); |
3583 | } |
3584 | else { |
3585 | NewP2 = ExtrapPoles; |
3586 | NewK2 = ExtrapNoeuds; |
3587 | NewM2 = ExtrapMults; |
3588 | } |
3589 | |
3590 | // noeuds plats du prolongement apres harmonisation des degres |
3591 | TColStd_Array1OfReal NewFK2(1, (CDegree+1)*2); |
3592 | BSplCLib::KnotSequence(NewK2,NewM2,NewFK2); |
3593 | |
3594 | |
3595 | //////////////////////////////////////////////////////////////////////// |
3596 | // |
3597 | // 2. concatenation C0 |
3598 | // |
3599 | //////////////////////////////////////////////////////////////////////// |
3600 | |
3601 | // ratio de reparametrisation |
3602 | Standard_Real Ratio=1, Delta; |
3603 | if ( (L1 > Precision::Confusion()) && (L2 > Precision::Confusion()) ) { |
3604 | Ratio = L2 / L1; |
3605 | } |
3606 | if ( (Ratio < 1.e-5) || (Ratio > 1.e5) ) Ratio = 1; |
3607 | |
3608 | if (After) { |
3609 | // on ne bouge pas la premiere BSpline |
3610 | Delta = Ratio*NewFK2(NewFK2.Lower()) - FlatKnots(FlatKnots.Upper()); |
3611 | } |
3612 | else { |
3613 | // on ne bouge pas la seconde BSpline |
3614 | Delta = Ratio*NewFK2(NewFK2.Upper()) - FlatKnots(FlatKnots.Lower()); |
3615 | } |
3616 | |
3617 | // resultat de la concatenation |
3618 | Standard_Integer NbP1 = NumPoles, NbP2 = CDegree+1; |
3619 | Standard_Integer NbK1 = FlatKnots.Length(), NbK2 = 2*(CDegree+1); |
3620 | TColStd_Array1OfReal NewPoles (1, (NbP1+ NbP2-1)*CDimension); |
3621 | TColStd_Array1OfReal NewFlats (1, NbK1+NbK2-CDegree-2); |
3622 | |
3623 | // les poles |
3624 | Standard_Integer indNP, indP, indEP; |
3625 | if (After) { |
3626 | |
3627 | for (ii=1; ii<=NbP1+NbP2-1; ii++) { |
3628 | |
3629 | for (jj=1; jj<=CDimension; jj++) { |
3630 | indNP = (ii-1)*CDimension+jj; |
3631 | indP = (ii-1)*CDimension+jj-1; |
3632 | indEP = (ii-NbP1)*CDimension+jj; |
3633 | if (ii<NbP1) NewPoles(indNP) = Padr[indP]; |
3634 | else NewPoles(indNP) = NewP2(indEP); |
3635 | } |
3636 | } |
3637 | } |
3638 | else { |
3639 | |
3640 | for (ii=1; ii<=NbP1+NbP2-1; ii++) { |
3641 | |
3642 | for (jj=1; jj<=CDimension; jj++) { |
3643 | indNP = (ii-1)*CDimension+jj; |
3644 | indEP = (ii-1)*CDimension+jj; |
3645 | indP = (ii-NbP2)*CDimension+jj-1; |
3646 | if (ii<NbP2) NewPoles(indNP) = NewP2(indEP); |
3647 | else NewPoles(indNP) = Padr[indP]; |
3648 | } |
3649 | } |
3650 | } |
3651 | |
3652 | // les noeuds plats |
3653 | if (After) { |
3654 | // on commence avec les noeuds de la surface initiale |
3655 | |
3656 | for (ii=1; ii<NbK1; ii++) { |
3657 | NewFlats(ii) = FlatKnots(FlatKnots.Lower()+ii-1); |
3658 | } |
3659 | // on continue avec les noeuds du prolongement reparametres |
3660 | |
3661 | for (ii=1; ii<=NbK2-CDegree-1; ii++) { |
3662 | NewFlats(NbK1+ii-1) = Ratio*NewFK2(NewFK2.Lower()+ii+CDegree) - Delta; |
3663 | } |
3664 | } |
3665 | else { |
3666 | // on commence avec les noeuds du prolongement reparametres |
3667 | |
3668 | for (ii=1; ii<NbK2-CDegree; ii++) { |
3669 | NewFlats(ii) = Ratio*NewFK2(NewFK2.Lower()+ii-1) - Delta; |
3670 | } |
3671 | // on continue avec les noeuds de la surface initiale |
3672 | |
3673 | for (ii=2; ii<=NbK1; ii++) { |
3674 | NewFlats(NbK2+ii-CDegree-2) = FlatKnots(FlatKnots.Lower()+ii-1); |
3675 | } |
3676 | } |
3677 | |
3678 | |
3679 | //////////////////////////////////////////////////////////////////////// |
3680 | // |
3681 | // 3. reduction de la multiplicite au noeud de raccord |
3682 | // |
3683 | //////////////////////////////////////////////////////////////////////// |
3684 | |
3685 | // nombre de noeuds distincts |
3686 | Standard_Integer KLength = 1; |
3687 | |
3688 | for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) { |
3689 | if (NewFlats(ii) != NewFlats(ii-1)) KLength++; |
3690 | } |
3691 | |
3692 | // noeuds plats --> noeuds + multiplicites |
3693 | TColStd_Array1OfReal NewKnots (1, KLength); |
3694 | TColStd_Array1OfInteger NewMults (1, KLength); |
3695 | NewMults.Init(1); |
3696 | jj = 1; |
3697 | NewKnots(jj) = NewFlats(1); |
3698 | |
3699 | for (ii=2; ii<=NbK1+NbK2-CDegree-2;ii++) { |
3700 | if (NewFlats(ii) == NewFlats(ii-1)) NewMults(jj)++; |
3701 | else { |
3702 | jj++; |
3703 | NewKnots(jj) = NewFlats(ii); |
3704 | } |
3705 | } |
3706 | |
3707 | // reduction de la multiplicite au second ou a l'avant-dernier noeud |
3708 | Standard_Integer Index = 2, M = CDegree; |
3709 | if (After) Index = KLength-1; |
3710 | TColStd_Array1OfReal ResultPoles (1, (NbP1+ NbP2-1)*CDimension); |
3711 | TColStd_Array1OfReal ResultKnots (1, KLength); |
3712 | TColStd_Array1OfInteger ResultMults (1, KLength); |
3713 | Standard_Real Tol = 1.e-6; |
3714 | Standard_Boolean Ok = Standard_True; |
3715 | |
3716 | while ( (M>CDegree-Continuity) && Ok) { |
3717 | Ok = RemoveKnot(Index, M-1, CDegree, Standard_False, CDimension, |
3718 | NewPoles, NewKnots, NewMults, |
3719 | ResultPoles, ResultKnots, ResultMults, Tol); |
3720 | if (Ok) M--; |
3721 | } |
3722 | |
3723 | if (M == CDegree) { |
3724 | // le nombre de poles de la concatenation |
3725 | NbPolesResult = NbP1 + NbP2 - 1; |
3726 | // les poles de la concatenation |
3727 | Standard_Integer PLength = NbPolesResult*CDimension; |
3728 | |
3729 | for (jj=1; jj<=PLength; jj++) { |
3730 | PRadr[jj-1] = NewPoles(jj); |
3731 | } |
3732 | |
3733 | // les noeuds plats de la concatenation |
3734 | Standard_Integer ideb = 0; |
3735 | |
3736 | for (jj=0; jj<NewKnots.Length(); jj++) { |
3737 | for (ii=0; ii<NewMults(jj+1); ii++) { |
3738 | KRadr[ideb+ii] = NewKnots(jj+1); |
3739 | } |
3740 | ideb += NewMults(jj+1); |
3741 | } |
3742 | NbKnotsResult = ideb; |
3743 | } |
3744 | |
3745 | else { |
3746 | // le nombre de poles du resultat |
3747 | NbPolesResult = NbP1 + NbP2 - 1 - CDegree + M; |
3748 | // les poles du resultat |
3749 | Standard_Integer PLength = NbPolesResult*CDimension; |
3750 | |
3751 | for (jj=0; jj<PLength; jj++) { |
3752 | PRadr[jj] = ResultPoles(jj+1); |
3753 | } |
3754 | |
3755 | // les noeuds plats du resultat |
3756 | Standard_Integer ideb = 0; |
3757 | |
3758 | for (jj=0; jj<ResultKnots.Length(); jj++) { |
3759 | for (ii=0; ii<ResultMults(jj+1); ii++) { |
3760 | KRadr[ideb+ii] = ResultKnots(jj+1); |
3761 | } |
3762 | ideb += ResultMults(jj+1); |
3763 | } |
3764 | NbKnotsResult = ideb; |
3765 | } |
3766 | } |
3767 | |
3768 | //======================================================================= |
3769 | //function : Resolution |
3770 | //purpose : |
3771 | // d |
3772 | // Soit C(t) = SUM Ci Bi(t) une courbe Bspline de degre d |
3773 | // i = 1,n |
3774 | // dont les noeuds sont tj pour j = 1,n+d+1 |
3775 | // |
3776 | // |
3777 | // ' C1 - Ci-1 d-1 |
3778 | // Alors C (t) = SUM d * --------- Bi (t) |
3779 | // i = 2,n ti+d - ti |
3780 | // |
3781 | // d-1 |
3782 | // pour la base de BSpline Bi (t) de degre d-1. |
3783 | // |
3784 | // Par suite un majorant de la norme de la derivee de C est : |
3785 | // |
3786 | // |
3787 | // | Ci - Ci-1 | |
3788 | // d * Max | --------- | |
3789 | // i = 2,n | ti+d - ti | |
3790 | // |
3791 | // N(t) |
3792 | // Dans le cas rationel on pose C(t) = ----- |
3793 | // D(t) |
3794 | // |
3795 | // |
3796 | // D(t) = SUM Di Bi(t) |
3797 | // i=1,n |
3798 | // |
3799 | // N(t) = SUM Di * Ci Bi(t) |
3800 | // i =1,n |
3801 | // |
3802 | // N'(t) - D'(t) C(t) |
3803 | // C'(t) = ----------------------- |
3804 | // D(t) |
3805 | // |
3806 | // |
3807 | // N'(t) - D'(t) C(t) = |
3808 | // |
3809 | // Di * (Ci - C(t)) - Di-1 * (Ci-1 - C(t)) d-1 |
3810 | // SUM d * ---------------------------------------- * Bi (t) = |
3811 | // i=2,n ti+d - ti |
3812 | // |
3813 | // |
3814 | // Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj) d-1 |
3815 | // SUM SUM d * ----------------------------------- * Betaj(t) * Bi (t) |
3816 | //i=2,n j=1,n ti+d - ti |
3817 | // |
3818 | // |
3819 | // |
3820 | // Dj Bj(t) |
3821 | // Betaj(t) = -------- |
3822 | // D(t) |
3823 | // |
3824 | // les Betaj(t) forment une partition >= 0 de l'unite dont le support |
3825 | // est tj, tj+d+1. Par suite si Rj = {j-d, ...., j+d+d+1} |
3826 | // obtient un majorant de la derivee de C en prenant : |
3827 | // |
3828 | // |
3829 | // |
3830 | // |
3831 | // |
3832 | // |
3833 | // Di * (Ci - Cj) - Di-1 * (Ci-1 - Cj) |
3834 | // Max Max d * ----------------------------------- |
3835 | // j=1,n i dans Rj ti+d - ti |
3836 | // |
3837 | // -------------------------------------------------------- |
3838 | // |
3839 | // Min Di |
3840 | // i =1,n |
3841 | // |
3842 | // |
3843 | //======================================================================= |
3844 | |
3845 | void BSplCLib::Resolution( Standard_Real& Poles, |
3846 | const Standard_Integer ArrayDimension, |
3847 | const Standard_Integer NumPoles, |
3848 | const TColStd_Array1OfReal& Weights, |
3849 | const TColStd_Array1OfReal& FlatKnots, |
3850 | const Standard_Integer Degree, |
3851 | const Standard_Real Tolerance3D, |
3852 | Standard_Real& UTolerance) |
3853 | { |
3854 | Standard_Integer ii,num_poles,ii_index,jj_index,ii_inDim; |
3855 | Standard_Integer lower,upper,ii_minus,jj,ii_miDim; |
3856 | Standard_Integer Deg1 = Degree + 1; |
3857 | Standard_Integer Deg2 = (Degree << 1) + 1; |
3858 | Standard_Real value,factor,W,min_weights,inverse; |
3859 | Standard_Real pa_ii_inDim_0, pa_ii_inDim_1, pa_ii_inDim_2, pa_ii_inDim_3; |
3860 | Standard_Real pa_ii_miDim_0, pa_ii_miDim_1, pa_ii_miDim_2, pa_ii_miDim_3; |
3861 | Standard_Real wg_ii_index, wg_ii_minus; |
3862 | Standard_Real *PA,max_derivative; |
3863 | const Standard_Real * FK = &FlatKnots(FlatKnots.Lower()); |
3864 | PA = &Poles; |
3865 | max_derivative = 0.0e0; |
3866 | num_poles = FlatKnots.Length() - Deg1; |
3867 | switch (ArrayDimension) { |
3868 | case 2 : { |
3869 | if (&Weights != NULL) { |
3870 | const Standard_Real * WG = &Weights(Weights.Lower()); |
3871 | min_weights = WG[0]; |
3872 | |
3873 | for (ii = 1 ; ii < NumPoles ; ii++) { |
3874 | W = WG[ii]; |
3875 | if (W < min_weights) min_weights = W; |
3876 | } |
3877 | |
3878 | for (ii = 1 ; ii < num_poles ; ii++) { |
3879 | ii_index = ii % NumPoles; |
3880 | ii_inDim = ii_index << 1; |
3881 | ii_minus = (ii - 1) % NumPoles; |
3882 | ii_miDim = ii_minus << 1; |
3883 | pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++; |
3884 | pa_ii_inDim_1 = PA[ii_inDim]; |
3885 | pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++; |
3886 | pa_ii_miDim_1 = PA[ii_miDim]; |
3887 | wg_ii_index = WG[ii_index]; |
3888 | wg_ii_minus = WG[ii_minus]; |
3889 | inverse = FK[ii + Degree] - FK[ii]; |
3890 | inverse = 1.0e0 / inverse; |
3891 | lower = ii - Deg1; |
3892 | if (lower < 0) lower = 0; |
3893 | upper = Deg2 + ii; |
3894 | if (upper > num_poles) upper = num_poles; |
3895 | |
3896 | for (jj = lower ; jj < upper ; jj++) { |
3897 | jj_index = jj % NumPoles; |
3898 | jj_index = jj_index << 1; |
3899 | value = 0.0e0; |
3900 | factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) - |
3901 | ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus)); |
3902 | if (factor < 0) factor = - factor; |
3903 | value += factor; jj_index++; |
3904 | factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) - |
3905 | ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus)); |
3906 | if (factor < 0) factor = - factor; |
3907 | value += factor; |
3908 | value *= inverse; |
3909 | if (max_derivative < value) max_derivative = value; |
3910 | } |
3911 | } |
3912 | max_derivative /= min_weights; |
3913 | } |
3914 | else { |
3915 | |
3916 | for (ii = 1 ; ii < num_poles ; ii++) { |
3917 | ii_index = ii % NumPoles; |
3918 | ii_index = ii_index << 1; |
3919 | ii_minus = (ii - 1) % NumPoles; |
3920 | ii_minus = ii_minus << 1; |
3921 | inverse = FK[ii + Degree] - FK[ii]; |
3922 | inverse = 1.0e0 / inverse; |
3923 | value = 0.0e0; |
3924 | factor = PA[ii_index] - PA[ii_minus]; |
3925 | if (factor < 0) factor = - factor; |
3926 | value += factor; ii_index++; ii_minus++; |
3927 | factor = PA[ii_index] - PA[ii_minus]; |
3928 | if (factor < 0) factor = - factor; |
3929 | value += factor; |
3930 | value *= inverse; |
3931 | if (max_derivative < value) max_derivative = value; |
3932 | } |
3933 | } |
3934 | break; |
3935 | } |
3936 | case 3 : { |
3937 | if (&Weights != NULL) { |
3938 | const Standard_Real * WG = &Weights(Weights.Lower()); |
3939 | min_weights = WG[0]; |
3940 | |
3941 | for (ii = 1 ; ii < NumPoles ; ii++) { |
3942 | W = WG[ii]; |
3943 | if (W < min_weights) min_weights = W; |
3944 | } |
3945 | |
3946 | for (ii = 1 ; ii < num_poles ; ii++) { |
3947 | ii_index = ii % NumPoles; |
3948 | ii_inDim = (ii_index << 1) + ii_index; |
3949 | ii_minus = (ii - 1) % NumPoles; |
3950 | ii_miDim = (ii_minus << 1) + ii_minus; |
3951 | pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++; |
3952 | pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++; |
3953 | pa_ii_inDim_2 = PA[ii_inDim]; |
3954 | pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++; |
3955 | pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++; |
3956 | pa_ii_miDim_2 = PA[ii_miDim]; |
3957 | wg_ii_index = WG[ii_index]; |
3958 | wg_ii_minus = WG[ii_minus]; |
3959 | inverse = FK[ii + Degree] - FK[ii]; |
3960 | inverse = 1.0e0 / inverse; |
3961 | lower = ii - Deg1; |
3962 | if (lower < 0) lower = 0; |
3963 | upper = Deg2 + ii; |
3964 | if (upper > num_poles) upper = num_poles; |
3965 | |
3966 | for (jj = lower ; jj < upper ; jj++) { |
3967 | jj_index = jj % NumPoles; |
3968 | jj_index = (jj_index << 1) + jj_index; |
3969 | value = 0.0e0; |
3970 | factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) - |
3971 | ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus)); |
3972 | if (factor < 0) factor = - factor; |
3973 | value += factor; jj_index++; |
3974 | factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) - |
3975 | ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus)); |
3976 | if (factor < 0) factor = - factor; |
3977 | value += factor; jj_index++; |
3978 | factor = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) - |
3979 | ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus)); |
3980 | if (factor < 0) factor = - factor; |
3981 | value += factor; |
3982 | value *= inverse; |
3983 | if (max_derivative < value) max_derivative = value; |
3984 | } |
3985 | } |
3986 | max_derivative /= min_weights; |
3987 | } |
3988 | else { |
3989 | |
3990 | for (ii = 1 ; ii < num_poles ; ii++) { |
3991 | ii_index = ii % NumPoles; |
3992 | ii_index = (ii_index << 1) + ii_index; |
3993 | ii_minus = (ii - 1) % NumPoles; |
3994 | ii_minus = (ii_minus << 1) + ii_minus; |
3995 | inverse = FK[ii + Degree] - FK[ii]; |
3996 | inverse = 1.0e0 / inverse; |
3997 | value = 0.0e0; |
3998 | factor = PA[ii_index] - PA[ii_minus]; |
3999 | if (factor < 0) factor = - factor; |
4000 | value += factor; ii_index++; ii_minus++; |
4001 | factor = PA[ii_index] - PA[ii_minus]; |
4002 | if (factor < 0) factor = - factor; |
4003 | value += factor; ii_index++; ii_minus++; |
4004 | factor = PA[ii_index] - PA[ii_minus]; |
4005 | if (factor < 0) factor = - factor; |
4006 | value += factor; |
4007 | value *= inverse; |
4008 | if (max_derivative < value) max_derivative = value; |
4009 | } |
4010 | } |
4011 | break; |
4012 | } |
4013 | case 4 : { |
4014 | if (&Weights != NULL) { |
4015 | const Standard_Real * WG = &Weights(Weights.Lower()); |
4016 | min_weights = WG[0]; |
4017 | |
4018 | for (ii = 1 ; ii < NumPoles ; ii++) { |
4019 | W = WG[ii]; |
4020 | if (W < min_weights) min_weights = W; |
4021 | } |
4022 | |
4023 | for (ii = 1 ; ii < num_poles ; ii++) { |
4024 | ii_index = ii % NumPoles; |
4025 | ii_inDim = ii_index << 2; |
4026 | ii_minus = (ii - 1) % NumPoles; |
4027 | ii_miDim = ii_minus << 2; |
4028 | pa_ii_inDim_0 = PA[ii_inDim]; ii_inDim++; |
4029 | pa_ii_inDim_1 = PA[ii_inDim]; ii_inDim++; |
4030 | pa_ii_inDim_2 = PA[ii_inDim]; ii_inDim++; |
4031 | pa_ii_inDim_3 = PA[ii_inDim]; |
4032 | pa_ii_miDim_0 = PA[ii_miDim]; ii_miDim++; |
4033 | pa_ii_miDim_1 = PA[ii_miDim]; ii_miDim++; |
4034 | pa_ii_miDim_2 = PA[ii_miDim]; ii_miDim++; |
4035 | pa_ii_miDim_3 = PA[ii_miDim]; |
4036 | wg_ii_index = WG[ii_index]; |
4037 | wg_ii_minus = WG[ii_minus]; |
4038 | inverse = FK[ii + Degree] - FK[ii]; |
4039 | inverse = 1.0e0 / inverse; |
4040 | lower = ii - Deg1; |
4041 | if (lower < 0) lower = 0; |
4042 | upper = Deg2 + ii; |
4043 | if (upper > num_poles) upper = num_poles; |
4044 | |
4045 | for (jj = lower ; jj < upper ; jj++) { |
4046 | jj_index = jj % NumPoles; |
4047 | jj_index = jj_index << 2; |
4048 | value = 0.0e0; |
4049 | factor = (((PA[jj_index] - pa_ii_inDim_0) * wg_ii_index) - |
4050 | ((PA[jj_index] - pa_ii_miDim_0) * wg_ii_minus)); |
4051 | if (factor < 0) factor = - factor; |
4052 | value += factor; jj_index++; |
4053 | factor = (((PA[jj_index] - pa_ii_inDim_1) * wg_ii_index) - |
4054 | ((PA[jj_index] - pa_ii_miDim_1) * wg_ii_minus)); |
4055 | if (factor < 0) factor = - factor; |
4056 | value += factor; jj_index++; |
4057 | factor = (((PA[jj_index] - pa_ii_inDim_2) * wg_ii_index) - |
4058 | ((PA[jj_index] - pa_ii_miDim_2) * wg_ii_minus)); |
4059 | if (factor < 0) factor = - factor; |
4060 | value += factor; jj_index++; |
4061 | factor = (((PA[jj_index] - pa_ii_inDim_3) * wg_ii_index) - |
4062 | ((PA[jj_index] - pa_ii_miDim_3) * wg_ii_minus)); |
4063 | if (factor < 0) factor = - factor; |
4064 | value += factor; |
4065 | value *= inverse; |
4066 | if (max_derivative < value) max_derivative = value; |
4067 | } |
4068 | } |
4069 | max_derivative /= min_weights; |
4070 | } |
4071 | else { |
4072 | |
4073 | for (ii = 1 ; ii < num_poles ; ii++) { |
4074 | ii_index = ii % NumPoles; |
4075 | ii_index = ii_index << 2; |
4076 | ii_minus = (ii - 1) % NumPoles; |
4077 | ii_minus = ii_minus << 2; |
4078 | inverse = FK[ii + Degree] - FK[ii]; |
4079 | inverse = 1.0e0 / inverse; |
4080 | value = 0.0e0; |
4081 | factor = PA[ii_index] - PA[ii_minus]; |
4082 | if (factor < 0) factor = - factor; |
4083 | value += factor; ii_index++; ii_minus++; |
4084 | factor = PA[ii_index] - PA[ii_minus]; |
4085 | if (factor < 0) factor = - factor; |
4086 | value += factor; ii_index++; ii_minus++; |
4087 | factor = PA[ii_index] - PA[ii_minus]; |
4088 | if (factor < 0) factor = - factor; |
4089 | value += factor; ii_index++; ii_minus++; |
4090 | factor = PA[ii_index] - PA[ii_minus]; |
4091 | if (factor < 0) factor = - factor; |
4092 | value += factor; |
4093 | value *= inverse; |
4094 | if (max_derivative < value) max_derivative = value; |
4095 | } |
4096 | } |
4097 | break; |
4098 | } |
4099 | default : { |
4100 | Standard_Integer kk; |
4101 | if (&Weights != NULL) { |
4102 | const Standard_Real * WG = &Weights(Weights.Lower()); |
4103 | min_weights = WG[0]; |
4104 | |
4105 | for (ii = 1 ; ii < NumPoles ; ii++) { |
4106 | W = WG[ii]; |
4107 | if (W < min_weights) min_weights = W; |
4108 | } |
4109 | |
4110 | for (ii = 1 ; ii < num_poles ; ii++) { |
4111 | ii_index = ii % NumPoles; |
4112 | ii_inDim = ii_index * ArrayDimension; |
4113 | ii_minus = (ii - 1) % NumPoles; |
4114 | ii_miDim = ii_minus * ArrayDimension; |
4115 | wg_ii_index = WG[ii_index]; |
4116 | wg_ii_minus = WG[ii_minus]; |
4117 | inverse = FK[ii + Degree] - FK[ii]; |
4118 | inverse = 1.0e0 / inverse; |
4119 | lower = ii - Deg1; |
4120 | if (lower < 0) lower = 0; |
4121 | upper = Deg2 + ii; |
4122 | if (upper > num_poles) upper = num_poles; |
4123 | |
4124 | for (jj = lower ; jj < upper ; jj++) { |
4125 | jj_index = jj % NumPoles; |
4126 | jj_index *= ArrayDimension; |
4127 | value = 0.0e0; |
4128 | |
4129 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
4130 | factor = (((PA[jj_index + kk] - PA[ii_inDim + kk]) * wg_ii_index) - |
4131 | ((PA[jj_index + kk] - PA[ii_miDim + kk]) * wg_ii_minus)); |
4132 | if (factor < 0) factor = - factor; |
4133 | value += factor; |
4134 | } |
4135 | value *= inverse; |
4136 | if (max_derivative < value) max_derivative = value; |
4137 | } |
4138 | } |
4139 | max_derivative /= min_weights; |
4140 | } |
4141 | else { |
4142 | |
4143 | for (ii = 1 ; ii < num_poles ; ii++) { |
4144 | ii_index = ii % NumPoles; |
4145 | ii_index *= ArrayDimension; |
4146 | ii_minus = (ii - 1) % NumPoles; |
4147 | ii_minus *= ArrayDimension; |
4148 | inverse = FK[ii + Degree] - FK[ii]; |
4149 | inverse = 1.0e0 / inverse; |
4150 | value = 0.0e0; |
4151 | |
4152 | for (kk = 0 ; kk < ArrayDimension ; kk++) { |
4153 | factor = PA[ii_index + kk] - PA[ii_minus + kk]; |
4154 | if (factor < 0) factor = - factor; |
4155 | value += factor; |
4156 | } |
4157 | value *= inverse; |
4158 | if (max_derivative < value) max_derivative = value; |
4159 | } |
4160 | } |
4161 | } |
4162 | } |
4163 | max_derivative *= Degree; |
4164 | if (max_derivative > RealSmall()) |
4165 | UTolerance = Tolerance3D / max_derivative; |
4166 | else |
4167 | UTolerance = Tolerance3D / RealSmall(); |
4168 | } |
4169 | |
4170 | //======================================================================= |
4171 | // function: FlatBezierKnots |
4172 | // purpose : |
4173 | //======================================================================= |
4174 | |
4175 | const Standard_Real& BSplCLib::FlatBezierKnots (const Standard_Integer Degree) |
4176 | { |
4177 | if ( Degree < 1 || Degree > MaxDegree() || MaxDegree() != 25 ) |
4178 | Standard_OutOfRange::Raise ("Bezier curve degree greater than maximal supported"); |
4179 | |
4180 | // array of flat knots for bezier curve of maximum 25 degree |
4181 | static Standard_Real knots[52] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, |
4182 | 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }; |
4183 | return knots[25-Degree]; |
4184 | } |