1// Aseprite Document Library
2// Copyright (c) 2018-2022 Igara Studio S.A.
3// Copyright (c) 2001-2018 David Capello
4//
5// This file is released under the terms of the MIT license.
6// Read LICENSE.txt for more information.
7
8#ifdef HAVE_CONFIG_H
9#include "config.h"
10#endif
11
12#include "doc/algo.h"
13
14#include "base/debug.h"
15
16#include <algorithm>
17#include <cmath>
18#include <utility>
19#include <vector>
20
21namespace doc {
22
23void algo_line_perfect(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// Special version of the perfect line algorithm specially done for
68// kLineBrushType so the whole line looks continuous without holes.
69//
70// TOOD in a future we should convert lines into scanlines and render
71// scanlines instead of drawing the brush on each pixel, that
72// would fix all cases
73void algo_line_perfect_with_fix_for_line_brush(int x1, int y1, int x2, int y2, void* data, AlgoPixel proc)
74{
75 bool yaxis;
76
77 if (ABS(y2-y1) > ABS(x2-x1)) {
78 std::swap(x1, y1);
79 std::swap(x2, y2);
80 yaxis = true;
81 }
82 else
83 yaxis = false;
84
85 const int w = ABS(x2-x1)+1;
86 const int h = ABS(y2-y1)+1;
87 const int dx = SGN(x2-x1);
88 const int dy = SGN(y2-y1);
89
90 int e = 0;
91 int y = y1;
92
93 x2 += dx;
94
95 for (int x=x1; x!=x2; x+=dx) {
96 if (yaxis)
97 proc(y, x, data);
98 else
99 proc(x, y, data);
100
101 e += h;
102 if (e >= w) {
103 y += dy;
104 e -= w;
105 if (x+dx != x2) {
106 if (yaxis)
107 proc(y, x, data);
108 else
109 proc(x, y, data);
110 }
111 }
112 }
113}
114
115// Line code based on Alois Zingl work released under the
116// MIT license http://members.chello.at/easyfilter/bresenham.html
117void algo_line_continuous(int x0, int y0, int x1, int y1, void* data, AlgoPixel proc)
118{
119 int dx = ABS(x1-x0), sx = (x0 < x1 ? 1: -1);
120 int dy = -ABS(y1-y0), sy = (y0 < y1 ? 1: -1);
121 int err = dx+dy, e2; // error value e_xy
122
123 for (;;) {
124 proc(x0, y0, data);
125 e2 = 2*err;
126 if (e2 >= dy) { // e_xy+e_x > 0
127 if (x0 == x1)
128 break;
129 err += dy;
130 x0 += sx;
131 }
132 if (e2 <= dx) { // e_xy+e_y < 0
133 if (y0 == y1)
134 break;
135 err += dx;
136 y0 += sy;
137 }
138 }
139}
140
141// Special version of the continuous line algorithm specially done for
142// kLineBrushType so the whole line looks continuous without holes.
143void algo_line_continuous_with_fix_for_line_brush(int x0, int y0, int x1, int y1, void* data, AlgoPixel proc)
144{
145 int dx = ABS(x1-x0), sx = (x0 < x1 ? 1: -1);
146 int dy = -ABS(y1-y0), sy = (y0 < y1 ? 1: -1);
147 int err = dx+dy, e2; // error value e_xy
148 bool x_changed;
149
150 for (;;) {
151 x_changed = false;
152
153 proc(x0, y0, data);
154 e2 = 2*err;
155 if (e2 >= dy) { // e_xy+e_x > 0
156 if (x0 == x1)
157 break;
158 err += dy;
159 x0 += sx;
160 x_changed = true;
161 }
162 if (e2 <= dx) { // e_xy+e_y < 0
163 if (y0 == y1)
164 break;
165 err += dx;
166 if (x_changed)
167 proc(x0, y0, data);
168 y0 += sy;
169 }
170 }
171}
172
173static int adjust_ellipse_args(int& x0, int& y0, int& x1, int& y1,
174 int& hPixels, int& vPixels)
175{
176 // hPixels : straight horizontal pixels added to mid region of the ellipse.
177 hPixels = std::max(hPixels, 0);
178 // vPixels : straight vertical pixels added to mid region of the ellipse.
179 vPixels = std::max(vPixels, 0);
180
181 // Conditioning swapped points
182 if (x0 > x1)
183 std::swap(x0, x1);
184 if (y0 > y1)
185 std::swap(y0, y1);
186 int w = x1 - x0 + 1;
187 int h = y1 - y0 + 1;
188
189 // hDiameter is the horizontal diameter of a circunference
190 // without the addition of straight pixels.
191 int hDiameter = w - hPixels;
192 // vDiameter is the vertical diameter of a circunference
193 // without the addition of straight pixels.
194 int vDiameter = h - vPixels;
195
196 // Manual adjustment
197 if (w == 8 || w == 12 || w == 22)
198 hPixels++;
199 if (h == 8 || h == 12 || h == 22)
200 vPixels++;
201
202 hPixels = (hDiameter > 5 ? hPixels : 0);
203 vPixels = (vDiameter > 5 ? vPixels : 0);
204
205 if ((hDiameter % 2 == 0) && (hDiameter > 5))
206 hPixels--;
207 if ((vDiameter % 2 == 0) && (vDiameter > 5))
208 vPixels--;
209
210 x1 -= hPixels;
211 y1 -= vPixels;
212
213 return h;
214}
215
216// Ellipse code based on Alois Zingl work released under the MIT
217// license http://members.chello.at/easyfilter/bresenham.html
218//
219// Adapted for Aseprite by David Capello
220
221void algo_ellipse(int x0, int y0, int x1, int y1,
222 int hPixels, int vPixels,
223 void* data, AlgoPixel proc)
224{
225 int h = adjust_ellipse_args(x0, y0, x1, y1, hPixels, vPixels);
226
227 long a = abs(x1-x0);
228 long b = abs(y1-y0); // diameter
229 long b1 = b&1;
230 double dx = 4*(1.0-a)*b*b; // error increment
231 double dy = 4*(b1+1)*a*a; // error increment
232 double err = dx + dy + b1*a*a; // error of 1.step
233 double e2;
234
235 y0 += (b+1)/2;
236 y1 = y0-b1; // starting pixel
237 a = 8*a*a;
238 b1 = 8*b*b;
239
240 int initialY0 = y0;
241 int initialY1 = y1;
242 int initialX0 = x0;
243 int initialX1 = x1 + hPixels;
244 do {
245 proc(x1 + hPixels, y0 + vPixels, data); // I. Quadrant
246 proc(x0, y0 + vPixels, data); // II. Quadrant
247 proc(x0, y1, data); // III. Quadrant
248 proc(x1 + hPixels, y1, data); // IV. Quadrant
249
250 e2 = 2*err;
251 if (e2 <= dy) { y0++; y1--; err += dy += a; } // y step
252 if (e2 >= dx || 2*err > dy) { x0++; x1--; err += dx += b1; } // x step
253 } while (x0 <= x1);
254
255 while (y0 + vPixels - y1 + 1 <= h) { // too early stop of flat ellipses a=1
256 proc(x0 - 1, y0 + vPixels, data); // -> finish tip of ellipse
257 proc(x1 + 1 + hPixels, y0++ + vPixels, data);
258 proc(x0 - 1, y1, data);
259 proc(x1 + 1 + hPixels, y1--, data);
260 }
261
262 // Extra horizontal straight pixels
263 if (hPixels > 0) {
264 for (int i = x0; i < x1 + hPixels + 1; i++) {
265 proc(i, y1 + 1, data);
266 proc(i, y0 + vPixels - 1, data);
267 }
268 }
269 // Extra vertical straight pixels
270 if (vPixels > 0) {
271 for (int i = initialY1 + 1; i < initialY0 + vPixels; i++) {
272 proc(initialX0, i, data);
273 proc(initialX1, i, data);
274 }
275 }
276}
277
278void algo_ellipsefill(int x0, int y0, int x1, int y1,
279 int hPixels, int vPixels,
280 void* data, AlgoHLine proc)
281{
282 int h = adjust_ellipse_args(x0, y0, x1, y1, hPixels, vPixels);
283
284 long a = abs(x1-x0), b = abs(y1-y0), b1 = b&1; // diameter
285 double dx = 4*(1.0-a)*b*b, dy = 4*(b1+1)*a*a; // error increment
286 double err = dx+dy+b1*a*a, e2; // error of 1.step
287
288 y0 += (b+1)/2; y1 = y0-b1; // starting pixel
289 a = 8*a*a; b1 = 8*b*b;
290
291 int initialY0 = y0;
292 int initialY1 = y1;
293 int initialX0 = x0;
294 int initialX1 = x1 + hPixels;
295
296 do {
297 proc(x0, y0 + vPixels, x1 + hPixels, data);
298 proc(x0, y1, x1 + hPixels, data);
299 e2 = 2*err;
300 if (e2 <= dy) { y0++; y1--; err += dy += a; } // y step
301 if (e2 >= dx || 2*err > dy) { x0++; x1--; err += dx += b1; } // x step
302 } while (x0 <= x1);
303
304 while (y0 + vPixels - y1 + 1 < h) { // too early stop of flat ellipses a=1
305 proc(x0-1, ++y0 + vPixels, x0-1, data); // -> finish tip of ellipse
306 proc(x1+1 + hPixels, y0 + vPixels, x1+1 + hPixels, data);
307 proc(x0-1, --y1, x0-1, data);
308 proc(x1+1 + hPixels, y1, x1+1 + hPixels, data);
309 }
310
311 if (vPixels > 0) {
312 for (int i = initialY1 + 1; i < initialY0 + vPixels; i++)
313 proc(initialX0, i, initialX1, data);
314 }
315}
316
317static void draw_quad_rational_bezier_seg(int x0, int y0,
318 int x1, int y1,
319 int x2, int y2,
320 double w,
321 void* data,
322 AlgoPixel proc)
323{ // plot a limited rational Bezier segment, squared weight
324 int sx = x2-x1; // relative values for checks
325 int sy = y2-y1;
326 int dx = x0-x2;
327 int dy = y0-y2;
328 int xx = x0-x1;
329 int yy = y0-y1;
330 double xy = xx*sy + yy*sx;
331 double cur = xx*sy - yy*sx; // curvature
332 double err;
333
334 ASSERT(xx*sx <= 0.0 && yy*sy <= 0.0); // sign of gradient must not change
335
336 if (cur != 0.0 && w > 0.0) { // no straight line
337 if (sx*sx+sy*sy > xx*xx+yy*yy) { // begin with shorter part
338 // swap P0 P2
339 x2 = x0;
340 x0 -= dx;
341 y2 = y0;
342 y0 -= dy;
343 cur = -cur;
344 }
345 xx = 2.0*(4.0*w*sx*xx+dx*dx); // differences 2nd degree
346 yy = 2.0*(4.0*w*sy*yy+dy*dy);
347 sx = x0 < x2 ? 1 : -1; // x step direction
348 sy = y0 < y2 ? 1 : -1; // y step direction
349 xy = -2.0*sx*sy*(2.0*w*xy+dx*dy);
350
351 if (cur*sx*sy < 0.0) { // negated curvature?
352 xx = -xx; yy = -yy; xy = -xy; cur = -cur;
353 }
354 dx = 4.0*w*(x1-x0)*sy*cur+xx/2.0+xy; // differences 1st degree
355 dy = 4.0*w*(y0-y1)*sx*cur+yy/2.0+xy;
356
357 if (w < 0.5 && (dy > xy || dx < xy)) { // flat ellipse, algorithm fails
358 cur = (w+1.0)/2.0;
359 w = std::sqrt(w);
360 xy = 1.0/(w+1.0);
361
362 sx = std::floor((x0+2.0*w*x1+x2)*xy/2.0+0.5); // subdivide curve in half
363 sy = std::floor((y0+2.0*w*y1+y2)*xy/2.0+0.5);
364
365 dx = std::floor((w*x1+x0)*xy+0.5);
366 dy = std::floor((y1*w+y0)*xy+0.5);
367 draw_quad_rational_bezier_seg(x0, y0, dx, dy, sx, sy, cur, data, proc); // plot separately
368
369 dx = std::floor((w*x1+x2)*xy+0.5);
370 dy = std::floor((y1*w+y2)*xy+0.5);
371 draw_quad_rational_bezier_seg(sx, sy, dx, dy, x2, y2, cur, data, proc);
372 return;
373 }
374 err = dx+dy-xy; // error 1.step
375 do {
376 // plot curve
377 proc(x0, y0, data);
378
379 if (x0 == x2 && y0 == y2)
380 return; // last pixel -> curve finished
381
382 x1 = 2*err > dy;
383 y1 = 2*(err+yy) < -dy; // save value for test of x step
384
385 if (2*err < dx || y1) {
386 // y step
387 y0 += sy;
388 dy += xy;
389 err += dx += xx;
390 }
391 if (2*err > dx || x1) {
392 // x step
393 x0 += sx;
394 dx += xy;
395 err += dy += yy;
396 }
397 } while (dy <= xy && dx >= xy); // gradient negates -> algorithm fails
398 }
399 algo_line_continuous(x0, y0, x2, y2, data, proc); // plot remaining needle to end
400}
401
402static void draw_rotated_ellipse_rect(int x0, int y0, int x1, int y1, double zd, void* data, AlgoPixel proc)
403{
404 int xd = x1-x0;
405 int yd = y1-y0;
406 double w = xd*yd;
407
408 if (zd == 0)
409 return algo_ellipse(x0, y0, x1, y1, 0, 0, data, proc);
410
411 if (w != 0.0)
412 w = (w-zd) / (w+w); // squared weight of P1
413
414 w = std::clamp(w, 0.0, 1.0);
415
416 xd = std::floor(w*xd + 0.5);
417 yd = std::floor(w*yd + 0.5);
418
419 draw_quad_rational_bezier_seg(x0, y0+yd, x0, y0, x0+xd, y0, 1.0-w, data, proc);
420 draw_quad_rational_bezier_seg(x0, y0+yd, x0, y1, x1-xd, y1, w, data, proc);
421 draw_quad_rational_bezier_seg(x1, y1-yd, x1, y1, x1-xd, y1, 1.0-w, data, proc);
422 draw_quad_rational_bezier_seg(x1, y1-yd, x1, y0, x0+xd, y0, w, data, proc);
423}
424
425void draw_rotated_ellipse(int cx, int cy, int a, int b, double angle, void* data, AlgoPixel proc)
426{
427 double xd = a*a;
428 double yd = b*b;
429 double s = std::sin(angle);
430 double zd = (xd-yd)*s; // ellipse rotation
431 xd = std::sqrt(xd-zd*s); // surrounding rectangle
432 yd = std::sqrt(yd+zd*s);
433
434 a = std::floor(xd+0.5);
435 b = std::floor(yd+0.5);
436 zd = zd*a*b / (xd*yd);
437 zd = 4*zd*std::cos(angle);
438
439 draw_rotated_ellipse_rect(cx-a, cy-b, cx+a, cy+b, zd, data, proc);
440}
441
442void fill_rotated_ellipse(int cx, int cy, int a, int b, double angle, void* data, AlgoHLine proc)
443{
444 struct Rows {
445 int y0;
446 std::vector<std::pair<int, int> > row;
447 Rows(int y0, int nrows)
448 : y0(y0), row(nrows, std::make_pair(1, -1)) { }
449 void update(int x, int y) {
450 int i = std::clamp(y-y0, 0, int(row.size()-1));
451 auto& r = row[i];
452 if (r.first > r.second) {
453 r.first = r.second = x;
454 }
455 else {
456 r.first = std::min(r.first, x);
457 r.second = std::max(r.second, x);
458 }
459 }
460 };
461
462 double xd = a*a;
463 double yd = b*b;
464 double s = std::sin(angle);
465 double zd = (xd-yd)*s;
466 xd = std::sqrt(xd-zd*s);
467 yd = std::sqrt(yd+zd*s);
468
469 a = std::floor(xd+0.5);
470 b = std::floor(yd+0.5);
471 zd = zd*a*b / (xd*yd);
472 zd = 4*zd*std::cos(angle);
473
474 Rows rows(cy-b, 2*b+1);
475
476 draw_rotated_ellipse_rect(
477 cx-a, cy-b, cx+a, cy+b, zd,
478 &rows,
479 [](int x, int y, void* data){
480 Rows* rows = (Rows*)data;
481 rows->update(x, y);
482 });
483
484 int y = rows.y0;
485 for (const auto& r : rows.row) {
486 if (r.first <= r.second)
487 proc(r.first, y, r.second, data);
488 ++y;
489 }
490}
491
492// Algorightm from Allegro (allegro/src/spline.c)
493// Adapted for Aseprite by David Capello.
494void algo_spline(double x0, double y0, double x1, double y1,
495 double x2, double y2, double x3, double y3,
496 void *data, AlgoLine proc)
497{
498 int npts;
499 int out_x1, out_x2;
500 int out_y1, out_y2;
501
502 /* Derivatives of x(t) and y(t). */
503 double x, dx, ddx, dddx;
504 double y, dy, ddy, dddy;
505 int i;
506
507 /* Temp variables used in the setup. */
508 double dt, dt2, dt3;
509 double xdt2_term, xdt3_term;
510 double ydt2_term, ydt3_term;
511
512#define MAX_POINTS 64
513#undef DIST
514#define DIST(x, y) (sqrt((x) * (x) + (y) * (y)))
515 npts = (int)(sqrt(DIST(x1-x0, y1-y0) +
516 DIST(x2-x1, y2-y1) +
517 DIST(x3-x2, y3-y2)) * 1.2);
518 if (npts > MAX_POINTS)
519 npts = MAX_POINTS;
520 else if (npts < 4)
521 npts = 4;
522
523 dt = 1.0 / (npts-1);
524 dt2 = (dt * dt);
525 dt3 = (dt2 * dt);
526
527 xdt2_term = 3 * (x2 - 2*x1 + x0);
528 ydt2_term = 3 * (y2 - 2*y1 + y0);
529 xdt3_term = x3 + 3 * (-x2 + x1) - x0;
530 ydt3_term = y3 + 3 * (-y2 + y1) - y0;
531
532 xdt2_term = dt2 * xdt2_term;
533 ydt2_term = dt2 * ydt2_term;
534 xdt3_term = dt3 * xdt3_term;
535 ydt3_term = dt3 * ydt3_term;
536
537 dddx = 6*xdt3_term;
538 dddy = 6*ydt3_term;
539 ddx = -6*xdt3_term + 2*xdt2_term;
540 ddy = -6*ydt3_term + 2*ydt2_term;
541 dx = xdt3_term - xdt2_term + 3 * dt * (x1 - x0);
542 dy = ydt3_term - ydt2_term + dt * 3 * (y1 - y0);
543 x = x0;
544 y = y0;
545
546 out_x1 = (int)x0;
547 out_y1 = (int)y0;
548
549 x += .5;
550 y += .5;
551 for (i=1; i<npts; i++) {
552 ddx += dddx;
553 ddy += dddy;
554 dx += ddx;
555 dy += ddy;
556 x += dx;
557 y += dy;
558
559 out_x2 = (int)x;
560 out_y2 = (int)y;
561
562 proc(out_x1, out_y1, out_x2, out_y2, data);
563
564 out_x1 = out_x2;
565 out_y1 = out_y2;
566 }
567}
568
569double algo_spline_get_y(double x0, double y0, double x1, double y1,
570 double x2, double y2, double x3, double y3,
571 double in_x)
572{
573 int npts;
574 double out_x, old_x;
575 double out_y, old_y;
576
577 /* Derivatives of x(t) and y(t). */
578 double x, dx, ddx, dddx;
579 double y, dy, ddy, dddy;
580 int i;
581
582 /* Temp variables used in the setup. */
583 double dt, dt2, dt3;
584 double xdt2_term, xdt3_term;
585 double ydt2_term, ydt3_term;
586
587#define MAX_POINTS 64
588#undef DIST
589#define DIST(x, y) (sqrt ((x) * (x) + (y) * (y)))
590 npts = (int) (sqrt (DIST(x1-x0, y1-y0) +
591 DIST(x2-x1, y2-y1) +
592 DIST(x3-x2, y3-y2)) * 1.2);
593 if (npts > MAX_POINTS)
594 npts = MAX_POINTS;
595 else if (npts < 4)
596 npts = 4;
597
598 dt = 1.0 / (npts-1);
599 dt2 = (dt * dt);
600 dt3 = (dt2 * dt);
601
602 xdt2_term = 3 * (x2 - 2*x1 + x0);
603 ydt2_term = 3 * (y2 - 2*y1 + y0);
604 xdt3_term = x3 + 3 * (-x2 + x1) - x0;
605 ydt3_term = y3 + 3 * (-y2 + y1) - y0;
606
607 xdt2_term = dt2 * xdt2_term;
608 ydt2_term = dt2 * ydt2_term;
609 xdt3_term = dt3 * xdt3_term;
610 ydt3_term = dt3 * ydt3_term;
611
612 dddx = 6*xdt3_term;
613 dddy = 6*ydt3_term;
614 ddx = -6*xdt3_term + 2*xdt2_term;
615 ddy = -6*ydt3_term + 2*ydt2_term;
616 dx = xdt3_term - xdt2_term + 3 * dt * (x1 - x0);
617 dy = ydt3_term - ydt2_term + dt * 3 * (y1 - y0);
618 x = x0;
619 y = y0;
620
621 old_x = x0;
622 out_y = old_y = y0;
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 out_y = y;
636 if (out_x > in_x) {
637 out_y = old_y + (out_y-old_y) * (in_x-old_x) / (out_x-old_x);
638 break;
639 }
640 old_x = out_x;
641 old_y = out_y;
642 }
643
644 return out_y;
645}
646
647double algo_spline_get_tan(double x0, double y0, double x1, double y1,
648 double x2, double y2, double x3, double y3,
649 double in_x)
650{
651 double out_x, old_x, old_dx, old_dy;
652 int npts;
653
654 /* Derivatives of x(t) and y(t). */
655 double x, dx, ddx, dddx;
656 double y, dy, ddy, dddy;
657 int i;
658
659 /* Temp variables used in the setup. */
660 double dt, dt2, dt3;
661 double xdt2_term, xdt3_term;
662 double ydt2_term, ydt3_term;
663
664#define MAX_POINTS 64
665#undef DIST
666#define DIST(x, y) (sqrt ((x) * (x) + (y) * (y)))
667 npts = (int) (sqrt (DIST(x1-x0, y1-y0) +
668 DIST(x2-x1, y2-y1) +
669 DIST(x3-x2, y3-y2)) * 1.2);
670 if (npts > MAX_POINTS)
671 npts = MAX_POINTS;
672 else if (npts < 4)
673 npts = 4;
674
675 dt = 1.0 / (npts-1);
676 dt2 = (dt * dt);
677 dt3 = (dt2 * dt);
678
679 xdt2_term = 3 * (x2 - 2*x1 + x0);
680 ydt2_term = 3 * (y2 - 2*y1 + y0);
681 xdt3_term = x3 + 3 * (-x2 + x1) - x0;
682 ydt3_term = y3 + 3 * (-y2 + y1) - y0;
683
684 xdt2_term = dt2 * xdt2_term;
685 ydt2_term = dt2 * ydt2_term;
686 xdt3_term = dt3 * xdt3_term;
687 ydt3_term = dt3 * ydt3_term;
688
689 dddx = 6*xdt3_term;
690 dddy = 6*ydt3_term;
691 ddx = -6*xdt3_term + 2*xdt2_term;
692 ddy = -6*ydt3_term + 2*ydt2_term;
693 dx = xdt3_term - xdt2_term + 3 * dt * (x1 - x0);
694 dy = ydt3_term - ydt2_term + dt * 3 * (y1 - y0);
695 x = x0;
696 y = y0;
697
698 old_x = x0;
699 old_dx = dx;
700 old_dy = dy;
701
702 x += .5;
703 y += .5;
704 for (i=1; i<npts; i++) {
705 ddx += dddx;
706 ddy += dddy;
707 dx += ddx;
708 dy += ddy;
709 x += dx;
710 y += dy;
711
712 out_x = x;
713 if (out_x > in_x) {
714 dx = old_dx + (dx-old_dx) * (in_x-old_x) / (out_x-old_x);
715 dy = old_dy + (dy-old_dy) * (in_x-old_x) / (out_x-old_x);
716 break;
717 }
718 old_x = out_x;
719 old_dx = dx;
720 old_dy = dy;
721 }
722
723 return dy / dx;
724}
725
726} // namespace doc
727