C++ Piecewise类代码示例

    virtual void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save) {
        cairo_set_source_rgba (cr, 0., 0., 0., 1);
        cairo_set_line_width (cr, 1);
        
        Piecewise<SBasis> pws;
        
        //pws = Geom::cos(Linear(0,100)) + 3;
        pws = Geom::sqrt(Linear(0,100));
        //pws = log(Interval(1,8));
        //Piecewise<SBasis> l(Linear(-100,100));
        //Piecewise<SBasis> one(Linear(1,1));
        //pws = Geom::reciprocal(l*l + one)*l + l;
        //pws = xlogx(Interval(0.5,3));
        //pws = Geom::reciprocal(pws);
        //pws = -integral(Geom::reciprocal(Linear(1,2)))*Piecewise<SBasis>(Linear(1,2));
        
        pws = -pws*width/4 + width/2;
        pws.scaleDomain(width/2);
        pws.offsetDomain(width/4);

        cairo_pw(cr, pws);

        cairo_stroke(cr);
        cairo_set_source_rgba (cr, 0., 0., .5, 1.);
        cairo_horiz(cr, 500, pws.cuts);
        cairo_stroke(cr);

        *notify << "total pieces: " << pws.size();

        Toy::draw(cr, notify, width, height, save);
    }

/** Reparameterise M to have unit speed.
 \param M the Element.
 \param tol the maximum error allowed.
 \param order the maximum degree to use for approximation
 \relates Piecewise, D2
*/
Piecewise<D2<SBasis> >
Geom::arc_length_parametrization(D2<SBasis> const &M,
                           unsigned order,
                           double tol){
    Piecewise<D2<SBasis> > u;
    u.push_cut(0);

    Piecewise<SBasis> s = arcLengthSb(Piecewise<D2<SBasis> >(M),tol);
    for (unsigned i=0; i < s.size();i++){
        double t0=s.cuts[i],t1=s.cuts[i+1];
        if ( are_near(s(t0),s(t1)) ) {
            continue;
        }
        D2<SBasis> sub_M = compose(M,Linear(t0,t1));
        D2<SBasis> sub_u;
        for (unsigned dim=0;dim<2;dim++){
            SBasis sub_s = s.segs[i];
            sub_s-=sub_s.at0();
            sub_s/=sub_s.at1();
            sub_u[dim]=compose_inverse(sub_M[dim],sub_s, order, tol);
        }
        u.push(sub_u,s(t1));
    }
    return u;
}

    void draw(cairo_t *cr,
	      std::ostringstream *notify,
	      int width, int height, bool save) {
    
        D2<SBasis> B1 = b1_handle.asBezier();
        D2<SBasis> B2 = b2_handle.asBezier();
        Piecewise<D2<SBasis> >B;
        B.concat(Piecewise<D2<SBasis> >(B1));
        B.concat(Piecewise<D2<SBasis> >(B2));

// testing fuse_nearby_ends
        std::vector< Piecewise<D2<SBasis> > > pieces;
        pieces = fuse_nearby_ends(split_at_discontinuities(B),9);
        Piecewise<D2<SBasis> > C;
        for (unsigned i=0; i<pieces.size(); i++){
            C.concat(pieces[i]);
        }
// testing fuse_nearby_ends

        cairo_set_line_width (cr, .5);
        cairo_set_source_rgba (cr, 0., 0.5, 0., 1);
        //cairo_d2_sb(cr, B1);
        cairo_pw_d2_sb(cr, C);
        //cairo_pw_d2_sb(cr, B);
        cairo_stroke(cr);

        Piecewise<D2<SBasis> > uniform_B = arc_length_parametrization(B);
        cairo_set_source_rgba (cr, 0., 0., 0.9, 1);
        dot_plot(cr,uniform_B);
        cairo_stroke(cr);
        *notify << "pieces = " << uniform_B.size() << ";\n";

        Toy::draw(cr, notify, width, height, save);
    }        

Geom::Piecewise<Geom::D2<Geom::SBasis> >
LPEOffset::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwd2_in)
{
    using namespace Geom;

    Piecewise<D2<SBasis> > output;

    double t = nearest_point(offset_pt, pwd2_in);
    Point A = pwd2_in.valueAt(t);
    double offset = L2(A - offset_pt);

    Piecewise<D2<SBasis> > der = unitVector(derivative(pwd2_in));
    Piecewise<D2<SBasis> > n   = rot90(der);

    output  = pwd2_in + n * offset;
    append_half_circle(output, pwd2_in.lastValue(), n.lastValue() * offset);
    output.continuousConcat(reverse(pwd2_in - n * offset));
    append_half_circle(output, pwd2_in.firstValue(), -n.firstValue() * offset);

    // TODO: here we should remove self-overlaps by applying the "union" boolop
    //       but we'd need to convert the path to a Shape, which is currently
    //       broken in 2geom, so we return the unaltered path

    return output;
}

Geom::Piecewise<Geom::D2<Geom::SBasis> >
LPECopyRotate::doEffect_pwd2 (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwd2_in)
{
    using namespace Geom;

    // I first suspected the minus sign to be a bug in 2geom but it is
    // likely due to SVG's choice of coordinate system orientation (max)
    start_pos = origin + dir * Rotate(-deg_to_rad(starting_angle)) * dist_angle_handle;
    double rotation_angle_end = rotation_angle;
    if(copiesTo360){
        rotation_angle_end = 360.0/(double)num_copies;
    }
    rot_pos = origin + dir * Rotate(-deg_to_rad(starting_angle + rotation_angle_end)) * dist_angle_handle;

    A = pwd2_in.firstValue();
    B = pwd2_in.lastValue();
    dir = unit_vector(B - A);

    Piecewise<D2<SBasis> > output;

    Affine pre = Translate(-origin) * Rotate(-deg_to_rad(starting_angle));
    for (int i = 0; i < num_copies; ++i) {
        // I first suspected the minus sign to be a bug in 2geom but it is
        // likely due to SVG's choice of coordinate system orientation (max)
        Rotate rot(-deg_to_rad(rotation_angle_end * i));
        Affine t = pre * rot * Translate(origin);
        output.concat(pwd2_in * t);
    }

    return output;
}

 void draw(cairo_t *cr, int t) {
     if(t - last_time > path.domain().max()) add_section(random_d2());
     if(t - last_time - length > path.cuts[1]) {
         Piecewise<D2<SBasis> > new_path;
         new_path.push_cut(0);
         for(unsigned i = 1; i < path.size(); i++) {
             new_path.push(path[i], path.cuts[i+1] - path.cuts[1]);
         }
         last_time = t - length;
         path = new_path;
     }
     cairo_set_source_rgb(cr, red, green, blue);
     Piecewise<D2<SBasis> > port = portion(path, std::max(t - last_time - length, 0.), t - last_time);
     cairo_pw_d2_sb(cr, port);
     cairo_stroke(cr);
     
     double d = 4;
     cairo_set_dash(cr, &d, 1, 0);
     for(unsigned i = 1; i < path.size(); i++) {
         if(path[i].at0() != path[i-1].at1()) {
             draw_line_seg(cr, path[i].at0(), path[i-1].at1());
         }
     }
     cairo_stroke(cr);
     cairo_set_dash(cr, &d, 0, 0);
     
     cairo_set_source_rgb(cr, 0., 0., 1.);
     dot_plot(cr, path, std::max(t - last_time - length, 0.), t - last_time);
 }

/** Compute the cosine of a function.
 \param f function
 \param tol maximum error
 \param order maximum degree polynomial to use
*/
Piecewise<SBasis> cos(Piecewise<SBasis> const &f, double tol, int order){
    Piecewise<SBasis> result;
    for (unsigned i=0; i<f.size(); i++){
        Piecewise<SBasis> cosfi = cos(f.segs[i],tol,order);
        cosfi.setDomain(Interval(f.cuts[i],f.cuts[i+1]));
        result.concat(cosfi);
    }
    return result;
}

/** Reparameterise M to have unit speed.
 \param M the Element.
 \param tol the maximum error allowed.
 \param order the maximum degree to use for approximation
 \relates Piecewise
*/
Piecewise<D2<SBasis> >
Geom::arc_length_parametrization(Piecewise<D2<SBasis> > const &M,
                                 unsigned order,
                                 double tol){
    Piecewise<D2<SBasis> > result;
    for (unsigned i=0; i<M.size(); i++) {
        result.concat( arc_length_parametrization(M[i],order,tol) );
    }
    return result;
}

/** returns a function giving the curvature at each point in M.
 \param M the Element.
 \param tol the maximum error allowed.
 \relates Piecewise
 \todo claimed incomplete.  Check.
*/
Piecewise<SBasis> 
Geom::curvature(Piecewise<D2<SBasis> > const &V, double tol){
    Piecewise<SBasis> result;
    Piecewise<D2<SBasis> > VV = cutAtRoots(V);
    result.cuts.push_back(VV.cuts.front());
    for (unsigned i=0; i<VV.size(); i++){
        Piecewise<SBasis> curv_seg;
        curv_seg = curvature(VV.segs[i],tol);
        curv_seg.setDomain(Interval(VV.cuts[i],VV.cuts[i+1]));
        result.concat(curv_seg);
    }
    return result;
}

/** Return a Piecewise<D2<SBasis> > which points in the same direction as V_in, but has unit_length.
 \param V_in the original path.
 \param tol the maximum error allowed.
 \param 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)

 \relates Piecewise
*/
Piecewise<D2<SBasis> >
Geom::unitVector(Piecewise<D2<SBasis> > const &V, double tol, unsigned order){
    Piecewise<D2<SBasis> > result;
    Piecewise<D2<SBasis> > VV = cutAtRoots(V);
    result.cuts.push_back(VV.cuts.front());
    for (unsigned i=0; i<VV.size(); i++){
        Piecewise<D2<SBasis> > unit_seg;
        unit_seg = unitVector(VV.segs[i],tol, order);
        unit_seg.setDomain(Interval(VV.cuts[i],VV.cuts[i+1]));
        result.concat(unit_seg);   
    }
    return result;
}

Piecewise<SBasis> convole(SBasisOf<double> const &f, Interval dom_f, 
                          SBasisOf<double> const &g, Interval dom_g,
                          bool f_cst_ends = false){

    if ( dom_f.extent() < dom_g.extent() ) return convole(g, dom_g, f, dom_f);

    Piecewise<SBasis> result;

    SBasisOf<SBasisOf<double> > u,v;
    u.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,1))));
    v.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,0)),
                                            SBasisOf<double>(LinearOf<double>(1,1))));
    SBasisOf<SBasisOf<double> > v_u = (v - u)*(dom_f.extent()/dom_g.extent());
    v_u += SBasisOf<SBasisOf<double> >(SBasisOf<double>(-dom_g.min()/dom_g.extent()));
    SBasisOf<SBasisOf<double> > gg = multi_compose(g,v_u);
    SBasisOf<SBasisOf<double> > ff = SBasisOf<SBasisOf<double> >(f);
    SBasisOf<SBasisOf<double> > hh = integral(ff*gg,0);
    
    result.cuts.push_back(dom_f.min()+dom_g.min());
    //Note: we know dom_f.extent() >= dom_g.extent()!!
    //double rho = dom_f.extent()/dom_g.extent();
    double t0 = dom_g.min()/dom_f.extent();
    double t1 = dom_g.max()/dom_f.extent();
    double t2 = t0+1;
    double t3 = t1+1;
    SBasisOf<double> a,b,t;
    SBasis seg;
    a = SBasisOf<double>(LinearOf<double>(0,0)); 
    b = SBasisOf<double>(LinearOf<double>(0,t1-t0)); 
    t = SBasisOf<double>(LinearOf<double>(t0,t1)); 
    seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
    result.push(seg,dom_f.min() + dom_g.max());
    if (dom_f.extent() > dom_g.extent()){
        a = SBasisOf<double>(LinearOf<double>(0,t2-t1)); 
        b = SBasisOf<double>(LinearOf<double>(t1-t0,1)); 
        t = SBasisOf<double>(LinearOf<double>(t1,t2)); 
        seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
        result.push(seg,dom_f.max() + dom_g.min());
    }
    a = SBasisOf<double>(LinearOf<double>(t2-t1,1.)); 
    b = SBasisOf<double>(LinearOf<double>(1.,1.)); 
    t = SBasisOf<double>(LinearOf<double>(t2,t3)); 
    seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
    result.push(seg,dom_f.max() + dom_g.max());
    result*=dom_f.extent();
    
    //--constant ends correction--------------
    if (f_cst_ends){
        SBasis ig = toSBasis(integraaal(g))*dom_g.extent();
        ig -= ig.at0();
        Piecewise<SBasis> cor;
        cor.cuts.push_back(dom_f.min()+dom_g.min());
        cor.push(reverse(ig)*f.at0(),dom_f.min()+dom_g.max());
        cor.push(Linear(0),dom_f.max()+dom_g.min());
        cor.push(ig*f.at1(),dom_f.max()+dom_g.max());
        result+=cor;
    }
    //----------------------------------------
    return result;
}

void
plot3d(cairo_t *cr, 
       Piecewise<SBasis> const &x, 
       Piecewise<SBasis> const &y, 
       Piecewise<SBasis> const &z, Frame frame){
    
    Piecewise<SBasis> xx = partition(x,y.cuts);
    Piecewise<SBasis> xxx = partition(xx,z.cuts);
    Piecewise<SBasis> yyy = partition(y,xxx.cuts);
    Piecewise<SBasis> zzz = partition(z,xxx.cuts);
    
    for (unsigned i=0; i<xxx.size(); i++){
        plot3d(cr, xxx[i], yyy[i], zzz[i], frame);
    }
}

static void dot_plot(cairo_t *cr, Piecewise<D2<SBasis> > const &M, double min, double max, double space=10){
    for( double t = min; t < max; t += space) {
        Point pos = M(t), perp = M.valueAndDerivatives(t, 2)[1].cw() * 3;
        draw_line_seg(cr, pos + perp, pos - perp);
    }
    cairo_stroke(cr);
}

 void add_section(const D2<SBasis> x) {
     Piecewise<D2<SBasis> > new_path(path);
     D2<SBasis> seg(x);
     seg[0][0][0] = path.segs.back()[0][0][1];
     seg[1][0][0] = path.segs.back()[1][0][1];
     new_path.push(seg, path.domain().max()+1);
     path = arc_length_parametrization(new_path);
 }

void cairo_pw(cairo_t *cr, Piecewise<SBasis> p, double hscale=1., double vscale=1.) {
    for(unsigned i = 0; i < p.size(); i++) {
        D2<SBasis> B;
        B[0] = Linear(150+p.cuts[i]*hscale, 150+p.cuts[i+1]*hscale);
        B[1] = Linear(450) - p[i]*vscale;
        cairo_d2_sb(cr, B);
    }
}