projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023953: Avoid unjustified use of $CASROOT and exporting to GIF in tests
[occt.git] / src / BSplSLib / BSplSLib.cxx
CommitLineData
b311480e 1// Created on: 1991-08-26
2// Created by: JCV
3// Copyright (c) 1991-1999 Matra Datavision
4// Copyright (c) 1999-2012 OPEN CASCADE SAS
5//
6// The content of this file is subject to the Open CASCADE Technology Public
7// License Version 6.5 (the "License"). You may not use the content of this file
8// except in compliance with the License. Please obtain a copy of the License
9// at http://www.opencascade.org and read it completely before using this file.
10//
11// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13//
14// The Original Code and all software distributed under the License is
15// distributed on an "AS IS" basis, without warranty of any kind, and the
16// Initial Developer hereby disclaims all such warranties, including without
17// limitation, any warranties of merchantability, fitness for a particular
18// purpose or non-infringement. Please see the License for the specific terms
19// and conditions governing the rights and limitations under the License.
20
7fd59977 21
22// Modifed RLE Aug 93 - Complete rewrite
23// xab 21-Mar-95 implemented cache mecanism
24// pmn 25-09-96 Interpolation
25// jct 25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative.
26// pmn 07-10-96 : Correction de DN dans le cas rationnal.
27// pmn 06-02-97 : Correction des poids dans RationalDerivative. (PRO700)
28
7fd59977 29#include <BSplSLib.ixx>
30#include <PLib.hxx>
41194117 31#include <PLib_LocalArray.hxx>
7fd59977 32#include <BSplCLib.hxx>
33#include <TColgp_Array2OfXYZ.hxx>
34#include <TColgp_Array1OfXYZ.hxx>
35#include <TColStd_HArray1OfInteger.hxx>
36#include <Standard_NotImplemented.hxx>
37#include <Standard_ConstructionError.hxx>
38#include <math_Matrix.hxx>
39
40// for null derivatives
41194117 41static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0};
7fd59977 42#ifdef WNT
43#define M_SQRT2 1.4142135623730950488016887
44#endif
45
46//=======================================================================
47//struct : BSplCLib_DataContainer
48//purpose: Auxiliary structure providing buffers for poles and knots used in
49// evaluation of bspline (allocated in the stack)
50//=======================================================================
51
41194117 52struct BSplSLib_DataContainer
7fd59977 53{
41194117 54 BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree)
7fd59977 55 {
41194117
K		 56 Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() ||
57 VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25,
58 "BSplSLib: bspline degree is greater than maximum supported");
7fd59977 59 }
60
61 Standard_Real poles[4*(25+1)*(25+1)];
62 Standard_Real knots1[2*25];
63 Standard_Real knots2[2*25];
64 Standard_Real ders[48];
65};
66
41194117 67typedef PLib_LocalArray BSplSLib_LocalArray;
7fd59977 68
69//**************************************************************************
70// Evaluation methods
71//**************************************************************************
72
73//=======================================================================
74//function : RationalDerivative
75//purpose : computes the rational derivatives when whe have the
76// the derivatives of the homogeneous numerator and the
77// the derivatives of the denominator
78//=======================================================================
79
80void BSplSLib::RationalDerivative(const Standard_Integer UDeg,
81 const Standard_Integer VDeg,
82 const Standard_Integer N,
83 const Standard_Integer M,
84 Standard_Real& HDerivatives,
85 Standard_Real& RDerivatives,
86 const Standard_Boolean All)
87{
88 //
89 // if All is True all derivatives are computed. if Not only
90 // the requested N, M is computed
91 //
92 // Numerator(u,v)
93 // let f(u,v) be a rational function = ------------------
94 // Denominator(u,v)
95 //
96 //
97 // Let (N,M) the order of the derivatives we want : then since
98 // we have :
99 //
100 // Numerator = f * Denominator
101 //
102 // we derive :
103 //
104 // (N,M) 1 ( (N M) (p q) (N -p M-q) )
105 // f = ------------ ( Numerator - SUM SUM a * f * Denominator )
106 // (0,0) ( p<N q<M p q )
107 // Denominator
108 //
109 // with :
110 //
111 // ( N ) ( M )
112 // a = ( ) ( )
113 // p q ( p ) ( q )
114 //
115 //
116 // HDerivatives is an array where derivatives are stored in the following form
117 // Numerator is assumee to have 3 functions that is a vector of dimension
118 // 3
119 //
120 // (0,0) (0,0) (0, DegV) (0, DegV)
121 // Numerator Denominator ... Numerator Denominator
122 //
123 // (1,0) (1,0) (1, DegV) (1, DegV)
124 // Numerator Denominator ... Numerator Denominator
125 //
126 // ...........................................................
127 //
128 //
129 // (DegU,0) (DegU,0) (DegU, DegV) (DegU, DegV)
130 // Numerator Denominator ... Numerator Denominator
131 //
132 //
133 Standard_Integer ii,jj,pp,qq,index,index1,index2;
134 Standard_Integer M1,M3,M4,N1,iiM1,iiM3,jjM1,ppM1,ppM3;
135 Standard_Integer MinN,MinN1,MinM,MinM1;
136 Standard_Integer index_u,index_u1,index_v,index_v1,index_w;
137
138 M1 = M + 1;
139 N1 = N + 1;
140 ii = N1 * M1;
141 M3 = (M1 << 1) + M1;
142 M4 = (VDeg + 1) << 2;
143
144 BSplSLib_LocalArray StoreDerivatives (All ? 0 : ii * 3);
145 Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives);
146 BSplSLib_LocalArray StoreW (ii);
147 Standard_Real *HomogeneousArray = &HDerivatives;
148 Standard_Real denominator,Pii,Pip,Pjq;
149
150 denominator = 1.0e0 / HomogeneousArray[3];
151 index_u = 0;
152 index_u1 = 0;
153 if (UDeg < N) MinN = UDeg;
154 else MinN = N;
155 if (VDeg < M) MinM = VDeg;
156 else MinM = M;
157 MinN1 = MinN + 1;
158 MinM1 = MinM + 1;
159 iiM1 = - M1;
160
161 for (ii = 0 ; ii < MinN1 ; ii++) {
162 iiM1 += M1;
163 index_v = index_u;
164 index_v1 = index_u1;
165 index_w = iiM1;
166
167 for (jj = 0 ; jj < MinM1 ; jj++) {
174178b9 168 RArray[index_v++] = HomogeneousArray[index_v1++];
169 RArray[index_v++] = HomogeneousArray[index_v1++];
170 RArray[index_v++] = HomogeneousArray[index_v1++];
171 StoreW[index_w++] = HomogeneousArray[index_v1++];
7fd59977 172 }
173
174 for (jj = MinM1 ; jj < M1 ; jj++) {
174178b9 175 RArray[index_v++] = 0.;
176 RArray[index_v++] = 0.;
177 RArray[index_v++] = 0.;
178 StoreW[index_w++] = 0.;
7fd59977 179 }
180 index_u1 += M4;
181 index_u += M3;
182 }
183 index_v = MinN1 * M3;
184 index_w = MinN1 * M1;
185
186 for (ii = MinN1 ; ii < N1 ; ii++) {
187
188 for (jj = 0 ; jj < M1 ; jj++) {
174178b9 189 RArray[index_v++] = 0.0e0;
190 RArray[index_v++] = 0.0e0;
191 RArray[index_v++] = 0.0e0;
192 StoreW[index_w++] = 0.0e0;
7fd59977 193 }
194 }
195
0d969553 196 // --------------- Calculation ----------------
7fd59977 197
198 iiM1 = - M1;
199 iiM3 = - M3;
7fd59977 200
201 for (ii = 0 ; ii <= N ; ii++) {
202 iiM1 += M1;
203 iiM3 += M3;
204 index1 = iiM3 - 3;
205 jjM1 = iiM1;
206
207 for (jj = 0 ; jj <= M ; jj++) {
208 jjM1 ++;
209 ppM1 = - M1;
210 ppM3 = - M3;
211 index1 += 3;
212
213 for (pp = 0 ; pp < ii ; pp++) {
214 ppM1 += M1;
215 ppM3 += M3;
216 index = ppM3;
217 index2 = jjM1 - ppM1;
218 Pip = PLib::Bin(ii,pp);
219
220 for (qq = 0 ; qq <= jj ; qq++) {
221 index2--;
222 Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2];
223 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
224 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
225 RArray[index1] -= Pjq * RArray[index]; index++;
226 index1 -= 2;
227 }
228 }
229 index = iiM3;
230 index2 = jj + 1;
231 Pii = PLib::Bin(ii,ii);
232
233 for (qq = 0 ; qq < jj ; qq++) {
234 index2--;
235 Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2];
236 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
237 RArray[index1] -= Pjq * RArray[index]; index++; index1++;
238 RArray[index1] -= Pjq * RArray[index]; index++;
239 index1 -= 2;
240 }
241 RArray[index1] *= denominator; index1++;
242 RArray[index1] *= denominator; index1++;
243 RArray[index1] *= denominator;
244 index1 -= 2;
245 }
246 }
247 if (!All) {
248 RArray = &RDerivatives;
249 index = N * M1 + M;
250 index = (index << 1) + index;
251 RArray[0] = StoreDerivatives[index]; index++;
252 RArray[1] = StoreDerivatives[index]; index++;
253 RArray[2] = StoreDerivatives[index];
254 }
255}
256
257//=======================================================================
258//function : PrepareEval
259//purpose :
260//=======================================================================
261
262//
263// PrepareEval :
264//
0d969553 265// Prepare all data for computing points :
7fd59977 266// local arrays of knots
267// local array of poles (multiplied by the weights if rational)
268//
269// The first direction to compute (smaller degree) is returned
270// and the poles are stored according to this direction.
271
272static Standard_Boolean PrepareEval
273(const Standard_Real U,
274 const Standard_Real V,
275 const Standard_Integer Uindex,
276 const Standard_Integer Vindex,
277 const Standard_Integer UDegree,
278 const Standard_Integer VDegree,
279 const Standard_Boolean URat,
280 const Standard_Boolean VRat,
281 const Standard_Boolean UPer,
282 const Standard_Boolean VPer,
283 const TColgp_Array2OfPnt& Poles,
284 const TColStd_Array2OfReal& Weights,
285 const TColStd_Array1OfReal& UKnots,
286 const TColStd_Array1OfReal& VKnots,
287 const TColStd_Array1OfInteger& UMults,
288 const TColStd_Array1OfInteger& VMults,
289 Standard_Real& u1, // first parameter to use
290 Standard_Real& u2, // second parameter to use
291 Standard_Integer& d1, // first degree
292 Standard_Integer& d2, // second degree
293 Standard_Boolean& rational,
294 BSplSLib_DataContainer& dc)
295{
296 rational = URat || VRat;
297 Standard_Integer uindex = Uindex;
298 Standard_Integer vindex = Vindex;
299 Standard_Integer UKLower = UKnots.Lower();
300 Standard_Integer UKUpper = UKnots.Upper();
301 Standard_Integer VKLower = VKnots.Lower();
302 Standard_Integer VKUpper = VKnots.Upper();
303 if (UDegree <= VDegree) {
304 // compute the indices
305 if (uindex < UKLower || uindex > UKUpper)
306 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1);
307 else u1 = U;
308 if (vindex < VKLower || vindex > VKUpper)
309 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2);
310 else u2 = V;
311 // get the knots
312 d1 = UDegree;
313 d2 = VDegree;
314 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1);
315 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2);
316 if (&UMults == NULL) uindex -= UKLower + UDegree;
317 else uindex = BSplCLib::PoleIndex
318 (UDegree,uindex,UPer,UMults);
319 if (&VMults == NULL) vindex -= VKLower + VDegree;
320 else vindex = BSplCLib::PoleIndex
321 (VDegree,vindex,VPer,VMults);
322 // get the poles
323// Standard_Integer i,j,k,ip,jp;
324 Standard_Integer i,j,ip,jp;
325 Standard_Real w, *pole = dc.poles;
326 d1 = UDegree;
327 d2 = VDegree;
328 Standard_Integer PLowerRow = Poles.LowerRow();
329 Standard_Integer PUpperRow = Poles.UpperRow();
330 Standard_Integer PLowerCol = Poles.LowerCol();
331 Standard_Integer PUpperCol = Poles.UpperCol();
332 if (rational) { // verify if locally non rational
333 rational = Standard_False;
334 ip = PLowerRow + uindex;
335 jp = PLowerCol + vindex;
336 w = Weights.Value(ip,jp);
337 Standard_Real eps = Epsilon(w);
338 Standard_Real dw;
339
340 for (i = 0; i <= UDegree && !rational; i++) {
341 jp = PLowerCol + vindex;
342
343 for (j = 0; j <= VDegree && !rational; j++) {
344 dw = Weights.Value(ip,jp) - w;
345 if (dw < 0) dw = - dw;
346 rational = dw > eps;
347 jp++;
348 if (jp > PUpperCol) jp = PLowerCol;
349 }
350 ip++;
351 if (ip > PUpperRow) ip = PLowerRow;
352 }
353 }
354 // copy the poles
355 ip = PLowerRow + uindex;
356 if (rational) {
357
358 for (i = 0; i <= d1; i++) {
359 jp = PLowerCol + vindex;
360
361 for (j = 0; j <= d2; j++) {
362 const gp_Pnt& P = Poles .Value(ip,jp);
363 pole[3] = w = Weights.Value(ip,jp);
364 pole[0] = P.X() * w;
365 pole[1] = P.Y() * w;
366 pole[2] = P.Z() * w;
367 pole += 4;
368 jp++;
369 if (jp > PUpperCol) jp = PLowerCol;
370 }
371 ip++;
372 if (ip > PUpperRow) ip = PLowerRow;
373 }
374 }
375 else {
376
377 for (i = 0; i <= d1; i++) {
378 jp = PLowerCol + vindex;
379
380 for (j = 0; j <= d2; j++) {
381 const gp_Pnt& P = Poles.Value(ip,jp);
382 pole[0] = P.X();
383 pole[1] = P.Y();
384 pole[2] = P.Z();
385 pole += 3;
386 jp++;
387 if (jp > PUpperCol) jp = PLowerCol;
388 }
389 ip++;
390 if (ip > PUpperRow) ip = PLowerRow;
391 }
392 }
393 return Standard_True;
394 }
395 else {
396 // compute the indices
397 if (uindex < UKLower || uindex > UKUpper)
398 BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2);
399 else u2 = U;
400 if (vindex < VKLower || vindex > VKUpper)
401 BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1);
402 else u1 = V;
403 // get the knots
404 d2 = UDegree;
405 d1 = VDegree;
406 BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2);
407 BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1);
408 if (&UMults == NULL) uindex -= UKLower + UDegree;
409 else uindex = BSplCLib::PoleIndex
410 (UDegree,uindex,UPer,UMults);
411 if (&VMults == NULL) vindex -= VKLower + VDegree;
412 else vindex = BSplCLib::PoleIndex
413 (VDegree,vindex,VPer,VMults);
414 // get the poles
415// Standard_Integer i,j,k,ip,jp;
416 Standard_Integer i,j,ip,jp;
417 Standard_Real w, *pole = dc.poles;
418 d1 = VDegree;
419 d2 = UDegree;
420 Standard_Integer PLowerRow = Poles.LowerRow();
421 Standard_Integer PUpperRow = Poles.UpperRow();
422 Standard_Integer PLowerCol = Poles.LowerCol();
423 Standard_Integer PUpperCol = Poles.UpperCol();
424 if (rational) { // verify if locally non rational
425 rational = Standard_False;
426 ip = PLowerRow + uindex;
427 jp = PLowerCol + vindex;
428 w = Weights.Value(ip,jp);
429 Standard_Real eps = Epsilon(w);
430 Standard_Real dw;
431
432 for (i = 0; i <= UDegree && !rational; i++) {
433 jp = PLowerCol + vindex;
434
435 for (j = 0; j <= VDegree && !rational; j++) {
436 dw = Weights.Value(ip,jp) - w;
437 if (dw < 0) dw = - dw;
438 rational = dw > eps;
439 jp++;
440 if (jp > PUpperCol) jp = PLowerCol;
441 }
442 ip++;
443 if (ip > PUpperRow) ip = PLowerRow;
444 }
445 }
446 // copy the poles
447 jp = PLowerCol + vindex;
448 if (rational) {
449
450 for (i = 0; i <= d1; i++) {
451 ip = PLowerRow + uindex;
452
453 for (j = 0; j <= d2; j++) {
454 const gp_Pnt& P = Poles .Value(ip,jp);
455 pole[3] = w = Weights.Value(ip,jp);
456 pole[0] = P.X() * w;
457 pole[1] = P.Y() * w;
458 pole[2] = P.Z() * w;
459 pole += 4;
460 ip++;
461 if (ip > PUpperRow) ip = PLowerRow;
462 }
463 jp++;
464 if (jp > PUpperCol) jp = PLowerCol;
465 }
466 }
467 else {
468
469 for (i = 0; i <= d1; i++) {
470 ip = PLowerRow + uindex;
471
472 for (j = 0; j <= d2; j++) {
473 const gp_Pnt& P = Poles.Value(ip,jp);
474 pole[0] = P.X();
475 pole[1] = P.Y();
476 pole[2] = P.Z();
477 pole += 3;
478 ip++;
479 if (ip > PUpperRow) ip = PLowerRow;
480 }
481 jp++;
482 if (jp > PUpperCol) jp = PLowerCol;
483 }
484 }
485 return Standard_False;
486 }
487}
488
489//=======================================================================
490//function : D0
491//purpose :
492//=======================================================================
493
494void BSplSLib::D0
495(const Standard_Real U,
496 const Standard_Real V,
497 const Standard_Integer UIndex,
498 const Standard_Integer VIndex,
499 const TColgp_Array2OfPnt& Poles,
500 const TColStd_Array2OfReal& Weights,
501 const TColStd_Array1OfReal& UKnots,
502 const TColStd_Array1OfReal& VKnots,
503 const TColStd_Array1OfInteger& UMults,
504 const TColStd_Array1OfInteger& VMults,
505 const Standard_Integer UDegree,
506 const Standard_Integer VDegree,
507 const Standard_Boolean URat,
508 const Standard_Boolean VRat,
509 const Standard_Boolean UPer,
510 const Standard_Boolean VPer,
511 gp_Pnt& P)
512{
513// Standard_Integer k ;
514 Standard_Real W ;
515 HomogeneousD0(U,
516 V,
517 UIndex,
518 VIndex,
519 Poles,
520 Weights,
521 UKnots,
522 VKnots,
523 UMults,
524 VMults,
525 UDegree,
526 VDegree,
527 URat,
528 VRat,
529 UPer,
530 VPer,
531 W,
532 P) ;
533 P.SetX(P.X() / W);
534 P.SetY(P.Y() / W);
535 P.SetZ(P.Z() / W);
536}
537
538//=======================================================================
539//function : D0
540//purpose :
541//=======================================================================
542
543void BSplSLib::HomogeneousD0
544(const Standard_Real U,
545 const Standard_Real V,
546 const Standard_Integer UIndex,
547 const Standard_Integer VIndex,
548 const TColgp_Array2OfPnt& Poles,
549 const TColStd_Array2OfReal& Weights,
550 const TColStd_Array1OfReal& UKnots,
551 const TColStd_Array1OfReal& VKnots,
552 const TColStd_Array1OfInteger& UMults,
553 const TColStd_Array1OfInteger& VMults,
554 const Standard_Integer UDegree,
555 const Standard_Integer VDegree,
556 const Standard_Boolean URat,
557 const Standard_Boolean VRat,
558 const Standard_Boolean UPer,
559 const Standard_Boolean VPer,
560 Standard_Real & W,
561 gp_Pnt& P)
562{
563 Standard_Boolean rational;
564// Standard_Integer k,dim;
565 Standard_Integer dim;
566 Standard_Real u1,u2;
567 Standard_Integer d1,d2;
568 W = 1.0e0 ;
569
570 BSplSLib_DataContainer dc (UDegree, VDegree);
571 PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
572 Poles,Weights,UKnots,VKnots,UMults,VMults,
573 u1,u2,d1,d2,rational,dc);
574 if (rational) {
575 dim = 4;
576 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
577 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
578 W = dc.poles[3];
579 P.SetX(dc.poles[0]);
580 P.SetY(dc.poles[1]);
581 P.SetZ(dc.poles[2]);
582 }
583 else {
584 dim = 3;
585 BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles);
586 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles);
587 P.SetX(dc.poles[0]);
588 P.SetY(dc.poles[1]);
589 P.SetZ(dc.poles[2]);
590 }
591}
592
593//=======================================================================
594//function : D1
595//purpose :
596//=======================================================================
597
598void BSplSLib::D1
599(const Standard_Real U,
600 const Standard_Real V,
601 const Standard_Integer UIndex,
602 const Standard_Integer VIndex,
603 const TColgp_Array2OfPnt& Poles,
604 const TColStd_Array2OfReal& Weights,
605 const TColStd_Array1OfReal& UKnots,
606 const TColStd_Array1OfReal& VKnots,
607 const TColStd_Array1OfInteger& UMults,
608 const TColStd_Array1OfInteger& VMults,
609 const Standard_Integer UDegree,
610 const Standard_Integer VDegree,
611 const Standard_Boolean URat,
612 const Standard_Boolean VRat,
613 const Standard_Boolean UPer,
614 const Standard_Boolean VPer,
615 gp_Pnt& P,
616 gp_Vec& Vu,
617 gp_Vec& Vv)
618{
619 Standard_Boolean rational;
620// Standard_Integer k,dim,dim2;
621 Standard_Integer dim,dim2;
622 Standard_Real u1,u2;
623 Standard_Integer d1,d2;
624 Standard_Real *result, *resVu, *resVv;
625 BSplSLib_DataContainer dc (UDegree, VDegree);
626 if (PrepareEval
627 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
628 Poles,Weights,UKnots,VKnots,UMults,VMults,
629 u1,u2,d1,d2,rational,dc)) {
630 if (rational) {
631 dim = 4;
632 dim2 = (d2 + 1) << 2;
633 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
634 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
635 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
636 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
637 result = dc.ders;
638 resVu = result + 6;
639 resVv = result + 3;
640 }
641 else {
642 dim = 3;
643 dim2 = d2 + 1;
644 dim2 = (dim2 << 1) + dim2;
645 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
646 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
647 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
648 result = dc.poles;
649 resVu = result + dim2;
650 resVv = result + 3;
651 }
652 }
653 else {
654 if (rational) {
655 dim = 4;
656 dim2 = (d2 + 1) << 2;
657 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
658 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
659 BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2));
660 BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders);
661 result = dc.ders;
662 resVu = result + 3;
663 resVv = result + 6;
664 }
665 else {
666 dim = 3;
667 dim2 = d2 + 1;
668 dim2 = (dim2 << 1) + dim2;
669 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles);
670 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles);
671 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2));
672 result = dc.poles;
673 resVu = result + 3;
674 resVv = result + dim2;
675 }
676 }
677
678 P .SetX(result[0]);
679 Vu.SetX(resVu [0]);
680 Vv.SetX(resVv [0]);
681
682 P .SetY(result[1]);
683 Vu.SetY(resVu [1]);
684 Vv.SetY(resVv [1]);
685
686 P .SetZ(result[2]);
687 Vu.SetZ(resVu [2]);
688 Vv.SetZ(resVv [2]);
689}
690
691//=======================================================================
692//function : D1
693//purpose :
694//=======================================================================
695
696void BSplSLib::HomogeneousD1
697(const Standard_Real U,
698 const Standard_Real V,
699 const Standard_Integer UIndex,
700 const Standard_Integer VIndex,
701 const TColgp_Array2OfPnt& Poles,
702 const TColStd_Array2OfReal& Weights,
703 const TColStd_Array1OfReal& UKnots,
704 const TColStd_Array1OfReal& VKnots,
705 const TColStd_Array1OfInteger& UMults,
706 const TColStd_Array1OfInteger& VMults,
707 const Standard_Integer UDegree,
708 const Standard_Integer VDegree,
709 const Standard_Boolean URat,
710 const Standard_Boolean VRat,
711 const Standard_Boolean UPer,
712 const Standard_Boolean VPer,
713 gp_Pnt& N,
714 gp_Vec& Nu,
715 gp_Vec& Nv,
716 Standard_Real& D,
717 Standard_Real& Du,
718 Standard_Real& Dv)
719{
720 Standard_Boolean rational;
721// Standard_Integer k,dim;
722 Standard_Integer dim;
723 Standard_Real u1,u2;
724 Standard_Integer d1,d2;
725
726 D = 1.0e0 ;
727 Du = 0.0e0 ;
728 Dv = 0.0e0 ;
729 BSplSLib_DataContainer dc (UDegree, VDegree);
730 Standard_Boolean ufirst = PrepareEval
731 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
732 Poles,Weights,UKnots,VKnots,UMults,VMults,
733 u1,u2,d1,d2,rational,dc);
734 dim = rational ? 4 : 3;
735
736 BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles);
737 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles);
738 BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1)));
739
740 Standard_Real *result, *resVu, *resVv;
741 result = dc.poles;
742 resVu = result + (ufirst ? dim*(d2+1) : dim);
743 resVv = result + (ufirst ? dim : dim*(d2+1));
744
745 N .SetX(result[0]);
746 Nu.SetX(resVu [0]);
747 Nv.SetX(resVv [0]);
748
749 N .SetY(result[1]);
750 Nu.SetY(resVu [1]);
751 Nv.SetY(resVv [1]);
752
753 N .SetZ(result[2]);
754 Nu.SetZ(resVu [2]);
755 Nv.SetZ(resVv [2]);
756
757 if (rational) {
758 D = result[3];
759 Du = resVu [3];
760 Dv = resVv [3];
761 }
762}
763
764//=======================================================================
765//function : D2
766//purpose :
767//=======================================================================
768
769void BSplSLib::D2
770(const Standard_Real U,
771 const Standard_Real V,
772 const Standard_Integer UIndex,
773 const Standard_Integer VIndex,
774 const TColgp_Array2OfPnt& Poles,
775 const TColStd_Array2OfReal& Weights,
776 const TColStd_Array1OfReal& UKnots,
777 const TColStd_Array1OfReal& VKnots,
778 const TColStd_Array1OfInteger& UMults,
779 const TColStd_Array1OfInteger& VMults,
780 const Standard_Integer UDegree,
781 const Standard_Integer VDegree,
782 const Standard_Boolean URat,
783 const Standard_Boolean VRat,
784 const Standard_Boolean UPer,
785 const Standard_Boolean VPer,
786 gp_Pnt& P,
787 gp_Vec& Vu,
788 gp_Vec& Vv,
789 gp_Vec& Vuu,
790 gp_Vec& Vvv,
791 gp_Vec& Vuv)
792{
793 Standard_Boolean rational;
794// Standard_Integer k,dim,dim2;
795 Standard_Integer dim,dim2;
796 Standard_Real u1,u2;
797 Standard_Integer d1,d2;
798 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv;
799 BSplSLib_DataContainer dc (UDegree, VDegree);
800 if (PrepareEval
801 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
802 Poles,Weights,UKnots,VKnots,UMults,VMults,
803 u1,u2,d1,d2,rational,dc)) {
804 if (rational) {
805 dim = 4;
806 dim2 = (d2 + 1) << 2;
807 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
808 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
809 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
810 if (d1 > 1)
811 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
812 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
813 result = dc.ders;
814 resVu = result + 9;
815 resVv = result + 3;
816 resVuu = result + 18;
817 resVvv = result + 6;
818 resVuv = result + 12;
819 }
820 else {
821 dim = 3;
822 dim2 = d2 + 1;
823 dim2 = (dim2 << 1) + dim2;
824 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
825 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
826 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
827 if (d1 > 1)
828 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
829 result = dc.poles;
830 resVu = result + dim2;
831 resVv = result + 3;
832 if (UDegree <= 1) resVuu = BSplSLib_zero;
833 else resVuu = result + (dim2 << 1);
834 if (VDegree <= 1) resVvv = BSplSLib_zero;
835 else resVvv = result + 6;
836 resVuv = result + (d2 << 1) + d2 + 6;
837 }
838 }
839 else {
840 if (rational) {
841 dim = 4;
842 dim2 = (d2 + 1) << 2;
843 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
844 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
845 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
846 if (d1 > 1)
847 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
848 BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders);
849 result = dc.ders;
850 resVu = result + 3;
851 resVv = result + 9;
852 resVuu = result + 6;
853 resVvv = result + 18;
854 resVuv = result + 12;
855 }
856 else {
857 dim = 3;
858 dim2 = d2 + 1;
859 dim2 = (dim2 << 1) + dim2;
860 BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles);
861 BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles);
862 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2));
863 if (d1 > 1)
864 BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
865 result = dc.poles;
866 resVu = result + 3;
867 resVv = result + dim2;
868 if (UDegree <= 1) resVuu = BSplSLib_zero;
869 else resVuu = result + 6;
870 if (VDegree <= 1) resVvv = BSplSLib_zero;
871 else resVvv = result + (dim2 << 1);
872 resVuv = result + (d2 << 1) + d2 + 6;
873 }
874 }
875
876 P .SetX(result[0]);
877 Vu .SetX(resVu [0]);
878 Vv .SetX(resVv [0]);
879 Vuu.SetX(resVuu[0]);
880 Vvv.SetX(resVvv[0]);
881 Vuv.SetX(resVuv[0]);
882
883 P .SetY(result[1]);
884 Vu .SetY(resVu [1]);
885 Vv .SetY(resVv [1]);
886 Vuu.SetY(resVuu[1]);
887 Vvv.SetY(resVvv[1]);
888 Vuv.SetY(resVuv[1]);
889
890 P .SetZ(result[2]);
891 Vu .SetZ(resVu [2]);
892 Vv .SetZ(resVv [2]);
893 Vuu.SetZ(resVuu[2]);
894 Vvv.SetZ(resVvv[2]);
895 Vuv.SetZ(resVuv[2]);
896}
897
898//=======================================================================
899//function : D3
900//purpose :
901//=======================================================================
902
903void BSplSLib::D3
904(const Standard_Real U,
905 const Standard_Real V,
906 const Standard_Integer UIndex,
907 const Standard_Integer VIndex,
908 const TColgp_Array2OfPnt& Poles,
909 const TColStd_Array2OfReal& Weights,
910 const TColStd_Array1OfReal& UKnots,
911 const TColStd_Array1OfReal& VKnots,
912 const TColStd_Array1OfInteger& UMults,
913 const TColStd_Array1OfInteger& VMults,
914 const Standard_Integer UDegree,
915 const Standard_Integer VDegree,
916 const Standard_Boolean URat,
917 const Standard_Boolean VRat,
918 const Standard_Boolean UPer,
919 const Standard_Boolean VPer,
920 gp_Pnt& P,
921 gp_Vec& Vu,
922 gp_Vec& Vv,
923 gp_Vec& Vuu,
924 gp_Vec& Vvv,
925 gp_Vec& Vuv,
926 gp_Vec& Vuuu,
927 gp_Vec& Vvvv,
928 gp_Vec& Vuuv,
929 gp_Vec& Vuvv)
930{
931 Standard_Boolean rational;
932// Standard_Integer k,dim,dim2;
933 Standard_Integer dim,dim2;
934 Standard_Real u1,u2;
935 Standard_Integer d1,d2;
936 Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv,
937 *resVuuu, *resVvvv, *resVuuv, *resVuvv;
938 BSplSLib_DataContainer dc (UDegree, VDegree);
939 if (PrepareEval
940 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
941 Poles,Weights,UKnots,VKnots,UMults,VMults,
942 u1,u2,d1,d2,rational,dc)) {
943 if (rational) {
944 dim = 4;
945 dim2 = (d2 + 1) << 2;
946 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
947 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
948 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
949 if (d1 > 1)
950 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
951 if (d1 > 2)
952 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
953 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
954 result = dc.ders;
955 resVu = result + 12;
956 resVv = result + 3;
957 resVuu = result + 24;
958 resVvv = result + 6;
959 resVuv = result + 15;
960 resVuuu = result + 36;
961 resVvvv = result + 9;
962 resVuuv = result + 27;
963 resVuvv = result + 18;
964 }
965 else {
966 dim = 3;
967 dim2 = (d2 + 1);
968 dim2 = (dim2 << 1) + dim2;
969 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
970 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
971 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
972 if (d1 > 1)
973 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
974 if (d1 > 2)
975 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
976 result = dc.poles;
977 resVu = result + dim2;
978 resVv = result + 3;
979 if (UDegree <= 1) {
980 resVuu = BSplSLib_zero;
981 resVuuv = BSplSLib_zero;
982 }
983 else {
984 resVuu = result + (dim2 << 1);
985 resVuuv = result + (dim2 << 1) + 3;
986 }
987 if (VDegree <= 1) {
988 resVvv = BSplSLib_zero;
989 resVuvv = BSplSLib_zero;
990 }
991 else {
992 resVvv = result + 6;
993 resVuvv = result + dim2 + 6;
994 }
995 resVuv = result + (d2 << 1) + d2 + 6;
996 if (UDegree <= 2) resVuuu = BSplSLib_zero;
997 else resVuuu = result + (dim2 << 1) + dim2;
998 if (VDegree <= 2) resVvvv = BSplSLib_zero;
999 else resVvvv = result + 9;
1000 }
1001 }
1002 else {
1003 if (rational) {
1004 dim = 4;
1005 dim2 = (d2 + 1) << 2;
1006 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1007 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1008 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1009 if (d1 > 1)
1010 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1011 if (d1 > 2)
1012 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1013 BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders);
1014 result = dc.ders;
1015 resVu = result + 3;
1016 resVv = result + 12;
1017 resVuu = result + 6;
1018 resVvv = result + 24;
1019 resVuv = result + 15;
1020 resVuuu = result + 9;
1021 resVvvv = result + 36;
1022 resVuuv = result + 18;
1023 resVuvv = result + 27;
1024 }
1025 else {
1026 dim = 3;
1027 dim2 = (d2 + 1);
1028 dim2 = (dim2 << 1) + dim2;
1029 BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles);
1030 BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles);
1031 BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2));
1032 if (d1 > 1)
1033 BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1)));
1034 if (d1 > 2)
1035 BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2));
1036 result = dc.poles;
1037 resVu = result + 3;
1038 resVv = result + dim2;
1039 if (UDegree <= 1) {
1040 resVuu = BSplSLib_zero;
1041 resVuuv = BSplSLib_zero;
1042 }
1043 else {
1044 resVuu = result + 6;
1045 resVuuv = result + dim2 + 6;
1046 }
1047 if (VDegree <= 1) {
1048 resVvv = BSplSLib_zero;
1049 resVuvv = BSplSLib_zero;
1050 }
1051 else {
1052 resVvv = result + (dim2 << 1);
1053 resVuvv = result + (dim2 << 1) + 3;
1054 }
1055 resVuv = result + (d2 << 1) + d2 + 6;
1056 if (UDegree <= 2) resVuuu = BSplSLib_zero;
1057 else resVuuu = result + 9;
1058 if (VDegree <= 2) resVvvv = BSplSLib_zero;
1059 else resVvvv = result + (dim2 << 1) + dim2;
1060 }
1061 }
1062
1063 P .SetX(result [0]);
1064 Vu .SetX(resVu [0]);
1065 Vv .SetX(resVv [0]);
1066 Vuu .SetX(resVuu [0]);
1067 Vvv .SetX(resVvv [0]);
1068 Vuv .SetX(resVuv [0]);
1069 Vuuu.SetX(resVuuu[0]);
1070 Vvvv.SetX(resVvvv[0]);
1071 Vuuv.SetX(resVuuv[0]);
1072 Vuvv.SetX(resVuvv[0]);
1073
1074 P .SetY(result [1]);
1075 Vu .SetY(resVu [1]);
1076 Vv .SetY(resVv [1]);
1077 Vuu .SetY(resVuu [1]);
1078 Vvv .SetY(resVvv [1]);
1079 Vuv .SetY(resVuv [1]);
1080 Vuuu.SetY(resVuuu[1]);
1081 Vvvv.SetY(resVvvv[1]);
1082 Vuuv.SetY(resVuuv[1]);
1083 Vuvv.SetY(resVuvv[1]);
1084
1085 P .SetZ(result [2]);
1086 Vu .SetZ(resVu [2]);
1087 Vv .SetZ(resVv [2]);
1088 Vuu .SetZ(resVuu [2]);
1089 Vvv .SetZ(resVvv [2]);
1090 Vuv .SetZ(resVuv [2]);
1091 Vuuu.SetZ(resVuuu[2]);
1092 Vvvv.SetZ(resVvvv[2]);
1093 Vuuv.SetZ(resVuuv[2]);
1094 Vuvv.SetZ(resVuvv[2]);
1095}
1096
1097//=======================================================================
1098//function : DN
1099//purpose :
1100//=======================================================================
1101
1102void BSplSLib::DN
1103(const Standard_Real U,
1104 const Standard_Real V,
1105 const Standard_Integer Nu,
1106 const Standard_Integer Nv,
1107 const Standard_Integer UIndex,
1108 const Standard_Integer VIndex,
1109 const TColgp_Array2OfPnt& Poles,
1110 const TColStd_Array2OfReal& Weights,
1111 const TColStd_Array1OfReal& UKnots,
1112 const TColStd_Array1OfReal& VKnots,
1113 const TColStd_Array1OfInteger& UMults,
1114 const TColStd_Array1OfInteger& VMults,
1115 const Standard_Integer UDegree,
1116 const Standard_Integer VDegree,
1117 const Standard_Boolean URat,
1118 const Standard_Boolean VRat,
1119 const Standard_Boolean UPer,
1120 const Standard_Boolean VPer,
1121 gp_Vec& Vn)
1122{
1123 Standard_Boolean rational;
1124 Standard_Integer k,dim;
1125 Standard_Real u1,u2;
1126 Standard_Integer d1,d2;
1127
1128 BSplSLib_DataContainer dc (UDegree, VDegree);
1129 Standard_Boolean ufirst = PrepareEval
1130 (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer,
1131 Poles,Weights,UKnots,VKnots,UMults,VMults,
1132 u1,u2,d1,d2,rational,dc);
1133 dim = rational ? 4 : 3;
1134
1135 if (!rational) {
1136 if ((Nu > UDegree) || (Nv > VDegree)) {
1137 Vn.SetX(0.);
1138 Vn.SetY(0.);
1139 Vn.SetZ(0.);
1140 return;
1141 }
1142 }
1143
1144 Standard_Integer n1 = ufirst ? Nu : Nv;
1145 Standard_Integer n2 = ufirst ? Nv : Nu;
1146
1147 BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles);
1148
1149 for (k = 0; k <= Min(n1,d1); k++)
1150 BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1)));
1151
1152 Standard_Real *result;
1153 if (rational) {
1154 BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False);
1155 result = dc.ders; // because Standard_False ci-dessus.
1156
1157 }
1158 else {
1159 result = dc.poles + (n1 * (d2+1) + n2) * dim ;
1160 }
1161
1162 Vn.SetX(result[0]);
1163 Vn.SetY(result[1]);
1164 Vn.SetZ(result[2]);
1165}
1166
1167//
1168// Surface modifications
1169//
1170// a surface is processed as a curve of curves.
1171// i.e : a curve of parameter u where the current point is the set of poles
1172// of the iso.
1173//
1174
1175//=======================================================================
1176//function : Iso
1177//purpose :
1178//=======================================================================
1179
1180void BSplSLib::Iso(const Standard_Real Param,
1181 const Standard_Boolean IsU,
1182 const TColgp_Array2OfPnt& Poles,
1183 const TColStd_Array2OfReal& Weights,
1184 const TColStd_Array1OfReal& Knots,
1185 const TColStd_Array1OfInteger& Mults,
1186 const Standard_Integer Degree,
1187 const Standard_Boolean Periodic,
1188 TColgp_Array1OfPnt& CPoles,
1189 TColStd_Array1OfReal& CWeights)
1190{
1191 Standard_Integer index = 0;
1192 Standard_Real u = Param;
1193 Standard_Boolean rational = &Weights != NULL;
1194 Standard_Integer dim = rational ? 4 : 3;
1195
1196 // compute local knots
1197
1198 BSplSLib_LocalArray locknots1 (2*Degree);
1199 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
1200 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1);
1201 if (&Mults == NULL)
1202 index -= Knots.Lower() + Degree;
1203 else
1204 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
1205
1206
1207 // copy the local poles
1208
1209// Standard_Integer f1,l1,f2,l2,i,j,k;
1210 Standard_Integer f1,l1,f2,l2,i,j;
1211
1212 if (IsU) {
1213 f1 = Poles.LowerRow();
1214 l1 = Poles.UpperRow();
1215 f2 = Poles.LowerCol();
1216 l2 = Poles.UpperCol();
1217 }
1218 else {
1219 f1 = Poles.LowerCol();
1220 l1 = Poles.UpperCol();
1221 f2 = Poles.LowerRow();
1222 l2 = Poles.UpperRow();
1223 }
1224
1225 BSplSLib_LocalArray locpoles ((Degree+1) * (l2-f2+1) * dim);
1226
1227 Standard_Real w, *pole = locpoles;
1228 index += f1;
1229
1230 for (i = 0; i <= Degree; i++) {
1231
1232 for (j = f2; j <= l2; j++) {
1233
1234 const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index);
1235 if (rational) {
1236 pole[3] = w = IsU ? Weights(index,j) : Weights(j,index);
1237 pole[0] = P.X() * w;
1238 pole[1] = P.Y() * w;
1239 pole[2] = P.Z() * w;
1240 }
1241 else {
1242 pole[0] = P.X();
1243 pole[1] = P.Y();
1244 pole[2] = P.Z();
1245 }
1246 pole += dim;
1247 }
1248 index++;
1249 if (index > l1) index = f1;
1250 }
1251
1252 // compute the iso
1253 BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles);
1254
1255 // get the result
1256 pole = locpoles;
1257
1258 for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) {
1259 gp_Pnt& P = CPoles(i);
1260 if (rational) {
1261 CWeights(i) = w = pole[3];
1262 P.SetX( pole[0] / w);
1263 P.SetY( pole[1] / w);
1264 P.SetZ( pole[2] / w);
1265 }
1266 else {
1267 P.SetX( pole[0]);
1268 P.SetY( pole[1]);
1269 P.SetZ( pole[2]);
1270 }
1271 pole += dim;
1272 }
1273
1274 // if the input is not rational but weights are wanted
1275 if (!rational && (&CWeights != NULL)) {
1276
1277 for (i = CWeights.Lower(); i <= CWeights.Upper(); i++)
1278 CWeights(i) = 1.;
1279 }
1280}
1281
1282//=======================================================================
1283//function : Reverse
1284//purpose :
1285//=======================================================================
1286
1287void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles,
1288 const Standard_Integer Last,
1289 const Standard_Boolean UDirection)
1290{
1291 Standard_Integer i,j, l = Last;
1292 if ( UDirection) {
1293 l = Poles.LowerRow() +
1294 (l - Poles.LowerRow())%(Poles.ColLength());
1295 TColgp_Array2OfPnt temp(0, Poles.ColLength()-1,
1296 Poles.LowerCol(), Poles.UpperCol());
1297
1298 for (i = Poles.LowerRow(); i <= l; i++) {
1299
1300 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1301 temp(l-i,j) = Poles(i,j);
1302 }
1303 }
1304
1305 for (i = l+1; i <= Poles.UpperRow(); i++) {
1306
1307 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1308 temp(l+Poles.ColLength()-i,j) = Poles(i,j);
1309 }
1310 }
1311
1312 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1313
1314 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1315 Poles(i,j) = temp (i-Poles.LowerRow(),j);
1316 }
1317 }
1318 }
1319 else {
1320 l = Poles.LowerCol() +
1321 (l - Poles.LowerCol())%(Poles.RowLength());
1322 TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(),
1323 0, Poles.RowLength()-1);
1324
1325 for (j = Poles.LowerCol(); j <= l; j++) {
1326
1327 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1328 temp(i,l-j) = Poles(i,j);
1329 }
1330 }
1331
1332 for (j = l+1; j <= Poles.UpperCol(); j++) {
1333
1334 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1335 temp(i,l+Poles.RowLength()-j) = Poles(i,j);
1336 }
1337 }
1338
1339 for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) {
1340
1341 for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) {
1342 Poles(i,j) = temp (i,j-Poles.LowerCol());
1343 }
1344 }
1345 }
1346}
1347
1348//=======================================================================
1349//function : Reverse
1350//purpose :
1351//=======================================================================
1352
1353void BSplSLib::Reverse( TColStd_Array2OfReal& Weights,
1354 const Standard_Integer Last,
1355 const Standard_Boolean UDirection)
1356{
1357 Standard_Integer i,j, l = Last;
1358 if ( UDirection) {
1359 l = Weights.LowerRow() +
1360 (l - Weights.LowerRow())%(Weights.ColLength());
1361 TColStd_Array2OfReal temp(0, Weights.ColLength()-1,
1362 Weights.LowerCol(), Weights.UpperCol());
1363
1364 for (i = Weights.LowerRow(); i <= l; i++) {
1365
1366 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1367 temp(l-i,j) = Weights(i,j);
1368 }
1369 }
1370
1371 for (i = l+1; i <= Weights.UpperRow(); i++) {
1372
1373 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1374 temp(l+Weights.ColLength()-i,j) = Weights(i,j);
1375 }
1376 }
1377
1378 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1379
1380 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1381 Weights(i,j) = temp (i-Weights.LowerRow(),j);
1382 }
1383 }
1384 }
1385 else {
1386 l = Weights.LowerCol() +
1387 (l - Weights.LowerCol())%(Weights.RowLength());
1388 TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(),
1389 0, Weights.RowLength()-1);
1390
1391 for (j = Weights.LowerCol(); j <= l; j++) {
1392
1393 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1394 temp(i,l-j) = Weights(i,j);
1395 }
1396 }
1397
1398 for (j = l+1; j <= Weights.UpperCol(); j++) {
1399
1400 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1401 temp(i,l+Weights.RowLength()-j) = Weights(i,j);
1402 }
1403 }
1404
1405 for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) {
1406
1407 for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) {
1408 Weights(i,j) = temp (i,j-Weights.LowerCol());
1409 }
1410 }
1411 }
1412}
1413
1414//=======================================================================
1415//function : IsRational
1416//purpose :
1417//=======================================================================
1418
1419Standard_Boolean BSplSLib::IsRational
1420(const TColStd_Array2OfReal& Weights,
1421 const Standard_Integer I1,
1422 const Standard_Integer I2,
1423 const Standard_Integer J1,
1424 const Standard_Integer J2,
1425 const Standard_Real Epsi)
1426{
1427 Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2));
1428 Standard_Integer i,j;
1429 Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength();
1430 Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength();
1431
1432 for (i = I1 - fi; i < I2 - fi; i++) {
1433
1434 for (j = J1 - fj; j < J2 - fj; j++) {
1435 if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps)
1436 return Standard_True;
1437 }
1438 }
1439 return Standard_False;
1440}
1441
1442//=======================================================================
1443//function : SetPoles
1444//purpose :
1445//=======================================================================
1446
1447void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1448 TColStd_Array1OfReal& FP,
1449 const Standard_Boolean UDirection)
1450{
1451 Standard_Integer i, j, l = FP.Lower();
1452 Standard_Integer PLowerRow = Poles.LowerRow();
1453 Standard_Integer PUpperRow = Poles.UpperRow();
1454 Standard_Integer PLowerCol = Poles.LowerCol();
1455 Standard_Integer PUpperCol = Poles.UpperCol();
1456 if (UDirection) {
1457
1458 for ( i = PLowerRow; i <= PUpperRow; i++) {
1459
1460 for ( j = PLowerCol; j <= PUpperCol; j++) {
1461 const gp_Pnt& P = Poles.Value(i,j);
1462 FP(l) = P.X(); l++;
1463 FP(l) = P.Y(); l++;
1464 FP(l) = P.Z(); l++;
1465 }
1466 }
1467 }
1468 else {
1469
1470 for ( j = PLowerCol; j <= PUpperCol; j++) {
1471
1472 for ( i = PLowerRow; i <= PUpperRow; i++) {
1473 const gp_Pnt& P = Poles.Value(i,j);
1474 FP(l) = P.X(); l++;
1475 FP(l) = P.Y(); l++;
1476 FP(l) = P.Z(); l++;
1477 }
1478 }
1479 }
1480}
1481
1482//=======================================================================
1483//function : SetPoles
1484//purpose :
1485//=======================================================================
1486
1487void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles,
1488 const TColStd_Array2OfReal& Weights,
1489 TColStd_Array1OfReal& FP,
1490 const Standard_Boolean UDirection)
1491{
1492 Standard_Integer i, j, l = FP.Lower();
1493 Standard_Integer PLowerRow = Poles.LowerRow();
1494 Standard_Integer PUpperRow = Poles.UpperRow();
1495 Standard_Integer PLowerCol = Poles.LowerCol();
1496 Standard_Integer PUpperCol = Poles.UpperCol();
1497 if (UDirection) {
1498
1499 for ( i = PLowerRow; i <= PUpperRow; i++) {
1500
1501 for ( j = PLowerCol; j <= PUpperCol; j++) {
1502 const gp_Pnt& P = Poles .Value(i,j);
1503 Standard_Real w = Weights.Value(i,j);
1504 FP(l) = P.X() * w; l++;
1505 FP(l) = P.Y() * w; l++;
1506 FP(l) = P.Z() * w; l++;
1507 FP(l) = w; l++;
1508 }
1509 }
1510 }
1511 else {
1512
1513 for ( j = PLowerCol; j <= PUpperCol; j++) {
1514
1515 for ( i = PLowerRow; i <= PUpperRow; i++) {
1516 const gp_Pnt& P = Poles .Value(i,j);
1517 Standard_Real w = Weights.Value(i,j);
1518 FP(l) = P.X() * w; l++;
1519 FP(l) = P.Y() * w; l++;
1520 FP(l) = P.Z() * w; l++;
1521 FP(l) = w; l++;
1522 }
1523 }
1524 }
1525}
1526
1527//=======================================================================
1528//function : GetPoles
1529//purpose :
1530//=======================================================================
1531
1532void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1533 TColgp_Array2OfPnt& Poles,
1534 const Standard_Boolean UDirection)
1535{
1536 Standard_Integer i, j, l = FP.Lower();
1537 Standard_Integer PLowerRow = Poles.LowerRow();
1538 Standard_Integer PUpperRow = Poles.UpperRow();
1539 Standard_Integer PLowerCol = Poles.LowerCol();
1540 Standard_Integer PUpperCol = Poles.UpperCol();
1541 if (UDirection) {
1542
1543 for ( i = PLowerRow; i <= PUpperRow; i++) {
1544
1545 for ( j = PLowerCol; j <= PUpperCol; j++) {
1546 gp_Pnt& P = Poles.ChangeValue(i,j);
1547 P.SetX(FP(l)); l++;
1548 P.SetY(FP(l)); l++;
1549 P.SetZ(FP(l)); l++;
1550 }
1551 }
1552 }
1553 else {
1554
1555 for ( j = PLowerCol; j <= PUpperCol; j++) {
1556
1557 for ( i = PLowerRow; i <= PUpperRow; i++) {
1558 gp_Pnt& P = Poles.ChangeValue(i,j);
1559 P.SetX(FP(l)); l++;
1560 P.SetY(FP(l)); l++;
1561 P.SetZ(FP(l)); l++;
1562 }
1563 }
1564 }
1565}
1566
1567//=======================================================================
1568//function : GetPoles
1569//purpose :
1570//=======================================================================
1571
1572void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP,
1573 TColgp_Array2OfPnt& Poles,
1574 TColStd_Array2OfReal& Weights,
1575 const Standard_Boolean UDirection)
1576{
1577 Standard_Integer i, j, l = FP.Lower();
1578 Standard_Integer PLowerRow = Poles.LowerRow();
1579 Standard_Integer PUpperRow = Poles.UpperRow();
1580 Standard_Integer PLowerCol = Poles.LowerCol();
1581 Standard_Integer PUpperCol = Poles.UpperCol();
1582 if (UDirection) {
1583
1584 for ( i = PLowerRow; i <= PUpperRow; i++) {
1585
1586 for ( j = PLowerCol; j <= PUpperCol; j++) {
1587 Standard_Real w = FP( l + 3);
1588 Weights(i,j) = w;
1589 gp_Pnt& P = Poles.ChangeValue(i,j);
1590 P.SetX(FP(l) / w); l++;
1591 P.SetY(FP(l) / w); l++;
1592 P.SetZ(FP(l) / w); l++;
1593 l++;
1594 }
1595 }
1596 }
1597 else {
1598
1599 for ( j = PLowerCol; j <= PUpperCol; j++) {
1600
1601 for ( i = PLowerRow; i <= PUpperRow; i++) {
1602 Standard_Real w = FP( l + 3);
1603 Weights(i,j) = w;
1604 gp_Pnt& P = Poles.ChangeValue(i,j);
1605 P.SetX(FP(l) / w); l++;
1606 P.SetY(FP(l) / w); l++;
1607 P.SetZ(FP(l) / w); l++;
1608 l++;
1609 }
1610 }
1611 }
1612}
1613
1614//=======================================================================
1615//function : InsertKnots
1616//purpose :
1617//=======================================================================
1618
1619void BSplSLib::InsertKnots(const Standard_Boolean UDirection,
1620 const Standard_Integer Degree,
1621 const Standard_Boolean Periodic,
1622 const TColgp_Array2OfPnt& Poles,
1623 const TColStd_Array2OfReal& Weights,
1624 const TColStd_Array1OfReal& Knots,
1625 const TColStd_Array1OfInteger& Mults,
1626 const TColStd_Array1OfReal& AddKnots,
1627 const TColStd_Array1OfInteger& AddMults,
1628 TColgp_Array2OfPnt& NewPoles,
1629 TColStd_Array2OfReal& NewWeights,
1630 TColStd_Array1OfReal& NewKnots,
1631 TColStd_Array1OfInteger& NewMults,
1632 const Standard_Real Epsilon,
1633 const Standard_Boolean Add )
1634{
1635 Standard_Boolean rational = &Weights != NULL;
1636 Standard_Integer dim = 3;
1637 if (rational) dim++;
1638
1639 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1640 TColStd_Array1OfReal
1641 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1642
1643 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1644 else SetPoles(Poles,poles,UDirection);
1645
1646 if (UDirection) {
1647 dim *= Poles.RowLength();
1648 }
1649 else {
1650 dim *= Poles.ColLength();
1651 }
1652 BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
1653 AddKnots,AddMults,newpoles,NewKnots,NewMults,
1654 Epsilon,Add);
1655
1656 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1657 else GetPoles(newpoles,NewPoles,UDirection);
1658}
1659
1660//=======================================================================
1661//function : RemoveKnot
1662//purpose :
1663//=======================================================================
1664
1665Standard_Boolean BSplSLib::RemoveKnot
1666(const Standard_Boolean UDirection,
1667 const Standard_Integer Index,
1668 const Standard_Integer Mult,
1669 const Standard_Integer Degree,
1670 const Standard_Boolean Periodic,
1671 const TColgp_Array2OfPnt& Poles,
1672 const TColStd_Array2OfReal& Weights,
1673 const TColStd_Array1OfReal& Knots,
1674 const TColStd_Array1OfInteger& Mults,
1675 TColgp_Array2OfPnt& NewPoles,
1676 TColStd_Array2OfReal& NewWeights,
1677 TColStd_Array1OfReal& NewKnots,
1678 TColStd_Array1OfInteger& NewMults,
1679 const Standard_Real Tolerance)
1680{
1681 Standard_Boolean rational = &Weights != NULL;
1682 Standard_Integer dim = 3;
1683 if (rational) dim++;
1684
1685 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1686 TColStd_Array1OfReal
1687 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1688
1689 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1690 else SetPoles(Poles,poles,UDirection);
1691
1692 if (UDirection) {
1693 dim *= Poles.RowLength();
1694 }
1695 else {
1696 dim *= Poles.ColLength();
1697 }
1698
1699 if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim,
1700 poles,Knots,Mults,newpoles,NewKnots,NewMults,
1701 Tolerance))
1702 return Standard_False;
1703
1704 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1705 else GetPoles(newpoles,NewPoles,UDirection);
1706 return Standard_True;
1707}
1708
1709//=======================================================================
1710//function : IncreaseDegree
1711//purpose :
1712//=======================================================================
1713
1714void BSplSLib::IncreaseDegree
1715(const Standard_Boolean UDirection,
1716 const Standard_Integer Degree,
1717 const Standard_Integer NewDegree,
1718 const Standard_Boolean Periodic,
1719 const TColgp_Array2OfPnt& Poles,
1720 const TColStd_Array2OfReal& Weights,
1721 const TColStd_Array1OfReal& Knots,
1722 const TColStd_Array1OfInteger& Mults,
1723 TColgp_Array2OfPnt& NewPoles,
1724 TColStd_Array2OfReal& NewWeights,
1725 TColStd_Array1OfReal& NewKnots,
1726 TColStd_Array1OfInteger& NewMults)
1727{
1728 Standard_Boolean rational = &Weights != NULL;
1729 Standard_Integer dim = 3;
1730 if (rational) dim++;
1731
1732 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1733 TColStd_Array1OfReal
1734 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1735
1736 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1737 else SetPoles(Poles,poles,UDirection);
1738
1739 if (UDirection) {
1740 dim *= Poles.RowLength();
1741 }
1742 else {
1743 dim *= Poles.ColLength();
1744 }
1745
1746 BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
1747 newpoles,NewKnots,NewMults);
1748
1749 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1750 else GetPoles(newpoles,NewPoles,UDirection);
1751}
1752
1753//=======================================================================
1754//function : Unperiodize
1755//purpose :
1756//=======================================================================
1757
1758void BSplSLib::Unperiodize
1759(const Standard_Boolean UDirection,
1760 const Standard_Integer Degree,
1761 const TColStd_Array1OfInteger& Mults,
1762 const TColStd_Array1OfReal& Knots,
1763 const TColgp_Array2OfPnt& Poles,
1764 const TColStd_Array2OfReal& Weights,
1765 TColStd_Array1OfInteger& NewMults,
1766 TColStd_Array1OfReal& NewKnots,
1767 TColgp_Array2OfPnt& NewPoles,
1768 TColStd_Array2OfReal& NewWeights)
1769{
1770 Standard_Boolean rational = &Weights != NULL;
1771 Standard_Integer dim = 3;
1772 if (rational) dim++;
1773
1774 TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength());
1775 TColStd_Array1OfReal
1776 newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength());
1777
1778 if (rational) SetPoles(Poles,Weights,poles,UDirection);
1779 else SetPoles(Poles,poles,UDirection);
1780
1781 if (UDirection) {
1782 dim *= Poles.RowLength();
1783 }
1784 else {
1785 dim *= Poles.ColLength();
1786 }
1787
1788 BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles,
1789 NewMults, NewKnots, newpoles);
1790
1791 if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection);
1792 else GetPoles(newpoles,NewPoles,UDirection);
1793}
1794
1795//=======================================================================
1796//function : BuildCache
1797//purpose : Stores theTaylor expansion normalized between 0,1 in the
1798// Cache : in case of a rational function the Poles are
1799// stored in homogeneous form
1800//=======================================================================
1801
1802void BSplSLib::BuildCache
1803(const Standard_Real U,
1804 const Standard_Real V,
1805 const Standard_Real USpanDomain,
1806 const Standard_Real VSpanDomain,
1807 const Standard_Boolean UPeriodic,
1808 const Standard_Boolean VPeriodic,
1809 const Standard_Integer UDegree,
1810 const Standard_Integer VDegree,
1811 const Standard_Integer UIndex,
1812 const Standard_Integer VIndex,
1813 const TColStd_Array1OfReal& UFlatKnots,
1814 const TColStd_Array1OfReal& VFlatKnots,
1815 const TColgp_Array2OfPnt& Poles,
1816 const TColStd_Array2OfReal& Weights,
1817 TColgp_Array2OfPnt& CachePoles,
1818 TColStd_Array2OfReal& CacheWeights)
1819{
1820 Standard_Boolean rational,rational_u,rational_v,flag_u_or_v;
1821 Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index;
1822 Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2;
1823 if (&Weights != NULL)
1824 rational_u = rational_v = Standard_True;
1825 else
1826 rational_u = rational_v = Standard_False;
1827 BSplSLib_DataContainer dc (UDegree, VDegree);
1828 flag_u_or_v =
1829 PrepareEval (U,
1830 V,
1831 UIndex,
1832 VIndex,
1833 UDegree,
1834 VDegree,
1835 rational_u,
1836 rational_v,
1837 UPeriodic,
1838 VPeriodic,
1839 Poles,
1840 Weights,
1841 UFlatKnots,
1842 VFlatKnots,
1843 (BSplCLib::NoMults()),
1844 (BSplCLib::NoMults()),
1845 u1,
1846 u2,
1847 d1,
1848 d2,
1849 rational,
1850 dc);
1851 d1p1 = d1 + 1;
1852 d2p1 = d2 + 1;
1853 if (rational) {
1854 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles);
1855
1856 for (kk = 0; kk <= d1 ; kk++)
1857 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1));
1858 if (flag_u_or_v) {
1859 min_degree_domain = USpanDomain ;
1860 max_degree_domain = VSpanDomain ;
1861 }
1862 else {
1863 min_degree_domain = VSpanDomain ;
1864 max_degree_domain = USpanDomain ;
1865 }
1866 factor[0] = 1.0e0 ;
1867
1868 for (ii = 0 ; ii <= d2 ; ii++) {
1869 iii = ii + 1;
1870 factor[1] = 1.0e0 ;
1871
1872 for (jj = 0 ; jj <= d1 ; jj++) {
1873 jjj = jj + 1;
1874 Index = jj * d2p1 + ii ;
1875 Index = Index << 2;
1876 gp_Pnt& P = CachePoles(iii,jjj);
1877 f = factor[0] * factor[1];
1878 P.SetX( f * dc.poles[Index]); Index++;
1879 P.SetY( f * dc.poles[Index]); Index++;
1880 P.SetZ( f * dc.poles[Index]); Index++;
1881 CacheWeights(iii ,jjj) = f * dc.poles[Index] ;
1882 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1883 }
1884 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1885 }
1886 }
1887 else {
1888 BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles);
1889
1890 for (kk = 0; kk <= d1 ; kk++)
1891 BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1));
1892 if (flag_u_or_v) {
1893 min_degree_domain = USpanDomain ;
1894 max_degree_domain = VSpanDomain ;
1895 }
1896 else {
1897 min_degree_domain = VSpanDomain ;
1898 max_degree_domain = USpanDomain ;
1899 }
1900 factor[0] = 1.0e0 ;
1901
1902 for (ii = 0 ; ii <= d2 ; ii++) {
1903 iii = ii + 1;
1904 factor[1] = 1.0e0 ;
1905
1906 for (jj = 0 ; jj <= d1 ; jj++) {
1907 jjj = jj + 1;
1908 Index = jj * d2p1 + ii ;
1909 Index = (Index << 1) + Index;
1910 gp_Pnt& P = CachePoles(iii,jjj);
1911 f = factor[0] * factor[1];
1912 P.SetX( f * dc.poles[Index]); Index++;
1913 P.SetY( f * dc.poles[Index]); Index++;
1914 P.SetZ( f * dc.poles[Index]);
1915 factor[1] *= min_degree_domain / (Standard_Real) (jjj) ;
1916 }
1917 factor[0] *= max_degree_domain / (Standard_Real) (iii) ;
1918 }
1919 if (&Weights != NULL) {
1920 //
1921 // means that PrepareEval did found out that the surface was
1922 // locally polynomial but since the surface is constructed
1923 // with some weights we need to set the weight polynomial to constant
1924 //
1925
1926 for (ii = 1 ; ii <= d2p1 ; ii++) {
1927
1928 for (jj = 1 ; jj <= d1p1 ; jj++) {
1929 CacheWeights(ii,jj) = 0.0e0 ;
1930 }
1931 }
1932 CacheWeights(1,1) = 1.0e0 ;
1933 }
1934 }
1935}
1936
1937//=======================================================================
1938//function : CacheD0
1939//purpose : Evaluates the polynomial cache of the Bspline Curve
1940//
1941//=======================================================================
1942
1943void BSplSLib::CacheD0(const Standard_Real UParameter,
1944 const Standard_Real VParameter,
1945 const Standard_Integer UDegree,
1946 const Standard_Integer VDegree,
1947 const Standard_Real UCacheParameter,
1948 const Standard_Real VCacheParameter,
1949 const Standard_Real USpanLenght,
1950 const Standard_Real VSpanLenght,
1951 const TColgp_Array2OfPnt& PolesArray,
1952 const TColStd_Array2OfReal& WeightsArray,
1953 gp_Pnt& aPoint)
1954{
1955 //
1956 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
1957 // form
1958 // the SpanLenght is the normalizing factor so that the polynomial is between
1959 // 0 and 1
1960 Standard_Integer
1961// ii,
1962 dimension,
1963 min_degree,
1964 max_degree ;
1965
1966 Standard_Real
1967 new_parameter[2] ,
1968 inverse ;
1969
1970 Standard_Real *
1971 PArray = (Standard_Real *)
1972 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
1973 Standard_Real *
1974 myPoint = (Standard_Real *) &aPoint ;
1975 if (UDegree <= VDegree) {
1976 min_degree = UDegree ;
1977 max_degree = VDegree ;
1978 new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ;
1979 new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ;
1980 dimension = 3 * (UDegree + 1) ;
1981 }
1982 else {
1983 min_degree = VDegree ;
1984 max_degree = UDegree ;
1985 new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ;
1986 new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ;
1987 dimension = 3 * (VDegree + 1) ;
1988 }
1989 BSplSLib_LocalArray locpoles(dimension);
1990
1991 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
1992 max_degree,
1993 dimension,
1994 max_degree*dimension,
1995 PArray[0],
1996 locpoles[0]) ;
1997
1998 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
1999 min_degree,
2000 3,
2001 (min_degree << 1) + min_degree,
2002 locpoles[0],
2003 myPoint[0]) ;
2004 if (&WeightsArray != NULL) {
2005 dimension = min_degree + 1 ;
2006 Standard_Real *
2007 WArray = (Standard_Real *)
2008 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2009 PLib::NoDerivativeEvalPolynomial(new_parameter[0],
2010 max_degree,
2011 dimension,
2012 max_degree*dimension,
2013 WArray[0],
2014 locpoles[0]) ;
2015
2016 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2017 min_degree,
2018 1,
2019 min_degree,
2020 locpoles[0],
2021 inverse) ;
2022 inverse = 1.0e0/ inverse ;
2023
2024 myPoint[0] *= inverse ;
2025 myPoint[1] *= inverse ;
2026 myPoint[2] *= inverse ;
2027 }
2028}
2029
2030//=======================================================================
2031//function : CacheD1
2032//purpose : Evaluates the polynomial cache of the Bspline Curve
2033//
2034//=======================================================================
2035
2036void BSplSLib::CacheD1(const Standard_Real UParameter,
2037 const Standard_Real VParameter,
2038 const Standard_Integer UDegree,
2039 const Standard_Integer VDegree,
2040 const Standard_Real UCacheParameter,
2041 const Standard_Real VCacheParameter,
2042 const Standard_Real USpanLenght,
2043 const Standard_Real VSpanLenght,
2044 const TColgp_Array2OfPnt& PolesArray,
2045 const TColStd_Array2OfReal& WeightsArray,
2046 gp_Pnt& aPoint,
2047 gp_Vec& aVecU,
2048 gp_Vec& aVecV)
2049{
2050 //
2051 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2052 // form
2053 // the SpanLenght is the normalizing factor so that the polynomial is between
2054 // 0 and 1
2055 Standard_Integer
2056// ii,
2057// jj,
2058// kk,
2059 dimension,
2060 min_degree,
2061 max_degree ;
2062
2063 Standard_Real
2064 inverse_min,
2065 inverse_max,
2066 new_parameter[2] ;
2067
2068 Standard_Real *
2069 PArray = (Standard_Real *)
2070 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2071 Standard_Real local_poles_array[2][2][3],
2072 local_poles_and_weights_array[2][2][4],
2073 local_weights_array[2][2] ;
2074 Standard_Real * my_vec_min,
2075 * my_vec_max,
2076 * my_point ;
2077 my_point = (Standard_Real *) &aPoint ;
2078 //
2079 // initialize in case of rational evaluation
2080 // because RationalDerivative will use all
2081 // the coefficients
2082 //
2083 //
2084 if (&WeightsArray != NULL) {
2085
2086 local_poles_array [0][0][0] = 0.0e0 ;
2087 local_poles_array [0][0][1] = 0.0e0 ;
2088 local_poles_array [0][0][2] = 0.0e0 ;
2089 local_weights_array [0][0] = 0.0e0 ;
2090 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2091 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2092 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2093 local_poles_and_weights_array[0][0][3] = 0.0e0 ;
2094
2095 local_poles_array [0][1][0] = 0.0e0 ;
2096 local_poles_array [0][1][1] = 0.0e0 ;
2097 local_poles_array [0][1][2] = 0.0e0 ;
2098 local_weights_array [0][1] = 0.0e0 ;
2099 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2100 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2101 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2102 local_poles_and_weights_array[0][1][3] = 0.0e0 ;
2103
2104 local_poles_array [1][0][0] = 0.0e0 ;
2105 local_poles_array [1][0][1] = 0.0e0 ;
2106 local_poles_array [1][0][2] = 0.0e0 ;
2107 local_weights_array [1][0] = 0.0e0 ;
2108 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2109 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2110 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2111 local_poles_and_weights_array[1][0][3] = 0.0e0 ;
2112
2113 local_poles_array [1][1][0] = 0.0e0 ;
2114 local_poles_array [1][1][1] = 0.0e0 ;
2115 local_poles_array [1][1][2] = 0.0e0 ;
2116 local_weights_array [1][1] = 0.0e0 ;
2117 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2118 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2119 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2120 local_poles_and_weights_array[1][1][3] = 0.0e0 ;
2121 }
2122
2123 if (UDegree <= VDegree) {
2124 min_degree = UDegree ;
2125 max_degree = VDegree ;
2126 inverse_min = 1.0e0/USpanLenght ;
2127 inverse_max = 1.0e0/VSpanLenght ;
2128 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2129 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2130
2131 dimension = 3 * (UDegree + 1) ;
2132 my_vec_min = (Standard_Real *) &aVecU ;
2133 my_vec_max = (Standard_Real *) &aVecV ;
2134 }
2135 else {
2136 min_degree = VDegree ;
2137 max_degree = UDegree ;
2138 inverse_min = 1.0e0/VSpanLenght ;
2139 inverse_max = 1.0e0/USpanLenght ;
2140 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2141 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2142 dimension = 3 * (VDegree + 1) ;
2143 my_vec_min = (Standard_Real *) &aVecV ;
2144 my_vec_max = (Standard_Real *) &aVecU ;
2145 }
2146
2147 BSplSLib_LocalArray locpoles (2 * dimension);
2148
2149 PLib::EvalPolynomial(new_parameter[0],
2150 1,
2151 max_degree,
2152 dimension,
2153 PArray[0],
2154 locpoles[0]) ;
2155
2156 PLib::EvalPolynomial(new_parameter[1],
2157 1,
2158 min_degree,
2159 3,
2160 locpoles[0],
2161 local_poles_array[0][0][0]) ;
2162 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2163 min_degree,
2164 3,
2165 (min_degree << 1) + min_degree,
2166 locpoles[dimension],
2167 local_poles_array[1][0][0]) ;
2168
2169 if (&WeightsArray != NULL) {
2170 dimension = min_degree + 1 ;
2171 Standard_Real *
2172 WArray = (Standard_Real *)
2173 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2174 PLib::EvalPolynomial(new_parameter[0],
2175 1,
2176 max_degree,
2177 dimension,
2178 WArray[0],
2179 locpoles[0]) ;
2180
2181 PLib::EvalPolynomial(new_parameter[1],
2182 1,
2183 min_degree,
2184 1,
2185 locpoles[0],
2186 local_weights_array[0][0]) ;
2187 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2188 min_degree,
2189 1,
2190 min_degree,
2191 locpoles[dimension],
2192 local_weights_array[1][0]) ;
2193
2194 local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ;
2195 local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ;
2196 local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ;
2197 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ;
2198
2199 local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ;
2200 local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ;
2201 local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ;
2202 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ;
2203
2204 local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ;
2205 local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ;
2206 local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ;
2207 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ;
2208
2209 local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ;
2210 local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ;
2211 local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ;
2212 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ;
2213
2214 BSplSLib::RationalDerivative(1,
2215 1,
2216 1,
2217 1,
2218 local_poles_and_weights_array[0][0][0],
2219 local_poles_array[0][0][0]) ;
2220 }
2221
2222 my_point [0] = local_poles_array [0][0][0] ;
2223 my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ;
2224 my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ;
2225
2226 my_point [1] = local_poles_array [0][0][1] ;
2227 my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ;
2228 my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ;
2229
2230 my_point [2] = local_poles_array [0][0][2] ;
2231 my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ;
2232 my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ;
2233}
2234
2235//=======================================================================
2236//function : CacheD2
2237//purpose : Evaluates the polynomial cache of the Bspline Curve
2238//
2239//=======================================================================
2240
2241void BSplSLib::CacheD2(const Standard_Real UParameter,
2242 const Standard_Real VParameter,
2243 const Standard_Integer UDegree,
2244 const Standard_Integer VDegree,
2245 const Standard_Real UCacheParameter,
2246 const Standard_Real VCacheParameter,
2247 const Standard_Real USpanLenght,
2248 const Standard_Real VSpanLenght,
2249 const TColgp_Array2OfPnt& PolesArray,
2250 const TColStd_Array2OfReal& WeightsArray,
2251 gp_Pnt& aPoint,
2252 gp_Vec& aVecU,
2253 gp_Vec& aVecV,
2254 gp_Vec& aVecUU,
2255 gp_Vec& aVecUV,
2256 gp_Vec& aVecVV)
2257{
2258 //
2259 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
2260 // form
2261 // the SpanLenght is the normalizing factor so that the polynomial is between
2262 // 0 and 1
2263 Standard_Integer
2264 ii,
2265// jj,
2266 kk,
2267 index,
2268 dimension,
2269 min_degree,
2270 max_degree ;
2271
2272 Standard_Real
2273 inverse_min,
2274 inverse_max,
2275 new_parameter[2] ;
2276
2277 Standard_Real *
2278 PArray = (Standard_Real *)
2279 &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ;
2280 Standard_Real local_poles_array[3][3][3],
2281 local_poles_and_weights_array[3][3][4],
2282 local_weights_array[3][3] ;
2283 Standard_Real * my_vec_min,
2284 * my_vec_max,
2285 * my_vec_min_min,
2286 * my_vec_max_max,
2287 * my_vec_min_max,
2288 * my_point ;
2289 my_point = (Standard_Real *) &aPoint ;
2290
2291 //
2292 // initialize in case the min and max degree are less than 2
2293 //
2294 local_poles_array[0][0][0] = 0.0e0 ;
2295 local_poles_array[0][0][1] = 0.0e0 ;
2296 local_poles_array[0][0][2] = 0.0e0 ;
2297 local_poles_array[0][1][0] = 0.0e0 ;
2298 local_poles_array[0][1][1] = 0.0e0 ;
2299 local_poles_array[0][1][2] = 0.0e0 ;
2300 local_poles_array[0][2][0] = 0.0e0 ;
2301 local_poles_array[0][2][1] = 0.0e0 ;
2302 local_poles_array[0][2][2] = 0.0e0 ;
2303
2304 local_poles_array[1][0][0] = 0.0e0 ;
2305 local_poles_array[1][0][1] = 0.0e0 ;
2306 local_poles_array[1][0][2] = 0.0e0 ;
2307 local_poles_array[1][1][0] = 0.0e0 ;
2308 local_poles_array[1][1][1] = 0.0e0 ;
2309 local_poles_array[1][1][2] = 0.0e0 ;
2310 local_poles_array[1][2][0] = 0.0e0 ;
2311 local_poles_array[1][2][1] = 0.0e0 ;
2312 local_poles_array[1][2][2] = 0.0e0 ;
2313
2314 local_poles_array[2][0][0] = 0.0e0 ;
2315 local_poles_array[2][0][1] = 0.0e0 ;
2316 local_poles_array[2][0][2] = 0.0e0 ;
2317 local_poles_array[2][1][0] = 0.0e0 ;
2318 local_poles_array[2][1][1] = 0.0e0 ;
2319 local_poles_array[2][1][2] = 0.0e0 ;
2320 local_poles_array[2][2][0] = 0.0e0 ;
2321 local_poles_array[2][2][1] = 0.0e0 ;
2322 local_poles_array[2][2][2] = 0.0e0 ;
2323 //
2324 // initialize in case of rational evaluation
2325 // because RationalDerivative will use all
2326 // the coefficients
2327 //
2328 //
2329 if (&WeightsArray != NULL) {
2330
2331 local_poles_and_weights_array[0][0][0] = 0.0e0 ;
2332 local_poles_and_weights_array[0][0][1] = 0.0e0 ;
2333 local_poles_and_weights_array[0][0][2] = 0.0e0 ;
2334 local_poles_and_weights_array[0][1][0] = 0.0e0 ;
2335 local_poles_and_weights_array[0][1][1] = 0.0e0 ;
2336 local_poles_and_weights_array[0][1][2] = 0.0e0 ;
2337 local_poles_and_weights_array[0][2][0] = 0.0e0 ;
2338 local_poles_and_weights_array[0][2][1] = 0.0e0 ;
2339 local_poles_and_weights_array[0][2][2] = 0.0e0 ;
2340
2341 local_poles_and_weights_array[1][0][0] = 0.0e0 ;
2342 local_poles_and_weights_array[1][0][1] = 0.0e0 ;
2343 local_poles_and_weights_array[1][0][2] = 0.0e0 ;
2344 local_poles_and_weights_array[1][1][0] = 0.0e0 ;
2345 local_poles_and_weights_array[1][1][1] = 0.0e0 ;
2346 local_poles_and_weights_array[1][1][2] = 0.0e0 ;
2347 local_poles_and_weights_array[1][2][0] = 0.0e0 ;
2348 local_poles_and_weights_array[1][2][1] = 0.0e0 ;
2349 local_poles_and_weights_array[1][2][2] = 0.0e0 ;
2350
2351 local_poles_and_weights_array[2][0][0] = 0.0e0 ;
2352 local_poles_and_weights_array[2][0][1] = 0.0e0 ;
2353 local_poles_and_weights_array[2][0][2] = 0.0e0 ;
2354 local_poles_and_weights_array[2][1][0] = 0.0e0 ;
2355 local_poles_and_weights_array[2][1][1] = 0.0e0 ;
2356 local_poles_and_weights_array[2][1][2] = 0.0e0 ;
2357 local_poles_and_weights_array[2][2][0] = 0.0e0 ;
2358 local_poles_and_weights_array[2][2][1] = 0.0e0 ;
2359 local_poles_and_weights_array[2][2][2] = 0.0e0 ;
2360
2361 local_poles_and_weights_array[0][0][3] =
2362 local_weights_array[0][0] = 0.0e0 ;
2363 local_poles_and_weights_array[0][1][3] =
2364 local_weights_array[0][1] = 0.0e0 ;
2365 local_poles_and_weights_array[0][2][3] =
2366 local_weights_array[0][2] = 0.0e0 ;
2367 local_poles_and_weights_array[1][0][3] =
2368 local_weights_array[1][0] = 0.0e0 ;
2369 local_poles_and_weights_array[1][1][3] =
2370 local_weights_array[1][1] = 0.0e0 ;
2371 local_poles_and_weights_array[1][2][3] =
2372 local_weights_array[1][2] = 0.0e0 ;
2373 local_poles_and_weights_array[2][0][3] =
2374 local_weights_array[2][0] = 0.0e0 ;
2375 local_poles_and_weights_array[2][1][3] =
2376 local_weights_array[2][1] = 0.0e0 ;
2377 local_poles_and_weights_array[2][2][3] =
2378 local_weights_array[2][2] = 0.0e0 ;
2379 }
2380
2381 if (UDegree <= VDegree) {
2382 min_degree = UDegree ;
2383 max_degree = VDegree ;
2384 inverse_min = 1.0e0/USpanLenght ;
2385 inverse_max = 1.0e0/VSpanLenght ;
2386 new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ;
2387 new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ;
2388
2389 dimension = 3 * (UDegree + 1) ;
2390 my_vec_min = (Standard_Real *) &aVecU ;
2391 my_vec_max = (Standard_Real *) &aVecV ;
2392 my_vec_min_min = (Standard_Real *) &aVecUU ;
2393 my_vec_min_max = (Standard_Real *) &aVecUV ;
2394 my_vec_max_max = (Standard_Real *) &aVecVV ;
2395 }
2396 else {
2397 min_degree = VDegree ;
2398 max_degree = UDegree ;
2399 inverse_min = 1.0e0/VSpanLenght ;
2400 inverse_max = 1.0e0/USpanLenght ;
2401 new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ;
2402 new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ;
2403 dimension = 3 * (VDegree + 1) ;
2404 my_vec_min = (Standard_Real *) &aVecV ;
2405 my_vec_max = (Standard_Real *) &aVecU ;
2406 my_vec_min_min = (Standard_Real *) &aVecVV ;
2407 my_vec_min_max = (Standard_Real *) &aVecUV ;
2408 my_vec_max_max = (Standard_Real *) &aVecUU ;
2409 }
2410
2411 BSplSLib_LocalArray locpoles (3 * dimension);
2412
2413 //
2414 // initialize in case min or max degree are less than 2
2415 //
2416 Standard_Integer MinIndMax = 2;
2417 if ( max_degree < 2) MinIndMax = max_degree;
2418 Standard_Integer MinIndMin = 2;
2419 if ( min_degree < 2) MinIndMin = min_degree;
2420
2421 index = MinIndMax * dimension ;
2422
2423 for (ii = MinIndMax ; ii < 3 ; ii++) {
2424
2425 for (kk = 0 ; kk < dimension ; kk++) {
2426 locpoles[index] = 0.0e0 ;
2427 index += 1 ;
2428 }
2429 }
2430
2431 PLib::EvalPolynomial(new_parameter[0],
2432 MinIndMax,
2433 max_degree,
2434 dimension,
2435 PArray[0],
2436 locpoles[0]) ;
2437
2438 PLib::EvalPolynomial(new_parameter[1],
2439 MinIndMin,
2440 min_degree,
2441 3,
2442 locpoles[0],
2443 local_poles_array[0][0][0]) ;
2444 PLib::EvalPolynomial(new_parameter[1],
2445 1,
2446 min_degree,
2447 3,
2448 locpoles[dimension],
2449 local_poles_array[1][0][0]) ;
2450
2451 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2452 min_degree,
2453 3,
2454 (min_degree << 1) + min_degree,
2455 locpoles[dimension + dimension],
2456 local_poles_array[2][0][0]) ;
2457
2458 if (&WeightsArray != NULL) {
2459 dimension = min_degree + 1 ;
2460 Standard_Real *
2461 WArray = (Standard_Real *)
2462 &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ;
2463 PLib::EvalPolynomial(new_parameter[0],
2464 MinIndMax,
2465 max_degree,
2466 dimension,
2467 WArray[0],
2468 locpoles[0]) ;
2469
2470 PLib::EvalPolynomial(new_parameter[1],
2471 MinIndMin,
2472 min_degree,
2473 1,
2474 locpoles[0],
2475 local_weights_array[0][0]) ;
2476 PLib::EvalPolynomial(new_parameter[1],
2477 1,
2478 min_degree,
2479 1,
2480 locpoles[dimension],
2481 local_weights_array[1][0]) ;
2482 PLib::NoDerivativeEvalPolynomial(new_parameter[1],
2483 min_degree,
2484 1,
2485 min_degree,
2486 locpoles[dimension + dimension],
2487 local_weights_array[2][0]) ;
2488
2489
2490 local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0];
2491 local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1];
2492 local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2];
2493 local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0];
2494 local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1];
2495 local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2];
2496 local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0];
2497 local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1];
2498 local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2];
2499
2500 local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0];
2501 local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1];
2502 local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2];
2503 local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0];
2504 local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1];
2505 local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2];
2506 local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0];
2507 local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1];
2508 local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2];
2509
2510 local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0];
2511 local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1];
2512 local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2];
2513 local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0];
2514 local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1];
2515 local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2];
2516 local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0];
2517 local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1];
2518 local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2];
2519
2520
2521 local_poles_and_weights_array[0][0][3] = local_weights_array[0][0];
2522 local_poles_and_weights_array[0][1][3] = local_weights_array[0][1];
2523 local_poles_and_weights_array[0][2][3] = local_weights_array[0][2];
2524 local_poles_and_weights_array[1][0][3] = local_weights_array[1][0];
2525 local_poles_and_weights_array[1][1][3] = local_weights_array[1][1];
2526 local_poles_and_weights_array[1][2][3] = local_weights_array[1][2];
2527 local_poles_and_weights_array[2][0][3] = local_weights_array[2][0];
2528 local_poles_and_weights_array[2][1][3] = local_weights_array[2][1];
2529 local_poles_and_weights_array[2][2][3] = local_weights_array[2][2];
2530
2531 BSplSLib::RationalDerivative(2,
2532 2,
2533 2,
2534 2,
2535 local_poles_and_weights_array[0][0][0],
2536 local_poles_array[0][0][0]) ;
2537 }
2538
2539
2540 Standard_Real minmin = inverse_min * inverse_min;
2541 Standard_Real minmax = inverse_min * inverse_max;
2542 Standard_Real maxmax = inverse_max * inverse_max;
2543
2544 my_point [0] = local_poles_array [0][0][0] ;
2545 my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ;
2546 my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ;
2547 my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ;
2548 my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ;
2549 my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ;
2550
2551 my_point [1] = local_poles_array [0][0][1] ;
2552 my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ;
2553 my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ;
2554 my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ;
2555 my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ;
2556 my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ;
2557
2558 my_point [2] = local_poles_array [0][0][2] ;
2559 my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ;
2560 my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ;
2561 my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ;
2562 my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ;
2563 my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ;
2564}
2565
2566//=======================================================================
2567//function : MovePoint
2568//purpose : Find the new poles which allows an old point (with a
2569// given u and v as parameters) to reach a new position
2570//=======================================================================
2571
2572void BSplSLib::MovePoint (const Standard_Real U,
2573 const Standard_Real V,
2574 const gp_Vec& Displ,
2575 const Standard_Integer UIndex1,
2576 const Standard_Integer UIndex2,
2577 const Standard_Integer VIndex1,
2578 const Standard_Integer VIndex2,
2579 const Standard_Integer UDegree,
2580 const Standard_Integer VDegree,
2581 const Standard_Boolean Rational,
2582 const TColgp_Array2OfPnt& Poles,
2583 const TColStd_Array2OfReal& Weights,
2584 const TColStd_Array1OfReal& UFlatKnots,
2585 const TColStd_Array1OfReal& VFlatKnots,
2586 Standard_Integer& UFirstIndex,
2587 Standard_Integer& ULastIndex,
2588 Standard_Integer& VFirstIndex,
2589 Standard_Integer& VLastIndex,
2590 TColgp_Array2OfPnt& NewPoles)
2591{
2592 // calculate the UBSplineBasis in the parameter U
2593 Standard_Integer UFirstNonZeroBsplineIndex;
2594 math_Matrix UBSplineBasis(1, 1,
2595 1, UDegree+1);
2596 Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1,
2597 0,
2598 UDegree+1,
2599 UFlatKnots,
2600 U,
2601 UFirstNonZeroBsplineIndex,
2602 UBSplineBasis);
2603 // calculate the VBSplineBasis in the parameter V
2604 Standard_Integer VFirstNonZeroBsplineIndex;
2605 math_Matrix VBSplineBasis(1, 1,
2606 1, VDegree+1);
2607 Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1,
2608 0,
2609 VDegree+1,
2610 VFlatKnots,
2611 V,
2612 VFirstNonZeroBsplineIndex,
2613 VBSplineBasis);
2614 if (ErrorCod1 || ErrorCod2) {
2615 UFirstIndex = 0;
2616 ULastIndex = 0;
2617 VFirstIndex = 0;
2618 VLastIndex = 0;
2619 return;
2620 }
2621
2622 // find the span which is predominant for parameter U
2623 UFirstIndex = UFirstNonZeroBsplineIndex;
2624 ULastIndex = UFirstNonZeroBsplineIndex + UDegree ;
2625 if (UFirstIndex < UIndex1) UFirstIndex = UIndex1;
2626 if (ULastIndex > UIndex2) ULastIndex = UIndex2;
2627
2628 Standard_Real maxValue = 0.0;
2629 Standard_Integer i, ukk1=0, ukk2;
2630
2631 for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) {
2632 if (UBSplineBasis(1,i) > maxValue) {
2633 ukk1 = i + UFirstNonZeroBsplineIndex - 1;
2634 maxValue = UBSplineBasis(1,i);
2635 }
2636 }
2637
2638 // find a ukk2 if symetriy
2639 ukk2 = ukk1;
2640 i = ukk1 - UFirstNonZeroBsplineIndex + 2;
2641 if ((ukk1+1) <= ULastIndex) {
2642 if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2643 ukk2 = ukk1+1;
2644 }
2645 }
2646
2647 // find the span which is predominant for parameter V
2648 VFirstIndex = VFirstNonZeroBsplineIndex;
2649 VLastIndex = VFirstNonZeroBsplineIndex + VDegree ;
2650
2651 if (VFirstIndex < VIndex1) VFirstIndex = VIndex1;
2652 if (VLastIndex > VIndex2) VLastIndex = VIndex2;
2653
2654 maxValue = 0.0;
2655 Standard_Integer j, vkk1=0, vkk2;
2656
2657 for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) {
2658 if (VBSplineBasis(1,j) > maxValue) {
2659 vkk1 = j + VFirstNonZeroBsplineIndex - 1;
2660 maxValue = VBSplineBasis(1,j);
2661 }
2662 }
2663
2664 // find a vkk2 if symetriy
2665 vkk2 = vkk1;
2666 j = vkk1 - VFirstNonZeroBsplineIndex + 2;
2667 if ((vkk1+1) <= VLastIndex) {
2668 if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
2669 vkk2 = vkk1+1;
2670 }
2671 }
2672
2673 // compute the vector of displacement
2674 Standard_Real D1 = 0.0;
2675 Standard_Real D2 = 0.0;
2676 Standard_Real hN, Coef, DvalU, DvalV;
2677
2678 Standard_Integer ii, jj;
2679
2680 for (i = 1; i <= UDegree+1; i++) {
2681 ii = i + UFirstNonZeroBsplineIndex - 1;
2682 if (ii < ukk1) {
2683 DvalU = ukk1-ii;
2684 }
2685 else if (ii > ukk2) {
2686 DvalU = ii - ukk2;
2687 }
2688 else {
2689 DvalU = 0.0;
2690 }
2691
2692 for (j = 1; j <= VDegree+1; j++) {
2693 jj = j + VFirstNonZeroBsplineIndex - 1;
2694 if (Rational) {
2695 hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j);
2696 D2 += hN;
2697 }
2698 else {
2699 hN = UBSplineBasis(1, i)*VBSplineBasis(1,j);
2700 }
2701 if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) {
2702 if (jj < vkk1) {
2703 DvalV = vkk1-jj;
2704 }
2705 else if (jj > vkk2) {
2706 DvalV = jj - vkk2;
2707 }
2708 else {
2709 DvalV = 0.0;
2710 }
2711 D1 += 1./(DvalU + DvalV + 1.) * hN;
2712 }
2713 }
2714 }
2715
2716 if (Rational) {
2717 Coef = D2/D1;
2718 }
2719 else {
2720 Coef = 1./D1;
2721 }
2722
2723 // compute the new poles
2724
2725 for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) {
2726 if (i < ukk1) {
2727 DvalU = ukk1-i;
2728 }
2729 else if (i > ukk2) {
2730 DvalU = i - ukk2;
2731 }
2732 else {
2733 DvalU = 0.0;
2734 }
2735
2736 for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) {
2737 if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) {
2738 if (j < vkk1) {
2739 DvalV = vkk1-j;
2740 }
2741 else if (j > vkk2) {
2742 DvalV = j - vkk2;
2743 }
2744 else {
2745 DvalV = 0.0;
2746 }
2747 NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ);
2748 }
2749 else {
2750 NewPoles(i,j) = Poles(i,j);
2751 }
2752 }
2753 }
2754}
2755
2756//=======================================================================
0d969553
Y		 2757// function : Resolution
2758// purpose : this computes an estimate for the maximum of the
7fd59977 2759// partial derivatives both in U and in V
2760//
2761//
0d969553
Y		 2762// The calculation resembles at the calculation of curves with
2763// additional index for the control point. Let Si,j be the
2764// control points for ls surface and Di,j the weights.
2765// The checking of upper bounds for the partial derivatives
2766// will be omitted and Su is the next upper bound in the polynomial case :
7fd59977 2767//
2768//
2769//
2770// | Si,j - Si-1,j |
2771// d * Max | ------------- |
2772// i = 2,n | ti+d - ti |
2773// i=1.m
2774//
2775//
0d969553 2776// and in the rational case :
7fd59977 2777//
2778//
2779//
2780// Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j)
2781// Max Max d * -----------------------------------------------
2782// k=1,n i dans Rj ti+d - ti
2783// j=1,m
2784// ----------------------------------------------------------------------
2785//
2786// Min Di,j
2787// i=1,n
2788// j=1,m
2789//
2790//
2791//
0d969553 2792// with Rj = {j-d, ...., j+d+d+1}.
7fd59977 2793//
2794//
2795//=======================================================================
2796
2797void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles,
2798 const TColStd_Array2OfReal& Weights,
2799 const TColStd_Array1OfReal& UKnots,
2800 const TColStd_Array1OfReal& VKnots,
2801 const TColStd_Array1OfInteger& UMults,
2802 const TColStd_Array1OfInteger& VMults,
2803 const Standard_Integer UDegree,
2804 const Standard_Integer VDegree,
2805 const Standard_Boolean URational,
2806 const Standard_Boolean VRational,
2807 const Standard_Boolean UPeriodic,
2808 const Standard_Boolean VPeriodic,
2809 const Standard_Real Tolerance3D,
2810 Standard_Real& UTolerance,
2811 Standard_Real& VTolerance)
2812{
2813 Standard_Real Wij,Wmj,Wji,Wjm;
2814 Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp;
2815 Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp;
2816 Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp;
2817 Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2];
2818
2819 max_derivative[0] = max_derivative[1] = 0.0e0 ;
2820
2821 Standard_Integer PRowLength, PColLength;
2822 Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index;
2823 Standard_Integer ii_minus,upper[2],lower[2],poles_length[2];
2824 Standard_Integer num_poles[2],num_flat_knots[2];
2825
2826 num_flat_knots[0] =
2827 BSplCLib::KnotSequenceLength(UMults,
2828 UDegree,
2829 UPeriodic) ;
2830 num_flat_knots[1] =
2831 BSplCLib::KnotSequenceLength(VMults,
2832 VDegree,
2833 VPeriodic) ;
2834 TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ;
2835 TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ;
2836 BSplCLib::KnotSequence(UKnots,
2837 UMults,
2838 UDegree,
2839 UPeriodic,
2840 flat_knots_in_u) ;
2841 BSplCLib::KnotSequence(VKnots,
2842 VMults,
2843 VDegree,
2844 VPeriodic,
2845 flat_knots_in_v) ;
2846 PRowLength = Poles.RowLength();
2847 PColLength = Poles.ColLength();
2848 if (URational || VRational) {
2849 Standard_Integer Wsize = PRowLength * PColLength;
2850 const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol());
2851 min_weights = WG[0];
2852
2853 for (ii = 1 ; ii < Wsize ; ii++) {
2854 min = WG[ii];
2855 if (min_weights > min) min_weights = min;
2856 }
2857 }
2858 Standard_Integer UD1 = UDegree + 1;
2859 Standard_Integer VD1 = VDegree + 1;
2860 num_poles[0] = num_flat_knots[0] - UD1;
2861 num_poles[1] = num_flat_knots[1] - VD1;
2862 poles_length[0] = PColLength;
2863 poles_length[1] = PRowLength;
2864 if(URational) {
2865 Standard_Integer UD2 = UDegree << 1;
2866 Standard_Integer VD2 = VDegree << 1;
2867
2868 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2869 ii_index = (ii - 1) % poles_length[0] + 1 ;
2870 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2871 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2872 inverse = 1.0e0 / inverse ;
2873 lower[0] = ii - UD1;
2874 if (lower[0] < 1) lower[0] = 1;
2875 upper[0] = ii + UD2 + 1;
2876 if (upper[0] > num_poles[0]) upper[0] = num_poles[0];
2877
2878 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2879 jj_index = (jj - 1) % poles_length[1] + 1 ;
2880 lower[1] = jj - VD1;
2881 if (lower[1] < 1) lower[1] = 1;
2882 upper[1] = jj + VD2 + 1;
2883 if (upper[1] > num_poles[1]) upper[1] = num_poles[1];
2884 const gp_Pnt& Pij = Poles .Value(ii_index,jj_index);
2885 Wij = Weights.Value(ii_index,jj_index);
2886 const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index);
2887 Wmj = Weights.Value(ii_minus,jj_index);
2888 Xij = Pij.X();
2889 Yij = Pij.Y();
2890 Zij = Pij.Z();
2891 Xmj = Pmj.X();
2892 Ymj = Pmj.Y();
2893 Zmj = Pmj.Z();
2894
2895 for (pp = lower[0] ; pp <= upper[0] ; pp++) {
2896 pp_index = (pp - 1) % poles_length[0] + 1 ;
2897
2898 for (qq = lower[1] ; qq <= upper[1] ; qq++) {
2899 value = 0.0e0 ;
2900 qq_index = (qq - 1) % poles_length[1] + 1 ;
2901 const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index);
2902 Xpq = Ppq.X();
2903 Ypq = Ppq.Y();
2904 Zpq = Ppq.Z();
2905 factor = (Xpq - Xij) * Wij;
2906 factor -= (Xpq - Xmj) * Wmj;
2907 if (factor < 0) factor = - factor;
2908 value += factor ;
2909 factor = (Ypq - Yij) * Wij;
2910 factor -= (Ypq - Ymj) * Wmj;
2911 if (factor < 0) factor = - factor;
2912 value += factor ;
2913 factor = (Zpq - Zij) * Wij;
2914 factor -= (Zpq - Zmj) * Wmj;
2915 if (factor < 0) factor = - factor;
2916 value += factor ;
2917 value *= inverse ;
2918 if (max_derivative[0] < value) max_derivative[0] = value ;
2919 }
2920 }
2921 }
2922 }
2923 max_derivative[0] /= min_weights ;
2924 }
2925 else {
2926
2927 for (ii = 2 ; ii <= num_poles[0] ; ii++) {
2928 ii_index = (ii - 1) % poles_length[0] + 1 ;
2929 ii_minus = (ii - 2) % poles_length[0] + 1 ;
2930 inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ;
2931 inverse = 1.0e0 / inverse ;
2932
2933 for ( jj = 1 ; jj <= num_poles[1] ; jj++) {
2934 jj_index = (jj - 1) % poles_length[1] + 1 ;
2935 value = 0.0e0 ;
2936 const gp_Pnt& Pij = Poles.Value(ii_index,jj_index);
2937 const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index);
2938 factor = Pij.X() - Pmj.X();
2939 if (factor < 0) factor = - factor;
2940 value += factor;
2941 factor = Pij.Y() - Pmj.Y();
2942 if (factor < 0) factor = - factor;
2943 value += factor;
2944 factor = Pij.Z() - Pmj.Z();
2945 if (factor < 0) factor = - factor;
2946 value += factor;
2947 value *= inverse ;
2948 if (max_derivative[0] < value) max_derivative[0] = value ;
2949 }
2950 }
2951 }
2952 max_derivative[0] *= UDegree ;
2953 if(VRational) {
2954 Standard_Integer UD2 = UDegree << 1;
2955 Standard_Integer VD2 = VDegree << 1;
2956
2957 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
2958 ii_index = (ii - 1) % poles_length[1] + 1 ;
2959 ii_minus = (ii - 2) % poles_length[1] + 1 ;
2960 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
2961 inverse = 1.0e0 / inverse ;
2962 lower[0] = ii - VD1;
2963 if (lower[0] < 1) lower[0] = 1;
2964 upper[0] = ii + VD2 + 1;
2965 if (upper[0] > num_poles[1]) upper[0] = num_poles[1];
2966
2967 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
2968 jj_index = (jj - 1) % poles_length[0] + 1 ;
2969 lower[1] = jj - UD1;
2970 if (lower[1] < 1) lower[1] = 1;
2971 upper[1] = jj + UD2 + 1;
2972 if (upper[1] > num_poles[0]) upper[1] = num_poles[0];
2973 const gp_Pnt& Pji = Poles .Value(jj_index,ii_index);
2974 Wji = Weights.Value(jj_index,ii_index);
2975 const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus);
2976 Wjm = Weights.Value(jj_index,ii_minus);
2977 Xji = Pji.X();
2978 Yji = Pji.Y();
2979 Zji = Pji.Z();
2980 Xjm = Pjm.X();
2981 Yjm = Pjm.Y();
2982 Zjm = Pjm.Z();
2983
2984 for (pp = lower[1] ; pp <= upper[1] ; pp++) {
2985 pp_index = (pp - 1) % poles_length[1] + 1 ;
2986
2987 for (qq = lower[0] ; qq <= upper[0] ; qq++) {
2988 value = 0.0e0 ;
2989 qq_index = (qq - 1) % poles_length[0] + 1 ;
2990 const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index);
2991 Xqp = Pqp.X();
2992 Yqp = Pqp.Y();
2993 Zqp = Pqp.Z();
2994 factor = (Xqp - Xji) * Wji;
2995 factor -= (Xqp - Xjm) * Wjm;
2996 if (factor < 0) factor = - factor;
2997 value += factor ;
2998 factor = (Yqp - Yji) * Wji;
2999 factor -= (Yqp - Yjm) * Wjm;
3000 if (factor < 0) factor = - factor;
3001 value += factor ;
3002 factor = (Zqp - Zji) * Wji;
3003 factor -= (Zqp - Zjm) * Wjm;
3004 if (factor < 0) factor = - factor;
3005 value += factor ;
3006 value *= inverse ;
3007 if (max_derivative[1] < value) max_derivative[1] = value ;
3008 }
3009 }
3010 }
3011 }
3012 max_derivative[1] /= min_weights ;
3013 }
3014 else {
3015
3016 for (ii = 2 ; ii <= num_poles[1] ; ii++) {
3017 ii_index = (ii - 1) % poles_length[1] + 1 ;
3018 ii_minus = (ii - 2) % poles_length[1] + 1 ;
3019 inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ;
3020 inverse = 1.0e0 / inverse ;
3021
3022 for ( jj = 1 ; jj <= num_poles[0] ; jj++) {
3023 jj_index = (jj - 1) % poles_length[0] + 1 ;
3024 value = 0.0e0 ;
3025 const gp_Pnt& Pji = Poles.Value(jj_index,ii_index);
3026 const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus);
3027 factor = Pji.X() - Pjm.X() ;
3028 if (factor < 0) factor = - factor;
3029 value += factor;
3030 factor = Pji.Y() - Pjm.Y() ;
3031 if (factor < 0) factor = - factor;
3032 value += factor;
3033 factor = Pji.Z() - Pjm.Z() ;
3034 if (factor < 0) factor = - factor;
3035 value += factor;
3036 value *= inverse ;
3037 if (max_derivative[1] < value) max_derivative[1] = value ;
3038 }
3039 }
3040 }
3041 max_derivative[1] *= VDegree ;
3042 max_derivative[0] *= M_SQRT2 ;
3043 max_derivative[1] *= M_SQRT2 ;
3044 if(max_derivative[0] && max_derivative[1]) {
3045 UTolerance = Tolerance3D / max_derivative[0] ;
3046 VTolerance = Tolerance3D / max_derivative[1] ;
3047 }
3048 else {
3049 UTolerance=VTolerance=0.0;
3050#ifdef DEB
3051 cout<<"ElSLib.cxx : maxderivative = 0.0 "<<endl;
3052#endif
3053 }
3054}
3055
3056//=======================================================================
3057//function : Interpolate
3058//purpose :
3059//=======================================================================
3060
3061void BSplSLib::Interpolate(const Standard_Integer UDegree,
3062 const Standard_Integer VDegree,
3063 const TColStd_Array1OfReal& UFlatKnots,
3064 const TColStd_Array1OfReal& VFlatKnots,
3065 const TColStd_Array1OfReal& UParameters,
3066 const TColStd_Array1OfReal& VParameters,
3067 TColgp_Array2OfPnt& Poles,
3068 TColStd_Array2OfReal& Weights,
3069 Standard_Integer& InversionProblem)
3070{
3071 Standard_Integer ii, jj, ll, kk, dimension;
3072 Standard_Integer ULength = UParameters.Length();
3073 Standard_Integer VLength = VParameters.Length();
3074 Standard_Real * poles_array;
3075
0d969553 3076 // extraction of iso u
7fd59977 3077 dimension = 4*ULength;
3078 TColStd_Array2OfReal Points(1, VLength,
3079 1, dimension);
3080
3081 Handle(TColStd_HArray1OfInteger) ContactOrder =
3082 new (TColStd_HArray1OfInteger)(1, VLength);
3083 ContactOrder->Init(0);
3084
3085 for (ii=1; ii <= VLength; ii++) {
3086
3087 for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) {
3088 Points(ii,ll) = Poles(jj, ii).X();
3089 Points(ii,ll+1) = Poles(jj, ii).Y();
3090 Points(ii,ll+2) = Poles(jj, ii).Z();
3091 Points(ii,ll+3) = Weights(jj,ii) ;
3092 }
3093 }
3094
0d969553 3095 // interpolation of iso u
7fd59977 3096 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3097 BSplCLib::Interpolate(VDegree,
3098 VFlatKnots,
3099 VParameters,
3100 ContactOrder->Array1(),
3101 dimension,
3102 poles_array[0],
3103 InversionProblem) ;
3104 if (InversionProblem != 0) return;
3105
0d969553 3106 // extraction of iso v
7fd59977 3107
3108 dimension = VLength*4;
3109 TColStd_Array2OfReal IsoPoles(1, ULength,
3110 1, dimension);
3111
3112 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3113 ContactOrder->Init(0);
3114 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3115
3116 for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) {
3117
3118 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3119 IsoPoles (ii,ll) = Points(jj, kk);
3120 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3121 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3122 IsoPoles (ii,ll+3) = Points(jj, kk+3);
3123 }
3124 }
0d969553 3125 // interpolation of iso v
7fd59977 3126 BSplCLib::Interpolate(UDegree,
3127 UFlatKnots,
3128 UParameters,
3129 ContactOrder->Array1(),
3130 dimension,
3131 poles_array[0],
3132 InversionProblem);
3133
0d969553 3134 // return results
7fd59977 3135
3136 for (ii=1; ii <= ULength; ii++) {
3137
3138 for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) {
3139 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3140 Poles.SetValue(ii, jj, Pnt);
3141 Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ;
3142 }
3143 }
3144}
3145
3146//=======================================================================
3147//function : Interpolate
3148//purpose :
3149//=======================================================================
3150
3151void BSplSLib::Interpolate(const Standard_Integer UDegree,
3152 const Standard_Integer VDegree,
3153 const TColStd_Array1OfReal& UFlatKnots,
3154 const TColStd_Array1OfReal& VFlatKnots,
3155 const TColStd_Array1OfReal& UParameters,
3156 const TColStd_Array1OfReal& VParameters,
3157 TColgp_Array2OfPnt& Poles,
3158 Standard_Integer& InversionProblem)
3159{
3160 Standard_Integer ii, jj, ll, kk, dimension;
3161 Standard_Integer ULength = UParameters.Length();
3162 Standard_Integer VLength = VParameters.Length();
3163 Standard_Real * poles_array;
3164
0d969553 3165 // extraction of iso u
7fd59977 3166 dimension = 3*ULength;
3167 TColStd_Array2OfReal Points(1, VLength,
3168 1, dimension);
3169
3170 Handle(TColStd_HArray1OfInteger) ContactOrder =
3171 new (TColStd_HArray1OfInteger)(1, VLength);
3172 ContactOrder->Init(0);
3173
3174 for (ii=1; ii <= VLength; ii++) {
3175
3176 for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) {
3177 Points(ii,ll) = Poles(jj, ii).X();
3178 Points(ii,ll+1) = Poles(jj, ii).Y();
3179 Points(ii,ll+2) = Poles(jj, ii).Z();
3180 }
3181 }
3182
0d969553 3183 // interpolation of iso u
7fd59977 3184 poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ;
3185 BSplCLib::Interpolate(VDegree,
3186 VFlatKnots,
3187 VParameters,
3188 ContactOrder->Array1(),
3189 dimension,
3190 poles_array[0],
3191 InversionProblem) ;
3192 if (InversionProblem != 0) return;
3193
0d969553 3194 // extraction of iso v
7fd59977 3195
3196 dimension = VLength*3;
3197 TColStd_Array2OfReal IsoPoles(1, ULength,
3198 1, dimension);
3199
3200 ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength);
3201 ContactOrder->Init(0);
3202 poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ;
3203
3204 for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) {
3205
3206 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3207 IsoPoles (ii,ll) = Points(jj, kk);
3208 IsoPoles (ii,ll+1) = Points(jj, kk+1);
3209 IsoPoles (ii,ll+2) = Points(jj, kk+2);
3210 }
3211 }
0d969553 3212 // interpolation of iso v
7fd59977 3213 BSplCLib::Interpolate(UDegree,
3214 UFlatKnots,
3215 UParameters,
3216 ContactOrder->Array1(),
3217 dimension,
3218 poles_array[0],
3219 InversionProblem);
3220
0d969553 3221 // return results
7fd59977 3222
3223 for (ii=1; ii <= ULength; ii++) {
3224
3225 for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) {
3226 gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2));
3227 Poles.SetValue(ii, jj, Pnt);
3228 }
3229 }
3230}
3231
3232//=======================================================================
3233//function : FunctionMultiply
3234//purpose :
3235//=======================================================================
3236
3237void BSplSLib::FunctionMultiply
3238(const BSplSLib_EvaluatorFunction& Function,
3239 const Standard_Integer UBSplineDegree,
3240 const Standard_Integer VBSplineDegree,
3241 const TColStd_Array1OfReal& UBSplineKnots,
3242 const TColStd_Array1OfReal& VBSplineKnots,
3243 const TColStd_Array1OfInteger & UMults,
3244 const TColStd_Array1OfInteger & VMults,
3245 const TColgp_Array2OfPnt& Poles,
3246 const TColStd_Array2OfReal& Weights,
3247 const TColStd_Array1OfReal& UFlatKnots,
3248 const TColStd_Array1OfReal& VFlatKnots,
3249 const Standard_Integer UNewDegree,
3250 const Standard_Integer VNewDegree,
3251 TColgp_Array2OfPnt& NewNumerator,
3252 TColStd_Array2OfReal& NewDenominator,
3253 Standard_Integer& Status)
3254{
3255 Standard_Integer num_uparameters,
3256// ii,jj,kk,
3257 ii,jj,
3258 error_code,
3259 num_vparameters ;
3260 Standard_Real result ;
3261
3262 num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ;
3263 num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ;
3264 TColStd_Array1OfReal UParameters(1,num_uparameters) ;
3265 TColStd_Array1OfReal VParameters(1,num_vparameters) ;
3266
3267 if ((NewNumerator.ColLength() == num_uparameters) &&
3268 (NewNumerator.RowLength() == num_vparameters) &&
3269 (NewDenominator.ColLength() == num_uparameters) &&
3270 (NewDenominator.RowLength() == num_vparameters)) {
3271
3272
3273 BSplCLib::BuildSchoenbergPoints(UNewDegree,
3274 UFlatKnots,
3275 UParameters) ;
3276
3277 BSplCLib::BuildSchoenbergPoints(VNewDegree,
3278 VFlatKnots,
3279 VParameters) ;
3280
3281 for (ii = 1 ; ii <= num_uparameters ; ii++) {
3282
3283 for (jj = 1 ; jj <= num_vparameters ; jj++) {
3284 HomogeneousD0(UParameters(ii),
3285 VParameters(jj),
3286 0,
3287 0,
3288 Poles,
3289 Weights,
3290 UBSplineKnots,
3291 VBSplineKnots,
3292 UMults,
3293 VMults,
3294 UBSplineDegree,
3295 VBSplineDegree,
3296 Standard_True,
3297 Standard_True,
3298 Standard_False,
3299 Standard_False,
3300 NewDenominator(ii,jj),
3301 NewNumerator(ii,jj)) ;
3302
41194117 3303 Function.Evaluate (0,
7fd59977 3304 UParameters(ii),
3305 VParameters(jj),
3306 result,
3307 error_code) ;
3308 if (error_code) {
3309 Standard_ConstructionError::Raise();
3310 }
3311 gp_Pnt& P = NewNumerator(ii,jj);
3312 P.SetX(P.X() * result);
3313 P.SetY(P.Y() * result);
3314 P.SetZ(P.Z() * result);
3315 NewDenominator(ii,jj) *= result ;
3316 }
3317 }
3318 Interpolate(UNewDegree,
3319 VNewDegree,
3320 UFlatKnots,
3321 VFlatKnots,
3322 UParameters,
3323 VParameters,
3324 NewNumerator,
3325 NewDenominator,
3326 Status) ;
3327 }
3328 else {
3329 Standard_ConstructionError::Raise();
3330 }
3331}
Open CASCADE Technology repository