1 #include <2geom/d2.h>
2 /* One would think that we would include d2-sbasis.h, however,
3 * you cannot actually include it in anything - only d2 may import it.
4 * This is due to the trickinesses of template submatching. */
6 namespace Geom {
8 SBasis L2(D2<SBasis> const & a, unsigned k) { return sqrt(dot(a, a), k); }
10 D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b) {
11 return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));
12 }
14 D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b) {
15 return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));
16 }
18 D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms) {
19 return D2<SBasis>(truncate(a[X], terms), truncate(a[Y], terms));
20 }
22 unsigned sbasis_size(D2<SBasis> const & a) {
23 return std::max((unsigned) a[0].size(), (unsigned) a[1].size());
24 }
26 //TODO: Is this sensical? shouldn't it be like pythagorean or something?
27 double tail_error(D2<SBasis> const & a, unsigned tail) {
28 return std::max(a[0].tailError(tail), a[1].tailError(tail));
29 }
31 Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a) {
32 Piecewise<SBasis> x = partition(a[0], a[1].cuts), y = partition(a[1], a[0].cuts);
33 assert(x.size() == y.size());
34 Piecewise<D2<SBasis> > ret;
35 for(unsigned i = 0; i < x.size(); i++)
36 ret.push_seg(D2<SBasis>(x[i], y[i]));
37 ret.cuts.insert(ret.cuts.end(), x.cuts.begin(), x.cuts.end());
38 return ret;
39 }
41 D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a) {
42 D2<Piecewise<SBasis> > ret;
43 for(unsigned d = 0; d < 2; d++) {
44 for(unsigned i = 0; i < a.size(); i++)
45 ret[d].push_seg(a[i][d]);
46 ret[d].cuts.insert(ret[d].cuts.end(), a.cuts.begin(), a.cuts.end());
47 }
48 return ret;
49 }
51 Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &M){
52 Piecewise<D2<SBasis> > result;
53 if (M.empty()) return M;
54 result.push_cut(M.cuts[0]);
55 for (unsigned i=0; i<M.size(); i++){
56 result.push(rot90(M[i]),M.cuts[i+1]);
57 }
58 return result;
59 }
61 /** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
62 * @return \f[
63 * f(t) \rightarrow \left\{
64 * \begin{array}{c}
65 * a_1 \bullet b_1 \\
66 * a_2 \bullet b_2 \\
67 * \ldots \\
68 * a_n \bullet b_n \\
69 * \end{array}\right.
70 * \f]
71 * @relates Piecewise */
72 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b)
73 {
74 Piecewise<SBasis > result;
75 if (a.empty() || b.empty()) return result;
76 Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
77 Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
79 result.push_cut(aa.cuts.front());
80 for (unsigned i=0; i<aa.size(); i++){
81 result.push(dot(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
82 }
83 return result;
84 }
86 /** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
87 * @return \f[
88 * f(t) \rightarrow \left\{
89 * \begin{array}{c}
90 * a_1 \bullet b \\
91 * a_2 \bullet b \\
92 * \ldots \\
93 * a_n \bullet b \\
94 * \end{array}\right.
95 * \f]
96 * @relates Piecewise */
97 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Point const &b)
98 {
99 Piecewise<SBasis > result;
100 if (a.empty()) return result;
102 result.push_cut(a.cuts.front());
103 for (unsigned i = 0; i < a.size(); ++i){
104 result.push(dot(a.segs[i],b), a.cuts[i+1]);
105 }
106 return result;
107 }
110 Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a,
111 Piecewise<D2<SBasis> > const &b){
112 Piecewise<SBasis > result;
113 if (a.empty() || b.empty()) return result;
114 Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
115 Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
117 result.push_cut(aa.cuts.front());
118 for (unsigned i=0; i<a.size(); i++){
119 result.push(cross(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
120 }
121 return result;
122 }
124 Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Matrix const &m) {
125 Piecewise<D2<SBasis> > result;
126 if(a.empty()) return result;
127 result.push_cut(a.cuts[0]);
128 for (unsigned i = 0; i < a.size(); i++) {
129 result.push(a[i] * m, a.cuts[i+1]);
130 }
131 return result;
132 }
134 //if tol>0, only force continuity where the jump is smaller than tol.
135 Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double tol, bool closed)
136 {
137 if (f.size()==0) return f;
138 Piecewise<D2<SBasis> > result=f;
139 unsigned cur = (closed)? 0:1;
140 unsigned prev = (closed)? f.size()-1:0;
141 while(cur<f.size()){
142 Point pt0 = f.segs[prev].at1();
143 Point pt1 = f.segs[cur ].at0();
144 if (tol<=0 || L2sq(pt0-pt1)<tol*tol){
145 pt0 = (pt0+pt1)/2;
146 for (unsigned dim=0; dim<2; dim++){
147 SBasis &prev_sb=result.segs[prev][dim];
148 SBasis &cur_sb =result.segs[cur][dim];
149 Coord const c=pt0[dim];
150 if (prev_sb.empty()) {
151 prev_sb = SBasis(Linear(0.0, c));
152 } else {
153 prev_sb[0][1] = c;
154 }
155 if (cur_sb.empty()) {
156 cur_sb = SBasis(Linear(c, 0.0));
157 } else {
158 cur_sb[0][0] = c;
159 }
160 }
161 }
162 prev = cur++;
163 }
164 return result;
165 }
167 std::vector<Geom::Piecewise<Geom::D2<Geom::SBasis> > >
168 split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbin, double tol)
169 {
170 using namespace Geom;
171 std::vector<Piecewise<D2<SBasis> > > ret;
172 unsigned piece_start = 0;
173 for (unsigned i=0; i<pwsbin.segs.size(); i++){
174 if (i==(pwsbin.segs.size()-1) || L2(pwsbin.segs[i].at1()- pwsbin.segs[i+1].at0()) > tol){
175 Piecewise<D2<SBasis> > piece;
176 piece.cuts.push_back(pwsbin.cuts[piece_start]);
177 for (unsigned j = piece_start; j<i+1; j++){
178 piece.segs.push_back(pwsbin.segs[j]);
179 piece.cuts.push_back(pwsbin.cuts[j+1]);
180 }
181 ret.push_back(piece);
182 piece_start = i+1;
183 }
184 }
185 return ret;
186 }
188 static void set_first_point(Piecewise<D2<SBasis> > &f, Point a){
189 if ( f.empty() ){
190 f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(Linear(a[X]),Linear(a[Y]))));
191 return;
192 }
193 for (unsigned dim=0; dim<2; dim++){
194 if (f.segs.front()[dim].size() == 0){
195 f.segs.front()[dim] = SBasis(Linear(a[dim],0));
196 }else{
197 f.segs.front()[dim][0][0] = a[dim];
198 }
199 }
200 }
201 static void set_last_point(Piecewise<D2<SBasis> > &f, Point a){
202 if ( f.empty() ){
203 f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(Linear(a[X]),Linear(a[Y]))));
204 return;
205 }
206 for (unsigned dim=0; dim<2; dim++){
207 if (f.segs.back()[dim].size() == 0){
208 f.segs.back()[dim] = SBasis(Linear(0,a[dim]));
209 }else{
210 f.segs.back()[dim][0][1] = a[dim];
211 }
212 }
213 }
215 std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol){
217 if ( f.size()==0 ) return f;
218 std::vector<Piecewise<D2<SBasis> > > result;
219 std::vector<std::vector<unsigned> > pre_result;
220 for (unsigned i=0; i<f.size(); i++){
221 bool inserted = false;
222 Point a = f[i].firstValue();
223 Point b = f[i].lastValue();
224 for (unsigned j=0; j<pre_result.size(); j++){
225 Point aj = f.at(pre_result[j].back()).lastValue();
226 Point bj = f.at(pre_result[j].front()).firstValue();
227 if ( L2(a-aj) < tol ) {
228 pre_result[j].push_back(i);
229 inserted = true;
230 break;
231 }
232 if ( L2(b-bj) < tol ) {
233 pre_result[j].insert(pre_result[j].begin(),i);
234 inserted = true;
235 break;
236 }
237 }
238 if (!inserted) {
239 pre_result.push_back(std::vector<unsigned>());
240 pre_result.back().push_back(i);
241 }
242 }
243 for (unsigned i=0; i<pre_result.size(); i++){
244 Piecewise<D2<SBasis> > comp;
245 for (unsigned j=0; j<pre_result[i].size(); j++){
246 Piecewise<D2<SBasis> > new_comp = f.at(pre_result[i][j]);
247 if ( j>0 ){
248 set_first_point( new_comp, comp.segs.back().at1() );
249 }
250 comp.concat(new_comp);
251 }
252 if ( L2(comp.firstValue()-comp.lastValue()) < tol ){
253 //TODO: check sizes!!!
254 set_last_point( comp, comp.segs.front().at0() );
255 }
256 result.push_back(comp);
257 }
258 return result;
259 return f;
260 }
263 } // namespace Geom
266 /*
267 Local Variables:
268 mode:c++
269 c-file-style:"stroustrup"
270 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
271 indent-tabs-mode:nil
272 fill-column:99
273 End:
274 */
275 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :