Logo Search packages:      
Sourcecode: inkscape version File versions

Piecewise< D2< SBasis > > Geom::unitVector ( D2< SBasis > const &  V_in,
double  tol,
unsigned  order 
)

Return a Piecewise<D2<SBasis> > which points in the same direction as V_in, but has unit_length.

Parameters:
V_in the original path.
tol the maximum error allowed.
order the maximum degree to use for approximation
unitVector(x,y) is computed as (b,-a) where a and b are solutions of: ax+by=0 (eqn1) and a^2+b^2=1 (eqn2)

Definition at line 223 of file sbasis-geometric.cpp.

References Geom::D2< T >::at0(), Geom::D2< T >::at1(), Geom::Piecewise< T >::concat(), Geom::Piecewise< T >::setDomain(), and unit_vector().

Referenced by curvature(), and unitVector().

                                                                  {
    //TODO: Handle vanishing vectors...
    // -This approach is numerically bad. Find a stable way to rescale V_in to have non vanishing ends.
    // -This done, unitVector will have jumps at zeros: fill the gaps with arcs of circles.
    D2<SBasis> V = RescaleForNonVanishingEnds(V_in);

    if (V[0].empty() && V[1].empty())
        return Piecewise<D2<SBasis> >(D2<SBasis>(Linear(1),SBasis()));
    SBasis x = V[0], y = V[1];
    SBasis r_eqn1, r_eqn2;

    Point v0 = unit_vector(V.at0());
    Point v1 = unit_vector(V.at1());
    SBasis a = SBasis(order+1, Linear(0.));
    a[0] = Linear(-v0[1],-v1[1]);
    SBasis b = SBasis(order+1, Linear(0.));
    b[0] = Linear( v0[0], v1[0]);

    r_eqn1 = -(a*x+b*y);
    r_eqn2 = Linear(1.)-(a*a+b*b);

    for (unsigned k=1; k<=order; k++){
        double r0  = (k<r_eqn1.size())? r_eqn1.at(k).at0() : 0;
        double r1  = (k<r_eqn1.size())? r_eqn1.at(k).at1() : 0;
        double rr0 = (k<r_eqn2.size())? r_eqn2.at(k).at0() : 0;
        double rr1 = (k<r_eqn2.size())? r_eqn2.at(k).at1() : 0;
        double a0,a1,b0,b1;// coeffs in a[k] and b[k]

        //the equations to solve at this point are:
        // a0*x(0)+b0*y(0)=r0 & 2*a0*a(0)+2*b0*b(0)=rr0
        //and
        // a1*x(1)+b1*y(1)=r1 & 2*a1*a(1)+2*b1*b(1)=rr1
        a0 = r0/dot(v0,V.at0())*v0[0]-rr0/2*v0[1];
        b0 = r0/dot(v0,V.at0())*v0[1]+rr0/2*v0[0];
        a1 = r1/dot(v1,V.at1())*v1[0]-rr1/2*v1[1];
        b1 = r1/dot(v1,V.at1())*v1[1]+rr1/2*v1[0];

        a[k] = Linear(a0,a1);
        b[k] = Linear(b0,b1);

        //TODO: use "incremental" rather than explicit formulas.
        r_eqn1 = -(a*x+b*y);
        r_eqn2 = Linear(1)-(a*a+b*b);
    }
    
    //our candidate is:
    D2<SBasis> unitV;
    unitV[0] =  b;
    unitV[1] = -a;

    //is it good?
    double rel_tol = std::max(1.,std::max(V_in[0].tailError(0),V_in[1].tailError(0)))*tol;
    if (r_eqn1.tailError(order)>rel_tol || r_eqn2.tailError(order)>tol){
        //if not: subdivide and concat results.
        Piecewise<D2<SBasis> > unitV0, unitV1;
        unitV0 = unitVector(compose(V,Linear(0,.5)),tol,order);
        unitV1 = unitVector(compose(V,Linear(.5,1)),tol,order);
        unitV0.setDomain(Interval(0.,.5));
        unitV1.setDomain(Interval(.5,1.));
        unitV0.concat(unitV1);
        return(unitV0);
    }else{
        //if yes: return it as pw.
        Piecewise<D2<SBasis> > result;
        result=(Piecewise<D2<SBasis> >)unitV;
        return result;
    }
}


Generated by  Doxygen 1.6.0   Back to index