0024002: Overall code and build procedure refactoring -- automatic
[occt.git] / src / AppParCurves / AppParCurves.cxx
CommitLineData
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
26void 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
62void 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;
68 Standard_Real u0, u1, y0, y1, xs, bj, bj1;;
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
103void 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
151void 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}