summaryrefslogtreecommitdiffstats
path: root/src/3rdparty/2geom/src/toys/curvature-curve.cpp
blob: 8f8f7bc85b2de0269391042a5bb821a21c3f70c3 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#include <2geom/d2.h>
#include <2geom/sbasis.h>
#include <2geom/sbasis-2d.h>
#include <2geom/bezier-to-sbasis.h>
#include <2geom/sbasis-geometric.h>

#include <toys/path-cairo.h>
#include <toys/toy-framework-2.h>

#include <time.h>
using std::vector;
using namespace Geom;
using namespace std;

//-----------------------------------------------

class CurvatureTester: public Toy {
    PointSetHandle curve_handle;
    Path current_curve;
    void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save, std::ostringstream *timer_stream) override {
        cairo_set_line_width (cr, 1);
        current_curve = Path();
        
        for(int base_i = 0; base_i < int(curve_handle.pts.size()/2) - 1; base_i++) {
            for(int i = 0; i < 2; i++) {
                Geom::Point center=curve_handle.pts[1+2*i+base_i*2];
                Geom::Point normal=center- curve_handle.pts[2*i+base_i*2];
                double radius = Geom::L2(normal);
                *notify<<"K="<<radius<<std::endl;
                if (fabs(radius)>1e-4){
                    
                    double ang = atan2(-normal);
                    cairo_arc(cr, center[0], center[1],fabs(radius), ang-0.3, ang+0.3);
                    cairo_set_source_rgba (cr, 0.75, 0.89, 1., 1);
                    cairo_stroke(cr);
                    
                    
                    //draw_handle(cr, center);
                }else{
                }
                {
                cairo_save(cr);
                double dashes[2] = {10, 10};
                cairo_set_dash(cr, dashes, 2, 0);
                cairo_move_to(cr, center);
                cairo_line_to(cr, center-normal);
                cairo_stroke(cr);
                cairo_restore(cr);
                }
            }
            cairo_set_source_rgba (cr, 0., 0, 0., 1);
                Geom::Point A = curve_handle.pts[0+base_i*2];
                Geom::Point B = curve_handle.pts[2+base_i*2];
            D2<SBasis> best_c = D2<SBasis>(SBasis(Linear(A[X],B[X])),SBasis(Linear(A[Y],B[Y])));
            double error = -1;
            for(int i = 0; i < 16; i++) {
                Geom::Point dA = curve_handle.pts[1+base_i*2]-A;
                Geom::Point dB = curve_handle.pts[3+base_i*2]-B;
                std::vector<D2<SBasis> > candidates = 
                    cubics_fitting_curvature(curve_handle.pts[0+base_i*2],curve_handle.pts[2+base_i*2],
                                             (i&2)?rot90(dA):-rot90(dA),
                                             (i&1)?rot90(dB):-rot90(dB),
                                             ((i&4)?-1:1)*L2sq(dA), ((i&8)?-1:1)*L2sq(dB));
	  
                if (candidates.empty()) {
                } else {
                    //TODO: I'm sure std algorithm could do that for me...
                    unsigned best = 0;
                    for (unsigned i=0; i<candidates.size(); i++){
                        Piecewise<SBasis> K = arcLengthSb(candidates[i]);
        
                        double l = Geom::length(candidates[i]);
                        //double l = K.segs.back().at1();//Geom::length(candidates[i]);
                        //printf("l = %g\n", l);
                        if ( l < error || error < 0 ){
                            error = l;
                            best = i;
                            best_c = candidates[best];
                        }
                    }
                }
            }
            if(error >= 0) {
                //cairo_d2_sb(cr, best_c);
                current_curve.append(best_c);
            }
        }
        
        cairo_path(cr, current_curve);
        cairo_stroke(cr);
        Toy::draw(cr, notify, width, height, save,timer_stream);
    }        
  
    void canvas_click(Geom::Point at, int button) override {
        std::cout << "clicked at:" << at << " with button " << button << std::endl;
        if(button == 1) {
            double dist;
            double t = current_curve.nearestTime(at, &dist).asFlatTime();
            if(dist > 5) {
                curve_handle.push_back(at);
                curve_handle.push_back(at+Point(100,100));
            } else {
                // split around t
                Piecewise<D2<SBasis> > pw =  current_curve.toPwSb();
                std::vector<Point > vnd = pw.valueAndDerivatives(t, 2);
                Point on_curve = current_curve(t);
                Point normal = rot90(vnd[1]);
                Piecewise<SBasis > K = curvature(pw);
                Point ps[2] = {on_curve, on_curve+unit_vector(normal)/K(t)};
                curve_handle.pts.insert(curve_handle.pts.begin()+2*(int(t)+1), ps, ps+2);
            }
        }
    }

public:
    CurvatureTester(){
        if(handles.empty()) {
            handles.push_back(&curve_handle);
            for(unsigned i = 0; i < 4; i++)
                curve_handle.push_back(150+uniform()*300,150+uniform()*300);
        }
    }
};

int main(int argc, char **argv) {
    std::cout << "testing unit_normal(multidim_sbasis) based offset." << std::endl;
    init(argc, argv, new CurvatureTester);
    return 0;
}

/*
  Local Variables:
  mode:c++
  c-file-style:"stroustrup"
  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
  indent-tabs-mode:nil
  fill-column:99
  End:
*/
//vim:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :