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