   ### 15.4.1345  src/lib/x-kit/draw/beta2-spline.pkg

`## beta2-spline.pkg`

`# Compiled by:`
`#     `src/lib/x-kit/draw/xkit-draw.sublib

`###             "When the blind beetle crawls over the surface of a globe,`
`###              he doesn't realize that the track he has covered is curved.`
`###              I was lucky enough to have spotted it."`
`###`
`###                                            -- Albert Einstein`

`stipulate`
`    package g2d =  geometry2d;                  # geometry2d    is from   `src/lib/std/2d/geometry2d.pkg
`herein`

`    package   beta2_spline`
`    : (weak)  Beta2_Spline                      # Beta2_Spline  is from   `src/lib/x-kit/draw/beta2-spline.api
`    {`
`        include package   geometry2d;           # geometry2d    is from   `src/lib/std/2d/geometry2d.pkg

`        fun round x`
`            =`
`            if (x > 0.0)    floor (x+0.5);`
`            else         -1*floor(-x+0.5);`
`            fi;`

`        fun add_seg ([], x0, y0, x1, y1)`
`                => `
`                [ { col => round x0, row => round y0 },`
`                  { col => round x1, row => round y1 }`
`                ];`

`            add_seg (l, x0, y0, _, _)`
`                =>`
`                { col => round x0, row => round y0 }  !  l;`
`        end;`

`        # is_flat:`
`        # Returns TRUE if the polygon determined by the four points`
`        # is flat enough. Flatness is measured by the maximum distance`
`        # of (x1, y1) and (x2, y2) from the line determined by (x0, y0)`
`        # and (x3, y3). In addition, check that p1, p2 are close to the `
`        # line segment. To do this, make sure they are roughly within `
`        # the circle with center (p0+p3)/2 and radius = ``|p3-p0|/2+flatness`
`        #`
`        fun is_flat { x0, y0, x1, y1, x2, y2, x3, y3 }`
`            =`
`            {   fun sqr x =   x * x;`
`                #`
`                dx =   x3 - x0;`
`                dy =   y3 - y0;`

`                midx  =   0.5*dx;`
`                midy  =   0.5*dy;`

`                dist2 =   sqr dy + sqr dx;`

`                flatness2 =   sqr 1.0 * dist2;`
`                halfd2 = 0.25*dist2;`

`                fun in_flat_range (x, y)`
`                    =`
`                    sqr (dy * x - dx * y) <= flatness2`
`                    and`
`                    {   d =   sqr (x-midx) + sqr (y-midy);`

`                        d <= halfd2`
`                        or `
`                        sqr (d-halfd2) <= flatness2;`
`                    };`

`                in_flat_range (x1-x0, y1-y0)`
`                and`
`                in_flat_range (x2-x0, y2-y0);`
`            };`

`        # bezier:`
`        # Recursively compute a Bezier cubic section. If the points`
`        # determine a polygon flat enough to be represented as a line`
`        # segment, the segment is added to the list. Otherwise, the`
`        # the curve is bisected and each part is recursively computed,`
`        # with the lists concatenated.`
`        #`
`        # From "The Beta2-split: A special case of the Beta-spline Curve and`
`        # Surface Representation." B. A. Barsky and A. D. DeRose. IEEE, 1985,`
`        # http://www.eecs.berkeley.edu/Pubs/TechRpts/1983/CSD-83-152.pdf`
`        # as adapted by Crispin Goswell for xps.`
`        #`
`        fun bezier (arg as { x0, y0, x1, y1, x2, y2, x3, y3 }, l)`
`            =`
`            if (is_flat arg)`
`                #`
`                add_seg (l, x0, y0, x3, y3);`
`            else`
`                mid_x =   (x0 + x3) / 8.0 + 3.0 * (x1 + x2) / 8.0;`
`                mid_y =   (y0 + y3) / 8.0 + 3.0 * (y1 + y2) / 8.0;`

`                l' =   bezier (`
`                           {   x0 => mid_x,                     y0 => mid_y,`
`                               x1 => (x1+x3) / 4.0 + x2 / 2.0,  y1 => (y1+y3) / 4.0 + y2 / 2.0,`
`                               x2 => (x2+x3) / 2.0,             y2 => (y2+y3) / 2.0,`
`                               x3,                        y3`
`                           },`
`                           l`
`                       );`

`                bezier (`
`                    {   x0,                        y0,`
`                        x1 => (x0+x1) / 2.0,             y1 => (y0+y1) / 2.0,`
`                        x2 => (x0+x2) / 4.0 + x1 / 2.0,  y2 => (y0+y2) / 4.0 + y1 / 2.0,`
`                        x3 => mid_x,                     y3 => mid_y`
`                    },`
`                    l'`
`                );`

`            fi;`

`        # curve:`
`        # Given four points [p0, p1, p2, p3], return a list of points corresponding to `
`        # to a Bezier cubic section, starting at p0, ending at p3, with p1, p2 as`
`        # control points.`
`        #`
`        fun curve`
`            ( { col=>x0, row=>y0 }, { col=>x1, row=>y1 },`
`              { col=>x2, row=>y2 }, { col=>x3, row=>y3 }`
`            )`
`            =`
`            bezier (`
`                {   x0 => float(x0), y0 => float(y0), `
`                    x1 => float(x1), y1 => float(y1), `
`                    x2 => float(x2), y2 => float(y2), `
`                    x3 => float(x3), y3 => float(y3)`
`                },`
`                []`
`            );`

`        # doSpline:`
`        # Given four points (p0, p1, p2, p3), return a list of points corresponding to `
`        # to the B-spline curve section, accumulating the results on the argument list.`
`        # We compute the curve by determining the corresponding Bezier control points,`
`        # and then use the Bezier routines above.`
`        #`
`        fun do_spline`
`            ( { col=>x0, row=>y0 },`
`              { col=>x1, row=>y1 },`
`              { col=>x2, row=>y2 },`
`              { col=>x3, row=>y3 },`
`              l`
`            )`
`            =`
`            {   x0 =   float(x0);   y0 = float(y0);`
`                x1 =   float(x1);   y1 = float(y1); `
`                x2 =   float(x2);   y2 = float(y2);`
`                x3 =   float(x3);   y3 = float(y3);`

`                bezier (`
`                    {   x0 => (x0 + 4.0*x1 + x2)/6.0,   y0 => (y0 + 4.0*y1 + y2)/6.0, `
`                        x1 => (2.0*x1 + x2)/3.0,        y1 => (2.0*y1 + y2)/3.0, `
`                        x2 => (x1 + 2.0*x2)/3.0,        y2 => (y1 + 2.0*y2)/3.0, `
`                        x3 => (x1 + 4.0*x2 + x3)/6.0,   y3 => (y1 + 4.0*y2 + y3)/6.0`
`                    },`
`                    l`
`                );`
`            };`

`        # loopSpline:`
`        # Given a list of spline control points, generate the corresponding`
`        # spline. Since we accumulate on the head of the list, we assume`
`        # the calling function has reversed the list of control points.`
`        # The loop continues as long as there are 4 control points left.`
`        #`
`        fun loop_spline (p3, p2, p1, p0 ! tl, l)`
`                => `
`                loop_spline (p2, p1, p0, tl, do_spline (p0, p1, p2, p3, l));`

`           loop_spline (_, _, _, [], l)`
`                =>`
`                l;`
`        end; `

`        # simpleBSpline:`
`        # Compute a simple B-spline with the given control points.`
`        #`
`        fun simple_bspline arg`
`            =`
`            case (reverse arg)`
`                #         `
`                (p3 ! p2 ! p1 ! tl) =>  loop_spline (p3, p2, p1, tl,[]);`
`                _                   =>  arg;`
`            esac;`

`        # bSpline:`
`        # Compute a B-spline using the given control points.`
`        # In addition, we constrain the resultant spline to connect the`
`        # first and last points by adding copies of these points.`
`        #`
`        fun b_spline (arg as (p0 ! _ ! _ ! _))`
`                =>`
`                loop_spline (pn, pn, pn, tl, [])`
`                where`
`                    my (pn, tl)`
`                        =`
`                        case (reverse (p0 ! p0 ! arg))   (pn ! tl) =>  (pn, tl);`
`                             /* */                           _         =>  raise exception DIE "Bug: Unsupported case in b_spline";     `
`                        esac;`
`                end;`

`            b_spline l`
`                =>`
`                l;`
`        end;`

`        # closed_bSpline:`
`        # Compute a closed B-spline. This is done by repeating the first`
`        # three points at the end of the list.`
`        # Note that the first and last points of the result are the same.`
`        #`
`        fun closed_bspline (arg as (p0 ! p1 ! p2 ! _))`
`                =>`
`                loop_spline (p2, p1, p0, reverse arg,[]);`

`            closed_bspline l`
`                =>`
`                l;`
`        end;`
`    };                                  # package spline `
`end;`   