9566e0a198fb5f5b153821ff73b60441c581f1c8
1 #include <2geom/sbasis-2d.h>
2 #include <2geom/sbasis-geometric.h>
4 namespace Geom{
6 SBasis extract_u(SBasis2d const &a, double u) {
7 SBasis sb;
8 double s = u*(1-u);
10 for(unsigned vi = 0; vi < a.vs; vi++) {
11 double sk = 1;
12 Linear bo(0,0);
13 for(unsigned ui = 0; ui < a.us; ui++) {
14 bo += (extract_u(a.index(ui, vi), u))*sk;
15 sk *= s;
16 }
17 sb.push_back(bo);
18 }
20 return sb;
21 }
23 SBasis extract_v(SBasis2d const &a, double v) {
24 SBasis sb;
25 double s = v*(1-v);
27 for(unsigned ui = 0; ui < a.us; ui++) {
28 double sk = 1;
29 Linear bo(0,0);
30 for(unsigned vi = 0; vi < a.vs; vi++) {
31 bo += (extract_v(a.index(ui, vi), v))*sk;
32 sk *= s;
33 }
34 sb.push_back(bo);
35 }
37 return sb;
38 }
40 SBasis compose(Linear2d const &a, D2<SBasis> const &p) {
41 D2<SBasis> omp(-p[X] + 1, -p[Y] + 1);
42 return multiply(omp[0], omp[1])*a[0] +
43 multiply(p[0], omp[1])*a[1] +
44 multiply(omp[0], p[1])*a[2] +
45 multiply(p[0], p[1])*a[3];
46 }
48 SBasis
49 compose(SBasis2d const &fg, D2<SBasis> const &p) {
50 SBasis B;
51 SBasis s[2];
52 SBasis ss[2];
53 for(unsigned dim = 0; dim < 2; dim++)
54 s[dim] = p[dim]*(Linear(1) - p[dim]);
55 ss[1] = Linear(1);
56 for(unsigned vi = 0; vi < fg.vs; vi++) {
57 ss[0] = ss[1];
58 for(unsigned ui = 0; ui < fg.us; ui++) {
59 unsigned i = ui + vi*fg.us;
60 B += ss[0]*compose(fg[i], p);
61 ss[0] *= s[0];
62 }
63 ss[1] *= s[1];
64 }
65 return B;
66 }
68 D2<SBasis>
69 compose_each(D2<SBasis2d> const &fg, D2<SBasis> const &p) {
70 return D2<SBasis>(compose(fg[X], p), compose(fg[Y], p));
71 }
73 SBasis2d partial_derivative(SBasis2d const &f, int dim) {
74 SBasis2d result;
75 for(unsigned i = 0; i < f.size(); i++) {
76 result.push_back(Linear2d(0,0,0,0));
77 }
78 result.us = f.us;
79 result.vs = f.vs;
81 for(unsigned i = 0; i < f.us; i++) {
82 for(unsigned j = 0; j < f.vs; j++) {
83 Linear2d lin = f.index(i,j);
84 Linear2d dlin(lin[1+dim]-lin[0], lin[1+2*dim]-lin[dim], lin[3-dim]-lin[2*(1-dim)], lin[3]-lin[2-dim]);
85 result[i+j*result.us] += dlin;
86 unsigned di = dim?j:i;
87 if (di>=1){
88 float motpi = dim?-1:1;
89 Linear2d ds_lin_low( lin[0], -motpi*lin[1], motpi*lin[2], -lin[3] );
90 result[(i+dim-1)+(j-dim)*result.us] += di*ds_lin_low;
92 Linear2d ds_lin_hi( lin[1+dim]-lin[0], lin[1+2*dim]-lin[dim], lin[3]-lin[2-dim], lin[3-dim]-lin[2-dim] );
93 result[i+j*result.us] += di*ds_lin_hi;
94 }
95 }
96 }
97 return result;
98 }
100 /**
101 * Finds a path which traces the 0 contour of f, traversing from A to B as a single d2<sbasis>.
102 * degmax specifies the degree (degree = 2*degmax-1, so a degmax of 2 generates a cubic fit).
103 * The algorithm is based on dividing out derivatives at each end point and does not use the curvature for fitting.
104 * It is less accurate than sb2d_cubic_solve, although this may be fixed in the future.
105 */
106 D2<SBasis>
107 sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigned degmax){
108 D2<SBasis>result(Linear(A[X],B[X]),Linear(A[Y],B[Y]));
109 //g_warning("check f(A)= %f = f(B) = %f =0!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y]));
111 SBasis2d dfdu = partial_derivative(f, 0);
112 SBasis2d dfdv = partial_derivative(f, 1);
113 Geom::Point dfA(dfdu.apply(A[X],A[Y]),dfdv.apply(A[X],A[Y]));
114 Geom::Point dfB(dfdu.apply(B[X],B[Y]),dfdv.apply(B[X],B[Y]));
115 Geom::Point nA = dfA/(dfA[X]*dfA[X]+dfA[Y]*dfA[Y]);
116 Geom::Point nB = dfB/(dfB[X]*dfB[X]+dfB[Y]*dfB[Y]);
118 double fact_k=1;
119 double sign = 1.;
120 for(unsigned k=1; k<degmax; k++){
121 // these two lines make the solutions worse!
122 //fact_k *= k;
123 //sign = -sign;
124 SBasis f_on_curve = compose(f,result);
125 Linear reste = f_on_curve[k];
126 double ax = -reste[0]/fact_k*nA[X];
127 double ay = -reste[0]/fact_k*nA[Y];
128 double bx = -sign*reste[1]/fact_k*nB[X];
129 double by = -sign*reste[1]/fact_k*nB[Y];
131 result[X].push_back(Linear(ax,bx));
132 result[Y].push_back(Linear(ay,by));
133 //sign *= 3;
134 }
135 return result;
136 }
138 /**
139 * Finds a path which traces the 0 contour of f, traversing from A to B as a single cubic d2<sbasis>.
140 * The algorithm is based on matching direction and curvature at each end point.
141 */
142 //TODO: handle the case when B is "behind" A for the natural orientation of the level set.
143 //TODO: more generally, there might be up to 4 solutions. Choose the best one!
144 D2<SBasis>
145 sb2d_cubic_solve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B){
146 D2<SBasis>result;//(Linear(A[X],B[X]),Linear(A[Y],B[Y]));
147 //g_warning("check 0 = %f = %f!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y]));
149 SBasis2d f_u = partial_derivative(f , 0);
150 SBasis2d f_v = partial_derivative(f , 1);
151 SBasis2d f_uu = partial_derivative(f_u, 0);
152 SBasis2d f_uv = partial_derivative(f_v, 0);
153 SBasis2d f_vv = partial_derivative(f_v, 1);
155 Geom::Point dfA(f_u.apply(A[X],A[Y]),f_v.apply(A[X],A[Y]));
156 Geom::Point dfB(f_u.apply(B[X],B[Y]),f_v.apply(B[X],B[Y]));
158 Geom::Point V0 = rot90(dfA);
159 Geom::Point V1 = rot90(dfB);
161 double D2fVV0 = f_uu.apply(A[X],A[Y])*V0[X]*V0[X]+
162 2*f_uv.apply(A[X],A[Y])*V0[X]*V0[Y]+
163 f_vv.apply(A[X],A[Y])*V0[Y]*V0[Y];
164 double D2fVV1 = f_uu.apply(B[X],B[Y])*V1[X]*V1[X]+
165 2*f_uv.apply(B[X],B[Y])*V1[X]*V1[Y]+
166 f_vv.apply(B[X],B[Y])*V1[Y]*V1[Y];
168 std::vector<D2<SBasis> > candidates = cubics_fitting_curvature(A,B,V0,V1,D2fVV0,D2fVV1);
169 if (candidates.size()==0) {
170 return D2<SBasis>(Linear(A[X],B[X]),Linear(A[Y],B[Y]));
171 }
172 //TODO: I'm sure std algorithm could do that for me...
173 double error = -1;
174 unsigned best = 0;
175 for (unsigned i=0; i<candidates.size(); i++){
176 Interval bounds = bounds_fast(compose(f,candidates[i]));
177 double new_error = (fabs(bounds.max())>fabs(bounds.min()) ? fabs(bounds.max()) : fabs(bounds.min()) );
178 if ( new_error < error || error < 0 ){
179 error = new_error;
180 best = i;
181 }
182 }
183 return candidates[best];
184 }
189 };
191 /*
192 Local Variables:
193 mode:c++
194 c-file-style:"stroustrup"
195 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
196 indent-tabs-mode:nil
197 fill-column:99
198 End:
199 */
200 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :