Task

Implement the   Xiaolin Wu's line algorithm   described in Wikipedia.



This algorithm draws anti-aliased lines.



Related task

   See   Bresenham's line algorithm   for aliased lines.


--
-- demo\rosetta\XiaolinWuLine.exw
-- ==============================
--
--  Resize the window to show lines at any angle
--
--  For education/comparision purposes only: see demo\pGUI\aaline.exw
--  for a much shorter version, but "wrong algorithm" for the RC task.
--  Also note this blends with BACK rather than the actual pixel,
--  whereas aaline.exw does it properly.
--
with javascript_semantics   -- not really fair: pwa/p2js uses OpenGL
-- and does not draw bresenham lines anyway/ever, plus the next line 
-- makes no difference whatsoever when running this in a browser.
constant USE_OPENGL = 0

constant TITLE = "Xiaolin Wu's line algorithm"

include pGUI.e

Ihandle dlg, canvas
cdCanvas cddbuffer, cdcanvas

bool wuline = true  -- space toggles, for comparison
                    -- when false, and with USE_OPENGL, lines are still smooth, 
                    --  but a bit thicker - and therefore less "ropey".
                    -- when false, but without USE_OPENGL, it draws bresenham 
                    --  lines (ie jagged, without anti-aliasing [desktop only]).

integer BACK = CD_PARCHMENT,
        LINE = CD_BLUE,
        {rB, gB, bB} = to_rgb(BACK),
        {rL, gL, bL} = to_rgb(LINE)

procedure plot(atom x, y, c, bool steep=false)
--  plot the pixel at (x, y) with brightness c (where 0 <= c <= 1)
    if steep then {x,y} = {y,x} end if
    atom C = 1-c
    c = rgb(rL*c+rB*C,gL*c+gB*C,bL*c+bB*C)
    cdCanvasPixel(cddbuffer, x, y, c)
end procedure

procedure plot2(atom x, y, f, xgap, bool steep)
    plot(x,y,(1-f)*xgap,steep)
    plot(x,y+1,(f)*xgap,steep)
end procedure

function fpart(atom x)
    return x-floor(x)   -- fractional part of x
end function

procedure wu_line(atom x0,y0,x1,y1)
    bool steep := abs(y1 - y0) > abs(x1 - x0)
    if steep then
        {x0, y0, x1, y1} = {y0, x0, y1, x1}
    end if
    if x0>x1 then
        {x0, x1, y0, y1} = {x1, x0, y1, y0}
    end if
    
    atom dx := x1 - x0,
         dy := y1 - y0,
         gradient := iff(dx=0? 1 : dy / dx)

    -- handle first endpoint
    atom xend := round(x0),
         yend := y0 + gradient * (xend - x0),
         xgap := 1 - fpart(x0 + 0.5),
         xpxl1 := xend, -- this will be used in the main loop
         ypxl1 := floor(yend)
    plot2(xpxl1, ypxl1, fpart(yend), xgap, steep)
    atom intery := yend + gradient -- first y-intersection for the main loop
    
    -- handle second endpoint
    xend := round(x1)
    yend := y1 + gradient * (xend - x1)
    xgap := fpart(x1 + 0.5)
    atom xpxl2 := xend, -- this will be used in the main loop
         ypxl2 := floor(yend)
    plot2(xpxl2, ypxl2, fpart(yend), xgap, steep)
    
    -- main loop
    for x = xpxl1+1 to xpxl2-1 do
        plot2(x, floor(intery), fpart(intery), 1, steep)
        intery += gradient
    end for
end procedure

procedure plot_4_points(integer cx, cy, x, y, atom f, angle1=0, angle2=360, angle=0)
integer x1 = cx+x, x2 = cx-x,
        y1 = cy+y, y2 = cy-y
    if angle<0 or angle>90.01 then ?9/0 end if
    if      angle >=angle1 and      angle <=angle2 then plot(x1, y1, f) end if  -- top right
    if (180-angle)>=angle1 and (180-angle)<=angle2 then plot(x2, y1, f) end if  -- top left
    if (180+angle)>=angle1 and (180+angle)<=angle2 then plot(x2, y2, f) end if  -- btm left
    if (360-angle)>=angle1 and (360-angle)<=angle2 then plot(x1, y2, f) end if  -- btm right
end procedure

procedure wu_ellipse(atom cx, cy, w, h, angle1=0, angle2=360)
--
-- (draws a circle when w=h) credit:
-- https://yellowsplash.wordpress.com/2009/10/23/fast-antialiased-circles-and-ellipses-from-xiaolin-wus-concepts/
--
    if w<=0 or h<=0 then return end if
    cx = round(cx)
    cy = round(cy)
    w = round(w)
    h = round(h)
    angle1 = mod(angle1,360)
    angle2 = mod(angle2,360)

    -- Match cdCanvasArc/Sector angles:
    angle1 = atan2((h/2)*sin(angle1*CD_DEG2RAD), (w/2)*cos(angle1*CD_DEG2RAD))*CD_RAD2DEG
    angle2 = atan2((h/2)*sin(angle2*CD_DEG2RAD), (w/2)*cos(angle2*CD_DEG2RAD))*CD_RAD2DEG

    if angle2<=angle1 then angle2 += 360 end if

    atom a := w/2, asq := a*a, 
         b := h/2, bsq := b*b,
         sqab = sqrt(asq+bsq),
         ffd = round(asq/sqab), -- forty-five-degree coord
         xj, yj, frc, flr, angle

    -- draw top right, and the 3 mirrors of it in horizontal fashion
    --  (ie 90 to 45 degrees for a circle)

    for xi=0 to ffd do
        yj := b*sqrt(1-xi*xi/asq)   -- real y value
        frc := fpart(yj)
        flr := floor(yj)
        angle := iff(xi=0?90:arctan(yj/xi)*CD_RAD2DEG)
        plot_4_points(cx, cy, xi, flr, 1-frc, angle1, angle2, angle)
        plot_4_points(cx, cy, xi, flr+1, frc, angle1, angle2, angle)
    end for

    -- switch from horizontal to vertial mode for the rest, ditto 3
    -- (ie 45..0 degrees for a circle)

    ffd = round(bsq/sqab)
    for yi=0 to ffd do
        xj := a*sqrt(1-yi*yi/bsq)   -- real x value
        frc := fpart(xj)
        flr := floor(xj)
        angle = iff(xj=0?0:arctan(yi/xj)*CD_RAD2DEG)
        plot_4_points (cx, cy, flr, yi, 1-frc, angle1, angle2, angle)
        plot_4_points (cx, cy, flr+1, yi, frc, angle1, angle2, angle)
    end for

end procedure

function redraw_cb(Ihandle /*ih*/)
    integer {w, h} = sq_sub(IupGetIntInt(canvas, "DRAWSIZE"),10)
    cdCanvasActivate(cddbuffer)
    cdCanvasClear(cddbuffer)
    if platform()=JS then
        cdCanvasSetLineWidth(cddbuffer,iff(wuline?1:4))
    end if
if wuline then
    wu_line(0,0,200,200)
    wu_line(w,0,200,200)
    wu_line(0,h,200,200)
    wu_line(w,h,200,200)
else
    cdCanvasLine(cddbuffer,0,0,200,200)
    cdCanvasLine(cddbuffer,w,0,200,200)
    cdCanvasLine(cddbuffer,0,h,200,200)
    cdCanvasLine(cddbuffer,w,h,200,200)
end if

if wuline then
    wu_ellipse(200,200,200,200)
--  cdCanvasSector(cddbuffer, 200, 200, 200, 200, 0, 360) 
    wu_ellipse(200,200,300,100)
--  wu_ellipse(200,200,300,100,15,85)
--  cdCanvasArc(cddbuffer, 205, 205, 300, 100, 15, 85) 
else
    cdCanvasArc(cddbuffer, 200, 200, 200, 200, 0, 360) 
--  cdCanvasSector(cddbuffer, 200, 200, 200, 200, 0, 360) 
    cdCanvasArc(cddbuffer, 200, 200, 300, 100, 0, 360) 
end if
--test - it works (much better) if you draw the polygon /after/ the lines!!
--  cdCanvasBegin(cddbuffer,CD_FILL)
--  cdCanvasVertex(cddbuffer,w,h)
--  cdCanvasVertex(cddbuffer,0,h)
--  cdCanvasVertex(cddbuffer,200,200)
--  cdCanvasEnd(cddbuffer)
--/test
    cdCanvasFlush(cddbuffer)
    if USE_OPENGL then
        if platform()!=JS then
            IupGLSwapBuffers(canvas)
        end if
    end if
    return IUP_DEFAULT
end function

function map_cb(Ihandle ih)
    if USE_OPENGL then
        IupGLMakeCurrent(canvas)
        if platform()=JS then
            cdcanvas = cdCreateCanvas(CD_IUP, canvas)
        else
            atom res = IupGetDouble(NULL, "SCREENDPI")/25.4
            cdcanvas = cdCreateCanvas(CD_GL, "10x10 %g", {res})
        end if
        cddbuffer = cdcanvas
    else
        cdcanvas = cdCreateCanvas(CD_IUP, ih)
        cddbuffer = cdCreateCanvas(CD_DBUFFER, cdcanvas)
    end if
    cdCanvasSetBackground(cddbuffer, BACK)
    cdCanvasSetForeground(cddbuffer, LINE)
    return IUP_DEFAULT
end function

function canvas_resize_cb(Ihandle /*canvas*/)
    if USE_OPENGL then
        integer {canvas_width, canvas_height} = IupGetIntInt(canvas, "DRAWSIZE")
        atom res = IupGetDouble(NULL, "SCREENDPI")/25.4
        cdCanvasSetAttribute(cdcanvas, "SIZE", "%dx%d %g", {canvas_width, canvas_height, res})
    end if
    return IUP_DEFAULT
end function

procedure set_title()
    string title = TITLE
    if USE_OPENGL then
        title &= iff(wuline?" (wu_line)":" (opengl)")
    else
        title &= iff(wuline?" (anti-aliased)":" (bresenham)")
    end if
    IupSetStrAttribute(dlg,"TITLE",title)
end procedure

function key_cb(Ihandle /*ih*/, atom c)
    if c=K_ESC then return IUP_CLOSE end if
    if c=' ' then
        wuline = not wuline
        set_title()
        IupRedraw(canvas)
    end if
    return IUP_CONTINUE
end function

procedure main()
    IupOpen()
    if USE_OPENGL then
        canvas = IupGLCanvas()
        IupSetAttribute(canvas, "BUFFER", "DOUBLE")
    else
        canvas = IupCanvas()
    end if

    IupSetAttribute(canvas, "RASTERSIZE", "640x480")
    IupSetCallbacks(canvas, {"MAP_CB", Icallback("map_cb"),
                             "ACTION", Icallback("redraw_cb"),
                             "RESIZE_CB", Icallback("canvas_resize_cb")})
    dlg = IupDialog(canvas)
    set_title()
    IupSetCallback(dlg, "KEY_CB", Icallback("key_cb"))
    IupShow(dlg)
    IupSetAttribute(canvas, "RASTERSIZE", NULL)
    if platform()!=JS then
        IupMainLoop()
        IupClose()
    end if
end procedure

main()
