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