Inkscape
Vector Graphics Editor
Loading...
Searching...
No Matches
convole.cpp
Go to the documentation of this file.
1#include <2geom/d2.h>
2#include <2geom/sbasis.h>
3#include <2geom/bezier-to-sbasis.h>
4#include <2geom/sbasis-geometric.h>
5
6#include <toys/path-cairo.h>
7#include <toys/toy-framework-2.h>
8
9#include <2geom/orphan-code/sbasis-of.h>
10
11#include <vector>
12using std::vector;
13using namespace Geom;
14
15SBasis toSBasis(SBasisOf<double> const &f){
16 SBasis result(f.size(), Linear());
17 for (unsigned i=0; i<f.size(); i++){
18 result[i] = Linear(f[i][0],f[i][1]);
19 }
20 return result;
21}
22SBasisOf<double> toSBasisOfDouble(SBasis const &f){
23 SBasisOf<double> result;
24 for (auto i : f){
25 result.push_back(LinearOf<double>(i[0],i[1]));
26 }
27 return result;
28}
29
30
31
32SBasis integral(SBasis const &c) {
33 SBasis a;
34 a.resize(c.size() + 1, Linear(0,0));
35 a[0] = Linear(0,0);
36
37 for(unsigned k = 1; k < c.size() + 1; k++) {
38 double ahat = -c[k-1].tri()/(2*k);
39 a[k][0] = a[k][1] = ahat;
40 }
41 double aTri = 0;
42 for(int k = c.size()-1; k >= 0; k--) {
43 aTri = (c[k].hat() + (k+1)*aTri/2)/(2*k+1);
44 a[k][0] -= aTri/2;
45 a[k][1] += aTri/2;
46 }
47 a.normalize();
48 return a;
49}
50template<typename T>
51SBasisOf<T> integraaal(SBasisOf<T> const &c){
52 SBasisOf<T> a;
53 a.resize(c.size() + 1, LinearOf<T>(T(0.),T(0.)));
54
55 for(unsigned k = 1; k < c.size() + 1; k++) {
56 T ahat = (c[k-1][0]-c[k-1][1])/(2*k);
57 a[k] = LinearOf<T>(ahat);
58 }
59
60 T aTri = T(0.);
61 for(int k = c.size()-1; k >= 0; k--) {
62 aTri = (HatOf<T>(c[k]).d + (k+1)*aTri/2)/(2*k+1);
63 a[k][0] -= aTri/2;
64 a[k][1] += aTri/2;
65 }
66 //a.normalize();
67 return a;
68}
69
70SBasisOf<SBasisOf<double> > integral(SBasisOf<SBasisOf<double> > const &f, unsigned var){
71 //variable of f = 1, variable of f's coefficients = 0.
72 if (var == 1) return integraaal(f);
73 SBasisOf<SBasisOf<double> > result;
74 for(unsigned i = 0; i< f.size(); i++) {
75 //result.push_back(LinearOf<SBasisOf<double> >( integral(f[i][0]),integral(f[i][1])));
76 result.push_back(LinearOf<SBasisOf<double> >( integraaal(f[i][0]),integraaal(f[i][1])));
77 }
78 return result;
79}
80
81template <typename T>
82SBasisOf<T> multi_compose(SBasisOf<double> const &f, SBasisOf<T> const &g ){
83
84 //assert( f.input_dim() <= g.size() );
85
86 SBasisOf<T> u1 = g;
87 SBasisOf<T> u0 = -g + LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(1,1)));
88 SBasisOf<T> s = multiply(u0,u1);
89 SBasisOf<T> r;
90
91 for(int i = f.size()-1; i >= 0; i--) {
92 r = s*r + f[i][0]*u0 + f[i][1]*u1;
93 }
94 return r;
95}
96
97SBasisOf<double> compose(SBasisOf<SBasisOf<double> > const &f,
98 SBasisOf<double> const &x,
99 SBasisOf<double> const &y){
100 SBasisOf<double> y0 = -y + LinearOf<double>(1,1);
101 SBasisOf<double> s = multiply(y0,y);
102 SBasisOf<double> r;
103
104 for(int i = f.size()-1; i >= 0; i--) {
105 r = s*r + compose(f[i][0],x)*y0 + compose(f[i][1],x)*y;
106 }
107 return r;
108}
109
110Piecewise<SBasis> convole(SBasisOf<double> const &f, Interval dom_f,
111 SBasisOf<double> const &g, Interval dom_g,
112 bool f_cst_ends = false){
113
114 if ( dom_f.extent() < dom_g.extent() ) return convole(g, dom_g, f, dom_f);
115
116 Piecewise<SBasis> result;
117
118 SBasisOf<SBasisOf<double> > u,v;
119 u.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,1))));
120 v.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,0)),
121 SBasisOf<double>(LinearOf<double>(1,1))));
122 SBasisOf<SBasisOf<double> > v_u = (v - u)*(dom_f.extent()/dom_g.extent());
123 v_u += SBasisOf<SBasisOf<double> >(SBasisOf<double>(-dom_g.min()/dom_g.extent()));
124 SBasisOf<SBasisOf<double> > gg = multi_compose(g,v_u);
125 SBasisOf<SBasisOf<double> > ff = SBasisOf<SBasisOf<double> >(f);
126 SBasisOf<SBasisOf<double> > hh = integral(ff*gg,0);
127
128 result.cuts.push_back(dom_f.min()+dom_g.min());
129 //Note: we know dom_f.extent() >= dom_g.extent()!!
130 //double rho = dom_f.extent()/dom_g.extent();
131 double t0 = dom_g.min()/dom_f.extent();
132 double t1 = dom_g.max()/dom_f.extent();
133 double t2 = t0+1;
134 double t3 = t1+1;
135 SBasisOf<double> a,b,t;
136 SBasis seg;
137 a = SBasisOf<double>(LinearOf<double>(0,0));
138 b = SBasisOf<double>(LinearOf<double>(0,t1-t0));
139 t = SBasisOf<double>(LinearOf<double>(t0,t1));
140 seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
141 result.push(seg,dom_f.min() + dom_g.max());
142 if (dom_f.extent() > dom_g.extent()){
143 a = SBasisOf<double>(LinearOf<double>(0,t2-t1));
144 b = SBasisOf<double>(LinearOf<double>(t1-t0,1));
145 t = SBasisOf<double>(LinearOf<double>(t1,t2));
146 seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
147 result.push(seg,dom_f.max() + dom_g.min());
148 }
149 a = SBasisOf<double>(LinearOf<double>(t2-t1,1.));
150 b = SBasisOf<double>(LinearOf<double>(1.,1.));
151 t = SBasisOf<double>(LinearOf<double>(t2,t3));
152 seg = toSBasis(compose(hh,b,t)-compose(hh,a,t));
153 result.push(seg,dom_f.max() + dom_g.max());
154 result*=dom_f.extent();
155
156 //--constant ends correction--------------
157 if (f_cst_ends){
158 SBasis ig = toSBasis(integraaal(g))*dom_g.extent();
159 ig -= ig.at0();
160 Piecewise<SBasis> cor;
161 cor.cuts.push_back(dom_f.min()+dom_g.min());
162 cor.push(reverse(ig)*f.at0(),dom_f.min()+dom_g.max());
163 cor.push(Linear(0),dom_f.max()+dom_g.min());
164 cor.push(ig*f.at1(),dom_f.max()+dom_g.max());
165 result+=cor;
166 }
167 //----------------------------------------
168 return result;
169}
170
171/*static void dot_plot(cairo_t *cr, Piecewise<D2<SBasis> > const &M, double space=10){
172 //double dt=(M[0].cuts.back()-M[0].cuts.front())/space;
173 Piecewise<D2<SBasis> > Mperp = rot90(derivative(M)) * 2;
174 for( double t = M.cuts.front(); t < M.cuts.back(); t += space) {
175 Point pos = M(t), perp = Mperp(t);
176 draw_line_seg(cr, pos + perp, pos - perp);
177 }
178 cairo_pw_d2_sb(cr, M);
179 cairo_stroke(cr);
180 }*/
181
182static void plot_graph(cairo_t *cr, Piecewise<SBasis> const &f,
183 double x_scale=300,
184 double y_scale=100){
185 //double dt=(M[0].cuts.back()-M[0].cuts.front())/space;
186 D2<Piecewise<SBasis> > g;
187 g[X] = Piecewise<SBasis>( SBasis(Linear(100+f.cuts.front()*x_scale,
188 100+f.cuts.back()*x_scale)));
189 g[X].setDomain(f.domain());
190 g[Y] = -f*y_scale+400;
191 cairo_d2_pw_sb(cr, g);
192}
193
194struct Frame
195{
196 Geom::Point O;
197 Geom::Point x;
198 Geom::Point y;
199 Geom::Point z;
200};
201
202void
203plot3d(cairo_t *cr, SBasis const &x, SBasis const &y, SBasis const &z, Frame frame){
204 D2<SBasis> curve;
205 for (unsigned dim=0; dim<2; dim++){
206 curve[dim] = x*frame.x[dim] + y*frame.y[dim] + z*frame.z[dim];
207 curve[dim] += frame.O[dim];
208 }
209 cairo_d2_sb(cr, curve);
210}
211void
212plot3d(cairo_t *cr, SBasisOf<SBasisOf<double> > const &f, Frame frame, int NbRays=5){
213 for (int i=0; i<=NbRays; i++){
214 SBasisOf<double> var = LinearOf<double>(0,1);
215 SBasisOf<double> cst = LinearOf<double>(i*1./NbRays,i*1./NbRays);
216 SBasis f_on_seg = toSBasis(compose(f,var,cst));
217 plot3d(cr,Linear(0,1),Linear(i*1./NbRays),f_on_seg,frame);
218 f_on_seg = toSBasis(compose(f,cst,var));
219 plot3d(cr,Linear(i*1./NbRays),Linear(0,1),f_on_seg,frame);
220 }
221}
222
223
224#define SIZE 2
225
226class ConvolutionToy: public Toy {
227
228 PointHandle adjuster, adjuster2, adjuster3;
229
230public:
231 PointSetHandle b1_handle;
232 PointSetHandle b2_handle;
233 void draw(cairo_t *cr,
234 std::ostringstream *notify,
235 int width, int height, bool save, std::ostringstream *timer_stream) override {
236
237 D2<SBasis> B1 = b1_handle.asBezier();
238 D2<SBasis> B2 = b2_handle.asBezier();
239 D2<Piecewise<SBasis> >B;
240 B[X].concat(Piecewise<SBasis>(B1[X]));
241 B[X].concat(Piecewise<SBasis>(B2[X]));
242 B[Y].concat(Piecewise<SBasis>(B1[Y]));
243 B[Y].concat(Piecewise<SBasis>(B2[Y]));
244
245//-----------------------------------------------------
246#if 0
247 Frame frame;
248 frame.O = Point(50,400);
249 frame.x = Point(300,0);
250 frame.y = Point(120,-75);
251 frame.z = Point(0,-300);
252 SBasisOf<SBasisOf<double> > u,v,cst;
253 cst.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(1,1))));
254 u.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,1))));
255 v.push_back(LinearOf<SBasisOf<double> >(SBasisOf<double>(LinearOf<double>(0,0)),
256 SBasisOf<double>(LinearOf<double>(1,1))));
257 plot3d(cr, integral(v,1) ,frame);
258 plot3d(cr, v ,frame);
259 cairo_set_source_rgba (cr, 0., 0.5, 0., 1);
260 cairo_stroke(cr);
261 plot3d(cr, Linear(0,1), Linear(0), Linear(0), frame);
262 plot3d(cr, Linear(0), Linear(0,1), Linear(0), frame);
263 plot3d(cr, Linear(0), Linear(0), Linear(0,1), frame);
264 plot3d(cr, Linear(0,1), Linear(1), Linear(0), frame);
265 plot3d(cr, Linear(1), Linear(0,1), Linear(0), frame);
266 cairo_set_source_rgba (cr, 0., 0.0, 0.9, 1);
267 cairo_stroke(cr);
268#endif
269//-----------------------------------------------------
270
271 SBasisOf<double> smoother;
272 smoother.push_back(LinearOf<double>(0.));
273 smoother.push_back(LinearOf<double>(0.));
274 smoother.push_back(LinearOf<double>(30.));//could be less degree hungry!
275
276 adjuster.pos[X] = 400;
277 if(adjuster.pos[Y]>400) adjuster.pos[Y] = 400;
278 if(adjuster.pos[Y]<100) adjuster.pos[Y] = 100;
279 double scale=(400.-adjuster.pos[Y])/300+.01;
280 D2<Piecewise<SBasis> > smoothB1,smoothB2;
281 smoothB1[X] = convole(toSBasisOfDouble(B1[X]),Interval(0,4),smoother/scale,Interval(-scale/2,scale/2));
282 smoothB1[Y] = convole(toSBasisOfDouble(B1[Y]),Interval(0,4),smoother/scale,Interval(-scale/2,scale/2));
283 smoothB1[X].push(Linear(0.), 8+scale/2);
284 smoothB1[Y].push(Linear(0.), 8+scale/2);
285 smoothB2[X] = Piecewise<SBasis>(Linear(0));
286 smoothB2[X].setDomain(Interval(-scale/2,4-scale/2));
287 smoothB2[Y] = smoothB2[X];
288 smoothB2[X].concat(convole(toSBasisOfDouble(B2[X]),Interval(4,8),smoother/scale,Interval(-scale/2,scale/2)));
289 smoothB2[Y].concat(convole(toSBasisOfDouble(B2[Y]),Interval(4,8),smoother/scale,Interval(-scale/2,scale/2)));
290
291 Piecewise<D2<SBasis> > smoothB;
292 smoothB = sectionize(smoothB1)+sectionize(smoothB2);
293 //cairo_d2_sb(cr, B1);
294 //cairo_d2_sb(cr, B2);
295 cairo_pw_d2_sb(cr, smoothB);
296
297 cairo_move_to(cr,100,400);
298 cairo_line_to(cr,500,400);
299 cairo_set_line_width (cr, .5);
300 cairo_set_source_rgba (cr, 0., 0., 0., 1);
301 cairo_stroke(cr);
302
303 Piecewise<SBasis>bx = Piecewise<SBasis>(B1[X]);
304 bx.setDomain(Interval(0,4));
305 Piecewise<SBasis>smth = Piecewise<SBasis>(toSBasis(smoother));
306 smth.setDomain(Interval(-scale/2,scale/2));
307 cairo_d2_sb(cr, B1);
308 plot_graph(cr, bx, 100, 1);
309 plot_graph(cr, smth/scale, 100, 100);
310 plot_graph(cr, smoothB1[X],100,1);
311
312 //cairo_pw_d2_sb(cr, Piecewise<D2<SBasis> >(B1) );
313 //cairo_pw_d2_sb(cr, sectionize(smoothB1));
314
315 cairo_set_line_width (cr, .5);
316 cairo_set_source_rgba (cr, 0.7, 0.2, 0., 1);
317 cairo_stroke(cr);
318
319
320 Toy::draw(cr, notify, width, height, save,timer_stream);
321 }
322
323public:
324 ConvolutionToy(){
325 for(int i = 0; i < SIZE; i++) {
326 b1_handle.push_back(150+uniform()*300,150+uniform()*300);
327 b2_handle.push_back(150+uniform()*300,150+uniform()*300);
328 }
329 handles.push_back(&b1_handle);
330 handles.push_back(&b2_handle);
331
332 adjuster.pos = Geom::Point(400,100+300*uniform());
333 handles.push_back(&adjuster);
334
335 }
336};
337
338int main(int argc, char **argv) {
339 init(argc, argv, new ConvolutionToy);
340 return 0;
341}
342
343/*
344 Local Variables:
345 mode:c++
346 c-file-style:"stroustrup"
347 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
348 indent-tabs-mode:nil
349 fill-column:99
350 End:
351*/
352//vim:filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99:
scale
double scale
Definition aa.cpp:228
main
int main()
Definition bendpath-test.cpp:108
bezier-to-sbasis.h
Conversion between Bezier control points and SBasis curves.
SIZE
static int constexpr SIZE
Definition cielab.cpp:45
Geom::D2
Adaptor that creates 2D functions from 1D ones.
Definition d2.h:55
Geom::GenericInterval::extent
constexpr C extent() const
Definition generic-interval.h:99
Geom::GenericInterval::min
constexpr C min() const
Definition generic-interval.h:97
Geom::GenericInterval::max
constexpr C max() const
Definition generic-interval.h:98
Geom::Interval
Range of real numbers that is never empty.
Definition interval.h:59
Geom::Linear
Function that interpolates linearly between two values.
Definition linear.h:55
Geom::Piecewise
Function defined as discrete pieces.
Definition piecewise.h:71
Geom::Piecewise::domain
Interval domain() const
Definition piecewise.h:215
Geom::Piecewise::push
void push(const T &s, double to)
Convenience/implementation hiding function to add segment/cut pairs.
Definition piecewise.h:141
Geom::Piecewise::cuts
std::vector< double > cuts
Definition piecewise.h:75
Geom::Piecewise::setDomain
void setDomain(Interval dom)
Definition piecewise.h:218
Geom::Point
Two-dimensional point that doubles as a vector.
Definition point.h:66
Geom::SBasisOf::at1
T at1() const
Definition sbasis-of.h:99
Geom::SBasisOf::at0
T at0() const
Definition sbasis-of.h:96
Geom::SBasis
Polynomial in symmetric power basis.
Definition sbasis.h:70
Geom::SBasis::normalize
void normalize()
Definition sbasis.h:247
Geom::SBasis::resize
void resize(unsigned n)
Definition sbasis.h:98
Geom::SBasis::at0
Coord at0() const
Definition sbasis.h:212
PointHandle::pos
Geom::Point pos
Definition toy-framework-2.h:171
PointSetHandle::push_back
void push_back(double x, double y)
Definition toy-framework-2.h:188
PointSetHandle::asBezier
Geom::D2< Geom::SBasis > asBezier()
Definition toy-framework-2.cpp:913
Toy::handles
vector< Handle * > handles
Definition toy-framework-2.h:317
Toy::save
virtual void save(FILE *f)
Definition toy-framework-2.cpp:270
Toy::draw
virtual void draw(cairo_t *cr, std::ostringstream *notify, int w, int h, bool save, std::ostringstream *timing_stream)
Definition toy-framework-2.cpp:166
multi_compose
SBasisOf< T > multi_compose(SBasisOf< double > const &f, SBasisOf< T > const &g)
Definition convole.cpp:82
toSBasisOfDouble
SBasisOf< double > toSBasisOfDouble(SBasis const &f)
Definition convole.cpp:22
integraaal
SBasisOf< T > integraaal(SBasisOf< T > const &c)
Definition convole.cpp:51
plot3d
void plot3d(cairo_t *cr, SBasis const &x, SBasis const &y, SBasis const &z, Frame frame)
Definition convole.cpp:203
plot_graph
static void plot_graph(cairo_t *cr, Piecewise< SBasis > const &f, double x_scale=300, double y_scale=100)
Definition convole.cpp:182
convole
Piecewise< SBasis > convole(SBasisOf< double > const &f, Interval dom_f, SBasisOf< double > const &g, Interval dom_g, bool f_cst_ends=false)
Definition convole.cpp:110
result
Css & result
Definition ctrl-handle-styling.cpp:32
c
double c[8][4]
Definition cylinder3d.cpp:40
d2.h
Lifts one dimensional objects into 2D.
Geom::Y
@ Y
Definition coord.h:48
Geom::X
@ X
Definition coord.h:48
Geom::SL::v
size_t v
Definition multi-index.h:105
Geom
Various utility functions.
Definition affine.h:22
Geom::reverse
Bezier reverse(const Bezier &a)
Definition bezier.h:342
Geom::multiply
Bezier multiply(Bezier const &f, Bezier const &g)
Definition bezier.h:337
Geom::sectionize
Piecewise< D2< SBasis > > sectionize(D2< Piecewise< SBasis > > const &a)
Definition d2-sbasis.cpp:65
Geom::compose
D2< T > compose(D2< T > const &a, T const &b)
Definition d2.h:405
Geom::integral
Bezier integral(Bezier const &a)
Definition bezier.cpp:294
Geom::toSBasis
SBasis toSBasis(SBasisN< 1 > f)
Definition sbasisN.h:640
cairo_line_to
void cairo_line_to(cairo_t *cr, Geom::Point p1)
Definition cairo-helpers.cpp:19
cairo_t
struct _cairo cairo_t
Definition path-cairo.h:16
cairo_d2_pw_sb
void cairo_d2_pw_sb(cairo_t *cr, Geom::D2< Geom::Piecewise< Geom::SBasis > > const &p)
Definition cairo-helpers.cpp:127
cairo_move_to
void cairo_move_to(cairo_t *cr, Geom::Point p1)
Definition cairo-helpers.cpp:14
cairo_d2_sb
void cairo_d2_sb(cairo_t *cr, Geom::D2< Geom::SBasis > const &p)
Definition cairo-helpers.cpp:123
cairo_pw_d2_sb
void cairo_pw_d2_sb(cairo_t *cr, Geom::Piecewise< Geom::D2< Geom::SBasis > > const &p)
Definition cairo-helpers.cpp:131
sbasis-geometric.h
two-dimensional geometric operators.
sbasis-of.h
Defines S-power basis function class with coefficient in arbitrary ring.
sbasis.h
Polynomial in symmetric power basis (S-basis)
curve
Definition curve.h:24
height
double height
Definition convex_hull.cpp:107
width
double width
Definition convex_hull.cpp:106
uniform
double uniform()
Definition toy-framework-2.cpp:30
cairo_set_source_rgba
void cairo_set_source_rgba(cairo_t *cr, colour c)
Definition toy-framework-2.cpp:119
init
void init(int argc, char **argv, Toy *t, int width=600, int height=600)
Definition toy-framework-2.cpp:600