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