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