1 /**
2 * \file
3 * \brief \todo brief description
4 *
5 * Authors:
6 * Nathan Hurst <?@?.?>
7 * JFBarraud <?@?.?>
8 *
9 * Copyright 2006-2008 authors
10 *
11 * This library is free software; you can redistribute it and/or
12 * modify it either under the terms of the GNU Lesser General Public
13 * License version 2.1 as published by the Free Software Foundation
14 * (the "LGPL") or, at your option, under the terms of the Mozilla
15 * Public License Version 1.1 (the "MPL"). If you do not alter this
16 * notice, a recipient may use your version of this file under either
17 * the MPL or the LGPL.
18 *
19 * You should have received a copy of the LGPL along with this library
20 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 * You should have received a copy of the MPL along with this library
23 * in the file COPYING-MPL-1.1
24 *
25 * The contents of this file are subject to the Mozilla Public License
26 * Version 1.1 (the "License"); you may not use this file except in
27 * compliance with the License. You may obtain a copy of the License at
28 * http://www.mozilla.org/MPL/
29 *
30 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
31 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
32 * the specific language governing rights and limitations.
33 *
34 */
36 #ifndef SEEN_SBASIS_2D_H
37 #define SEEN_SBASIS_2D_H
38 #include <vector>
39 #include <cassert>
40 #include <algorithm>
41 #include <2geom/d2.h>
42 #include <2geom/sbasis.h>
43 #include <iostream>
45 namespace Geom{
47 class Linear2d{
48 public:
49 /*
50 u 0,1
51 v 0,2
52 */
53 double a[4];
54 Linear2d() {
55 a[0] = 0;
56 a[1] = 0;
57 a[2] = 0;
58 a[3] = 0;
59 }
60 Linear2d(double aa) {
61 for(unsigned i = 0 ; i < 4; i ++)
62 a[i] = aa;
63 }
64 Linear2d(double a00, double a01, double a10, double a11)
65 {
66 a[0] = a00;
67 a[1] = a01;
68 a[2] = a10;
69 a[3] = a11;
70 }
72 double operator[](const int i) const {
73 assert(i >= 0);
74 assert(i < 4);
75 return a[i];
76 }
77 double& operator[](const int i) {
78 assert(i >= 0);
79 assert(i < 4);
80 return a[i];
81 }
82 double apply(double u, double v) {
83 return (a[0]*(1-u)*(1-v) +
84 a[1]*u*(1-v) +
85 a[2]*(1-u)*v +
86 a[3]*u*v);
87 }
88 };
90 inline Linear extract_u(Linear2d const &a, double u) {
91 return Linear(a[0]*(1-u) +
92 a[1]*u,
93 a[2]*(1-u) +
94 a[3]*u);
95 }
96 inline Linear extract_v(Linear2d const &a, double v) {
97 return Linear(a[0]*(1-v) +
98 a[2]*v,
99 a[1]*(1-v) +
100 a[3]*v);
101 }
102 inline Linear2d operator-(Linear2d const &a) {
103 return Linear2d(-a.a[0], -a.a[1],
104 -a.a[2], -a.a[3]);
105 }
106 inline Linear2d operator+(Linear2d const & a, Linear2d const & b) {
107 return Linear2d(a[0] + b[0],
108 a[1] + b[1],
109 a[2] + b[2],
110 a[3] + b[3]);
111 }
112 inline Linear2d operator-(Linear2d const & a, Linear2d const & b) {
113 return Linear2d(a[0] - b[0],
114 a[1] - b[1],
115 a[2] - b[2],
116 a[3] - b[3]);
117 }
118 inline Linear2d& operator+=(Linear2d & a, Linear2d const & b) {
119 for(unsigned i = 0; i < 4; i++)
120 a[i] += b[i];
121 return a;
122 }
123 inline Linear2d& operator-=(Linear2d & a, Linear2d const & b) {
124 for(unsigned i = 0; i < 4; i++)
125 a[i] -= b[i];
126 return a;
127 }
128 inline Linear2d& operator*=(Linear2d & a, double b) {
129 for(unsigned i = 0; i < 4; i++)
130 a[i] *= b;
131 return a;
132 }
134 inline bool operator==(Linear2d const & a, Linear2d const & b) {
135 for(unsigned i = 0; i < 4; i++)
136 if(a[i] != b[i])
137 return false;
138 return true;
139 }
140 inline bool operator!=(Linear2d const & a, Linear2d const & b) {
141 for(unsigned i = 0; i < 4; i++)
142 if(a[i] == b[i])
143 return false;
144 return true;
145 }
146 inline Linear2d operator*(double const a, Linear2d const & b) {
147 return Linear2d(a*b[0], a*b[1],
148 a*b[2], a*b[3]);
149 }
151 class SBasis2d : public std::vector<Linear2d>{
152 public:
153 // vector in u,v
154 unsigned us, vs; // number of u terms, v terms
155 SBasis2d() {}
156 SBasis2d(Linear2d const & bo)
157 : us(1), vs(1) {
158 push_back(bo);
159 }
160 SBasis2d(SBasis2d const & a)
161 : std::vector<Linear2d>(a), us(a.us), vs(a.vs) {}
163 Linear2d& index(unsigned ui, unsigned vi) {
164 assert(ui < us);
165 assert(vi < vs);
166 return (*this)[ui + vi*us];
167 }
169 Linear2d index(unsigned ui, unsigned vi) const {
170 if(ui >= us)
171 return Linear2d(0);
172 if(vi >= vs)
173 return Linear2d(0);
174 return (*this)[ui + vi*us];
175 }
177 double apply(double u, double v) const {
178 double s = u*(1-u);
179 double t = v*(1-v);
180 Linear2d p;
181 double tk = 1;
182 // XXX rewrite as horner
183 for(unsigned vi = 0; vi < vs; vi++) {
184 double sk = 1;
185 for(unsigned ui = 0; ui < us; ui++) {
186 p += (sk*tk)*index(ui, vi);
187 sk *= s;
188 }
189 tk *= t;
190 }
191 return p.apply(u,v);
192 }
194 void clear() {
195 fill(begin(), end(), Linear2d(0));
196 }
198 void normalize(); // remove extra zeros
200 double tail_error(unsigned tail) const;
202 void truncate(unsigned k);
203 };
205 inline SBasis2d operator-(const SBasis2d& p) {
206 SBasis2d result;
207 result.reserve(p.size());
209 for(unsigned i = 0; i < p.size(); i++) {
210 result.push_back(-p[i]);
211 }
212 return result;
213 }
215 inline SBasis2d operator+(const SBasis2d& a, const SBasis2d& b) {
216 SBasis2d result;
217 result.us = std::max(a.us, b.us);
218 result.vs = std::max(a.vs, b.vs);
219 const unsigned out_size = result.us*result.vs;
220 result.resize(out_size);
222 for(unsigned vi = 0; vi < result.vs; vi++) {
223 for(unsigned ui = 0; ui < result.us; ui++) {
224 Linear2d bo;
225 if(ui < a.us && vi < a.vs)
226 bo += a.index(ui, vi);
227 if(ui < b.us && vi < b.vs)
228 bo += b.index(ui, vi);
229 result.index(ui, vi) = bo;
230 }
231 }
232 return result;
233 }
235 inline SBasis2d operator-(const SBasis2d& a, const SBasis2d& b) {
236 SBasis2d result;
237 result.us = std::max(a.us, b.us);
238 result.vs = std::max(a.vs, b.vs);
239 const unsigned out_size = result.us*result.vs;
240 result.resize(out_size);
242 for(unsigned vi = 0; vi < result.vs; vi++) {
243 for(unsigned ui = 0; ui < result.us; ui++) {
244 Linear2d bo;
245 if(ui < a.us && vi < a.vs)
246 bo += a.index(ui, vi);
247 if(ui < b.us && vi < b.vs)
248 bo -= b.index(ui, vi);
249 result.index(ui, vi) = bo;
250 }
251 }
252 return result;
253 }
256 inline SBasis2d& operator+=(SBasis2d& a, const Linear2d& b) {
257 if(a.size() < 1)
258 a.push_back(b);
259 else
260 a[0] += b;
261 return a;
262 }
264 inline SBasis2d& operator-=(SBasis2d& a, const Linear2d& b) {
265 if(a.size() < 1)
266 a.push_back(-b);
267 else
268 a[0] -= b;
269 return a;
270 }
272 inline SBasis2d& operator+=(SBasis2d& a, double b) {
273 if(a.size() < 1)
274 a.push_back(Linear2d(b));
275 else {
276 for(unsigned i = 0; i < 4; i++)
277 a[0] += double(b);
278 }
279 return a;
280 }
282 inline SBasis2d& operator-=(SBasis2d& a, double b) {
283 if(a.size() < 1)
284 a.push_back(Linear2d(-b));
285 else {
286 a[0] -= b;
287 }
288 return a;
289 }
291 inline SBasis2d& operator*=(SBasis2d& a, double b) {
292 for(unsigned i = 0; i < a.size(); i++)
293 a[i] *= b;
294 return a;
295 }
297 inline SBasis2d& operator/=(SBasis2d& a, double b) {
298 for(unsigned i = 0; i < a.size(); i++)
299 a[i] *= (1./b);
300 return a;
301 }
303 SBasis2d operator*(double k, SBasis2d const &a);
304 SBasis2d operator*(SBasis2d const &a, SBasis2d const &b);
306 SBasis2d shift(SBasis2d const &a, int sh);
308 SBasis2d shift(Linear2d const &a, int sh);
310 SBasis2d truncate(SBasis2d const &a, unsigned terms);
312 SBasis2d multiply(SBasis2d const &a, SBasis2d const &b);
314 SBasis2d integral(SBasis2d const &c);
316 SBasis2d partial_derivative(SBasis2d const &a, int dim);
318 SBasis2d sqrt(SBasis2d const &a, int k);
320 // return a kth order approx to 1/a)
321 SBasis2d reciprocal(Linear2d const &a, int k);
323 SBasis2d divide(SBasis2d const &a, SBasis2d const &b, int k);
325 // a(b(t))
326 SBasis2d compose(SBasis2d const &a, SBasis2d const &b);
327 SBasis2d compose(SBasis2d const &a, SBasis2d const &b, unsigned k);
328 SBasis2d inverse(SBasis2d const &a, int k);
330 // these two should probably be replaced with compose
331 SBasis extract_u(SBasis2d const &a, double u);
332 SBasis extract_v(SBasis2d const &a, double v);
334 SBasis compose(Linear2d const &a, D2<SBasis> const &p);
336 SBasis compose(SBasis2d const &fg, D2<SBasis> const &p);
338 D2<SBasis> compose_each(D2<SBasis2d> const &fg, D2<SBasis> const &p);
340 inline std::ostream &operator<< (std::ostream &out_file, const Linear2d &bo) {
341 out_file << "{" << bo[0] << ", " << bo[1] << "}, ";
342 out_file << "{" << bo[2] << ", " << bo[3] << "}";
343 return out_file;
344 }
346 inline std::ostream &operator<< (std::ostream &out_file, const SBasis2d & p) {
347 for(unsigned i = 0; i < p.size(); i++) {
348 out_file << p[i] << "s^" << i << " + ";
349 }
350 return out_file;
351 }
353 D2<SBasis>
354 sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigned degmax=2);
356 D2<SBasis>
357 sb2d_cubic_solve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B);
359 };
361 /*
362 Local Variables:
363 mode:c++
364 c-file-style:"stroustrup"
365 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
366 indent-tabs-mode:nil
367 fill-column:99
368 End:
369 */
370 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
371 #endif