' Display
SCREEN _NEWIMAGE(800, 600, 32)
_SCREENMOVE (_DESKTOPWIDTH \ 2 - _WIDTH \ 2) - 3, (_DESKTOPHEIGHT \ 2 - _HEIGHT \ 2) - 29
_TITLE "If these curves were smoother they'd steal your wife."
_DELAY 1

' Meta
start:
CLEAR
CLS
RANDOMIZE TIMER

' Data structures
TYPE Vector
    x AS DOUBLE
    y AS DOUBLE
END TYPE

' Object type
TYPE Object
    Elements AS INTEGER
    Shade AS _UNSIGNED LONG
END TYPE

' Object storage
DIM SHARED Shape(300) AS Object
DIM SHARED PointChain(300, 500) AS Vector
DIM SHARED TempChain(300, 500) AS Vector
DIM SHARED ShapeCount AS INTEGER
DIM SHARED SelectedShape AS INTEGER

' Initialize
ShapeCount = 0

' Main loop
DO
    IF (UserInput = -1) THEN GOTO start
    CALL Graphics
    _LIMIT 120
LOOP

END

FUNCTION UserInput
    TheReturn = 0
    ' Keyboard input
    kk = _KEYHIT
    SELECT CASE kk
        CASE 32
            DO: LOOP UNTIL _KEYHIT
            WHILE _MOUSEINPUT: WEND
            _KEYCLEAR
            CALL NewMouseShape(7.5, 150, 15)
            CLS
    END SELECT
    IF (kk) THEN
        _KEYCLEAR
    END IF
    UserInput = TheReturn
END FUNCTION

SUB Graphics
    LINE (0, 0)-(_WIDTH, _HEIGHT), _RGBA(0, 0, 0, 200), BF
    CALL cprintstring(16 * 17, "PRESS SPACE and then drag MOUSE 1 to draw a new shape.")
    FOR ShapeIndex = 1 TO ShapeCount
        FOR i = 1 TO Shape(ShapeIndex).Elements - 1
            CALL lineSmooth(PointChain(ShapeIndex, i).x, PointChain(ShapeIndex, i).y, PointChain(ShapeIndex, i + 1).x, PointChain(ShapeIndex, i + 1).y, Shape(ShapeIndex).Shade)
        NEXT
    NEXT
    _DISPLAY
END SUB

SUB NewMouseShape (rawresolution AS DOUBLE, targetpoints AS INTEGER, smoothiterations AS INTEGER)
    ShapeCount = ShapeCount + 1
    numpoints = 0
    xold = 999 ^ 999
    yold = 999 ^ 999
    DO
        DO WHILE _MOUSEINPUT
            x = _MOUSEX
            y = _MOUSEY
            IF (x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT) THEN
                IF _MOUSEBUTTON(1) THEN
                    x = x - (_WIDTH / 2)
                    y = -y + (_HEIGHT / 2)
                    delta = SQR((x - xold) ^ 2 + (y - yold) ^ 2)
                    IF (delta > rawresolution) AND (numpoints < targetpoints - 1) THEN
                        numpoints = numpoints + 1
                        PointChain(ShapeCount, numpoints).x = x
                        PointChain(ShapeCount, numpoints).y = y
                        CALL cpset(x, y, _RGB(0, 255, 255))
                        xold = x
                        yold = y
                    END IF
                END IF
            END IF
        LOOP
        _DISPLAY
    LOOP UNTIL NOT _MOUSEBUTTON(1) AND (numpoints > 1)

    DO WHILE (numpoints < targetpoints)
        rad2max = -1
        kmax = -1
        FOR k = 1 TO numpoints - 1
            xfac = PointChain(ShapeCount, k).x - PointChain(ShapeCount, k + 1).x
            yfac = PointChain(ShapeCount, k).y - PointChain(ShapeCount, k + 1).y
            rad2 = xfac ^ 2 + yfac ^ 2
            IF rad2 > rad2max THEN
                kmax = k
                rad2max = rad2
            END IF
        NEXT
        FOR j = numpoints TO kmax + 1 STEP -1
            PointChain(ShapeCount, j + 1).x = PointChain(ShapeCount, j).x
            PointChain(ShapeCount, j + 1).y = PointChain(ShapeCount, j).y
        NEXT
        PointChain(ShapeCount, kmax + 1).x = (1 / 2) * (PointChain(ShapeCount, kmax).x + PointChain(ShapeCount, kmax + 2).x)
        PointChain(ShapeCount, kmax + 1).y = (1 / 2) * (PointChain(ShapeCount, kmax).y + PointChain(ShapeCount, kmax + 2).y)
        numpoints = numpoints + 1
    LOOP

    FOR j = 1 TO smoothiterations
        FOR k = 2 TO numpoints - 1
            TempChain(ShapeCount, k).x = (1 / 2) * (PointChain(ShapeCount, k - 1).x + PointChain(ShapeCount, k + 1).x)
            TempChain(ShapeCount, k).y = (1 / 2) * (PointChain(ShapeCount, k - 1).y + PointChain(ShapeCount, k + 1).y)
        NEXT
        FOR k = 2 TO numpoints - 1
            PointChain(ShapeCount, k).x = TempChain(ShapeCount, k).x
            PointChain(ShapeCount, k).y = TempChain(ShapeCount, k).y
        NEXT
    NEXT

    Shape(ShapeCount).Elements = numpoints
    Shape(ShapeCount).Shade = _RGB(100 + INT(RND * 155), 100 + INT(RND * 155), 100 + INT(RND * 155))
    SelectedShape = ShapeCount
END SUB

SUB cpset (x1, y1, col AS _UNSIGNED LONG)
    PSET (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), col
END SUB

SUB cprintstring (y, a AS STRING)
    _PRINTSTRING (_WIDTH / 2 - (LEN(a) * 8) / 2, -y + _HEIGHT / 2), a
END SUB

SUB lineSmooth (x0, y0, x1, y1, c AS _UNSIGNED LONG)
    'translated from
    'https://en.wikipedia.org/w/index.php?title=Xiaolin_Wu%27s_line_algorithm&oldid=852445548

    DIM plX AS INTEGER, plY AS INTEGER, plI

    DIM steep AS _BYTE
    steep = ABS(y1 - y0) > ABS(x1 - x0)

    IF steep THEN
        SWAP x0, y0
        SWAP x1, y1
    END IF

    IF x0 > x1 THEN
        SWAP x0, x1
        SWAP y0, y1
    END IF

    DIM dx, dy, gradient
    dx = x1 - x0
    dy = y1 - y0
    gradient = dy / dx

    IF dx = 0 THEN
        gradient = 1
    END IF

    'handle first endpoint
    DIM xend, yend, xgap, xpxl1, ypxl1
    xend = _ROUND(x0)
    yend = y0 + gradient * (xend - x0)
    xgap = (1 - ((x0 + .5) - INT(x0 + .5)))
    xpxl1 = xend 'this will be used in the main loop
    ypxl1 = INT(yend)
    IF steep THEN
        plX = ypxl1
        plY = xpxl1
        plI = (1 - (yend - INT(yend))) * xgap
        GOSUB plot

        plX = ypxl1 + 1
        plY = xpxl1
        plI = (yend - INT(yend)) * xgap
        GOSUB plot
    ELSE
        plX = xpxl1
        plY = ypxl1
        plI = (1 - (yend - INT(yend))) * xgap
        GOSUB plot

        plX = xpxl1
        plY = ypxl1 + 1
        plI = (yend - INT(yend)) * xgap
        GOSUB plot
    END IF

    DIM intery
    intery = yend + gradient 'first y-intersection for the main loop

    'handle second endpoint
    DIM xpxl2, ypxl2
    xend = _ROUND(x1)
    yend = y1 + gradient * (xend - x1)
    xgap = ((x1 + .5) - INT(x1 + .5))
    xpxl2 = xend 'this will be used in the main loop
    ypxl2 = INT(yend)
    IF steep THEN
        plX = ypxl2
        plY = xpxl2
        plI = (1 - (yend - INT(yend))) * xgap
        GOSUB plot

        plX = ypxl2 + 1
        plY = xpxl2
        plI = (yend - INT(yend)) * xgap
        GOSUB plot
    ELSE
        plX = xpxl2
        plY = ypxl2
        plI = (1 - (yend - INT(yend))) * xgap
        GOSUB plot

        plX = xpxl2
        plY = ypxl2 + 1
        plI = (yend - INT(yend)) * xgap
        GOSUB plot
    END IF

    'main loop
    DIM x
    IF steep THEN
        FOR x = xpxl1 + 1 TO xpxl2 - 1
            plX = INT(intery)
            plY = x
            plI = (1 - (intery - INT(intery)))
            GOSUB plot

            plX = INT(intery) + 1
            plY = x
            plI = (intery - INT(intery))
            GOSUB plot

            intery = intery + gradient
        NEXT
    ELSE
        FOR x = xpxl1 + 1 TO xpxl2 - 1
            plX = x
            plY = INT(intery)
            plI = (1 - (intery - INT(intery)))
            GOSUB plot

            plX = x
            plY = INT(intery) + 1
            plI = (intery - INT(intery))
            GOSUB plot

            intery = intery + gradient
        NEXT
    END IF

    EXIT SUB

    plot:
    ' Change to regular PSET for standard coordinate orientation.
    CALL cpset(plX, plY, _RGB32(_RED32(c), _GREEN32(c), _BLUE32(c), plI * 255))
    RETURN
END SUB

