+ smooth_draw(x1, y1, x2, y2, x3, y3);
+}
+
+/*
+ Non-Parametric Smooth Curve Generation. Don Kelly 1984
+
+ P+-----+Q'= virtual
+ / / origin
+ / /
+ Q+-----+R
+
+ Let the starting point be P. the ending point R. and the tangent vertex Q.
+ A general point Z on the curve is then
+ Z = (P + R - Q) + (Q - P) sin t + (Q - R) cos t
+
+ Expanding the Cartesian coordinates around (P + R - Q) gives
+ [x y] = Z - (P + R - Q)
+ [a c] = Q - P
+ [b d] = Q - R
+ x = a*sin(t) + b*cos(t)
+ y = c*sin(t) + d*cos(t)
+
+ from which t can now be eliminated via
+ c*x - a*y = (c*b - a*d)*cos(t)
+ d*x - b*y = (a*d - c*b)*sin(t)
+
+ giving the Cartesian equation for the ellipse as
+ f(x, y) = (c*x - a*y)**2 + (d*x - b*y)**2 - (a*d - c*b)**2 = 0
+
+ or: f(x, y) = A*x**2 - 2*B*x*y + C*y**2 + B**2 - A*C = 0
+ where: A = c**2 + d**2, B = a*c + b*d, C = a**2 + b**2
+
+ The maximum y extent of the ellipse may now be derived as follows:
+ let df/dx = 0, 2*A*x = 2*B*y, x = y*B/A
+ f(x, y) == B**2 * y**2 / A - 2*B**2 * y**2 / A + C*y**2 + B**2 - A*C = 0
+ (A*C - B**2)*y = (A*C - B**2)*A
+ max x = sqrt(C), at y = B/sqrt(C)
+ max y = sqrt(A), at x = B/sqrt(A)
+
+ */
+
+
+/* x1,y1 = P, x2,y2 = Q, x3,y3=R,
+ * draw from P to Q to R if top=0
+ * or from P to (x,ymax) if top>0
+ * or from Q to (x,ymax) if top<0
+ */
+void smooth_line::init0(int x1,int y1, int x2,int y2, int x3,int y3, int top)
+{
+ int x0 = x1+x3-x2, y0 = y1+y3-y2; // Q'
+
+ int a = x2-x1, c = y2-y1;
+ int b = x2-x3, d = y2-y3;
+ A = c*c + d*d; C = a*a + b*b; B = a*c + b*d;
+
+ sx = top >= 0 ? x1 : x3;
+ sy = top >= 0 ? y1 : y3;
+ xs = x2 > sx || (x2==sx && (x1+x3-sx)>=x2) ? 1 : -1;
+ int64_t px = sx-x0, py = sy-y0;
+ dx = A*px - B*py; dy = C*py - B*px;
+ r = 0;
+
+ if( top ) {
+ double ymy = sqrt(A), ymx = B/ymy;
+ ex = x0 + rnd(ymx);
+ ey = y0 + rnd(ymy);
+ }
+ else {
+ ex = x3; ey = y3;
+ }
+
+ ys = a*b > 0 && (!top || top*xs*(b*c - a*d) > 0) ? -1 : 1;
+ if( ys < 0 ) {
+ double xmx = xs*sqrt(C), xmy = B/xmx;
+ xmxx = x0 + rnd(xmx);
+ xmxy = y0 + rnd(xmy);
+ }
+ else {
+ xmxx = ex; xmxy = ey;
+ }
+}
+
+/* x1,y1 = P, x2,y2 = Q, x3,y3=R,
+ * draw from (x,ymax) to P
+ */
+void smooth_line::init1(int x1,int y1, int x2,int y2, int x3,int y3)
+{
+ int x0 = x1+x3-x2, y0 = y1+y3-y2; // Q'
+
+ int a = x2-x1, c = y2-y1;
+ int b = x2-x3, d = y2-y3;
+ A = c*c + d*d; C = a*a + b*b; B = a*c + b*d;
+
+ double ymy = -sqrt(A), ymx = B/ymy;
+ int64_t px = rnd(ymx), py = rnd(ymy);
+ sx = x0 + px; ex = x1;
+ sy = y0 + py; ey = y1;
+ xs = x2 > x1 || (x2==x1 && x3>=x2) ? 1 : -1;
+ dx = A*px - B*py; dy = C*py - B*px;
+ r = 4 * (A*px*px - 2*B*px*py + C*py*py + B*B - A*C);
+
+ ys = a*b > 0 && xs*(b*c - a*d) < 0 ? -1 : 1;
+ if( ys < 0 ) {
+ double xmx = xs*sqrt(C), xmy = B/xmx;
+ xmxx = x0 + rnd(xmx);
+ xmxy = y0 + rnd(xmy);
+ }
+ else {
+ xs = -xs;
+ xmxx = ex; xmxy = ey;
+ }
+ if( xs > 0 )
+ vframe->draw_pixel(sx, sy);
+ while( xs*(sx-xmxx) < 0 && (xs*dx < 0 || rx() < 0) ) {
+ moveX(rx());
+ vframe->draw_pixel(sx, sy);
+ }
+}
+
+
+void VFrame::smooth_draw(int x1, int y1, int x2, int y2, int x3, int y3)
+{
+//printf("p smooth_draw( %d,%d, %d,%d, %d,%d )\n", x1,y1,x2,y2,x3,y3);
+ if( y1 > y3 ) { // make y3 >= y1
+ int xt = x1; x1 = x3; x3 = xt;
+ int yt = y1; y1 = y3; y3 = yt;
+ }
+ if( y1 > y2 && y3 > y2 ) {
+ smooth_line lt(this), rt(this); // Q on bottom
+ lt.init1(x1, y1, x2, y2, x3, y3);
+ rt.init1(x3, y3, x2, y2, x1, y1);
+ while( !lt.done || !rt.done ) {
+ lt.draw();
+ rt.draw();