diff options
Diffstat (limited to 'src/2geom/sbasis-2d.cpp')
| -rw-r--r-- | src/2geom/sbasis-2d.cpp | 117 |
1 files changed, 117 insertions, 0 deletions
diff --git a/src/2geom/sbasis-2d.cpp b/src/2geom/sbasis-2d.cpp index 79e99519a..9566e0a19 100644 --- a/src/2geom/sbasis-2d.cpp +++ b/src/2geom/sbasis-2d.cpp @@ -1,4 +1,5 @@ #include <2geom/sbasis-2d.h> +#include <2geom/sbasis-geometric.h> namespace Geom{ @@ -69,6 +70,122 @@ compose_each(D2<SBasis2d> const &fg, D2<SBasis> const &p) { return D2<SBasis>(compose(fg[X], p), compose(fg[Y], p)); } +SBasis2d partial_derivative(SBasis2d const &f, int dim) { + SBasis2d result; + for(unsigned i = 0; i < f.size(); i++) { + result.push_back(Linear2d(0,0,0,0)); + } + result.us = f.us; + result.vs = f.vs; + + for(unsigned i = 0; i < f.us; i++) { + for(unsigned j = 0; j < f.vs; j++) { + Linear2d lin = f.index(i,j); + 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]); + result[i+j*result.us] += dlin; + unsigned di = dim?j:i; + if (di>=1){ + float motpi = dim?-1:1; + Linear2d ds_lin_low( lin[0], -motpi*lin[1], motpi*lin[2], -lin[3] ); + result[(i+dim-1)+(j-dim)*result.us] += di*ds_lin_low; + + 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] ); + result[i+j*result.us] += di*ds_lin_hi; + } + } + } + return result; +} + +/** + * Finds a path which traces the 0 contour of f, traversing from A to B as a single d2<sbasis>. + * degmax specifies the degree (degree = 2*degmax-1, so a degmax of 2 generates a cubic fit). + * The algorithm is based on dividing out derivatives at each end point and does not use the curvature for fitting. + * It is less accurate than sb2d_cubic_solve, although this may be fixed in the future. + */ +D2<SBasis> +sb2dsolve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B, unsigned degmax){ + D2<SBasis>result(Linear(A[X],B[X]),Linear(A[Y],B[Y])); + //g_warning("check f(A)= %f = f(B) = %f =0!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y])); + + SBasis2d dfdu = partial_derivative(f, 0); + SBasis2d dfdv = partial_derivative(f, 1); + Geom::Point dfA(dfdu.apply(A[X],A[Y]),dfdv.apply(A[X],A[Y])); + Geom::Point dfB(dfdu.apply(B[X],B[Y]),dfdv.apply(B[X],B[Y])); + Geom::Point nA = dfA/(dfA[X]*dfA[X]+dfA[Y]*dfA[Y]); + Geom::Point nB = dfB/(dfB[X]*dfB[X]+dfB[Y]*dfB[Y]); + + double fact_k=1; + double sign = 1.; + for(unsigned k=1; k<degmax; k++){ + // these two lines make the solutions worse! + //fact_k *= k; + //sign = -sign; + SBasis f_on_curve = compose(f,result); + Linear reste = f_on_curve[k]; + double ax = -reste[0]/fact_k*nA[X]; + double ay = -reste[0]/fact_k*nA[Y]; + double bx = -sign*reste[1]/fact_k*nB[X]; + double by = -sign*reste[1]/fact_k*nB[Y]; + + result[X].push_back(Linear(ax,bx)); + result[Y].push_back(Linear(ay,by)); + //sign *= 3; + } + return result; +} + +/** + * Finds a path which traces the 0 contour of f, traversing from A to B as a single cubic d2<sbasis>. + * The algorithm is based on matching direction and curvature at each end point. + */ +//TODO: handle the case when B is "behind" A for the natural orientation of the level set. +//TODO: more generally, there might be up to 4 solutions. Choose the best one! +D2<SBasis> +sb2d_cubic_solve(SBasis2d const &f, Geom::Point const &A, Geom::Point const &B){ + D2<SBasis>result;//(Linear(A[X],B[X]),Linear(A[Y],B[Y])); + //g_warning("check 0 = %f = %f!", f.apply(A[X],A[Y]), f.apply(B[X],B[Y])); + + SBasis2d f_u = partial_derivative(f , 0); + SBasis2d f_v = partial_derivative(f , 1); + SBasis2d f_uu = partial_derivative(f_u, 0); + SBasis2d f_uv = partial_derivative(f_v, 0); + SBasis2d f_vv = partial_derivative(f_v, 1); + + Geom::Point dfA(f_u.apply(A[X],A[Y]),f_v.apply(A[X],A[Y])); + Geom::Point dfB(f_u.apply(B[X],B[Y]),f_v.apply(B[X],B[Y])); + + Geom::Point V0 = rot90(dfA); + Geom::Point V1 = rot90(dfB); + + double D2fVV0 = f_uu.apply(A[X],A[Y])*V0[X]*V0[X]+ + 2*f_uv.apply(A[X],A[Y])*V0[X]*V0[Y]+ + f_vv.apply(A[X],A[Y])*V0[Y]*V0[Y]; + double D2fVV1 = f_uu.apply(B[X],B[Y])*V1[X]*V1[X]+ + 2*f_uv.apply(B[X],B[Y])*V1[X]*V1[Y]+ + f_vv.apply(B[X],B[Y])*V1[Y]*V1[Y]; + + std::vector<D2<SBasis> > candidates = cubics_fitting_curvature(A,B,V0,V1,D2fVV0,D2fVV1); + if (candidates.size()==0) { + return D2<SBasis>(Linear(A[X],B[X]),Linear(A[Y],B[Y])); + } + //TODO: I'm sure std algorithm could do that for me... + double error = -1; + unsigned best = 0; + for (unsigned i=0; i<candidates.size(); i++){ + Interval bounds = bounds_fast(compose(f,candidates[i])); + double new_error = (fabs(bounds.max())>fabs(bounds.min()) ? fabs(bounds.max()) : fabs(bounds.min()) ); + if ( new_error < error || error < 0 ){ + error = new_error; + best = i; + } + } + return candidates[best]; +} + + + + }; /* |