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