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