0024624: Lost word in license statement in source files
[occt.git] / src / AppParCurves / AppParCurves.cxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #define No_Standard_RangeError
16 #define No_Standard_OutOfRange
17
18 #include <AppParCurves.ixx>
19 #include <BSplCLib.hxx>
20 #include <TColStd_Array1OfReal.hxx>
21 #include <gp_Pnt2d.hxx>
22 #include <gp_Vec2d.hxx>
23
24
25 void AppParCurves::BernsteinMatrix(const Standard_Integer NbPoles, 
26                                    const math_Vector& U,
27                                    math_Matrix& A) {
28
29   Standard_Integer i, j, id;
30   Standard_Real u0, u1, y0, y1, xs;
31   Standard_Integer first = U.Lower(), last = U.Upper();
32   math_Vector B(1, NbPoles-1);
33
34
35   for (i = first; i <= last; i++) {
36     B(1) = 1;
37     u0 = U(i);
38     u1 = 1.-u0;
39     
40     for (id = 2; id <= NbPoles-1; id++) {
41       y0 = B(1);
42       y1 = u0*y0;
43       B(1) = y0-y1;
44       for (j = 2; j <= id-1; j++) {
45         xs = y1;
46         y0 = B(j);
47         y1 = u0*y0;
48         B(j) = y0-y1+xs;
49       }
50       B(id) = y1;
51     }
52     A(i, 1) = u1*B(1);
53     A(i, NbPoles) = u0*B(NbPoles-1);
54     for (j = 2; j <= NbPoles-1; j++) {
55       A(i, j) = u1*B(j)+u0*B(j-1);
56     }
57   }
58 }
59
60
61 void AppParCurves::Bernstein(const Standard_Integer NbPoles, 
62                              const math_Vector& U,
63                              math_Matrix& A,
64                              math_Matrix& DA) {
65   
66   Standard_Integer i, j, id, Ndeg = NbPoles-1;
67   Standard_Real u0, u1, y0, y1, xs, bj, bj1;;
68   Standard_Integer first = U.Lower(), last = U.Upper();
69   math_Vector B(1, NbPoles-1);
70
71   
72   for (i = first; i <= last; i++) {
73     B(1) = 1;
74     u0 = U(i);
75     u1 = 1.-u0;
76     
77     for (id = 2; id <= NbPoles-1; id++) {
78       y0 = B(1);
79       y1 = u0*y0;
80       B(1) = y0-y1;
81       for (j = 2; j <= id-1; j++) {
82         xs = y1;
83         y0 = B(j);
84         y1 = u0*y0;
85         B(j) = y0-y1+xs;
86       }
87       B(id) = y1;
88     }
89     DA(i, 1) = -Ndeg*B(1);
90     DA(i, NbPoles) = Ndeg*B(NbPoles-1);
91     A(i, 1) = u1*B(1);
92     A(i, NbPoles) = u0*B(NbPoles-1);
93     for (j = 2; j <= NbPoles-1; j++) {
94       bj = B(j);  bj1 = B(j-1);
95       DA(i,j) = Ndeg*(bj1-bj);
96       A(i, j) = u1*bj+u0*bj1;
97     }
98   }
99 }
100
101
102 void AppParCurves::SecondDerivativeBernstein(const Standard_Real U,
103                                              math_Vector&        DDA) {
104 //  Standard_Real U1 = 1-U, Y0, Y1, Xs;
105   Standard_Real Y0, Y1, Xs;
106   Standard_Integer NbPoles = DDA.Length();
107   Standard_Integer id, j, N4, deg = NbPoles-1;
108   N4 = deg*(deg-1);
109   math_Vector B(1, deg-1);
110   B(1) = 1.;
111   
112   // Cas particulier si degre = 1:
113   if (deg == 1) {
114     DDA(1) = 0.0;
115     DDA(2) = 0.0;
116   }
117   else if (deg == 2) {
118     DDA(1) = 2.0;
119     DDA(2) = -4.0;
120     DDA(3) = 2.0;
121   }
122   else {
123
124     for (id = 2; id <= deg-1; id++) {
125       Y0 = B(1);
126       Y1 = U*Y0;
127       B(1) = Y0-Y1;
128       for (j = 2; j <= id-1; j++) {
129         Xs = Y1;
130         Y0 = B(j);
131         Y1 = U*Y0;
132         B(j) = Y0 - Y1 +Xs;
133       }
134       B(id) = Y1;
135     }
136     
137     DDA(1) = N4*B(1);
138     DDA(2) = N4*(-2*B(1) + B(2));
139     DDA(deg) = N4*(B(deg-2) - 2*B(deg-1));
140     DDA(deg+1) = N4*B(deg-1);
141     
142     for(j = 2; j <= deg-2; j++) {
143       DDA(j+1) = N4*(B(j-1)-2*B(j)+B(j+1));
144     }
145   }
146 }
147
148
149
150 void AppParCurves::SplineFunction(const Standard_Integer nbpoles,
151                                   const Standard_Integer deg,
152                                   const math_Vector&     Parameters,
153                                   const math_Vector&     flatknots,
154                                   math_Matrix&           A,
155                                   math_Matrix&           DA,
156                                   math_IntegerVector&    index)
157 {
158 //  Standard_Real U, NewU, co, diff, t1, t2;
159   Standard_Real U, NewU;
160 //  gp_Pnt2d Pt, P0;
161 //  gp_Vec2d V1;
162 //  Standard_Integer i, j, k, iter, in, ik, deg1 = deg+1;
163   Standard_Integer i, j, deg1 = deg+1;
164 //  Standard_Integer oldkindex, kindex, theindex, ttindex;
165   Standard_Integer oldkindex, kindex, theindex;
166   math_Vector locpoles(1 , deg1);
167   math_Vector locdpoles(1 , deg1);
168   Standard_Integer firstp = Parameters.Lower(), lastp = Parameters.Upper();
169
170   TColStd_Array1OfReal Aflatknots(flatknots.Lower(), flatknots.Upper());
171   for (i = flatknots.Lower(); i <= flatknots.Upper(); i++) {
172     Aflatknots(i) = flatknots(i);
173   }
174
175
176   oldkindex = 1;
177   
178   Standard_Integer pp, qq;
179   Standard_Real Saved, Inverse, LocalInverse, locqq, locdqq, val;
180
181   for (i = firstp; i <= lastp; i++) {
182     U = Parameters(i);
183     NewU = U;
184     kindex = oldkindex;
185     BSplCLib::LocateParameter(deg, Aflatknots, U, Standard_False, 
186                               deg1, nbpoles+1, kindex, NewU);
187     
188     oldkindex = kindex;
189     
190     // On stocke les index:
191     index(i) = kindex - deg-1;
192     
193     locpoles(1) = 1.0;
194     
195     for (qq = 2; qq <= deg; qq++) {
196       locpoles(qq) = 0.0;
197       for (pp = 1; pp <= qq-1; pp++) {
198         Inverse = 1.0 / (flatknots(kindex + pp)  - flatknots(kindex - qq + pp + 1)) ; 
199         Saved = (U - flatknots(kindex - qq + pp + 1)) * Inverse * locpoles(pp);
200         locpoles(pp) *= (flatknots(kindex + pp) - U) * Inverse ;
201         locpoles(pp) += locpoles(qq) ;
202         locpoles(qq) = Saved ;
203       }
204     }
205
206     qq = deg+1;
207     for (pp = 1; pp <= deg; pp++) {
208       locdpoles(pp)= locpoles(pp);
209     }
210
211     locqq = 0.0;
212     locdqq = 0.0;
213     for (pp = 1; pp <= deg; pp++) {
214       Inverse =  1.0 / (flatknots(kindex + pp)  - flatknots(kindex - qq + pp + 1)); 
215       Saved = (U - flatknots(kindex - qq + pp + 1)) * Inverse * locpoles(pp);
216       locpoles(pp) *= (flatknots(kindex + pp) - U) * Inverse;
217       locpoles(pp) += locqq;
218       locqq = Saved ;
219       LocalInverse = (Standard_Real) (deg) * Inverse;
220       Saved = LocalInverse * locdpoles(pp);
221       locdpoles(pp) *= - LocalInverse ;
222       locdpoles(pp) +=   locdqq;
223       locdqq = Saved ;
224     }
225
226     locpoles(qq) = locqq;
227     locdpoles(qq) = locdqq;
228
229     for (j = 1; j <= deg1; j++) {
230       val = locpoles(j);
231       theindex = j+oldkindex-deg1;
232       A(i, theindex) = val;
233       DA(i, theindex) = locdpoles(j);
234     }
235
236     for (j = 1; j < oldkindex-deg; j++)
237       A(i, j) = DA(i, j) = 0.0;
238     for (j = oldkindex+1; j <= nbpoles; j++)
239       A(i, j) = DA(i, j) = 0.0;
240     
241   }
242
243 }