1 // Aseprite Document Library
2 // Copyright (c) 2001-2018 David Capello
3 //
4 // This file is released under the terms of the MIT license.
5 // Read LICENSE.txt for more information.
6
7 #ifdef HAVE_CONFIG_H
8 #include "config.h"
9 #endif
10
11 #include "doc/algo.h"
12
13 #include "base/base.h"
14 #include "base/debug.h"
15
16 #include <algorithm>
17 #include <cmath>
18 #include <utility>
19 #include <vector>
20
21 namespace doc {
22
algo_line(int x1,int y1,int x2,int y2,void * data,AlgoPixel proc)23 void algo_line(int x1, int y1, int x2, int y2, void* data, AlgoPixel proc)
24 {
25 bool yaxis;
26
27 // If the height if the line is bigger than the width, we'll iterate
28 // over the y-axis.
29 if (ABS(y2-y1) > ABS(x2-x1)) {
30 std::swap(x1, y1);
31 std::swap(x2, y2);
32 yaxis = true;
33 }
34 else
35 yaxis = false;
36
37 const int w = ABS(x2-x1)+1;
38 const int h = ABS(y2-y1)+1;
39 const int dx = SGN(x2-x1);
40 const int dy = SGN(y2-y1);
41
42 int e = 0;
43 int y = y1;
44
45 // Move x2 one extra pixel to the dx direction so we can use
46 // operator!=() instead of operator<(). Here I prefer operator!=()
47 // instead of swapping x1 with x2 so the error always start from 0
48 // in the origin (x1,y1).
49 x2 += dx;
50
51 for (int x=x1; x!=x2; x+=dx) {
52 if (yaxis)
53 proc(y, x, data);
54 else
55 proc(x, y, data);
56
57 // The error advances "h/w" per each "x" step. As we're using a
58 // integer value for "e", we use "w" as the unit.
59 e += h;
60 if (e >= w) {
61 y += dy;
62 e -= w;
63 }
64 }
65 }
66
67 /* Additional helper functions for the ellipse-drawing helper function
68 below corresponding to routines in Bresenham's algorithm. */
69 namespace {
bresenham_ellipse_error(int rx,int ry,int x,int y)70 int bresenham_ellipse_error(int rx, int ry, int x, int y) {
71 return x*x*ry*ry + y*y*rx*rx - rx*rx*ry*ry;
72 }
73
74 // Initialize positions x and y for Bresenham's algorithm
bresenham_ellipse_init(int rx,int ry,int & x,int & y)75 void bresenham_ellipse_init(int rx, int ry, int& x, int& y) {
76 // Start at the fatter pole
77 if (rx > ry) {
78 x = 0;
79 y = ry;
80 }
81 else {
82 x = rx;
83 y = 0;
84 }
85 }
86
87 // Move to next pixel to draw, according to Bresenham's algorithm
bresenham_ellipse_step(int rx,int ry,int & x,int & y)88 void bresenham_ellipse_step(int rx, int ry, int& x, int& y) {
89 // Move towards the skinnier pole. Having 2 cases isn't needed, but it ensures
90 // swapping rx and ry is the same as rotating 90 degrees.
91 if (rx > ry) {
92 int ex = bresenham_ellipse_error(rx, ry, x, y-1);
93 int ey = bresenham_ellipse_error(rx, ry, x+1, y);
94 int exy = bresenham_ellipse_error(rx, ry, x+1, y-1);
95 if (ex + exy < 0) ++x;
96 if (ey + exy > 0) --y;
97 }
98 else {
99 int ex = bresenham_ellipse_error(rx, ry, x, y+1);
100 int ey = bresenham_ellipse_error(rx, ry, x-1, y);
101 int exy = bresenham_ellipse_error(rx, ry, x-1, y+1);
102 if (ex + exy > 0) --x;
103 if (ey + exy < 0) ++y;
104 }
105 }
106 }
107
108 /* Helper function for the ellipse drawing routines. Calculates the
109 points of an ellipse which fits onto the rectangle specified by x1,
110 y1, x2 and y2, and calls the specified routine for each one. The
111 output proc has the same format as for do_line, and if the width or
112 height of the ellipse is only 1 or 2 pixels, do_line will be
113 called.
114
115 Copyright (C) 2002 by Elias Pschernig (eliaspschernig@aon.at)
116 for Allegro 4.x.
117
118 Adapted for ASEPRITE by David A. Capello. */
algo_ellipse(int x1,int y1,int x2,int y2,void * data,AlgoPixel proc)119 void algo_ellipse(int x1, int y1, int x2, int y2, void *data, AlgoPixel proc)
120 {
121 int mx, my, rx, ry;
122 int x, y;
123
124 /* Cheap hack to get elllipses with integer diameter, by just offsetting
125 * some quadrants by one pixel. */
126 int mx2, my2;
127
128 mx = (x1 + x2) / 2;
129 mx2 = (x1 + x2 + 1) / 2;
130 my = (y1 + y2) / 2;
131 my2 = (y1 + y2 + 1) / 2;
132 rx = ABS(x1 - x2);
133 ry = ABS(y1 - y2);
134
135 if (rx == 1) { algo_line(x2, y1, x2, y2, data, proc); rx--; }
136 if (rx == 0) { algo_line(x1, y1, x1, y2, data, proc); return; }
137
138 if (ry == 1) { algo_line(x1, y2, x2, y2, data, proc); ry--; }
139 if (ry == 0) { algo_line(x1, y1, x2, y1, data, proc); return; }
140
141 rx /= 2;
142 ry /= 2;
143
144 /* Draw the 4 poles. */
145 proc(mx, my2 + ry, data);
146 proc(mx, my - ry, data);
147 proc(mx2 + rx, my, data);
148 proc(mx - rx, my, data);
149
150 /* For even diameter axis, double the poles. */
151 if (mx != mx2) {
152 proc(mx2, my2 + ry, data);
153 proc(mx2, my - ry, data);
154 }
155
156 if (my != my2) {
157 proc(mx2 + rx, my2, data);
158 proc(mx - rx, my2, data);
159 }
160
161 /* Initialize drawing position at a pole. */
162 bresenham_ellipse_init(rx, ry, x, y);
163
164 for (;;) {
165 /* Step to the next pixel to draw. */
166 bresenham_ellipse_step(rx, ry, x, y);
167
168 /* Edge conditions */
169 if (y == 0 && x < rx) ++y; // don't move to horizontal radius except at pole
170 if (x == 0 && y < ry) ++x; // don't move to vertical radius except at pole
171 if (y <= 0 || x <= 0) break; // stop before pole, since it's already drawn
172
173 /* Process pixel */
174 proc(mx2 + x, my - y, data);
175 proc(mx - x, my - y, data);
176 proc(mx2 + x, my2 + y, data);
177 proc(mx - x, my2 + y, data);
178 }
179 }
180
algo_ellipsefill(int x1,int y1,int x2,int y2,void * data,AlgoHLine proc)181 void algo_ellipsefill(int x1, int y1, int x2, int y2, void *data, AlgoHLine proc)
182 {
183 int mx, my, rx, ry;
184 int x, y;
185
186 /* Cheap hack to get elllipses with integer diameter, by just offsetting
187 * some quadrants by one pixel. */
188 int mx2, my2;
189
190 mx = (x1 + x2) / 2;
191 mx2 = (x1 + x2 + 1) / 2;
192 my = (y1 + y2) / 2;
193 my2 = (y1 + y2 + 1) / 2;
194 rx = ABS (x1 - x2);
195 ry = ABS (y1 - y2);
196
197 if (rx == 1) { int c; for (c=y1; c<=y2; c++) proc(x2, c, x2, data); rx--; }
198 if (rx == 0) { int c; for (c=y1; c<=y2; c++) proc(x1, c, x1, data); return; }
199
200 if (ry == 1) { proc(x1, y2, x2, data); ry--; }
201 if (ry == 0) { proc(x1, y1, x2, data); return; }
202
203 rx /= 2;
204 ry /= 2;
205
206 /* Draw the north and south poles (possibly 2 pixels) */
207 proc(mx, my2 + ry, mx2, data);
208 proc(mx, my - ry, mx2, data);
209
210 /* Draw the equator (possibly width 2) */
211 proc(mx - rx, my, mx2 + rx, data);
212 if (my != my2) proc(mx - rx, my2, mx2 + rx, data);
213
214 /* Initialize drawing position at a pole. */
215 bresenham_ellipse_init(rx, ry, x, y);
216
217 for (;;) {
218 /* Step to the next pixel to draw. */
219 bresenham_ellipse_step(rx, ry, x, y);
220
221 /* Edge conditions */
222 if (y == 0 && x < rx) ++y; // don't move to horizontal radius except at pole
223 if (x == 0 && y < ry) ++x; // don't move to vertical radius except at pole
224 if (y <= 0 || x <= 0) break; // stop before pole, since it's already drawn
225
226 /* Draw the north and south 'lines of latitude' */
227 proc(mx - x, my - y, mx2 + x, data);
228 proc(mx - x, my2 + y, mx2 + x, data);
229 }
230 }
231
232 // Rotated ellipse code based on Alois Zingl work released under the
233 // MIT license http://members.chello.at/easyfilter/bresenham.html
234 //
235 // Adapted for Aseprite by David Capello
236
draw_quad_rational_bezier_seg(int x0,int y0,int x1,int y1,int x2,int y2,double w,void * data,AlgoPixel proc)237 static void draw_quad_rational_bezier_seg(int x0, int y0,
238 int x1, int y1,
239 int x2, int y2,
240 double w,
241 void* data,
242 AlgoPixel proc)
243 { // plot a limited rational Bezier segment, squared weight
244 int sx = x2-x1; // relative values for checks
245 int sy = y2-y1;
246 int dx = x0-x2;
247 int dy = y0-y2;
248 int xx = x0-x1;
249 int yy = y0-y1;
250 double xy = xx*sy + yy*sx;
251 double cur = xx*sy - yy*sx; // curvature
252 double err;
253
254 ASSERT(xx*sx <= 0.0 && yy*sy <= 0.0); // sign of gradient must not change
255
256 if (cur != 0.0 && w > 0.0) { // no straight line
257 if (sx*sx+sy*sy > xx*xx+yy*yy) { // begin with shorter part
258 // swap P0 P2
259 x2 = x0;
260 x0 -= dx;
261 y2 = y0;
262 y0 -= dy;
263 cur = -cur;
264 }
265 xx = 2.0*(4.0*w*sx*xx+dx*dx); // differences 2nd degree
266 yy = 2.0*(4.0*w*sy*yy+dy*dy);
267 sx = x0 < x2 ? 1 : -1; // x step direction
268 sy = y0 < y2 ? 1 : -1; // y step direction
269 xy = -2.0*sx*sy*(2.0*w*xy+dx*dy);
270
271 if (cur*sx*sy < 0.0) { // negated curvature?
272 xx = -xx; yy = -yy; xy = -xy; cur = -cur;
273 }
274 dx = 4.0*w*(x1-x0)*sy*cur+xx/2.0+xy; // differences 1st degree
275 dy = 4.0*w*(y0-y1)*sx*cur+yy/2.0+xy;
276
277 if (w < 0.5 && (dy > xy || dx < xy)) { // flat ellipse, algorithm fails
278 cur = (w+1.0)/2.0;
279 w = std::sqrt(w);
280 xy = 1.0/(w+1.0);
281
282 sx = std::floor((x0+2.0*w*x1+x2)*xy/2.0+0.5); // subdivide curve in half
283 sy = std::floor((y0+2.0*w*y1+y2)*xy/2.0+0.5);
284
285 dx = std::floor((w*x1+x0)*xy+0.5);
286 dy = std::floor((y1*w+y0)*xy+0.5);
287 draw_quad_rational_bezier_seg(x0, y0, dx, dy, sx, sy, cur, data, proc); // plot separately
288
289 dx = std::floor((w*x1+x2)*xy+0.5);
290 dy = std::floor((y1*w+y2)*xy+0.5);
291 draw_quad_rational_bezier_seg(sx, sy, dx, dy, x2, y2, cur, data, proc);
292 return;
293 }
294 err = dx+dy-xy; // error 1.step
295 do {
296 // plot curve
297 proc(x0, y0, data);
298
299 if (x0 == x2 && y0 == y2)
300 return; // last pixel -> curve finished
301
302 x1 = 2*err > dy;
303 y1 = 2*(err+yy) < -dy; // save value for test of x step
304
305 if (2*err < dx || y1) {
306 // y step
307 y0 += sy;
308 dy += xy;
309 err += dx += xx;
310 }
311 if (2*err > dx || x1) {
312 // x step
313 x0 += sx;
314 dx += xy;
315 err += dy += yy;
316 }
317 } while (dy <= xy && dx >= xy); // gradient negates -> algorithm fails
318 }
319 algo_line(x0, y0, x2, y2, data, proc); // plot remaining needle to end
320 }
321
draw_rotated_ellipse_rect(int x0,int y0,int x1,int y1,double zd,void * data,AlgoPixel proc)322 static void draw_rotated_ellipse_rect(int x0, int y0, int x1, int y1, double zd, void* data, AlgoPixel proc)
323 {
324 int xd = x1-x0;
325 int yd = y1-y0;
326 double w = xd*yd;
327
328 if (zd == 0)
329 return algo_ellipse(x0, y0, x1, y1, data, proc);
330
331 if (w != 0.0)
332 w = (w-zd) / (w+w); // squared weight of P1
333
334 w = MID(0.0, w, 1.0);
335
336 xd = std::floor(w*xd + 0.5);
337 yd = std::floor(w*yd + 0.5);
338
339 draw_quad_rational_bezier_seg(x0, y0+yd, x0, y0, x0+xd, y0, 1.0-w, data, proc);
340 draw_quad_rational_bezier_seg(x0, y0+yd, x0, y1, x1-xd, y1, w, data, proc);
341 draw_quad_rational_bezier_seg(x1, y1-yd, x1, y1, x1-xd, y1, 1.0-w, data, proc);
342 draw_quad_rational_bezier_seg(x1, y1-yd, x1, y0, x0+xd, y0, w, data, proc);
343 }
344
draw_rotated_ellipse(int cx,int cy,int a,int b,double angle,void * data,AlgoPixel proc)345 void draw_rotated_ellipse(int cx, int cy, int a, int b, double angle, void* data, AlgoPixel proc)
346 {
347 double xd = a*a;
348 double yd = b*b;
349 double s = std::sin(angle);
350 double zd = (xd-yd)*s; // ellipse rotation
351 xd = std::sqrt(xd-zd*s); // surrounding rectangle
352 yd = std::sqrt(yd+zd*s);
353
354 a = std::floor(xd+0.5);
355 b = std::floor(yd+0.5);
356 zd = zd*a*b / (xd*yd);
357 zd = 4*zd*std::cos(angle);
358
359 draw_rotated_ellipse_rect(cx-a, cy-b, cx+a, cy+b, zd, data, proc);
360 }
361
fill_rotated_ellipse(int cx,int cy,int a,int b,double angle,void * data,AlgoHLine proc)362 void fill_rotated_ellipse(int cx, int cy, int a, int b, double angle, void* data, AlgoHLine proc)
363 {
364 struct Rows {
365 int y0;
366 std::vector<std::pair<int, int> > row;
367 Rows(int y0, int nrows)
368 : y0(y0), row(nrows, std::make_pair(1, -1)) { }
369 void update(int x, int y) {
370 int i = MID(0, y-y0, row.size()-1);
371 auto& r = row[i];
372 if (r.first > r.second) {
373 r.first = r.second = x;
374 }
375 else {
376 r.first = MIN(r.first, x);
377 r.second = MAX(r.second, x);
378 }
379 }
380 };
381
382 double xd = a*a;
383 double yd = b*b;
384 double s = std::sin(angle);
385 double zd = (xd-yd)*s;
386 xd = std::sqrt(xd-zd*s);
387 yd = std::sqrt(yd+zd*s);
388
389 a = std::floor(xd+0.5);
390 b = std::floor(yd+0.5);
391 zd = zd*a*b / (xd*yd);
392 zd = 4*zd*std::cos(angle);
393
394 Rows rows(cy-b, 2*b+1);
395
396 draw_rotated_ellipse_rect(
397 cx-a, cy-b, cx+a, cy+b, zd,
398 &rows,
399 [](int x, int y, void* data){
400 Rows* rows = (Rows*)data;
401 rows->update(x, y);
402 });
403
404 int y = rows.y0;
405 for (const auto& r : rows.row) {
406 if (r.first <= r.second)
407 proc(r.first, y, r.second, data);
408 ++y;
409 }
410 }
411
412 // Algorightm from Allegro (allegro/src/spline.c)
413 // Adapted for Aseprite by David Capello.
algo_spline(double x0,double y0,double x1,double y1,double x2,double y2,double x3,double y3,void * data,AlgoLine proc)414 void algo_spline(double x0, double y0, double x1, double y1,
415 double x2, double y2, double x3, double y3,
416 void *data, AlgoLine proc)
417 {
418 int npts;
419 int out_x1, out_x2;
420 int out_y1, out_y2;
421
422 /* Derivatives of x(t) and y(t). */
423 double x, dx, ddx, dddx;
424 double y, dy, ddy, dddy;
425 int i;
426
427 /* Temp variables used in the setup. */
428 double dt, dt2, dt3;
429 double xdt2_term, xdt3_term;
430 double ydt2_term, ydt3_term;
431
432 #define MAX_POINTS 64
433 #undef DIST
434 #define DIST(x, y) (sqrt((x) * (x) + (y) * (y)))
435 npts = (int)(sqrt(DIST(x1-x0, y1-y0) +
436 DIST(x2-x1, y2-y1) +
437 DIST(x3-x2, y3-y2)) * 1.2);
438 if (npts > MAX_POINTS)
439 npts = MAX_POINTS;
440 else if (npts < 4)
441 npts = 4;
442
443 dt = 1.0 / (npts-1);
444 dt2 = (dt * dt);
445 dt3 = (dt2 * dt);
446
447 xdt2_term = 3 * (x2 - 2*x1 + x0);
448 ydt2_term = 3 * (y2 - 2*y1 + y0);
449 xdt3_term = x3 + 3 * (-x2 + x1) - x0;
450 ydt3_term = y3 + 3 * (-y2 + y1) - y0;
451
452 xdt2_term = dt2 * xdt2_term;
453 ydt2_term = dt2 * ydt2_term;
454 xdt3_term = dt3 * xdt3_term;
455 ydt3_term = dt3 * ydt3_term;
456
457 dddx = 6*xdt3_term;
458 dddy = 6*ydt3_term;
459 ddx = -6*xdt3_term + 2*xdt2_term;
460 ddy = -6*ydt3_term + 2*ydt2_term;
461 dx = xdt3_term - xdt2_term + 3 * dt * (x1 - x0);
462 dy = ydt3_term - ydt2_term + dt * 3 * (y1 - y0);
463 x = x0;
464 y = y0;
465
466 out_x1 = (int)x0;
467 out_y1 = (int)y0;
468
469 x += .5;
470 y += .5;
471 for (i=1; i<npts; i++) {
472 ddx += dddx;
473 ddy += dddy;
474 dx += ddx;
475 dy += ddy;
476 x += dx;
477 y += dy;
478
479 out_x2 = (int)x;
480 out_y2 = (int)y;
481
482 proc(out_x1, out_y1, out_x2, out_y2, data);
483
484 out_x1 = out_x2;
485 out_y1 = out_y2;
486 }
487 }
488
algo_spline_get_y(double x0,double y0,double x1,double y1,double x2,double y2,double x3,double y3,double in_x)489 double algo_spline_get_y(double x0, double y0, double x1, double y1,
490 double x2, double y2, double x3, double y3,
491 double in_x)
492 {
493 int npts;
494 double out_x, old_x;
495 double out_y, old_y;
496
497 /* Derivatives of x(t) and y(t). */
498 double x, dx, ddx, dddx;
499 double y, dy, ddy, dddy;
500 int i;
501
502 /* Temp variables used in the setup. */
503 double dt, dt2, dt3;
504 double xdt2_term, xdt3_term;
505 double ydt2_term, ydt3_term;
506
507 #define MAX_POINTS 64
508 #undef DIST
509 #define DIST(x, y) (sqrt ((x) * (x) + (y) * (y)))
510 npts = (int) (sqrt (DIST(x1-x0, y1-y0) +
511 DIST(x2-x1, y2-y1) +
512 DIST(x3-x2, y3-y2)) * 1.2);
513 if (npts > MAX_POINTS)
514 npts = MAX_POINTS;
515 else if (npts < 4)
516 npts = 4;
517
518 dt = 1.0 / (npts-1);
519 dt2 = (dt * dt);
520 dt3 = (dt2 * dt);
521
522 xdt2_term = 3 * (x2 - 2*x1 + x0);
523 ydt2_term = 3 * (y2 - 2*y1 + y0);
524 xdt3_term = x3 + 3 * (-x2 + x1) - x0;
525 ydt3_term = y3 + 3 * (-y2 + y1) - y0;
526
527 xdt2_term = dt2 * xdt2_term;
528 ydt2_term = dt2 * ydt2_term;
529 xdt3_term = dt3 * xdt3_term;
530 ydt3_term = dt3 * ydt3_term;
531
532 dddx = 6*xdt3_term;
533 dddy = 6*ydt3_term;
534 ddx = -6*xdt3_term + 2*xdt2_term;
535 ddy = -6*ydt3_term + 2*ydt2_term;
536 dx = xdt3_term - xdt2_term + 3 * dt * (x1 - x0);
537 dy = ydt3_term - ydt2_term + dt * 3 * (y1 - y0);
538 x = x0;
539 y = y0;
540
541 out_x = old_x = x0;
542 out_y = old_y = y0;
543
544 x += .5;
545 y += .5;
546 for (i=1; i<npts; i++) {
547 ddx += dddx;
548 ddy += dddy;
549 dx += ddx;
550 dy += ddy;
551 x += dx;
552 y += dy;
553
554 out_x = x;
555 out_y = y;
556 if (out_x > in_x) {
557 out_y = old_y + (out_y-old_y) * (in_x-old_x) / (out_x-old_x);
558 break;
559 }
560 old_x = out_x;
561 old_y = out_y;
562 }
563
564 return out_y;
565 }
566
algo_spline_get_tan(double x0,double y0,double x1,double y1,double x2,double y2,double x3,double y3,double in_x)567 double algo_spline_get_tan(double x0, double y0, double x1, double y1,
568 double x2, double y2, double x3, double y3,
569 double in_x)
570 {
571 double out_x, old_x, old_dx, old_dy;
572 int npts;
573
574 /* Derivatives of x(t) and y(t). */
575 double x, dx, ddx, dddx;
576 double y, dy, ddy, dddy;
577 int i;
578
579 /* Temp variables used in the setup. */
580 double dt, dt2, dt3;
581 double xdt2_term, xdt3_term;
582 double ydt2_term, ydt3_term;
583
584 #define MAX_POINTS 64
585 #undef DIST
586 #define DIST(x, y) (sqrt ((x) * (x) + (y) * (y)))
587 npts = (int) (sqrt (DIST(x1-x0, y1-y0) +
588 DIST(x2-x1, y2-y1) +
589 DIST(x3-x2, y3-y2)) * 1.2);
590 if (npts > MAX_POINTS)
591 npts = MAX_POINTS;
592 else if (npts < 4)
593 npts = 4;
594
595 dt = 1.0 / (npts-1);
596 dt2 = (dt * dt);
597 dt3 = (dt2 * dt);
598
599 xdt2_term = 3 * (x2 - 2*x1 + x0);
600 ydt2_term = 3 * (y2 - 2*y1 + y0);
601 xdt3_term = x3 + 3 * (-x2 + x1) - x0;
602 ydt3_term = y3 + 3 * (-y2 + y1) - y0;
603
604 xdt2_term = dt2 * xdt2_term;
605 ydt2_term = dt2 * ydt2_term;
606 xdt3_term = dt3 * xdt3_term;
607 ydt3_term = dt3 * ydt3_term;
608
609 dddx = 6*xdt3_term;
610 dddy = 6*ydt3_term;
611 ddx = -6*xdt3_term + 2*xdt2_term;
612 ddy = -6*ydt3_term + 2*ydt2_term;
613 dx = xdt3_term - xdt2_term + 3 * dt * (x1 - x0);
614 dy = ydt3_term - ydt2_term + dt * 3 * (y1 - y0);
615 x = x0;
616 y = y0;
617
618 out_x = x0;
619
620 old_x = x0;
621 old_dx = dx;
622 old_dy = dy;
623
624 x += .5;
625 y += .5;
626 for (i=1; i<npts; i++) {
627 ddx += dddx;
628 ddy += dddy;
629 dx += ddx;
630 dy += ddy;
631 x += dx;
632 y += dy;
633
634 out_x = x;
635 if (out_x > in_x) {
636 dx = old_dx + (dx-old_dx) * (in_x-old_x) / (out_x-old_x);
637 dy = old_dy + (dy-old_dy) * (in_x-old_x) / (out_x-old_x);
638 break;
639 }
640 old_x = out_x;
641 old_dx = dx;
642 old_dy = dy;
643 }
644
645 return dy / dx;
646 }
647
648 } // namespace doc
649