Jump to content

Xiaolin Wu's line algorithm

From Rosetta Code
Task
Xiaolin Wu's line algorithm
You are encouraged to solve this task according to the task description, using any language you may know.
Task

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


This algorithm draws anti-aliased lines.


Related task



Action!

INCLUDE "H6:REALMATH.ACT"

REAL one

PROC PutBigPixel(INT x,y,col)
  IF x>=0 AND x<=79 AND y>=0 AND y<=47 THEN
    y==LSH 2
    IF col<0 THEN col=0
    ELSEIF col>15 THEN col=15 FI
    Color=col
    Plot(x,y)
    DrawTo(x,y+3)
  FI
RETURN

INT FUNC Abs(INT x)
  IF x<0 THEN RETURN (-x) FI
RETURN (x)

PROC Swap(INT POINTER a,b)
  INT tmp

  tmp=a^ a^=b^ b^=tmp
RETURN

PROC Line(INT x1,y1,x2,y2,col)
  INT x,y,dx,dy,c
  INT POINTER px,py
  REAL rx,ry,grad,rcol,tmp1,tmp2

  dx=Abs(x2-x1)
  dy=Abs(y2-y1)
  IF dy>dx THEN
    Swap(@x1,@y1)
    Swap(@x2,@y2)
    px=@y py=@x
  ELSE
    px=@x py=@y
  FI

  IF x1>x2 THEN
    Swap(@x1,@x2)
    Swap(@y1,@y2)
  FI

  x=x2-x1
  IF x=0 THEN x=1 FI
  IntToRealForNeg(x,rx)
  IntToRealForNeg(y2-y1,ry)
  RealDiv(ry,rx,grad)
  
  IntToRealForNeg(y1,ry)
  IntToReal(col,rcol)
  FOR x=x1 TO x2
  DO
    Frac(ry,tmp1)
    IF IsNegative(tmp1) THEN
      RealAdd(one,tmp1,tmp2)
      RealAssign(tmp2,tmp1)
    FI
    RealMult(rcol,tmp1,tmp2)
    c=Round(tmp2)
    y=Floor(ry)
    PutBigPixel(px^,py^,col-c)
    y==+1
    PutBigPixel(px^,py^,c)
    RealAdd(ry,grad,tmp1)
    RealAssign(tmp1,ry)
  OD
RETURN

PROC Main()
  BYTE CH=$02FC ;Internal hardware value for last key pressed
  REAL tmp,c
  BYTE i,x1,y1,x2,y2

  MathInit()
  IntToReal(1,one)

  Graphics(9)
  FOR i=0 TO 11
  DO
    Line(0,i*4,38,1+i*5,15)
    Line(40+i*4,0,41+i*6,47,15)
  OD

  DO UNTIL CH#$FF OD
  CH=$FF
RETURN
Output:

Screenshot from Atari 8-bit computer

ARM Assembly

Works with: as version Raspberry Pi
/* ARM assembly Raspberry PI  */
/*  program xiaolin1.s   */

/* REMARK 1 : this program use routines in a include file 
   see task Include a file language arm assembly 
   for the routine affichageMess displayerror 
   see at end oh this program the instruction include */

/* REMARK 2 : display use a FrameBuffer device : see raspberry pi FrameBuffer documentation
              this solution write directly on the screen of raspberry pi
              other solution is to use X11 windows but X11 has a function drawline !! */

/* REMARK 3 : this program do not respect the convention for use, save and restau registers
              in rhe routine call !!!!   */

/*******************************************/
/* Constantes    */
/*******************************************/
.equ STDOUT,              1     @ Linux output console
.equ EXIT,                1     @ Linux syscall
.equ WRITE,               4     @ Linux syscall
.equ OPEN,                5
.equ CLOSE,               6
.equ IOCTL,               0x36
.equ MMAP,                0xC0
.equ UNMAP,               0x5B
.equ O_RDWR,              0x0002    @ open for reading and writing
.equ MAP_SHARED,          0x01      @ Share changes.
.equ PROT_READ,           0x1       @ Page can be read.
.equ PROT_WRITE,          0x2       @ Page can be written.

/*******************************************/
/* Initialized data                        */
/*******************************************/
.data
szMessErreur:   .asciz "File open error.\n"
szMessErreur1:  .asciz "File close error.\n"
szMessErreur2:  .asciz "File mapping error.\n"
szMessDebutPgm: .asciz "Program start. \n"
szMessFinOK:    .asciz "Normal end program. \n"
szMessErrFix:   .asciz  "Read error info fix framebuffer  \n"
szMessErrVar:   .asciz  "Read error info var framebuffer  \n"
szRetourligne:  .asciz  "\n"
szParamNom:     .asciz "/dev/fb0"         @ FrameBuffer device name
szLigneVar:     .ascii "Variables info : "
sWidth:        .fill 11, 1, ' ' 
                .ascii " * "
sHeight:        .fill 11, 1, ' ' 
                .ascii " Bits par pixel : "
sBits:           .fill 11, 1, ' '
                .asciz  "\n"
/*************************************************/
szMessErr: .ascii	"Error code hexa : "
sHexa: .space 9,' '
         .ascii "  decimal :  "
sDeci: .space 15,' '
         .asciz "\n"
.align 4
/* codes fonction pour la récupération des données fixes et variables */
FBIOGET_FSCREENINFO: .int 0x4602  @ function code for read infos fixes Framebuffer
FBIOGET_VSCREENINFO: .int 0x4600  @ function code for read infos variables Framebuffer

/*******************************************/
/* UnInitialized data */
/*******************************************/ 
.bss
.align 4
fix_info: .skip FBFIXSCinfo_fin                @ memory reserve for structure FSCREENINFO
.align 4
var_info: .skip FBVARSCinfo_fin                @ memory reserve for structure VSCREENINFO
/**********************************************/
/* -- Code section                            */
/**********************************************/
.text
.global main

main:
    ldr r0,iAdrszMessDebutPgm
    bl affichageMess                  @ display message
    ldr r0,iAdrszParamNom             @ frameBuffer device name
    mov r1,#O_RDWR                    @ flags read/write
    mov r2,#0                         @ mode 
    mov r7,#OPEN                      @ open device FrameBuffer 
    svc 0 
    cmp r0,#0                         @ error ?
    ble erreur
    mov r10,r0                        @ save FD du device FrameBuffer in r10
    
    ldr r1,iAdrFBIOGET_VSCREENINFO    @ read variables datas of FrameBuffer
    ldr r1,[r1]                       @ load code function
    ldr r2,iAdrvar_info               @ structure memory address
    mov r7, #IOCTL                    @ call system
    swi 0 
    cmp r0,#0
    blt erreurVar
    ldr r2,iAdrvar_info
    ldr r0,[r2,#FBVARSCinfo_xres]     @ load screen width
    ldr r1,iAdrsWidth                 @ and convert in string for display
    bl conversion10S
    ldr r0,[r2,#FBVARSCinfo_yres]     @ load screen height 
    ldr r1,iAdrsHeight                @ and convert in string for display
    bl conversion10S
    ldr r0,[r2,#FBVARSCinfo_bits_per_pixel]  @ load bits by pixel  
    ldr r1,iAdrsBits                  @ and convert in string for display
    bl conversion10S
    ldr r0,iAdrszLigneVar             @ display result 
    bl affichageMess

    mov r0,r10                        @ FD du FB
    ldr r1,iAdrFBIOGET_FSCREENINFO    @ read fixes datas of FrameBuffe
    ldr r1,[r1]                       @ load code function
    ldr r2,iAdrfix_info               @ structure memory address
    mov r7, #IOCTL                    @ call system
    svc 0 
    cmp r0,#0                         @ error ?
    blt erreurFix
    ldr r0,iAdrfix_info

    ldr r1,iAdrfix_info               @ read size memory for datas
    ldr r1,[r1,#FBFIXSCinfo_smem_len] @ in octets
                                      @ datas mapping
    mov r0,#0
    ldr r2,iFlagsMmap
    mov r3,#MAP_SHARED
    mov r4,r10
    mov r5,#0
    mov r7, #MMAP                     @ 192 call system for mapping
    swi #0 
    cmp r0,#0                         @ error ?
    beq erreur2    
    mov r9,r0                         @ save mapping address in r9
    /*************************************/
    /* display draw                      */
    bl dessin
    /************************************/
    mov r0,r9                         @ mapping close
    ldr r1,iAdrfix_info
    ldr r1,[r1,#FBFIXSCinfo_smem_len] @ mapping memory size
    mov r7,#UNMAP                     @call system 91 for unmapping
    svc #0                            @ error ?
    cmp r0,#0
    blt erreur1    
                                      @ close device FrameBuffer
    mov r0,r10                        @ load FB du device
    mov r7, #CLOSE                    @ call system
    swi 0 
    ldr r0,iAdrszMessFinOK            @ display end message
    bl affichageMess
    mov r0,#0                         @ return code = OK
    b 100f
erreurFix:                            @ display read error datas fix
    ldr r1,iAdrszMessErrFix           @ message address
    bl   displayError                 @ call display
    mov r0,#1                         @ return code = error
    b 100f
erreurVar:                            @ display read error datas var
    ldr r1,iAdrszMessErrVar
    bl   displayError
    mov r0,#1
    b 100f
erreur:                               @ display open error 
    ldr r1,iAdrszMessErreur
    bl   displayError
    mov r0,#1
    b 100f
erreur1:                              @ display unmapped error
    ldr r1,iAdrszMessErreur1
    bl   displayError
    mov r0,#1
    b 100f
erreur2:                              @ display mapped error
    ldr r1,iAdrszMessErreur2
    bl   displayError
    mov r0,#1
    b 100f
100:                                  @ end program
    mov r7, #EXIT
    svc 0 
/************************************/
iAdrszParamNom:           .int szParamNom
iFlagsMmap:               .int PROT_READ|PROT_WRITE
iAdrszMessErreur:         .int szMessErreur
iAdrszMessErreur1:        .int szMessErreur1
iAdrszMessErreur2:        .int szMessErreur2
iAdrszMessDebutPgm:       .int szMessDebutPgm
iAdrszMessFinOK:          .int szMessFinOK
iAdrszMessErrFix:         .int szMessErrFix
iAdrszMessErrVar:         .int szMessErrVar
iAdrszLigneVar:           .int szLigneVar
iAdrvar_info:             .int var_info
iAdrfix_info:             .int fix_info
iAdrFBIOGET_FSCREENINFO:  .int FBIOGET_FSCREENINFO
iAdrFBIOGET_VSCREENINFO:  .int FBIOGET_VSCREENINFO
iAdrsWidth:               .int sWidth
iAdrsHeight:              .int sHeight
iAdrsBits:                .int sBits
/***************************************************/
/*   dessin                  */
/***************************************************/
/* r9 framebuffer memory address   */
dessin:
    push {r1-r12,lr}                  @ save registers
    mov r0,#255                       @ red
    mov r1,#255                       @ green
    mov r2,#255                       @ blue    3 bytes 255 = white
    bl codeRGB                        @ code color RGB  32 bits
    mov r1,r0                         @ background color
    ldr r0,iAdrfix_info               @ load memory mmap size 
    ldr r0,[r0,#FBFIXSCinfo_smem_len]    
    bl coloriageFond                  @
    /* draw line 1  */
    mov r0,#200                       @ X start line
    mov r1,#200                       @ Y start line
    mov r2,#200                       @ X end line
    mov r3,#100                       @ Y end line
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]     @ load screen width
    bl drawLine
    /* draw line 2  */
    mov r0,#200
    mov r1,#200
    mov r2,#200
    mov r3,#300
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine
    /* draw line 3  */
    mov r0,#200
    mov r1,#200
    mov r2,#100
    mov r3,#200
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine
    /* draw line 4  */
    mov r0,#200
    mov r1,#200
    mov r2,#300
    mov r3,#200
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine
    /* draw line 5  */
    mov r0,#200
    mov r1,#200
    mov r2,#100
    mov r3,#100
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine
    /* draw line 6  */
    mov r0,#200
    mov r1,#200
    mov r2,#100         
    mov r3,#300
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine
    /* draw line 7  */
    mov r0,#200
    mov r1,#200
    mov r2,#300
    mov r3,#300
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine
    /* draw line 8  */
    mov r0,#200
    mov r1,#200
    mov r2,#300
    mov r3,#100
    ldr r4,iAdrvar_info
    ldr r4,[r4,#FBVARSCinfo_xres]
    bl drawLine

100:
    pop {r1-r12,lr}                 @ restaur registers
    bx lr                           @ end function

/********************************************************/
/*   set background color                               */
/********************************************************/
/* r0 contains size screen memory  */
/* r1 contains rgb code color      */
/* r9 contains screen memory address */
coloriageFond:
    push {r2,lr}
    mov r2,#0                     @ counter 
1:                                @ begin loop
    str r1,[r9,r2]
    add r2,#4
    cmp r2,r0
    blt 1b
    pop {r2,lr}
    bx lr
/********************************************************/
/*   Xiaolin Wu  line algorithm                        */
/*  no floating point compute,  multiply value for 128  */
/*  for integer compute                                 */
/********************************************************/
/* r0  x1 start line */
/* r1  y1 start line */
/* r2  x2 end line */
/* r3  y2 end line */
/* r4  screen width */
drawLine:
    push {fp,lr}      @ save registers ( no other registers save )
    mov r5,r0         @ save x1
    mov r6,r1         @ save y1
    cmp r2,r5         @ compar x2,x1
    subgt r1,r2,r5
    suble r1,r5,r2    @ compute dx=abs val de x1-x2
    cmp r3,r6         @ compar y2,y1
    subgt r0,r3,r6
    suble r0,r6,r3    @ compute dy = abs val de y1-y2
    cmp r1,r0         @ compare dx , dy
    blt 5f            @ dx < dy
                      @ dx > dy
    cmp r2,r5         @ compare x2,x1
    movlt r8,r5       @ x2 < x1 
    movlt r5,r2       @ swap x2,x1
    movlt r2,r8
    movlt r8,r6       @ swap y2,y1
    movlt r6,r3
    movlt r3,r8
    lsl r0,#7         @ * by 128
    mov r7,r2         @ save x2
    mov r8,r3         @ save y2
    cmp r1,#0         @ divisor = 0 ?
    moveq r10,#128
    beq 1f
    bl division       @ gradient compute (* 128)
    mov r10,r2        @ r10 contient le gradient
1:
    @ display start points
    mov r0,#64
    bl colorPixel
    mov r3,r0              @ RGB color 
    mov r0,r5              @ x1
    mov r1,r6              @ y1
    mov r2,r4              @ screen witdh 
    bl aff_pixel_codeRGB32 @ display pixel
    add r1,#1              @ increment y1
    bl aff_pixel_codeRGB32
    @ display end points
    mov r0,r7              @ x2
    mov r1,r8              @ y2
    bl aff_pixel_codeRGB32
    add r1,#1              @ increment y2 
    bl aff_pixel_codeRGB32
    cmp r8,r6              @ compar y2,y1
    blt 3f                 @ y2 < y1
    mov r4,r5              @ x =  x1 
    lsl r5,r6,#7           @ compute y1 * 128
    add r5,r10             @ compute intery = (y1 * 128 + gradient * 128)
2:                         @ start loop draw line pixels
    lsr r1,r5,#7           @ intery / 128  = y
    lsl r8,r1,#7
    sub r6,r5,r8           @ reminder of intery /128 = brightness
    mov r0,r6
    bl colorPixel          @ compute rgb color brightness
    mov r3,r0              @ rgb color
    mov r0,r4              @ x 
    bl aff_pixel_codeRGB32 @ display pixel 1
    add r1,#1              @ increment y
    rsb r0,r6,#128         @ compute 128 - brightness
    bl colorPixel          @ compute new rgb color
    mov r3,r0
    mov r0,r4
    bl aff_pixel_codeRGB32 @ display pixel 2
    add r5,r10             @ add gradient to intery
    add r4,#1              @ increment x
    cmp r4,r7              @ x < x2
    ble 2b                 @ yes -> loop
    b 100f                 @ else end
3:                         @ y2 < y1  
    mov r4,r7              @ x = x2 
    mov r7,r5              @ save x1
    lsl r5,r8,#7           @ y = y1 * 128 
    add r5,r10             @ compute intery = (y1 * 128 + gradient * 128)
4:
    lsr r1,r5,#7           @ y = ent(intery / 128)
    lsl r8,r1,#7
    sub r8,r5,r8           @ brightness = remainder
    mov r0,r8
    bl colorPixel
    mov r3,r0
    mov r0,r4
    bl aff_pixel_codeRGB32
    add r1,#1
    rsb r0,r8,#128
    bl colorPixel
    mov r3,r0
    mov r0,r4
    bl aff_pixel_codeRGB32
    add r5,r10
    sub r4,#1             @ decrement x
    cmp r4,r7             @ x > x1
    bgt 4b                @ yes -> loop
    b 100f
5:                        @ dx < dy
    cmp r3,r6             @ compare y2,y1
    movlt r8,r5           @ y2 < y1 
    movlt r5,r2           @ swap x1,x2
    movlt r2,r8
    movlt r8,r6           @ swap y1,y2
    movlt r6,r3
    movlt r3,r8
    mov r8,r1             @ swap r0,r1 for routine division
    mov r1,r0
    lsl r0,r8,#7          @ dx * by 128
    mov r7,r2             @ save x2
    mov r8,r3             @ save y2
    cmp r1,#0             @ dy = zero ?
    moveq r10,#128
    beq 6f
    bl division           @  compute gradient * 128
    mov r10,r2            @  gradient -> r10
6:
    @ display start points
    mov r0,#64
    bl colorPixel
    mov r3,r0             @ color pixel
    mov r0,r5             @ x1
    mov r1,r6             @ y1
    mov r2,r4             @ screen width
    bl aff_pixel_codeRGB32
    add r1,#1
    bl aff_pixel_codeRGB32
    @ display end points
    mov r0,r7
    mov r1,r8
    bl aff_pixel_codeRGB32
    add r1,#1
    bl aff_pixel_codeRGB32
    cmp r5,r7                  @ x1 < x2 ?
    blt 8f
    mov r4,r6                  @  y = y1
    lsl r5,#7                  @ compute x1 * 128
    add r5,r10                 @ compute interx
7:
    lsr r1,r5,#7               @ compute x = ent ( interx / 128)
    lsl r3,r1,#7
    sub r6,r5,r3               @ brightness = remainder
    mov r0,r6
    bl colorPixel
    mov r3,r0
    mov r0,r1                  @ new x
    add r7,r0,#1
    mov r1,r4                  @ y
    bl aff_pixel_codeRGB32
    rsb r0,r6,#128
    bl colorPixel
    mov r3,r0
    mov r0,r7                  @ new x + 1
    mov r1,r4                  @ y
    bl aff_pixel_codeRGB32
    add r5,r10
    add r4,#1
    cmp r4,r8
    ble 7b
    b 100f
8:
    mov r4,r8                  @  y = y2
    lsl r5,#7                  @ compute x1 * 128
    add r5,r10                 @ compute interx
9:
    lsr r1,r5,#7               @ compute x
    lsl r3,r1,#7
    sub r8,r5,r3
    mov r0,r8
    bl colorPixel
    mov r3,r0
    mov r0,r1                  @ new x
    add r7,r0,#1
    mov r1,r4                  @ y
    bl aff_pixel_codeRGB32
    rsb r0,r8,#128
    bl colorPixel
    mov r3,r0
    mov r0,r7                  @ new x + 1
    mov r1,r4                  @ y
    bl aff_pixel_codeRGB32
    add r5,r10
    sub r4,#1
    cmp r4,r6
    bgt 9b
    b 100f
100:
    pop {fp,lr}
    bx lr
/********************************************************/
/*   brightness color pixel                              */
/********************************************************/
/* r0 % brightness ( 0 to 128)  */
colorPixel:
    push {r1,r2,lr}    /* save des  2 registres frame et retour */
    cmp r0,#0
    beq 100f
    cmp r0,#128
    mov r0,#127
    lsl r0,#1          @ red = brightness * 2 ( 2 to 254)
    mov r1,r0          @ green = red
    mov r2,r0          @ blue = red
    bl codeRGB         @ compute rgb code color 32 bits
100:
    pop {r1,r2,lr}
    bx lr 

/***************************************************/
/*   display pixels  32 bits                       */
/***************************************************/
/* r9 framebuffer memory address */
/* r0 = x */
/* r1 = y */
/* r2 screen width in pixels */
/* r3 code color RGB 32 bits  */
aff_pixel_codeRGB32:
    push {r0-r4,lr}       @  save registers
                          @ compute location pixel
    mul r4,r1,r2          @ compute y * screen width
    add r0,r0,r4          @ + x
    lsl r0,#2             @ * 4 octets
    str r3,[r9,r0]        @ store rgb code in mmap memory
    pop {r0-r4,lr}        @ restaur registers
    bx lr
/********************************************************/
/*   Code color RGB                                     */
/********************************************************/
/* r0 red r1 green  r2 blue */
/* r0 returns RGB code      */
codeRGB:
    lsl r0,#16               @ shift red color 16 bits
    lsl r1,#8                @ shift green color 8 bits
    eor r0,r1                @ or two colors
    eor r0,r2                @ or 3 colors in r0
    bx lr

/***************************************************/
/*      ROUTINES INCLUDE                 */
/***************************************************/
.include "./affichage.inc"

/***************************************************/
/*      DEFINITION DES STRUCTURES                 */
/***************************************************/
/* structure FSCREENINFO */    
/* voir explication détaillée : https://www.kernel.org/doc/Documentation/fb/api.txt */
    .struct  0
FBFIXSCinfo_id:          /* identification string eg "TT Builtin" */
    .struct FBFIXSCinfo_id + 16  
FBFIXSCinfo_smem_start:    /* Start of frame buffer mem */
    .struct FBFIXSCinfo_smem_start + 4   
FBFIXSCinfo_smem_len:       /* Length of frame buffer mem */
    .struct FBFIXSCinfo_smem_len + 4   
FBFIXSCinfo_type:    /* see FB_TYPE_*        */
    .struct FBFIXSCinfo_type + 4  
FBFIXSCinfo_type_aux:      /* Interleave for interleaved Planes */
    .struct FBFIXSCinfo_type_aux + 4  
FBFIXSCinfo_visual:    /* see FB_VISUAL_*        */
    .struct FBFIXSCinfo_visual + 4  
FBFIXSCinfo_xpanstep:    /* zero if no hardware panning  */
    .struct FBFIXSCinfo_xpanstep + 2      
FBFIXSCinfo_ypanstep:    /* zero if no hardware panning  */
    .struct FBFIXSCinfo_ypanstep + 2 
FBFIXSCinfo_ywrapstep:      /* zero if no hardware ywrap    */
    .struct FBFIXSCinfo_ywrapstep + 4 
FBFIXSCinfo_line_length:    /* length of a line in bytes    */
    .struct FBFIXSCinfo_line_length + 4 
FBFIXSCinfo_mmio_start:     /* Start of Memory Mapped I/O   */
    .struct FBFIXSCinfo_mmio_start + 4     
FBFIXSCinfo_mmio_len:        /* Length of Memory Mapped I/O  */
    .struct FBFIXSCinfo_mmio_len + 4 
FBFIXSCinfo_accel:     /* Indicate to driver which    specific chip/card we have    */
    .struct FBFIXSCinfo_accel + 4 
FBFIXSCinfo_capabilities:     /* see FB_CAP_*            */
    .struct FBFIXSCinfo_capabilities + 4 
FBFIXSCinfo_reserved:     /* Reserved for future compatibility */
    .struct FBFIXSCinfo_reserved + 8    
FBFIXSCinfo_fin:

/* structure VSCREENINFO */    
    .struct  0
FBVARSCinfo_xres:           /* visible resolution        */ 
    .struct FBVARSCinfo_xres + 4  
FBVARSCinfo_yres:          
    .struct FBVARSCinfo_yres + 4 
FBVARSCinfo_xres_virtual:          /* virtual resolution        */
    .struct FBVARSCinfo_xres_virtual + 4 
FBVARSCinfo_yres_virtual:          
    .struct FBVARSCinfo_yres_virtual + 4 
FBVARSCinfo_xoffset:          /* offset from virtual to visible resolution */
    .struct FBVARSCinfo_xoffset + 4 
FBVARSCinfo_yoffset:          
    .struct FBVARSCinfo_yoffset + 4 
FBVARSCinfo_bits_per_pixel:          /* bits par pixel */
    .struct FBVARSCinfo_bits_per_pixel + 4     
FBVARSCinfo_grayscale:          /* 0 = color, 1 = grayscale,  >1 = FOURCC    */
    .struct FBVARSCinfo_grayscale + 4 
FBVARSCinfo_red:          /* bitfield in fb mem if true color, */
    .struct FBVARSCinfo_red + 4 
FBVARSCinfo_green:          /* else only length is significant */
    .struct FBVARSCinfo_green + 4 
FBVARSCinfo_blue:          
    .struct FBVARSCinfo_blue + 4 
FBVARSCinfo_transp:          /* transparency            */
    .struct FBVARSCinfo_transp + 4     
FBVARSCinfo_nonstd:          /* != 0 Non standard pixel format */
    .struct FBVARSCinfo_nonstd + 4 
FBVARSCinfo_activate:          /* see FB_ACTIVATE_*        */
    .struct FBVARSCinfo_activate + 4     
FBVARSCinfo_height:              /* height of picture in mm    */
    .struct FBVARSCinfo_height + 4 
FBVARSCinfo_width:           /* width of picture in mm     */
    .struct FBVARSCinfo_width + 4 
FBVARSCinfo_accel_flags:          /* (OBSOLETE) see fb_info.flags */
    .struct FBVARSCinfo_accel_flags + 4 
/* Timing: All values in pixclocks, except pixclock (of course) */    
FBVARSCinfo_pixclock:          /* pixel clock in ps (pico seconds) */
    .struct FBVARSCinfo_pixclock + 4     
FBVARSCinfo_left_margin:          
    .struct FBVARSCinfo_left_margin + 4 
FBVARSCinfo_right_margin:          
    .struct FBVARSCinfo_right_margin + 4 
FBVARSCinfo_upper_margin:          
    .struct FBVARSCinfo_upper_margin + 4 
FBVARSCinfo_lower_margin:          
    .struct FBVARSCinfo_lower_margin + 4 
FBVARSCinfo_hsync_len:          /* length of horizontal sync    */
    .struct FBVARSCinfo_hsync_len + 4     
FBVARSCinfo_vsync_len:          /* length of vertical sync    */
    .struct FBVARSCinfo_vsync_len + 4 
FBVARSCinfo_sync:          /* see FB_SYNC_*        */
    .struct FBVARSCinfo_sync + 4 
FBVARSCinfo_vmode:          /* see FB_VMODE_*        */
    .struct FBVARSCinfo_vmode + 4 
FBVARSCinfo_rotate:          /* angle we rotate counter clockwise */
    .struct FBVARSCinfo_rotate + 4     
FBVARSCinfo_colorspace:          /* colorspace for FOURCC-based modes */
    .struct FBVARSCinfo_colorspace + 4     
FBVARSCinfo_reserved:          /* Reserved for future compatibility */
    .struct FBVARSCinfo_reserved + 16        
FBVARSCinfo_fin:

ATS

Translation of: ObjectIcon

The following program writes (to standard output) a Portable Grayscale Map without gamma correction. One can use Netpbm to do the gamma correction. Thus, for example:

$ patscc -g -O2 -std=gnu2x -DATS_MEMALLOC_GCBDW xiaolin_wu_line_algorithm.dats -lgc -lm
$ ./a.out | pnmgamma -lineartobt709 | pnmtopng > image.png

For mathematics, I simply make foreign language calls to C. There is some mathematics support in libats/libc, and ats2-xprelude provides even more. However, making the foreign language calls is trivial, as you can see below.

#include "share/atspre_staload.hats"
staload UN = "prelude/SATS/unsafe.sats"

%{^
#include <float.h>
#include <math.h>
%}

typedef color = double

typedef drawing_surface (u0 : int, v0 : int,
                         u1 : int, v1 : int) =
  [u0 <= u1; v0 <= v1]
  '{
    u0 = int u0,
    v0 = int v0,
    u1 = int u1,
    v1 = int v1,
    pixels = matrixref (color, (u1 - u0) + 1, (v1 - v0) + 1)
  }

fn
drawing_surface_make
          {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
          (u0 : int u0,
           v0 : int v0,
           u1 : int u1,
           v1 : int v1)
    : drawing_surface (u0, v0, u1, v1) =
  let
    val w = succ (u1 - u0) and h = succ (v1 - v0)
  in
    '{
      u0 = u0,
      v0 = v0,
      u1 = u1,
      v1 = v1,
      pixels = matrixref_make_elt (i2sz w, i2sz h, 0.0)
    }
  end

fn
drawing_surface_get_at
          {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
          (s : drawing_surface (u0, v0, u1, v1),
           x : int,
           y : int)
    : color =
  let
    val '{ u0 = u0, v0 = v0, u1 = u1, v1 = v1, pixels = pixels } = s
    and x = g1ofg0 x and y = g1ofg0 y
  in
    if (u0 <= x) * (x <= u1) * (v0 <= y) * (y <= v1) then
      (* Notice there are THREE values in the square brackets. The
         matrixref does not store its dimensions and so we have to
         specify a stride as the second value. The value must,
         however, be the CORRECT stride. This is checked at compile
         time. (Here ATS is striving to be efficient rather than
         convenient!) *)
      pixels[x - u0, succ (v1 - v0), v1 - y]
    else
      $extval (double, "nan (\"0\")")
  end

fn
drawing_surface_set_at
          {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
          (s : drawing_surface (u0, v0, u1, v1),
           x : int,
           y : int,
           c : color)
    : void =
  let
    val '{ u0 = u0, v0 = v0, u1 = u1, v1 = v1, pixels = pixels } = s
    and x = g1ofg0 x and y = g1ofg0 y
  in
    if (u0 <= x) * (x <= u1) * (v0 <= y) * (y <= v1) then
      pixels[x - u0, succ (v1 - v0), v1 - y] := c
  end

(* Indices outside the drawing_surface are allowed. They are handled
   by treating them as if you were trying to draw on the air. *)
overload [] with drawing_surface_get_at
overload [] with drawing_surface_set_at

fn
write_PGM {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
          (outf : FILEref,
           s    : drawing_surface (u0, v0, u1, v1))
    : void =
  (* Write a Portable Grayscale Map. *)
  let
    val '{ u0 = u0, v0 = v0, u1 = u1, v1 = v1, pixels = pixels } = s

    stadef w = (u1 - u0) + 1
    stadef h = (v1 - v0) + 1
    val w : int w = succ (u1 - u0)
    and h : int h = succ (v1 - v0)

    fun
    loop {x, y : int | 0 <= x; x <= w; 0 <= y; y <= h}
         .<h - y, w - x>.
         (x : int x,
          y : int y) : void =
      if y = h then
        ()
      else if x = w then
        loop (0, succ y)
      else
        let
          (* I do no gamma correction, but gamma correction can be
             done afterwards by running the output through "pnmgamma
             -lineartobt709" *)
          val intensity = 1.0 - pixels[x, h, y]
          val pgm_value : int =
            g0f2i ($extfcall (double, "rint", 65535 * intensity))
          val more_significant_byte = pgm_value / 256
          and less_significant_byte = pgm_value mod 256
          val msbyte = int2uchar0 more_significant_byte
          and lsbyte = int2uchar0 less_significant_byte
        in
          fprint_val<uchar> (outf, msbyte);
          fprint_val<uchar> (outf, lsbyte);
          loop (succ x, y)
        end
  in
    fprintln! (outf, "P5");
    fprintln! (outf, w, " ", h);
    fprintln! (outf, "65535");
    loop (0, 0)
  end

fn
ipart (x : double) : int =
  g0f2i ($extfcall (double, "floor", x))

fn
iround (x : double) : int =
  ipart (x + 0.5)

fn
fpart (x : double) : double =
  x - $extfcall (double, "floor", x)

fn
rfpart (x : double) : double =
  1.0 - fpart (x)

fn
plot {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
     (s : drawing_surface (u0, v0, u1, v1),
      x : int,
      y : int,
      c : color)
    : void =
  let
    (* One might prefer a more sophisticated function than mere
       addition. *)
    val combined_color = s[x, y] + c
  in
    s[x, y] := min (combined_color, 1.0)
  end

fn
_drawln {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
        (s     : drawing_surface (u0, v0, u1, v1),
         x0    : double,
         y0    : double,
         x1    : double,
         y1    : double,
         steep : bool)
    : void =
  let
    val dx = x1 - x0 and dy = y1 - y0
    val gradient = (if dx = 0.0 then 1.0 else dy / dx) : double

    (* Handle the first endpoint. *)
    val xend = iround x0
    val yend = y0 + (gradient * (g0i2f xend - x0))
    val xgap = rfpart (x0 + 0.5)
    val xpxl1 = xend and ypxl1 = ipart yend
    val () =
      if steep then
        begin
          plot (s, ypxl1, xpxl1, rfpart yend * xgap);
          plot (s, succ ypxl1, xpxl1, fpart yend * xgap)
        end
      else
        begin
          plot (s, xpxl1, ypxl1, rfpart yend * xgap);
          plot (s, xpxl1, succ ypxl1, fpart yend * xgap)
        end

    (* The first y-intersection. Notice it is a "var" (a variable)
       instead of a "val" (an immutable value). There is no need to
       box it as a "ref", the way one must typically do in an ML
       dialect. We could have done so, but the following treats the
       variable as an ordinary C automatic variable, and is more
       efficient. *)
    var intery : double = yend + gradient

    (* Handle the second endpoint. *)
    val xend = iround (x1)
    val yend = y1 + (gradient * (g0i2f xend - x1))
    val xgap = fpart (x1 + 0.5)
    val xpxl2 = xend and ypxl2 = ipart yend
    val () =
      if steep then
        begin
          plot (s, ypxl2, xpxl2, rfpart yend * xgap);
          plot (s, succ ypxl2, xpxl2, fpart yend * xgap)
        end
      else
        begin
          plot (s, xpxl2, ypxl2, rfpart yend * xgap);
          plot (s, xpxl2, succ ypxl2, fpart yend * xgap)
      end
  in
    (* Loop over the rest of the points. I use procedural "for"-loops
       instead of the more usual (for ATS) tail recursion. *)
    if steep then
      let
        var x : int
      in
        for (x := succ xpxl1; x <> xpxl2; x := succ x)
          begin
            plot (s, ipart intery, x, rfpart intery);
            plot (s, succ (ipart intery), x, fpart intery);
            intery := intery + gradient
          end
      end
    else
      let
        var x : int
      in
        for (x := succ xpxl1; x <> xpxl2; x := succ x)
          begin
            plot (s, x, ipart intery, rfpart intery);
            plot (s, x, succ (ipart intery), fpart intery);
            intery := intery + gradient
          end
      end
  end

fn
draw_line {u0, v0, u1, v1 : int | u0 <= u1; v0 <= v1}
          (s  : drawing_surface (u0, v0, u1, v1),
           x0 : double,
           y0 : double,
           x1 : double,
           y1 : double)
    : void =
  let
    val xdiff = abs (x1 - x0) and ydiff = abs (y1 - y0)
  in
    if ydiff <= xdiff then
      begin
        if x0 <= x1 then
          _drawln (s, x0, y0, x1, y1, false)
        else
          _drawln (s, x1, y1, x0, y0, false)
      end
    else
      begin
        if y0 <= y1 then
          _drawln (s, y0, x0, y1, x1, true)
        else
          _drawln (s, y1, x1, y0, x0, true)
      end
  end

implement
main0 () =
  let
    macdef M_PI = $extval (double, "M_PI")

    val u0 = 0 and v0 = 0
    and u1 = 639 and v1 = 479

    val s = drawing_surface_make (u0, v0, u1, v1)

    fun
    loop (theta : double) : void =
      if theta < 360.0 then
        let
          val cos_theta = $extfcall (double, "cos",
                                     theta * (M_PI / 180.0))
          and sin_theta = $extfcall (double, "sin",
                                     theta * (M_PI / 180.0))
          and x0 = 380.0 and y0 = 130.0
          val x1 = x0 + (cos_theta * 1000.0)
          and y1 = y0 + (sin_theta * 1000.0)
        in
          draw_line (s, x0, y0, x1, y1);
          loop (theta + 5.0)
        end
  in
    loop 0.0;
    write_PGM (stdout_ref, s)
  end
Output:

The following image has been gamma-corrected with pnmgamma (although the thumbnail image may look strange).

A starburst of straight lines, black on white.

AutoHotkey

Library: GDIP
#SingleInstance, Force
#NoEnv
SetBatchLines, -1

pToken := Gdip_Startup()
global pBitmap := Gdip_CreateBitmap(500, 500)
drawLine(100,50,400,400)
Gdip_SaveBitmapToFile(pBitmap, A_ScriptDir "\linetest.png")
Gdip_DisposeImage(pBitmap)
Gdip_Shutdown(pToken)
Run, % A_ScriptDir "\linetest.png"
ExitApp

plot(x, y, c) {
    A := DecToBase(255 * c, 16)
    Gdip_SetPixel(pBitmap, x, y, "0x" A "000000")
}
 
; integer part of x
ipart(x) {
    return x // 1
}
 
rnd(x) {
    return ipart(x + 0.5)
}
 
; fractional part of x
fpart(x) {
    if (x < 0)
        return 1 - (x - floor(x))
    return x - floor(x)
}
 
rfpart(x) {
    return 1 - fpart(x)
}
 
drawLine(x0,y0,x1,y1) {
    steep := abs(y1 - y0) > abs(x1 - x0)
 
    if (steep) {
        temp := x0, x0 := y0, y0 := temp
        temp := x1, x1 := y1, y1 := temp
    }
    if (x0 > x1 then) {
        temp := x0, x0 := x1, x1 := temp
        temp := y0, y0 := y1, y1 := temp
    }
 
    dx := x1 - x0
    dy := y1 - y0
    gradient := dy / dx
 
    ; handle first endpoint
    xend := rnd(x0)
    yend := y0 + gradient * (xend - x0)
    xgap := rfpart(x0 + 0.5)
    xpxl1 := xend ; this will be used in the main loop
    ypxl1 := ipart(yend)
    if (steep) {
        plot(ypxl1,   xpxl1, rfpart(yend) * xgap)
        plot(ypxl1+1, xpxl1,  fpart(yend) * xgap)
    }   
    else {
        plot(xpxl1, ypxl1  , rfpart(yend) * xgap)
        plot(xpxl1, ypxl1+1,  fpart(yend) * xgap)
    }
    intery := yend + gradient ; first y-intersection for the main loop
 
    ; handle second endpoint
    xend := rnd(x1)
    yend := y1 + gradient * (xend - x1)
    xgap := fpart(x1 + 0.5)
    xpxl2 := xend ;this will be used in the main loop
    ypxl2 := ipart(yend)
    if (steep) {
        plot(ypxl2  , xpxl2, rfpart(yend) * xgap)
        plot(ypxl2+1, xpxl2,  fpart(yend) * xgap)
    }
    else {
        plot(xpxl2, ypxl2,  rfpart(yend) * xgap)
        plot(xpxl2, ypxl2+1, fpart(yend) * xgap)
    }
 
    ; main loop
    while (x := xpxl1 + A_Index) < xpxl2 {
        if (steep) {
            plot(ipart(intery)  , x, rfpart(intery))
            plot(ipart(intery)+1, x,  fpart(intery))
        }
        else {
            plot(x, ipart (intery),  rfpart(intery))
            plot(x, ipart (intery)+1, fpart(intery))
        }
        intery := intery + gradient
    }
}

DecToBase(n, Base) {
    static U := A_IsUnicode ? "w" : "a"
    VarSetCapacity(S,65,0)
    DllCall("msvcrt\_i64to" U, "Int64",n, "Str",S, "Int",Base)
    return, S
}

BBC BASIC

      PROCdrawAntiAliasedLine(100, 100, 600, 400, 0, 0, 0)
      END
      
      DEF PROCdrawAntiAliasedLine(x1, y1, x2, y2, r%, g%, b%)
      LOCAL dx, dy, xend, yend, grad, yf, xgap, ix1%, iy1%, ix2%, iy2%, x%
      
      dx = x2 - x1
      dy = y2 - y1
      IF ABS(dx) < ABS(dy) THEN
        SWAP x1, y1
        SWAP x2, y2
        SWAP dx, dy
      ENDIF
      
      IF x2 < x1 THEN
        SWAP x1, x2
        SWAP y1, y2
      ENDIF
      
      grad = dy / dx
      
      xend = INT(x1 + 0.5)
      yend = y1 + grad * (xend - x1)
      xgap = xend + 0.5 - x1
      ix1% = xend
      iy1% = INT(yend)
      PROCplot(ix1%, iy1%, r%, b%, g%, (INT(yend) + 1 - yend) * xgap)
      PROCplot(ix1%, iy1% + 1, r%, b%, g%, (yend - INT(yend)) * xgap)
      yf = yend + grad
      
      xend = INT(x2 + 0.5)
      yend = y2 + grad * (xend - x2)
      xgap = x2 + 0.5 - xend
      ix2% = xend
      iy2% = INT(yend)
      PROCplot(ix2%, iy2%, r%, b%, g%, (INT(yend) + 1 - yend) * xgap)
      PROCplot(ix2%, iy2% + 1, r%, b%, g%, (yend - INT(yend)) * xgap)
      
      FOR x% = ix1% + 1 TO ix2% - 1
        PROCplot(x%, INT(yf), r%, b%, g%, INT(yf) + 1 - yf)
        PROCplot(x%, INT(yf) + 1, r%, b%, g%, yf - INT(yf))
        yf += grad
      NEXT
      ENDPROC
      
      DEF PROCplot(X%, Y%, R%, G%, B%, a)
      LOCAL C%
      C% = TINT(X%*2,Y%*2)
      COLOUR 1, R%*a + (C% AND 255)*(1-a), \
      \         G%*a + (C% >> 8 AND 255)*(1-a), \
      \         B%*a + (C% >> 16 AND 255)*(1-a)
      GCOL 1
      LINE X%*2, Y%*2, X%*2, Y%*2
      ENDPROC

C fast fixed-point

This is an implementation in C using fixed-point to speed things up significantly. Suitable for 32-bit+ architectures. For reference and comparison, the floating-point version is also included.

This implementation of plot() only draws white on a fixed canvas, but can easily be modified.

// Something to draw on
static uint8_t canvas[240][240];

// Paint pixel white
static void inline plot(int16_t x, int16_t y, uint16_t alpha) {
	canvas[y][x] = 255 - (((255 - canvas[y][x]) * (alpha & 0x1FF)) >> 8);
}

// Xiaolin Wu's line algorithm
// Coordinates are Q16 fixed point, ie 0x10000 == 1
void wuline(int32_t x0, int32_t y0, int32_t x1, int32_t y1) {
	bool steep = ((y1 > y0) ? (y1 - y0) : (y0 - y1)) > ((x1 > x0) ? (x1 - x0) : (x0 - x1));

	if(steep)   { int32_t z = x0; x0 = y0; y0 = z; z = x1; x1 = y1; y1 = z; }
	if(x0 > x1) { int32_t z = x0; x0 = x1; x1 = z; z = y0; y0 = y1; y1 = z; }

	int32_t dx = x1 - x0, dy = y1 - y0;
	int32_t gradient = ((dx >> 8) == 0) ? 0x10000 : (dy << 8) / (dx >> 8);

	// handle first endpoint
	int32_t xend = (x0 + 0x8000) & 0xFFFF0000;
	int32_t yend = y0 + ((gradient * (xend - x0)) >> 16);
	int32_t xgap = 0x10000 - ((x0 + 0x8000) & 0xFFFF);
	int16_t xpxl1 = xend >> 16; // this will be used in the main loop
	int16_t ypxl1 = yend >> 16;
	if(steep) {
		plot(ypxl1,     xpxl1,     0x100 - (((0x100 - ((yend >> 8) & 0xFF)) * xgap) >> 16));
		plot(ypxl1 + 1, xpxl1,     0x100 - ((         ((yend >> 8) & 0xFF)  * xgap) >> 16));
	} else {
		plot(xpxl1,     ypxl1,     0x100 - (((0x100 - ((yend >> 8) & 0xFF)) * xgap) >> 16));
		plot(xpxl1,     ypxl1 + 1, 0x100 - ((         ((yend >> 8) & 0xFF)  * xgap) >> 16));
	}
	int32_t intery = yend + gradient; // first y-intersection for the main loop

	// handle second endpoint
	xend = (x1 + 0x8000) & 0xFFFF0000;
	yend = y1 + ((gradient * (xend - x1)) >> 16);
	xgap = (x1 + 0x8000) & 0xFFFF;
	int16_t xpxl2 = xend >> 16; //this will be used in the main loop
	int16_t ypxl2 = yend >> 16;
	if(steep) {
		plot(ypxl2,     xpxl2,     0x100 - (((0x100 - ((yend >> 8) & 0xFF)) * xgap) >> 16));
		plot(ypxl2 + 1, xpxl2,     0x100 - ((         ((yend >> 8) & 0xFF)  * xgap) >> 16));
	} else {                       
		plot(xpxl2,     ypxl2,     0x100 - (((0x100 - ((yend >> 8) & 0xFF)) * xgap) >> 16));
		plot(xpxl2,     ypxl2 + 1, 0x100 - ((         ((yend >> 8) & 0xFF)  * xgap) >> 16));
	}

	// main loop
	if(steep) {
		for(int32_t x = xpxl1 + 1; x < xpxl2; x++) {
			plot((intery >> 16)    , x,          (intery >> 8) & 0xFF );
			plot((intery >> 16) + 1, x, 0x100 - ((intery >> 8) & 0xFF));
			intery += gradient;
		}
	} else {
		for(int32_t x = xpxl1 + 1; x < xpxl2; x++) {
			plot(x, (intery >> 16),              (intery >> 8) & 0xFF );
			plot(x, (intery >> 16) + 1, 0x100 - ((intery >> 8) & 0xFF));
			intery += gradient;
		}
	}
}

// Paint pixel white (floating point version, for reference only)
static void inline plotf(int16_t x, int16_t y, float alpha) {
	canvas[y][x] = 255 - ((255 - canvas[y][x]) * (1.0 - alpha));
}

// Xiaolin Wu's line algorithm (floating point version, for reference only)
void wulinef(float x0, float y0, float x1, float y1) {
	bool steep = fabs(y1 - y0) > fabs(x1 - x0);

	if(steep)   { float z = x0; x0 = y0; y0 = z; z = x1; x1 = y1; y1 = z; }
	if(x0 > x1) { float z = x0; x0 = x1; x1 = z; z = y0; y0 = y1; y1 = z; }

	float dx = x1 - x0, dy = y1 - y0;
	float gradient = (dx == 0.0) ? 1.0 : dy / dx;

	// handle first endpoint
	uint16_t xend = round(x0);
	float yend = y0 + gradient * ((float)xend - x0);
	float xgap = 1.0 - (x0 + 0.5 - floor(x0 + 0.5));
	int16_t xpxl1 = xend; // this will be used in the main loop
	int16_t ypxl1 = floor(yend);
	if(steep) {
		plotf(ypxl1,       xpxl1, (1.0 - (yend - floor(yend))) * xgap);
		plotf(ypxl1 + 1.0, xpxl1,        (yend - floor(yend))  * xgap);
	} else {
		plotf(xpxl1, ypxl1,       (1.0 - (yend - floor(yend))) * xgap);
		plotf(xpxl1, ypxl1 + 1.0,        (yend - floor(yend))  * xgap);
	}
	float intery = yend + gradient; // first y-intersection for the main loop

	// handle second endpoint
	xend = round(x1);
	yend = y1 + gradient * ((float)xend - x1);
	xgap = x1 + 0.5 - floor(x1 + 0.5);
	int16_t xpxl2 = xend; //this will be used in the main loop
	int16_t ypxl2 = floor(yend);
	if(steep) {
		plotf(ypxl2,       xpxl2, (1.0 - (yend - floor(yend))) * xgap);
		plotf(ypxl2 + 1.0, xpxl2,        (yend - floor(yend))  * xgap);
	} else {
		plotf(xpxl2, ypxl2,       (1.0 - (yend - floor(yend))) * xgap);
		plotf(xpxl2, ypxl2 + 1.0,        (yend - floor(yend))  * xgap);
	}

	// main loop
	if(steep) {
		for(uint16_t x = xpxl1 + 1; x < xpxl2; x++) {
			plotf(floor(intery),     x, (1.0 - (intery - floor(intery))));
			plotf(floor(intery) + 1, x,        (intery - floor(intery) ));
			intery += gradient;
		}
	} else {
		for(uint16_t x = xpxl1 + 1; x < xpxl2; x++) {
			plotf(x, floor(intery),     (1.0 - (intery - floor(intery))));
			plotf(x, floor(intery) + 1,        (intery - floor(intery) ));
			intery += gradient;
		}
	}
}

void wudemo() {

	// Clear the canvas
	memset(canvas, 0, sizeof(canvas));

	// Of course it doesn't make sense to use slow floating point trig. functions here
	// This is just for demo purposes
	static float wudemo_v;
	wudemo_v += 0.005;
	float x = sinf(wudemo_v) * 50;
	float y = cosf(wudemo_v) * 50;

	// Draw using fast fixed-point version
	wuline ((x + 125) * (1 << 16), (y + 125) * (1 << 16), (-x + 125) * (1 << 16), (-y + 125) * (1 << 16));

	// Draw using reference version for comparison
	wulinef(  x + 115,                y + 115,               -x + 115,                -y + 115             );

	// -- insert display code here --
	showme(canvas);

}

C

This implementation follows straightforwardly the pseudocode given on Wikipedia. (Further analysis of the code could give suggestions for improvements).

void draw_line_antialias(
        image img,
        unsigned int x0, unsigned int y0,
        unsigned int x1, unsigned int y1,
        color_component r,
        color_component g,
        color_component b );
inline void _dla_changebrightness(rgb_color_p from,
				  rgb_color_p to, float br)
{
  if ( br > 1.0 ) br = 1.0;
  /* linear... Maybe something more complex could give better look */
  to->red = br * (float)from->red;
  to->green = br * (float)from->green;
  to->blue = br * (float)from->blue;
}

#define plot_(X,Y,D) do{ rgb_color f_;				\
  f_.red = r; f_.green = g; f_.blue = b;			\
  _dla_plot(img, (X), (Y), &f_, (D)) ; }while(0)

inline void _dla_plot(image img, int x, int y, rgb_color_p col, float br)
{
  rgb_color oc;
  _dla_changebrightness(col, &oc, br);
  put_pixel_clip(img, x, y, oc.red, oc.green, oc.blue);
}

#define ipart_(X) ((int)(X))
#define round_(X) ((int)(((double)(X))+0.5))
#define fpart_(X) (((double)(X))-(double)ipart_(X))
#define rfpart_(X) (1.0-fpart_(X))

#define swap_(a, b) do{ __typeof__(a) tmp;  tmp = a; a = b; b = tmp; }while(0)
void draw_line_antialias(
  image img,
  unsigned int x1, unsigned int y1,
  unsigned int x2, unsigned int y2,
  color_component r,
  color_component g,
  color_component b )
{
  double dx = (double)x2 - (double)x1;
  double dy = (double)y2 - (double)y1;
  if ( fabs(dx) > fabs(dy) ) {
    if ( x2 < x1 ) {
      swap_(x1, x2);
      swap_(y1, y2);
    }
    double gradient = dy / dx;
    double xend = round_(x1);
    double yend = y1 + gradient*(xend - x1);
    double xgap = rfpart_(x1 + 0.5);
    int xpxl1 = xend;
    int ypxl1 = ipart_(yend);
    plot_(xpxl1, ypxl1, rfpart_(yend)*xgap);
    plot_(xpxl1, ypxl1+1, fpart_(yend)*xgap);
    double intery = yend + gradient;

    xend = round_(x2);
    yend = y2 + gradient*(xend - x2);
    xgap = fpart_(x2+0.5);
    int xpxl2 = xend;
    int ypxl2 = ipart_(yend);
    plot_(xpxl2, ypxl2, rfpart_(yend) * xgap);
    plot_(xpxl2, ypxl2 + 1, fpart_(yend) * xgap);

    int x;
    for(x=xpxl1+1; x < xpxl2; x++) {
      plot_(x, ipart_(intery), rfpart_(intery));
      plot_(x, ipart_(intery) + 1, fpart_(intery));
      intery += gradient;
    }
  } else {
    if ( y2 < y1 ) {
      swap_(x1, x2);
      swap_(y1, y2);
    }
    double gradient = dx / dy;
    double yend = round_(y1);
    double xend = x1 + gradient*(yend - y1);
    double ygap = rfpart_(y1 + 0.5);
    int ypxl1 = yend;
    int xpxl1 = ipart_(xend);
    plot_(xpxl1, ypxl1, rfpart_(xend)*ygap);
    plot_(xpxl1 + 1, ypxl1, fpart_(xend)*ygap);
    double interx = xend + gradient;

    yend = round_(y2);
    xend = x2 + gradient*(yend - y2);
    ygap = fpart_(y2+0.5);
    int ypxl2 = yend;
    int xpxl2 = ipart_(xend);
    plot_(xpxl2, ypxl2, rfpart_(xend) * ygap);
    plot_(xpxl2 + 1, ypxl2, fpart_(xend) * ygap);

    int y;
    for(y=ypxl1+1; y < ypxl2; y++) {
      plot_(ipart_(interx), y, rfpart_(interx));
      plot_(ipart_(interx) + 1, y, fpart_(interx));
      interx += gradient;
    }
  }
}
#undef swap_
#undef plot_
#undef ipart_
#undef fpart_
#undef round_
#undef rfpart_

C++

#include <functional>
#include <algorithm>
#include <utility>

void WuDrawLine(float x0, float y0, float x1, float y1,
                const std::function<void(int x, int y, float brightess)>& plot) {
    auto ipart = [](float x) -> int {return int(std::floor(x));};
    auto round = [](float x) -> float {return std::round(x);};
    auto fpart = [](float x) -> float {return x - std::floor(x);};
    auto rfpart = [=](float x) -> float {return 1 - fpart(x);};
        
    const bool steep = abs(y1 - y0) > abs(x1 - x0);
    if (steep) {
        std::swap(x0,y0);
        std::swap(x1,y1);
    }
    if (x0 > x1) {
        std::swap(x0,x1);
        std::swap(y0,y1);
    }
        
    const float dx = x1 - x0;
    const float dy = y1 - y0;
    const float gradient = (dx == 0) ? 1 : dy/dx;
        
    int xpx11;
    float intery;
    {
        const float xend = round(x0);
        const float yend = y0 + gradient * (xend - x0);
        const float xgap = rfpart(x0 + 0.5);
        xpx11 = int(xend);
        const int ypx11 = ipart(yend);
        if (steep) {
            plot(ypx11,     xpx11, rfpart(yend) * xgap);
            plot(ypx11 + 1, xpx11,  fpart(yend) * xgap);
        } else {
            plot(xpx11, ypx11,    rfpart(yend) * xgap);
            plot(xpx11, ypx11 + 1, fpart(yend) * xgap);
        }
        intery = yend + gradient;
    }
    
    int xpx12;
    {
        const float xend = round(x1);
        const float yend = y1 + gradient * (xend - x1);
        const float xgap = rfpart(x1 + 0.5);
        xpx12 = int(xend);
        const int ypx12 = ipart(yend);
        if (steep) {
            plot(ypx12,     xpx12, rfpart(yend) * xgap);
            plot(ypx12 + 1, xpx12,  fpart(yend) * xgap);
        } else {
            plot(xpx12, ypx12,    rfpart(yend) * xgap);
            plot(xpx12, ypx12 + 1, fpart(yend) * xgap);
        }
    }
        
    if (steep) {
        for (int x = xpx11 + 1; x < xpx12; x++) {
            plot(ipart(intery),     x, rfpart(intery));
            plot(ipart(intery) + 1, x,  fpart(intery));
            intery += gradient;
        }
    } else {
        for (int x = xpx11 + 1; x < xpx12; x++) {
            plot(x, ipart(intery),     rfpart(intery));
            plot(x, ipart(intery) + 1,  fpart(intery));
            intery += gradient;
        }
    }
}

C#

public class Line
    {
        private double x0, y0, x1, y1;
        private Color foreColor;
        private byte lineStyleMask;
        private int thickness;
        private float globalm;

        public Line(double x0, double y0, double x1, double y1, Color color, byte lineStyleMask, int thickness)
        {
            this.x0 = x0;
            this.y0 = y0;
            this.y1 = y1;
            this.x1 = x1;

            this.foreColor = color;

            this.lineStyleMask = lineStyleMask;

            this.thickness = thickness;

        }

        private void plot(Bitmap bitmap, double x, double y, double c)
        {
            int alpha = (int)(c * 255);
            if (alpha > 255) alpha = 255;
            if (alpha < 0) alpha = 0;
            Color color = Color.FromArgb(alpha, foreColor);
            if (BitmapDrawHelper.checkIfInside((int)x, (int)y, bitmap))
            {
                bitmap.SetPixel((int)x, (int)y, color);
            }
        }

        int ipart(double x) { return (int)x;}

        int round(double x) {return ipart(x+0.5);}
    
        double fpart(double x) {
            if(x<0) return (1-(x-Math.Floor(x)));
            return (x-Math.Floor(x));
        }
    
        double rfpart(double x) {
            return 1-fpart(x);
        }


        public void draw(Bitmap bitmap) {
            bool steep = Math.Abs(y1-y0)>Math.Abs(x1-x0);
             double temp;
            if(steep){
                temp=x0; x0=y0; y0=temp;
                temp=x1;x1=y1;y1=temp;
            }
            if(x0>x1){
                temp = x0;x0=x1;x1=temp;
                temp = y0;y0=y1;y1=temp;
            }

            double dx = x1-x0;
            double dy = y1-y0;
            double gradient = dy/dx;

            double xEnd = round(x0);
            double yEnd = y0+gradient*(xEnd-x0);
            double xGap = rfpart(x0+0.5);
            double xPixel1 = xEnd;
            double yPixel1 = ipart(yEnd);

            if(steep){
                plot(bitmap, yPixel1,   xPixel1, rfpart(yEnd)*xGap);
                plot(bitmap, yPixel1+1, xPixel1,  fpart(yEnd)*xGap);
            }else{
                plot(bitmap, xPixel1,yPixel1, rfpart(yEnd)*xGap);
                plot(bitmap, xPixel1, yPixel1+1, fpart(yEnd)*xGap);
            }
            double intery = yEnd+gradient;

            xEnd = round(x1);
            yEnd = y1+gradient*(xEnd-x1);
            xGap = fpart(x1+0.5);
            double xPixel2 = xEnd;
            double yPixel2 = ipart(yEnd);
            if(steep){
                plot(bitmap, yPixel2,   xPixel2, rfpart(yEnd)*xGap);
                plot(bitmap, yPixel2+1, xPixel2, fpart(yEnd)*xGap);
            }else{
                plot(bitmap, xPixel2, yPixel2, rfpart(yEnd)*xGap);
                plot(bitmap, xPixel2, yPixel2+1, fpart(yEnd)*xGap);
            }

            if(steep){
                for(int x=(int)(xPixel1+1);x<=xPixel2-1;x++){
                    plot(bitmap, ipart(intery), x, rfpart(intery));
                    plot(bitmap, ipart(intery)+1, x, fpart(intery));
                    intery+=gradient;
                }
            }else{
                for(int x=(int)(xPixel1+1);x<=xPixel2-1;x++){
                    plot(bitmap, x,ipart(intery), rfpart(intery));
                    plot(bitmap, x, ipart(intery)+1, fpart(intery));
                    intery+=gradient;
                }
            }
        }
    }

Common Lisp

Translation of: Scheme

The program outputs a transparency mask in Portable Gray Map format. It draws lines normal to a tractrix. They trace out a catenary.

;;; The program outputs a transparency mask in plain Portable Gray Map
;;; format.
;;; -------------------------------------------------------------------

(defstruct (drawing-surface
            (:constructor make-drawing-surface (u0 v0 u1 v1)))
  (u0 0 :type fixnum :read-only t)
  (v0 0 :type fixnum :read-only t)
  (u1 0 :type fixnum :read-only t)
  (v1 0 :type fixnum :read-only t)
  (pixels (make-array (* (- u1 u0 -1) (- v1 v0 -1))
                      :element-type 'single-float
                      :initial-element 0.0)))

;;; In calls to drawing-surface-ref and drawing-surface-set, indices
;;; outside the drawing_surface are allowed. Such indices are treated
;;; as if you were trying to draw on the air.

(defun drawing-surface-ref (s x y)
  (let ((u0 (drawing-surface-u0 s))
        (v0 (drawing-surface-v0 s))
        (u1 (drawing-surface-u1 s))
        (v1 (drawing-surface-v1 s)))
    (if (and (<= u0 x) (<= x u1) (<= v0 y) (<= y v1))
        (aref (drawing-surface-pixels s)
              (+ (- x u0) (* (- v1 y) (- u1 u0 -1))))
        0.0))) ;; The Scheme for this returns +nan.0

(defun drawing-surface-set (s x y opacity)
  (let ((u0 (drawing-surface-u0 s))
        (v0 (drawing-surface-v0 s))
        (u1 (drawing-surface-u1 s))
        (v1 (drawing-surface-v1 s)))
    (when (and (<= u0 x) (<= x u1) (<= v0 y) (<= y v1))
      (setf (aref (drawing-surface-pixels s)
                  (+ (- x u0) (* (- v1 y) (- u1 u0 -1))))
            opacity))))

(defun write-transparency-mask (s)
  ;; In the Scheme, I had the program write a Portable Arbitrary Map
  ;; with both a color and a transparency map. Here, by contrast, only
  ;; the transparency map will be output. It will be in plain Portable
  ;; Gray Map format, but representing opacities rather than
  ;; whitenesses. (Thus there will be no need for gamma corrections.)
  ;; See the pgm(5) manpage for a discussion of this use of PGM
  ;; format.
  (let* ((u0 (drawing-surface-u0 s))
         (v0 (drawing-surface-v0 s))
         (u1 (drawing-surface-u1 s))
         (v1 (drawing-surface-v1 s))
         (w (- u1 u0 -1))
         (h (- v1 v0 -1))
         (|(w * h) - 1| (1- (* w h)))
         (opacities (drawing-surface-pixels s)))
    ;; "format" is not standard in Scheme, although it is widely
    ;; implemented as an extension. However, in Common Lisp it is
    ;; standardized. So let us use it.
    (format t "P2~%")
    (format t "# transparency map~%")
    (format t "~a ~a~%" w h)
    (format t "255~%")
    (loop for i from 0 to |(w * h) - 1|
          do (let* ((opacity (aref opacities i))
                    (byteval (round (* 255 opacity))))
               ;; Using "plain" PGM format avoids the issue of how to
               ;; write raw bytes. OTOH it makes the output file large
               ;; and slow to work with. (In the R7RS Scheme,
               ;; "bytevectors" provided an obvious way to write
               ;; bytes.)
               (princ byteval)
               (terpri)))))

;;;-------------------------------------------------------------------

(defun ipart (x) (floor x))
(defun iround (x) (ipart (+ x 0.5)))
(defun fpart (x) (nth-value 1 (floor x)))
(defun rfpart (x) (- 1.0 (fpart x)))

(defun plot-shallow (s x y opacity)
  (let ((combined-opacity
          (+ opacity (drawing-surface-ref s x y))))
    (drawing-surface-set s x y (min combined-opacity 1.0))))

(defun plot-steep (s x y opacity)
  (plot-shallow s y x opacity))

(defun drawln% (s x0 y0 x1 y1 plot)
  (let* ((dx (- x1 x0))
         (dy (- y1 y0))
         (gradient (if (zerop dx) 1.0 (/ dy dx)))

         ;; Handle the first endpoint.
         (xend (iround x0))
         (yend (+ y0 (* gradient (- xend x0))))
         (xgap (rfpart (+ x0 0.5)))
         (xpxl1 xend)
         (ypxl1 (ipart yend))
         (_ (funcall plot s xpxl1 ypxl1 (* (rfpart yend) xgap)))
         (_ (funcall plot s xpxl1 (1+ ypxl1) (* (fpart yend) xgap)))

         (first-y-intersection (+ yend gradient))

         ;; Handle the second endpoint.
         (xend (iround x1))
         (yend (+ y1 (* gradient (- xend x1))))
         (xgap (fpart (+ x1 0.5)))
         (xpxl2 xend)
         (ypxl2 (ipart yend))
         (_ (funcall plot s xpxl2 ypxl2 (* (rfpart yend) xgap)))
         (_ (funcall plot s xpxl2 (1+ ypxl2) (* (fpart yend) xgap))))

    ;; Loop over the rest of the points.
    (do ((x (+ xpxl1 1) (1+ x))
         (intery first-y-intersection (+ intery gradient)))
        ((= x xpxl2))
      (funcall plot s x (ipart intery) (rfpart intery))
      (funcall plot s x (1+ (ipart intery)) (fpart intery)))))

(defun draw-line (s x0 y0 x1 y1)
  (let ((x0 (coerce x0 'single-float))
        (y0 (coerce y0 'single-float))
        (x1 (coerce x1 'single-float))
        (y1 (coerce y1 'single-float)))
    (let ((xdiff (abs (- x1 x0)))
          (ydiff (abs (- y1 y0))))
      (if (<= ydiff xdiff)
          (if (<= x0 x1)
              (drawln% s x0 y0 x1 y1 #'plot-shallow)
              (drawln% s x1 y1 x0 y0 #'plot-shallow))
          (if (<= y0 y1)
              (drawln% s y0 x0 y1 x1 #'plot-steep)
              (drawln% s y1 x1 y0 x0 #'plot-steep))))))

;;;-------------------------------------------------------------------
;;; Draw a catenary as the evolute of a tractrix. See
;;; https://en.wikipedia.org/w/index.php?title=Tractrix&oldid=1143719802#Properties
;;; See also https://archive.is/YfgXW

(defvar u0 -399)
(defvar v0 -199)
(defvar u1 400)
(defvar v1 600)

(defvar s (make-drawing-surface u0 v0 u1 v1))


(loop for i from -300 to 300 by 10
      for t_ = (/ i 100.0)          ; Independent parameter.
      for x = (- t_  (tanh t_))     ; Parametric tractrix coordinates.
      for y = (/ (cosh t_))         ;
      for u = y                     ; Parametric normal vector.
      for v = (* y (sinh t_))       ;
      for x0 = (* 100.0 (- x (* 10.0 u))) ; Scaled for plotting.
      for y0 = (* 100.0 (- y (* 10.0 v)))
      for x1 = (* 100.0 (+ x (* 10.0 u)))
      for y1 = (* 100.0 (+ y (* 10.0 v)))
      do (draw-line s x0 y0 x1 y1))

(write-transparency-mask s)

;;;-------------------------------------------------------------------

Here is a script to run the program and produce a PNG image.

#!/bin/sh

sbcl --script xiaolin_wu_line_algorithm.lisp > alpha.pgm
pamgradient black black darkblue darkblue 800 800 > bluegradient.pam
pamgradient red red magenta magenta 800 800 > redgradient.pam
pamcomp -alpha=alpha.pgm redgradient.pam bluegradient.pam | pamtopng > image.png
Output:
A catenary as the evolute of a tractrix. Reddish gradient on bluish gradient.

D

Translation of: Go

This performs the mixing of the colors, both in grey scale and RGB.

import std.math, std.algorithm, grayscale_image;

/// Plots anti-aliased line by Xiaolin Wu's line algorithm.
void aaLine(Color)(ref Image!Color img,
                   double x1, double y1,
                   double x2, double y2,
                   in Color color) pure nothrow @safe @nogc {
    // Straight translation of Wikipedia pseudocode.

    // std.math.round is not pure. **
    static double round(in double x) pure nothrow @safe @nogc {
        return floor(x + 0.5);
    }

    static double fpart(in double x) pure nothrow @safe @nogc {
        return x - x.floor;
    }

    static double rfpart(in double x) pure nothrow @safe @nogc {
        return 1 - fpart(x);
    }

    auto dx = x2 - x1;
    auto dy = y2 - y1;
    immutable ax = dx.abs;
    immutable ay = dy.abs;

    static Color mixColors(in Color c1, in Color c2, in double p)
    pure nothrow @safe @nogc {
        static if (is(Color == RGB))
            return Color(cast(ubyte)(c1.r * p + c2.r * (1 - p)),
                         cast(ubyte)(c1.g * p + c2.g * (1 - p)),
                         cast(ubyte)(c1.b * p + c2.b * (1 - p)));
        else
            // This doesn't work for every kind of Color.
            return Color(cast(ubyte)(c1 * p + c2 * (1 - p)));
    }

    // Plot function set here to handle the two cases of slope.
    void function(ref Image!Color, in int, in int, in double, in Color)
    pure nothrow @safe @nogc plot;

    if (ax < ay) {
        swap(x1, y1);
        swap(x2, y2);
        swap(dx, dy);
        //plot = (img, x, y, p, col) {
        plot = (ref img, x, y, p, col) {
            assert(p >= 0.0 && p <= 1.0);
            img[y, x] = mixColors(col, img[y, x], p);
        };
    } else {
        //plot = (img, x, y, p, col) {
        plot = (ref img, x, y, p, col) {
            assert(p >= 0.0 && p <= 1.0);
            img[x, y] = mixColors(col, img[x, y], p);
        };
    }

    if (x2 < x1) {
        swap(x1, x2);
        swap(y1, y2);
    }
    immutable gradient = dy / dx;

    // Handle first endpoint.
    auto xEnd = round(x1);
    auto yEnd = y1 + gradient * (xEnd - x1);
    auto xGap = rfpart(x1 + 0.5);
    // This will be used in the main loop.
    immutable xpxl1 = cast(int)xEnd;
    immutable ypxl1 = cast(int)yEnd.floor;
    plot(img, xpxl1, ypxl1, rfpart(yEnd) * xGap, color);
    plot(img, xpxl1, ypxl1 + 1, fpart(yEnd) * xGap, color);
    // First y-intersection for the main loop.
    auto yInter = yEnd + gradient;

    // Handle second endpoint.
    xEnd = round(x2);
    yEnd = y2 + gradient * (xEnd - x2);
    xGap = fpart(x2 + 0.5);
    // This will be used in the main loop.
    immutable xpxl2 = cast(int)xEnd;
    immutable ypxl2 = cast(int)yEnd.floor;
    plot(img, xpxl2, ypxl2, rfpart(yEnd) * xGap, color);
    plot(img, xpxl2, ypxl2 + 1, fpart(yEnd) * xGap, color);

    // Main loop.
    foreach (immutable x; xpxl1 + 1 .. xpxl2) {
        plot(img, x, cast(int)yInter.floor, rfpart(yInter), color);
        plot(img, x, cast(int)yInter.floor + 1, fpart(yInter), color);
        yInter += gradient;
    }
}

void main() {
    auto im1 = new Image!Gray(400, 300);
    im1.clear(Gray.white);
    im1.aaLine(7.4, 12.3, 307, 122.5, Gray.black);
    im1.aaLine(177.4, 12.3, 127, 222.5, Gray.black);
    im1.savePGM("xiaolin_lines1.pgm");

    auto im2 = new Image!RGB(400, 300);
    im2.clear(RGB(0, 255, 0));
    immutable red = RGB(255, 0, 0);
    im2.aaLine(7.4, 12.3, 307, 122.5, red);
    im2.aaLine(177.4, 12.3, 127, 222.5, red);
    im2.savePPM6("xiaolin_lines2.ppm");
}

Fortran

Translation of: Common Lisp

The program outputs a Portable Gray Map representing a transparency mask. The mask is full of straight lines, drawn by a solution of this Rosetta Code task. Some of the lines draw a piecewise approximation of an ellipse, and others draw the normals of the ellipse. The normals form an envelope that is the evolute of the ellipse.

program xiaolin_wu_line_algorithm
  use, intrinsic :: ieee_arithmetic
  implicit none

  type :: drawing_surface
     integer :: u0, v0, u1, v1
     real, allocatable :: pixels(:)
  end type drawing_surface

  interface
     subroutine point_plotter (s, x, y, opacity)
       import drawing_surface
       type(drawing_surface), intent(inout) :: s
       integer, intent(in) :: x, y
       real, intent(in) :: opacity
     end subroutine point_plotter
  end interface

  real, parameter :: pi = 4.0 * atan (1.0)

  integer, parameter :: u0 = -499
  integer, parameter :: v0 = -374
  integer, parameter :: u1 = 500
  integer, parameter :: v1 = 375

  real, parameter :: a = 200.0  ! Ellipse radius on x axis.
  real, parameter :: b = 350.0  ! Ellipse radius on y axis.

  type(drawing_surface) :: s
  integer :: i, step_size
  real :: t, c, d
  real :: x, y
  real :: xnext, ynext
  real :: u, v
  real :: rhs, ad, bc
  real :: x0, y0, x1, y1

  s = make_drawing_surface (u0, v0, u1, v1)

  ! Draw both an ellipse and the normals of that ellipse.
  step_size = 2
  do i = 0, 359, step_size
     ! Parametric representation of an ellipse.
     t = i * (pi / 180.0)
     c = cos (t)
     d = sin (t)
     x = a * c
     y = b * d

     ! Draw a piecewise linear approximation of the ellipse. The
     ! endpoints of the line segments will lie on the curve.
     xnext = a * cos ((i + step_size) * (pi / 180.0))
     ynext = b * sin ((i + step_size) * (pi / 180.0))
     call draw_line (s, x, y, xnext, ynext)

     ! The parametric equation of the normal:
     !
     !   (a * sin (t) * xnormal) - (b * cos (t) * ynormal)
     !        = (a**2 - b**2) * cos (t) * sin (t)
     !
     ! That is:
     !
     !   (a * d * xnormal) - (b * c * ynormal) = (a**2 - b**2) * c * d
     !
     rhs = (a**2 - b**2) * c * d
     ad = a * d
     bc = b * c
     if (abs (ad) < abs (bc)) then
        x0 = -1000.0
        y0 = ((ad * x0) - rhs) / bc
        x1 = 1000.0
        y1 = ((ad * x1) - rhs) / bc
     else
        y0 = -1000.0
        x0 = (rhs - (bc * y0)) / ad
        y1 = 1000.0
        x1 = (rhs - (bc * y1)) / ad
     end if

     call draw_line (s, x0, y0, x1, y1)
  end do

  call write_transparency_mask (s)

contains

  function make_drawing_surface (u0, v0, u1, v1) result (s)
    integer, intent(in) :: u0, v0, u1, v1
    type(drawing_surface) :: s

    integer :: w, h

    if (u1 < u0 .or. v1 < v0) error stop
    s%u0 = u0; s%v0 = v0
    s%u1 = u1; s%v1 = v1
    w = u1 - u0 + 1
    h = v1 - v0 + 1
    allocate (s%pixels(0:(w * h) - 1), source = 0.0)
  end function make_drawing_surface

  function drawing_surface_ref (s, x, y) result (c)
    type(drawing_surface), intent(in) :: s
    integer, intent(in) :: x, y
    real :: c

    ! In calls to drawing_surface_ref and drawing_surface_set, indices
    ! outside the drawing_surface are allowed. Such indices are
    ! treated as if you were trying to draw on the air.
    
    if (s%u0 <= x .and. x <= s%u1 .and. s%v0 <= y .and. y <= s%v1) then
       c = s%pixels((x - s%u0) + ((s%v1 - y) * (s%u1 - s%u0 + 1)))
    else
       c = ieee_value (s%pixels(0), ieee_quiet_nan) 
    end if
  end function drawing_surface_ref

  subroutine drawing_surface_set (s, x, y, c)
    type(drawing_surface), intent(inout) :: s
    integer, intent(in) :: x, y
    real, intent(in) :: c

    ! In calls to drawing_surface_ref and drawing_surface_set, indices
    ! outside the drawing_surface are allowed. Such indices are
    ! treated as if you were trying to draw on the air.

    if (s%u0 <= x .and. x <= s%u1 .and. s%v0 <= y .and. y <= s%v1) then
       s%pixels((x - s%u0) + ((s%v1 - y) * (s%u1 - s%u0 + 1))) = c
    end if
  end subroutine drawing_surface_set

  subroutine write_transparency_mask (s)
    type(drawing_surface), intent(in) :: s

    ! Write a transparency mask, in plain (ASCII or EBCDIC) Portable
    ! Gray Map format, but representing opacities rather than
    ! whitenesses. (Thus there will be no need for gamma corrections.)
    ! See the pgm(5) manpage for a discussion of this use of PGM
    ! format.

    integer :: w, h
    integer :: i

    w = s%u1 - s%u0 + 1
    h = s%v1 - s%v0 + 1

    write (*, '("P2")')
    write (*, '("# transparency mask")')
    write (*, '(I0, 1X, I0)') w, h
    write (*, '("255")')
    write (*, '(15I4)') (nint (255 * s%pixels(i)), i = 0, (w * h) - 1)
  end subroutine write_transparency_mask

  subroutine draw_line (s, x0, y0, x1, y1)
    type(drawing_surface), intent(inout) :: s
    real, intent(in) :: x0, y0, x1, y1

    real :: xdiff, ydiff

    xdiff = abs (x1 - x0)
    ydiff = abs (y1 - y0)
    if (ydiff <= xdiff) then
       if (x0 <= x1) then
          call drawln (s, x0, y0, x1, y1, plot_shallow)
       else
          call drawln (s, x1, y1, x0, y0, plot_shallow)
       end if
    else
       if (y0 <= y1) then
          call drawln (s, y0, x0, y1, x1, plot_steep)
       else
          call drawln (s, y1, x1, y0, x0, plot_steep)
       end if
    end if
  end subroutine draw_line

  subroutine drawln (s, x0, y0, x1, y1, plot)
    type(drawing_surface), intent(inout) :: s
    real, intent(in) :: x0, y0, x1, y1
    procedure(point_plotter) :: plot

    real :: dx, dy, gradient
    real :: yend, xgap
    real :: first_y_intersection, intery
    integer :: xend
    integer :: xpxl1, ypxl1
    integer :: xpxl2, ypxl2
    integer :: x

    dx = x1 - x0;  dy = y1 - y0
    if (dx == 0.0) then
       gradient = 1.0
    else
       gradient = dy / dx
    end if

    ! Handle the first endpoint.
    xend = iround (x0)
    yend = y0 + (gradient * (xend - x0))
    xgap = rfpart (x0 + 0.5)
    xpxl1 = xend
    ypxl1 = ipart (yend)
    call plot (s, xpxl1, ypxl1, rfpart (yend) * xgap)
    call plot (s, xpxl1, ypxl1 + 1, fpart (yend) * xgap)

    first_y_intersection = yend + gradient

    ! Handle the second endpoint.
    xend = iround (x1)
    yend = y1 + (gradient * (xend - x1))
    xgap = fpart (x1 + 0.5)
    xpxl2 = xend
    ypxl2 = ipart (yend)
    call plot (s, xpxl2, ypxl2, (rfpart (yend) * xgap))
    call plot (s, xpxl2, ypxl2 + 1, fpart (yend) * xgap)

    ! Loop over the rest of the points.
    intery = first_y_intersection
    do x = xpxl1 + 1, xpxl2 - 1
       call plot (s, x, ipart (intery), rfpart (intery))
       call plot (s, x, ipart (intery) + 1, fpart (intery))
       intery = intery + gradient
    end do
  end subroutine drawln

  subroutine plot_shallow (s, x, y, opacity)
    type(drawing_surface), intent(inout) :: s
    integer, intent(in) :: x, y
    real, intent(in) :: opacity

    real :: combined_opacity

    ! Let us simply add opacities, up to the maximum of 1.0. You might
    ! wish to do something different, of course.
    combined_opacity = opacity + drawing_surface_ref (s, x, y)
    call drawing_surface_set (s, x, y, min (combined_opacity, 1.0))
  end subroutine plot_shallow

  subroutine plot_steep (s, x, y, opacity)
    type(drawing_surface), intent(inout) :: s
    integer, intent(in) :: x, y
    real, intent(in) :: opacity
    call plot_shallow (s, y, x, opacity)
  end subroutine plot_steep

  elemental function ipart (x) result (i)
    real, intent(in) :: x
    integer :: i
    i = floor (x)
  end function ipart

  elemental function iround (x) result (i)
    real, intent(in) :: x
    integer :: i
    i = ipart (x + 0.5)
  end function iround

  elemental function fpart (x) result (y)
    real, intent(in) :: x
    real :: y
    y = modulo (x, 1.0)
  end function fpart

  elemental function rfpart (x) result (y)
    real, intent(in) :: x
    real :: y
    y = 1.0 - fpart (x)
  end function rfpart

end program xiaolin_wu_line_algorithm

Here is a shell script that runs the program and creates a PNG. The background of the image is one pattern, and the foreground is another. The foreground, however, is masked by the transparency mask, and so only all those straight lines we drew show up in the PNG.

#!/bin/sh

# Using the optimizer, even at low settings, avoids trampolines and
# executable stacks.
gfortran -std=f2018 -g -O1 xiaolin_wu_line_algorithm.f90

./a.out > alpha.pgm
ppmpat -anticamo -randomseed=36 1000 750 | pambrighten -value=-60 -saturation=50 > fg.pam
ppmpat -poles -randomseed=57 1000 750 | pambrighten -value=+200 -saturation=-80 > bg.pam
pamcomp -alpha=alpha.pgm fg.pam bg.pam | pamtopng > image.png
Output:
An ellipse and its normals (forming an evolute as their envelope), all of it in goofy colors.

FreeBASIC

This implementation follows the pseudocode given on Wikipedia. Only changed xend=round() in xend=ipart() to make it more in line with FreeBASIC's own line drawing routine. Rfpart give me some trouble so I changed if somewhat. The small functions where all converted into macro's

' version 21-06-2015
' compile with: fbc -s console or fbc -s gui
' Xiaolin Wu’s line-drawing algorithm
'shared var and macro's

Dim Shared As UInteger wu_color

#Macro ipart(x)
Int(x)             ' integer part
#EndMacro

#Macro round(x)
Int((x) + .5)      ' round off
#EndMacro

#Macro fpart(x)
Frac(x)            ' fractional part
#EndMacro

#Macro rfpart(x)
' 1 - Frac(x)    ' seems to give problems for very small x
IIf(1 - Frac(x) >= 1, 1, 1 - Frac(x))
#EndMacro

#Macro plot(x, y , c)
' use the alpha channel to set the amount of color
PSet(x,y), wu_color Or (Int(c * 255)) Shl 24
#EndMacro

Sub drawline(x0 As Single, y0 As Single, x1 As Single, y1 As Single,_
    col As UInteger = RGB(255,255,255))

    wu_color = col And &HFFFFFF ' strip off the alpha channel information

    Dim As Single gradient
    Dim As Single xend, yend, xgap, intery
    Dim As UInteger xpxl1, ypxl1, xpxl2, ypxl2, x
    Dim As Integer steep = Abs(y1 - y0) > Abs(x1 - x0) ' boolean

    If steep Then
        Swap x0, y0
        Swap x1, y1
    End If

    If x0 > x1 Then
        Swap x0, x1
        Swap y0, y1
    End If

    gradient = (y1 - y0) / (x1 - x0)

    ' first endpoint
    ' xend = round(x0)
    xend = ipart(x0)
    yend = y0 + gradient * (xend - x0)
    xgap = rfpart(x0 + .5)
    xpxl1 = xend              ' this will be used in the main loop
    ypxl1 = ipart(yend)
    If steep Then
        plot(ypxl1,   xpxl1, rfpart(yend) * xgap)
        plot(ypxl1+1, xpxl1,  fpart(yend) * xgap)
    Else
        plot(xpxl1, ypxl1,   rfpart(yend) * xgap)
        plot(xpxl1, ypxl1+1,  fpart(yend) * xgap)
    End If
    intery = yend + gradient  ' first y-intersecction for the main loop

    ' handle second endpoint
    ' xend = round(x1)
    xend = ipart(x1)
    yend = y1 + gradient * (xend - x1)
    xgap = fpart(x1 + .5)
    xpxl2 = xend              ' this will be used in the main loop
    ypxl2 = ipart(yend)
    If steep Then
        plot(ypxl2,   xpxl2, rfpart(yend) * xgap)
        plot(ypxl2+1, xpxl2,  fpart(yend) * xgap)
    Else
        plot(xpxl2, ypxl2,   rfpart(yend) * xgap)
        plot(xpxl2, ypxl2+1,  fpart(yend) * xgap)
    End If

    ' main loop
    If steep Then
        For x = xpxl1 + 1 To xpxl2 - 1
            plot(ipart(intery),   x, rfpart(intery))
            plot(ipart(intery)+1, x,  fpart(intery))
            intery = intery + gradient
        Next
    Else
        For x = xpxl1 + 1 To xpxl2 - 1
            plot(x, ipart(intery),   rfpart(intery))
            plot(x, ipart(intery)+1,  fpart(intery))
            intery = intery + gradient
        Next
    End If

End Sub

' ------=< MAIN >=------

#Define W_  600
#Define H_  600

#Include Once "fbgfx.bi"   ' needed setting the screen attributes
Dim As Integer i
Dim As String fname = __FILE__

ScreenRes W_, H_, 32,, FB.GFX_ALPHA_PRIMITIVES

Randomize Timer

For i = 0 To H_ Step H_\30
    drawline(0, 0, W_, i, Int(Rnd * &HFFFFFF))
Next

For i = 0 To W_ Step W_\30
    drawline(0, 0, i, H_, Int(Rnd * &HFFFFFF))
Next

i = InStr(fname,".bas")
fname = Left(fname, Len(fname)-i+1)
WindowTitle fname + "    hit any key to end program"

While Inkey <> "" : Wend
Sleep
End

Go

package raster

import "math"

func ipart(x float64) float64 {
    return math.Floor(x)
}

func round(x float64) float64 {
    return ipart(x + .5)
}

func fpart(x float64) float64 {
    return x - ipart(x)
}

func rfpart(x float64) float64 {
    return 1 - fpart(x)
}

// AaLine plots anti-aliased line by Xiaolin Wu's line algorithm.
func (g *Grmap) AaLine(x1, y1, x2, y2 float64) {
    // straight translation of WP pseudocode
    dx := x2 - x1
    dy := y2 - y1
    ax := dx
    if ax < 0 {
        ax = -ax
    }
    ay := dy
    if ay < 0 {
        ay = -ay
    }
    // plot function set here to handle the two cases of slope
    var plot func(int, int, float64)
    if ax < ay {
        x1, y1 = y1, x1
        x2, y2 = y2, x2
        dx, dy = dy, dx
        plot = func(x, y int, c float64) {
            g.SetPx(y, x, uint16(c*math.MaxUint16))
        }
    } else {
        plot = func(x, y int, c float64) {
            g.SetPx(x, y, uint16(c*math.MaxUint16))
        }
    }
    if x2 < x1 {
        x1, x2 = x2, x1
        y1, y2 = y2, y1
    }
    gradient := dy / dx

    // handle first endpoint
    xend := round(x1)
    yend := y1 + gradient*(xend-x1)
    xgap := rfpart(x1 + .5)
    xpxl1 := int(xend) // this will be used in the main loop
    ypxl1 := int(ipart(yend))
    plot(xpxl1, ypxl1, rfpart(yend)*xgap)
    plot(xpxl1, ypxl1+1, fpart(yend)*xgap)
    intery := yend + gradient // first y-intersection for the main loop

    // handle second endpoint
    xend = round(x2)
    yend = y2 + gradient*(xend-x2)
    xgap = fpart(x2 + 0.5)
    xpxl2 := int(xend) // this will be used in the main loop
    ypxl2 := int(ipart(yend))
    plot(xpxl2, ypxl2, rfpart(yend)*xgap)
    plot(xpxl2, ypxl2+1, fpart(yend)*xgap)

    // main loop
    for x := xpxl1 + 1; x <= xpxl2-1; x++ {
        plot(x, int(ipart(intery)), rfpart(intery))
        plot(x, int(ipart(intery))+1, fpart(intery))
        intery = intery + gradient
    }
}

Demonstration program:

package main

// Files required to build supporting package raster are found in:
// * This task (immediately above)
// * Bitmap
// * Grayscale image
// * Write a PPM file

import "raster"

func main() {
    g := raster.NewGrmap(400, 300)
    g.AaLine(7.4, 12.3, 307, 122.5)
    g.AaLine(177.4, 12.3, 127, 222.5)
    g.Bitmap().WritePpmFile("wu.ppm")
}

Haskell

Example makes use of JuicyPixels for serialization to PNG format and and primitive to abstract away memory-related operations. This is a fairly close translation of the algorithm as described on Wikipedia:

{-# LANGUAGE ScopedTypeVariables #-}

module Main (main) where

import Codec.Picture (writePng)
import Codec.Picture.Types (Image, MutableImage(..), Pixel, PixelRGB8(..), createMutableImage, unsafeFreezeImage, writePixel)
import Control.Monad (void)
import Control.Monad.Primitive (PrimMonad, PrimState)
import Data.Foldable (foldlM)

type MImage m px = MutableImage (PrimState m) px

-- | Create an image given a function to apply to an empty mutable image
withMutableImage
    :: (Pixel px, PrimMonad m)
    => Int                      -- ^ image width
    -> Int                      -- ^ image height
    -> px                       -- ^ background colour
    -> (MImage m px -> m ())    -- ^ function to apply to mutable image
    -> m (Image px)             -- ^ action
withMutableImage w h px f = createMutableImage w h px >>= \m -> f m >> unsafeFreezeImage m

-- | Plot a pixel at the given point in the given colour
plot
    :: (Pixel px, PrimMonad m)
    => MImage m px  -- ^ mutable image
    -> Int          -- ^ x-coordinate of point
    -> Int          -- ^ y-coordinate of point
    -> px           -- ^ colour
    -> m ()         -- ^ action
plot = writePixel

-- | Draw an antialiased line from first point to second point in given colour
drawAntialiasedLine
    :: forall px m . (Pixel px, PrimMonad m)
    => MImage m px      -- ^ mutable image
    -> Int              -- ^ x-coordinate of first point
    -> Int              -- ^ y-coordinate of first point
    -> Int              -- ^ x-coordinate of second point
    -> Int              -- ^ y-coordinate of second point
    -> (Double -> px)   -- ^ colour generator function
    -> m ()             -- ^ action
drawAntialiasedLine m p1x p1y p2x p2y colour = do
    let steep = abs (p2y - p1y) > abs (p2x - p1x)
        ((p3x, p4x), (p3y, p4y)) = swapIf steep ((p1x, p2x), (p1y, p2y))
        ((ax, ay), (bx, by)) = swapIf (p3x > p4x) ((p3x, p3y), (p4x, p4y))
        dx = bx - ax
        dy = by - ay
        gradient = if dx == 0 then 1.0 else fromIntegral dy / fromIntegral dx

    -- handle first endpoint
    let xpxl1 = ax -- round (fromIntegral ax)
        yend1 = fromIntegral ay + gradient * fromIntegral (xpxl1 - ax)
        xgap1 = rfpart (fromIntegral ax + 0.5)
    endpoint steep xpxl1 yend1 xgap1

    -- handle second endpoint
    let xpxl2 = bx -- round (fromIntegral bx)
        yend2 = fromIntegral by + gradient * fromIntegral (xpxl2 - bx)
        xgap2 = fpart (fromIntegral bx + 0.5)
    endpoint steep xpxl2 yend2 xgap2

    -- main loop
    let intery = yend1 + gradient
    void $ if steep
        then foldlM (\i x -> do
            plot m (ipart i) x (colour (rfpart i))
            plot m (ipart i + 1) x (colour (fpart i))
            pure $ i + gradient) intery [xpxl1 + 1..xpxl2 - 1]
        else foldlM (\i x -> do
            plot m x (ipart i) (colour (rfpart i))
            plot m x (ipart i + 1) (colour (fpart i))
            pure $ i + gradient) intery [xpxl1 + 1..xpxl2 - 1]

    where
        endpoint :: Bool -> Int -> Double -> Double -> m ()
        endpoint True xpxl yend xgap = do
            plot m ypxl xpxl (colour (rfpart yend * xgap))
            plot m (ypxl + 1) xpxl (colour (fpart yend * xgap))
            where ypxl = ipart yend
        endpoint False xpxl yend xgap = do
            plot m xpxl ypxl (colour (rfpart yend * xgap))
            plot m xpxl (ypxl + 1) (colour (fpart yend * xgap))
            where ypxl = ipart yend

swapIf :: Bool -> (a, a) -> (a, a)
swapIf False p = p
swapIf True (x, y) = (y, x)

ipart :: Double -> Int
ipart = truncate

fpart :: Double -> Double
fpart x
    | x > 0 = x - temp
    | otherwise = x - (temp + 1)
    where temp = fromIntegral (ipart x)

rfpart :: Double -> Double
rfpart x = 1 - fpart x

main :: IO ()
main = do
    -- We start and end the line with sufficient clearance from the edge of the
    -- image to be able to see the endpoints
    img <- withMutableImage 640 480 (PixelRGB8 0 0 0) $ \m@(MutableImage w h _) ->
            drawAntialiasedLine m 2 2 (w - 2) (h - 2)
            (\brightness -> let level = round (brightness * 255) in PixelRGB8 level level level)

    -- Write it out to a file on disc
    writePng "xiaolin-wu-algorithm.png" img

Building and running this program will generate an output PNG file named xiaolin-wu-algorithm.png showing a white antialiased diagonal line.

Icon

Translation of: ObjectIcon

Please be aware that the program below is written for classical "University of Arizona" Icon, and not for Unicon or Object Icon. It uses a graphics system designed for machines of past times.

I have taken the Object Icon program and made the minimum number of changes needed to get it running as an Arizona Icon program. As with the Object Icon, a window comes up and then you can draw lines on it by pressing the left mouse button. Pressing "q" or "Q" will quit the program. Instead of a PNG, the program writes a GIF. Instead of varying the opacity of a line's pixels, the program varies the shade of gray.

link graphics

$define YES 1
$define NO &null

procedure main ()
  local width, height
  local done, w, event
  local x1, y1, x2, y2, press_is_active

  width := 640
  height := 480

  w := WOpen ("size=" || width || "," || height, "bg=white") |
    stop ("failed to open a window")

  press_is_active := NO
  done := NO
  while /done do
  {
    if *(Pending (w)) ~= 0 then
    {
      event := Event (w)
      case event of
      {
        QuitEvents () : done := YES
        &lpress :
        {
          if /press_is_active then
          {
            x1 := &x;  y1 := &y
            press_is_active := YES
          }
          else
          {
            x2 := &x;  y2 := &y
            draw_line (w, x1, y1, x2, y2)
            press_is_active := NO
          }
        }
      }
    }
  }

  # GIF is the only format "regular" Icon is documented to be able to
  # write on all platforms. (Icon used to run on, for instance, 16-bit
  # MSDOS boxes.)
  WriteImage (w, "xiaolin_wu_line_algorithm_Arizona.gif") |
    stop ("failed to write a GIF")
end

procedure draw_line (w, x0, y0, x1, y1)
  local steep
  local dx, dy, gradient
  local xend, yend, xgap, intery
  local xpxl1, ypxl1
  local xpxl2, ypxl2
  local x

  x0 := real (x0)
  y0 := real (y0)
  x1 := real (x1)
  y1 := real (y1)

  # In Icon, comparisons DO NOT return boolean values! They either
  # SUCCEED or they FAIL. Thus the need for an "if-then-else" here.
  steep := (if abs (y1 - y0) > abs (x1 - x0) then YES else NO)

  if \steep then { x0 :=: y0;  x1 :=: y1 }
  if x0 > x1 then { x0 :=: x1;  y0 :=: y1 }
  dx := x1 - x0;  dy := y1 - y0
  gradient := (if dx = 0 then 1.0 else dy / dx)

  # Handle the first endpoint.
  xend := round (x0);  yend := y0 + (gradient * (xend - x0))
  xgap := rfpart (x0 + 0.5)
  xpxl1 := xend;  ypxl1 := ipart (yend)
  if \steep then
  {
    plot (w, ypxl1, xpxl1, rfpart (yend) * xgap)
    plot (w, ypxl1 + 1, xpxl1, fpart(yend) * xgap)
  }
  else
  {
    plot (w, xpxl1, ypxl1, rfpart (yend) * xgap)
    plot (w, xpxl1, ypxl1 + 1, fpart (yend) * xgap)
  }

  # The first y-intersection.
  intery := yend + gradient

  # Handle the second endpoint.
  xend := round (x1);  yend := y1 + (gradient * (xend - x1))
  xgap := fpart (x1 + 0.5)
  xpxl2 := xend;  ypxl2 := ipart (yend)
  if \steep then
  {
    plot (w, ypxl2, xpxl2, rfpart (yend) * xgap)
    plot (w, ypxl2 + 1, xpxl2, fpart (yend) * xgap)
  }
  else
  {
    plot (w, xpxl2, ypxl2, rfpart (yend) * xgap)
    plot (w, xpxl2, ypxl2 + 1, fpart (yend) * xgap)
  }

  if \steep then
    every x := xpxl1 + 1 to xpxl2 - 1 do
    {
      plot (w, ipart (intery), x, rfpart (intery))
      plot (w, ipart (intery) + 1, x, fpart (intery))
      intery := intery + gradient
    }
  else
    every x := xpxl1 + 1 to xpxl2 - 1 do
    {
      plot(w, x, ipart (intery), rfpart (intery))
      plot(w, x, ipart (intery) + 1, fpart (intery))
      intery := intery + gradient
    }

  return
end

procedure plot (w, x, y, c)
  # In Object Icon, one can vary the opacity. But not in "regular"
  # Icon. Here we vary, instead, the shade of gray.
  c := round ((1.0 - c) * 65535)
  Fg (w, c || "," || c || "," || c)
  DrawPoint (w, x, y)
  return
end

procedure ipart (x)
  local i
  i := integer (x)
  return (if i = x then i else if x < 0 then i - 1 else i)
end

procedure round (x)
  return ipart (x + 0.5)
end

procedure fpart (x)
  return x - ipart (x)
end

procedure rfpart (x)
  return 1.0 - fpart (x)
end

# Unlike Object Icon, "regular" Icon has no "abs" built-in.
procedure abs (x)
  return (if x < 0 then -x else x)
end
Output:

An example GIF:

A GIF of the window drawn by a run of the Icon program. An assortment of anti-aliased straight lines in black on white.

J

Solution:

load'gl2'
coinsert'jgl2'

drawpt=:4 :0"0 1
   glrgb <.(-.x)*255 255 255
   glpixel y
)

drawLine=:3 :0 NB. drawline x1,y1,x2,y2
   pts=. 2 2$y
   isreversed=. </ |d=. -~/pts
   r=. |.^:isreversed"1
   pts=. /:~ pts \:"1 |d
   gradient=. %~/ (\:|)d

   'x y'=. |:pts
   xend=. <.0.5+ x
   yend=. y + gradient* xend-x
   xgap=. -.1|x+0.5

   n=. i. >: -~/ xend
   'xlist ylist'=. (n*/~1,gradient) + ({.xend),({.yend)
   weights=. ((2&}.,~ xgap*2&{.)&.(_1&|.) (,.~-.) 1|ylist)
   weights (drawpt r)"1 2 (,:+&0 1)"1 xlist,.<.ylist
)

Example use:

   wd'pc win closeok; xywh 0 0 300 200;cc g isigraph; pas 0 0; pshow;' NB. J6 or earlier
   wd'pc win closeok; minwh 600 400;cc g isidraw flush; pshow;'        NB. J802 or later
   glpaint glclear ''
   glpaint drawLine 10 10 590 390

Java

Works with: Java version 8
import java.awt.*;
import static java.lang.Math.*;
import javax.swing.*;

public class XiaolinWu extends JPanel {

    public XiaolinWu() {
        Dimension dim = new Dimension(640, 640);
        setPreferredSize(dim);
        setBackground(Color.white);
    }

    void plot(Graphics2D g, double x, double y, double c) {
        g.setColor(new Color(0f, 0f, 0f, (float)c));
        g.fillOval((int) x, (int) y, 2, 2);
    }

    int ipart(double x) {
        return (int) x;
    }

    double fpart(double x) {
        return x - floor(x);
    }

    double rfpart(double x) {
        return 1.0 - fpart(x);
    }

    void drawLine(Graphics2D g, double x0, double y0, double x1, double y1) {

        boolean steep = abs(y1 - y0) > abs(x1 - x0);
        if (steep)
            drawLine(g, y0, x0, y1, x1);

        if (x0 > x1)
            drawLine(g, x1, y1, x0, y0);

        double dx = x1 - x0;
        double dy = y1 - y0;
        double gradient = dy / dx;

        // handle first endpoint
        double xend = round(x0);
        double yend = y0 + gradient * (xend - x0);
        double xgap = rfpart(x0 + 0.5);
        double xpxl1 = xend; // this will be used in the main loop
        double ypxl1 = ipart(yend);

        if (steep) {
            plot(g, ypxl1, xpxl1, rfpart(yend) * xgap);
            plot(g, ypxl1 + 1, xpxl1, fpart(yend) * xgap);
        } else {
            plot(g, xpxl1, ypxl1, rfpart(yend) * xgap);
            plot(g, xpxl1, ypxl1 + 1, fpart(yend) * xgap);
        }

        // first y-intersection for the main loop
        double intery = yend + gradient;

        // handle second endpoint
        xend = round(x1);
        yend = y1 + gradient * (xend - x1);
        xgap = fpart(x1 + 0.5);
        double xpxl2 = xend; // this will be used in the main loop
        double ypxl2 = ipart(yend);

        if (steep) {
            plot(g, ypxl2, xpxl2, rfpart(yend) * xgap);
            plot(g, ypxl2 + 1, xpxl2, fpart(yend) * xgap);
        } else {
            plot(g, xpxl2, ypxl2, rfpart(yend) * xgap);
            plot(g, xpxl2, ypxl2 + 1, fpart(yend) * xgap);
        }

        // main loop
        for (double x = xpxl1 + 1; x <= xpxl2 - 1; x++) {
            if (steep) {
                plot(g, ipart(intery), x, rfpart(intery));
                plot(g, ipart(intery) + 1, x, fpart(intery));
            } else {
                plot(g, x, ipart(intery), rfpart(intery));
                plot(g, x, ipart(intery) + 1, fpart(intery));
            }
            intery = intery + gradient;
        }
    }

    @Override
    public void paintComponent(Graphics gg) {
        super.paintComponent(gg);
        Graphics2D g = (Graphics2D) gg;

        drawLine(g, 550, 170, 50, 435);
    }

    public static void main(String[] args) {
        SwingUtilities.invokeLater(() -> {
            JFrame f = new JFrame();
            f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
            f.setTitle("Xiaolin Wu's line algorithm");
            f.setResizable(false);
            f.add(new XiaolinWu(), BorderLayout.CENTER);
            f.pack();
            f.setLocationRelativeTo(null);
            f.setVisible(true);
        });
    }
}

Julia

Works with: Julia version 0.6
using Images

fpart(x) = mod(x, one(x))
rfpart(x) = one(x) - fpart(x)

function drawline!(img::Matrix{Gray{N0f8}}, x0::Integer, y0::Integer, x1::Integer, y1::Integer)
    steep = abs(y1 - y0) > abs(x1 - x0)

    if steep
        x0, y0 = y0, x0
        x1, y1 = y1, x1
    end
    if x0 > x1
        x0, x1 = x1, x0
        y0, y1 = y1, y0
    end

    dx = x1 - x0
    dy = y1 - y0
    grad = dy / dx

    if iszero(dx)
        grad = oftype(grad, 1.0)
    end

    # handle first endpoint
    xend = round(Int, x0)
    yend = y0 + grad * (xend - x0)
    xgap = rfpart(x0 + 0.5)
    xpxl1 = xend
    ypxl1 = floor(Int, yend)

    if steep
        img[ypxl1,   xpxl1] = rfpart(yend) * xgap
        img[ypxl1+1, xpxl1] =  fpart(yend) * xgap
    else
        img[xpxl1, ypxl1  ] = rfpart(yend) * xgap
        img[xpxl1, ypxl1+1] =  fpart(yend) * xgap
    end
    intery = yend + grad # first y-intersection for the main loop

    # handle second endpoint
    xend = round(Int, x1)
    yend = y1 + grad * (xend - x1)
    xgap = fpart(x1 + 0.5)
    xpxl2 = xend
    ypxl2 = floor(Int, yend)
    if steep
        img[ypxl2,   xpxl2] = rfpart(yend) * xgap
        img[ypxl2+1, xpxl2] =  fpart(yend) * xgap
    else
        img[xpxl2, ypxl2  ] = rfpart(yend) * xgap
        img[xpxl2, ypxl2+1] =  fpart(yend) * xgap
    end

    # main loop
    if steep
        for x in xpxl1+1:xpxl2-1
            img[floor(Int, intery),   x] = rfpart(intery)
            img[floor(Int, intery)+1, x] =  fpart(intery)
            intery += grad
        end
    else
        for x in xpxl1+1:xpxl2-1
            img[x, floor(Int, intery)  ] = rfpart(intery)
            img[x, floor(Int, intery)+1] =  fpart(intery)
            intery += grad
        end
    end

    return img
end

img = fill(Gray(1.0N0f8), 250, 250);
drawline!(img, 8, 8, 192, 154)

Kotlin

Translation of: Java
// version 1.1.2

import java.awt.*
import javax.swing.*

class XiaolinWu: JPanel() {
    init {
        preferredSize = Dimension(640, 640)
        background = Color.white
    }

    private fun plot(g: Graphics2D, x: Double, y: Double, c: Double) {
        g.color = Color(0f, 0f, 0f, c.toFloat())
        g.fillOval(x.toInt(), y.toInt(), 2, 2)
    }

    private fun ipart(x: Double) = x.toInt()

    private fun fpart(x: Double) = x - Math.floor(x)

    private fun rfpart(x: Double) = 1.0 - fpart(x)

    private fun drawLine(g: Graphics2D, x0: Double, y0: Double, x1: Double, y1: Double) {
        val steep = Math.abs(y1 - y0) > Math.abs(x1 - x0)
        if (steep) drawLine(g, y0, x0, y1, x1)
        if (x0 > x1) drawLine(g, x1, y1, x0, y0)

        val dx = x1 - x0
        val dy = y1 - y0
        val gradient = dy / dx

        // handle first endpoint
        var xend = Math.round(x0).toDouble()
        var yend = y0 + gradient * (xend - x0)
        var xgap = rfpart(x0 + 0.5)
        val xpxl1 = xend  // this will be used in the main loop
        val ypxl1 = ipart(yend).toDouble()

        if (steep) {
            plot(g, ypxl1, xpxl1, rfpart(yend) * xgap)
            plot(g, ypxl1 + 1.0, xpxl1, fpart(yend) * xgap)
        }
        else {
            plot(g, xpxl1, ypxl1, rfpart(yend) * xgap)
            plot(g, xpxl1, ypxl1 + 1.0, fpart(yend) * xgap)
        }

        // first y-intersection for the main loop
        var intery = yend + gradient

        // handle second endpoint
        xend = Math.round(x1).toDouble()
        yend = y1 + gradient * (xend - x1)
        xgap = fpart(x1 + 0.5)
        val xpxl2 = xend  // this will be used in the main loop
        val ypxl2 = ipart(yend).toDouble()

        if (steep) {
            plot(g, ypxl2, xpxl2, rfpart(yend) * xgap)
            plot(g, ypxl2 + 1.0, xpxl2, fpart(yend) * xgap)
        }
        else {
            plot(g, xpxl2, ypxl2, rfpart(yend) * xgap)
            plot(g, xpxl2, ypxl2 + 1.0, fpart(yend) * xgap)
        }

        // main loop
        var x = xpxl1 + 1.0
        while (x <=  xpxl2 - 1) {
            if (steep) {
                plot(g, ipart(intery).toDouble(), x, rfpart(intery))
                plot(g, ipart(intery).toDouble() + 1.0, x, fpart(intery))
            }
            else {
                plot(g, x, ipart(intery).toDouble(), rfpart(intery))
                plot(g, x, ipart(intery).toDouble() + 1.0, fpart(intery))
            }
            intery += gradient
            x++
        }
    }

    override protected fun paintComponent(gg: Graphics) {
        super.paintComponent(gg)
        val g = gg as Graphics2D
        drawLine(g, 550.0, 170.0, 50.0, 435.0)
    }
}

fun main(args: Array<String>) {
    SwingUtilities.invokeLater {
        val f = JFrame()
        f.defaultCloseOperation = JFrame.EXIT_ON_CLOSE
        f.title = "Xiaolin Wu's line algorithm"
        f.isResizable = false
        f.add(XiaolinWu(), BorderLayout.CENTER)
        f.pack()
        f.setLocationRelativeTo(null)
        f.isVisible = true
    }
}

Liberty BASIC

NoMainWin
WindowWidth = 270
WindowHeight = 290
UpperLeftX=int((DisplayWidth-WindowWidth)/2)
UpperLeftY=int((DisplayHeight-WindowHeight)/2)

Global variablesInitialized : variablesInitialized = 0
Global BackColor$ : BackColor$ = "0 0 0"
'    BackColor$ = "255 255 255"
    'now, right click randomizes BG
Global size : size = 1'4
global mousepoints.mouseX0,  mousepoints.mouseY0, mousepoints.mouseX1, mousepoints.mouseY1

'StyleBits #main.gbox, 0, _WS_BORDER, 0, 0
GraphicBox #main.gbox, 0, 0, 253, 252

Open "Click Twice to Form Line" For Window As #main
Print #main, "TrapClose quit"
Print #main.gbox, "Down; Color Black"
Print #main.gbox, "Down; fill ";BackColor$
Print #main.gbox, "When leftButtonUp gBoxClick"
Print #main.gbox, "When rightButtonUp RandomBG"
Print #main.gbox, "Size "; size

result = drawAntiAliasedLine(126.5, 0, 126.5, 252, "255 0 0")
result = drawAntiAliasedLine(0, 126, 253, 126, "255 0 0")
result = drawAntiAliasedLine(0, 0, 253, 252, "255 0 0")
result = drawAntiAliasedLine(253, 0, 0, 252, "255 0 0")
Wait


    Sub quit handle$
        Close #main
        End
    End Sub

sub RandomBG handle$, MouseX, MouseY
    BackColor$ = int(rnd(1)*256);" ";int(rnd(1)*256);" ";int(rnd(1)*256)
    Print #main.gbox, "CLS; fill ";BackColor$
    variablesInitialized = 0
end sub

    Sub gBoxClick handle$, MouseX, MouseY
        'We will use the mousepoints "struct" to hold the values
        'that way they are retained between subroutine calls
        If variablesInitialized = 0 Then
            Print #main.gbox, "CLS; fill ";BackColor$
            mousepoints.mouseX0 = MouseX
            mousepoints.mouseY0 = MouseY
            variablesInitialized = 1
        Else
            If variablesInitialized = 1 Then
                mousepoints.mouseX1 = MouseX
                mousepoints.mouseY1 = MouseY
                variablesInitialized = 0
                result = drawAntiAliasedLine(mousepoints.mouseX0, mousepoints.mouseY0, mousepoints.mouseX1, mousepoints.mouseY1, "255 0 0")
            End If
        End If
    End Sub

    Function Swap(Byref a,Byref b)
        aTemp = b
        b = a
        a = aTemp
    End Function

    Function RoundtoInt(val)
        RoundtoInt = Int(val + 0.5)
    End Function

    Function PlotAntiAliased(x, y, RGB$, b, steep)

        RGB$ = Int(Val(Word$(BackColor$, 1))*(1-b) + Val(Word$(RGB$, 1)) * b) ; " " ; _
               Int(Val(Word$(BackColor$, 2))*(1-b) + Val(Word$(RGB$, 3)) * b) ; " " ; _
               Int(Val(Word$(BackColor$, 3))*(1-b) + Val(Word$(RGB$, 2)) * b)

        if steep then 'x and y reversed
            Print #main.gbox, "Down; Color " + RGB$ + "; Set " + str$(y) + " " + str$(x)
        else
            Print #main.gbox, "Down; Color " + RGB$ + "; Set " + str$(x) + " " + str$(y)
        end if
    End Function

    Function fracPart(x)
        fracPart = (x Mod 1)
    End function

    Function invFracPart(x)
        invFracPart = (1 - fracPart(x))
    End Function

    Function drawAntiAliasedLine(x1, y1, x2, y2, RGB$)
        If (x2 - x1)=0 Or (y2 - y1)=0 Then
            Print #main.gbox, "Down; Color " + RGB$
            result = BresenhamLine(x1, y1, x2, y2)
            Exit Function
        End If
        steep = abs(x2 - x1) < abs(y2 - y1)
        if steep then   'x and y should be reversed
            result = Swap(x1, y1)
            result = Swap(x2, y2)
        end if

        If (x2 < x1) Then
            result = Swap(x1, x2)
            result = Swap(y1, y2)
        End If
        dx = (x2 - x1)
        dy = (y2 - y1)
        grad = (dy/ dx)
        'Handle the First EndPoint
        xend = RoundtoInt(x1)
        yend = y1 + grad * (xend - x1)
        xgap = invFracPart(x1 + 0.5)
        ix1 = xend
        iy1 = Int(yend)
        result = PlotAntiAliased(ix1, iy1, RGB$, invFracPart(yend) * xgap, steep )
        result = PlotAntiAliased(ix1, (iy1 + size), RGB$, fracPart(yend) * xgap, steep )
        yf = (yend + grad)
        'Handle the Second EndPoint
        xend = RoundtoInt(x2)
        yend = y2 + grad * (xend - x2)
        xgap = fracPart(x2 + 0.5)
        ix2 = xend
        iy2 = Int(yend)
        result = PlotAntiAliased(ix2, iy2, RGB$, invFracPart(yend) * xgap, steep )
        result = PlotAntiAliased(ix2, (iy2 + size), RGB$, fracPart(yend) * xgap, steep )
        For x = ix1 + 1 To ix2 - 1
            result = PlotAntiAliased(x, Int(yf), RGB$, invFracPart(yf), steep )
            result = PlotAntiAliased(x, (Int(yf) + size), RGB$, fracPart(yf), steep )
            yf = (yf + grad)
        Next x
    End Function


    Function BresenhamLine(x0, y0, x1, y1)
        dx = Abs(x1 - x0)
        dy = Abs(y1 - y0)
        sx = ((x1 > x0) + Not(x0 < x1))
        sy = ((y1 > y0) + Not(y0 < y1))
        errornum = (dx - dy)
        Do While 1
            Print #main.gbox, "Set " + str$(x0) + " " + str$(y0)
            If (x0 = x1) And (y0 = y1) Then Exit Do
            errornum2 = (2 * errornum)
            If errornum2 > (-1 * dy) Then
                errornum = (errornum - dy)
                x0 = (x0 + sx)
            End If
            If errornum2 < dx Then
                errornum = (errornum + dx)
                y0 = (y0 + sy)
            End If
        Loop
    End Function

Mathematica /Wolfram Language

ClearAll[ReverseFractionalPart, ReplacePixelWithAlpha, DrawEndPoint, DrawLine]
ReverseFractionalPart[x_] := 1 - FractionalPart[x]
ReplacePixelWithAlpha[img_Image, pos_ -> colvals : {_, _, _}, 
  alpha_] := Module[{vals,},
  vals = PixelValue[img, pos];
  vals = (1 - alpha) vals + alpha colvals;
  ReplacePixelValue[img, pos -> vals]
  ]
DrawEndPoint[img_Image, pt : {x_, y_}, grad_, p_] := 
 Module[{xend, yend, xgap, px, py, i},
  xend = Round[x];
  yend = y + grad (xend - x);
  xgap = ReverseFractionalPart[x + 0.5];
  {px, py} = Floor[{xend, yend}];
  i = ReplacePixelWithAlpha[img, p[{x, py}] -> {1, 1, 1}, ReverseFractionalPart[yend] xgap];
  i = ReplacePixelWithAlpha[i, p[{x, py + 1}] -> {1, 1, 1}, FractionalPart[yend] xgap];
  {px, i}
  ]
DrawLine[img_Image, p1 : {_, _}, p2 : {_, _}] := 
 Module[{x1, x2, y1, y2, steep, p, grad, intery, xend, yend, x, y, 
   xstart, ystart, dx, dy, i},
  {x1, y1} = p1;
  {x2, y2} = p2;
  dx = x2 - x1;
  dy = y2 - y1;
  steep = Abs[dx] < Abs[dy];
  p = If[steep, Reverse[#], #] &;
  If[steep,
   {x1, y1, x2, y2, dx, dy} = {y1, x1, y2, x2, dy, dx}
   ];
  If[x2 < x1,
   {x1, x2, y1, y2} = {x2, x1, y2, y1}
   ];
  grad = dy/dx;
  intery = y1 + ReverseFractionalPart[x1] grad;
  {xstart, i} = DrawEndPoint[img, p[p1], grad, p];
  xstart += 1;
  {xend, i} = DrawEndPoint[i, p[p2], grad, p];
  Do[
   y = Floor[intery];
   i = ReplacePixelWithAlpha[i, p[{x, y}] -> {1, 1, 1}, ReverseFractionalPart[intery]];
   i = ReplacePixelWithAlpha[i, p[{x, y + 1}] -> {1, 1, 1}, FractionalPart[intery]];
   intery += grad
   ,
   {x, xstart, xend}
   ];
  i
  ]
image = ConstantImage[Black, {100, 100}];
Fold[DrawLine[#1, {20, 10}, #2] &, image, AngleVector[{20, 10}, {75, #}] & /@ Subdivide[0, Pi/2, 10]]

MATLAB

Translation of: Julia
clear all;close all;clc;
% Example usage:
img = ones(250, 250);
img = drawline(img, 8, 8, 192, 154);
imshow(img); % Display the image

function img = drawline(img, x0, y0, x1, y1)
    function f = fpart(x)
        f = mod(x, 1);
    end

    function rf = rfpart(x)
        rf = 1 - fpart(x);
    end

    steep = abs(y1 - y0) > abs(x1 - x0);

    if steep
        [x0, y0] = deal(y0, x0);
        [x1, y1] = deal(y1, x1);
    end
    if x0 > x1
        [x0, x1] = deal(x1, x0);
        [y0, y1] = deal(y1, y0);
    end

    dx = x1 - x0;
    dy = y1 - y0;
    grad = dy / dx;

    if dx == 0
        grad = 1.0;
    end

    % handle first endpoint
    xend = round(x0);
    yend = y0 + grad * (xend - x0);
    xgap = rfpart(x0 + 0.5);
    xpxl1 = xend;
    ypxl1 = floor(yend);

    if steep
        img(ypxl1,   xpxl1) = rfpart(yend) * xgap;
        img(ypxl1+1, xpxl1) = fpart(yend) * xgap;
    else
        img(xpxl1, ypxl1  ) = rfpart(yend) * xgap;
        img(xpxl1, ypxl1+1) = fpart(yend) * xgap;
    end
    intery = yend + grad; % first y-intersection for the main loop

    % handle second endpoint
    xend = round(x1);
    yend = y1 + grad * (xend - x1);
    xgap = fpart(x1 + 0.5);
    xpxl2 = xend;
    ypxl2 = floor(yend);
    if steep
        img(ypxl2,   xpxl2) = rfpart(yend) * xgap;
        img(ypxl2+1, xpxl2) = fpart(yend) * xgap;
    else
        img(xpxl2, ypxl2  ) = rfpart(yend) * xgap;
        img(xpxl2, ypxl2+1) = fpart(yend) * xgap;
    end

    % main loop
    if steep
        for x = (xpxl1+1):(xpxl2-1)
            img(floor(intery),   x) = rfpart(intery);
            img(floor(intery)+1, x) = fpart(intery);
            intery = intery + grad;
        end
    else
        for x = (xpxl1+1):(xpxl2-1)
            img(x, floor(intery)  ) = rfpart(intery);
            img(x, floor(intery)+1) = fpart(intery);
            intery = intery + grad;
        end
    end
end


Modula-2

Translation of: Fortran

The program outputs a transparency map in Portable Gray Map format. It draws lines normal to a parabola. The envelope formed is a semicubic parabola.

MODULE Xiaolin_Wu_Task;

(* The program is for ISO Modula-2. To compile with GNU Modula-2
   (gm2), use the "-fiso" option. *)

IMPORT RealMath;
IMPORT SRawIO;
IMPORT STextIO;
IMPORT SWholeIO;
IMPORT SYSTEM;

CONST MaxDrawingSurfaceIndex = 1999;
CONST MaxDrawingSurfaceSize =
      (MaxDrawingSurfaceIndex + 1) * (MaxDrawingSurfaceIndex + 1);

TYPE DrawingSurfaceIndex = [0 .. MaxDrawingSurfaceIndex];
TYPE PixelsIndex = [0 .. MaxDrawingSurfaceSize - 1];
TYPE DrawingSurface =
     RECORD
       u0, v0, u1, v1 : INTEGER;
       pixels : ARRAY PixelsIndex OF REAL;
     END;
TYPE PointPlotter = PROCEDURE (VAR DrawingSurface,
                               INTEGER, INTEGER, REAL);

PROCEDURE InitializeDrawingSurface (VAR s : DrawingSurface;
                                    u0, v0, u1, v1 : INTEGER);
  VAR i : PixelsIndex;
BEGIN
  s.u0 := u0; s.v0 := v0;
  s.u1 := u1; s.v1 := v1;
  FOR i := 0 TO MaxDrawingSurfaceSize - 1 DO
    s.pixels[i] := 0.0
  END
END InitializeDrawingSurface;

PROCEDURE DrawingSurfaceRef (VAR s : DrawingSurface;
                             x, y : DrawingSurfaceIndex) : REAL;
  VAR c : REAL;
BEGIN
  IF (s.u0 <= x) AND (x <= s.u1) AND (s.v0 <= y) AND (y <= s.v1) THEN
    c := s.pixels[(x - s.u0) + ((s.v1 - y) * (s.u1 - s.u0 + 1))]
  ELSE
    (* (x,y) is outside the drawing surface. Return a somewhat
       arbitrary value. "Not a number" would be better. *)
    c := 0.0
  END;
  RETURN c
END DrawingSurfaceRef;

PROCEDURE DrawingSurfaceSet (VAR s : DrawingSurface;
                             x, y : DrawingSurfaceIndex;
                             c : REAL);
BEGIN
  (* Store the value only if (x,y) is within the drawing surface. *)
  IF (s.u0 <= x) AND (x <= s.u1) AND (s.v0 <= y) AND (y <= s.v1) THEN
    s.pixels[(x - s.u0) + ((s.v1 - y) * (s.u1 - s.u0 + 1))] := c
  END
END DrawingSurfaceSet;

PROCEDURE WriteTransparencyMask (VAR s : DrawingSurface);
  VAR w, h : INTEGER;
      i : DrawingSurfaceIndex;
      byteval : [0 .. 255];
      byte : SYSTEM.LOC;
BEGIN
  (* Send to standard output a transparency map in raw Portable Gray
     Map format. *)
  w := s.u1 - s.u0 + 1;
  h := s.v1 - s.v0 + 1;
  STextIO.WriteString ('P5');
  STextIO.WriteLn;
  STextIO.WriteString ('# transparency mask');
  STextIO.WriteLn;
  SWholeIO.WriteCard (VAL (CARDINAL, w), 0);
  STextIO.WriteString (' ');
  SWholeIO.WriteCard (VAL (CARDINAL, h), 0);
  STextIO.WriteLn;
  STextIO.WriteString ('255');
  STextIO.WriteLn;
  FOR i := 0 TO (w * h) - 1 DO
    byteval := RealMath.round (255.0 * s.pixels[i]);
    byte := SYSTEM.CAST (SYSTEM.LOC, byteval);
    SRawIO.Write (byte)
  END
END WriteTransparencyMask;

PROCEDURE ipart (x : REAL) : INTEGER;
  VAR i : INTEGER;
BEGIN
  i := VAL (INTEGER, x);
  IF x < VAL (REAL, i) THEN
    i := i - 1;
  END;
  RETURN i
END ipart;

PROCEDURE iround (x : REAL) : INTEGER;
BEGIN
  RETURN ipart (x + 0.5)
END iround;

PROCEDURE fpart (x : REAL) : REAL;
BEGIN
  RETURN x - VAL (REAL, ipart (x))
END fpart;

PROCEDURE rfpart (x : REAL) : REAL;
BEGIN
  RETURN 1.0 - fpart (x)
END rfpart;

PROCEDURE PlotShallow (VAR s : DrawingSurface;
                       x, y : INTEGER;
                       opacity : REAL);
  VAR combined_opacity : REAL;
BEGIN
  (* Let us simply add opacities, up to the maximum of 1.0. You might,
     of course, wish to do something different. *)
  combined_opacity := opacity + DrawingSurfaceRef (s, x, y);
  IF combined_opacity > 1.0 THEN
    combined_opacity := 1.0
  END;
  DrawingSurfaceSet (s, x, y, combined_opacity)
END PlotShallow;

PROCEDURE PlotSteep (VAR s : DrawingSurface;
                     x, y : INTEGER;
                     opacity : REAL);
BEGIN
  PlotShallow (s, y, x, opacity)
END PlotSteep;

PROCEDURE drawln (VAR s : DrawingSurface;
                  x0, y0, x1, y1 : REAL;
                  plot : PointPlotter);
  VAR dx, dy, gradient : REAL;
      yend, xgap : REAL;
      first_y_intersection, intery : REAL;
      xend : INTEGER;
      xpxl1, ypxl1 : INTEGER;
      xpxl2, ypxl2 : INTEGER;
      x : INTEGER;
BEGIN
  dx := x1 - x0;  dy := y1 - y0;
  IF dx = 0.0 THEN
    gradient := 1.0
  ELSE
    gradient := dy / dx
  END;

  (* Handle the first endpoint. *)
  xend := iround (x0);
  yend := y0 + (gradient * (VAL (REAL, xend) - x0));
  xgap := rfpart (x0 + 0.5);
  xpxl1 := xend;
  ypxl1 := ipart (yend);
  plot (s, xpxl1, ypxl1, rfpart (yend) * xgap);
  plot (s, xpxl1, ypxl1 + 1, fpart (yend) * xgap);

  first_y_intersection := yend + gradient;

  (* Handle the second endpoint. *)
  xend := iround (x1);
  yend := y1 + (gradient * (VAL (REAL, xend) - x1));
  xgap := fpart (x1 + 0.5);
  xpxl2 := xend;
  ypxl2 := ipart (yend);
  plot (s, xpxl2, ypxl2, (rfpart (yend) * xgap));
  plot (s, xpxl2, ypxl2 + 1, fpart (yend) * xgap);

  (* Loop over the rest of the points. *)
  intery := first_y_intersection;
  FOR x := xpxl1 + 1 TO xpxl2 - 1 DO
    plot (s, x, ipart (intery), rfpart (intery));
    plot (s, x, ipart (intery) + 1, fpart (intery));
    intery := intery + gradient
  END
END drawln;

PROCEDURE DrawLine (VAR s : DrawingSurface;
                    x0, y0, x1, y1 : REAL);
  VAR xdiff, ydiff : REAL;
BEGIN
  xdiff := ABS (x1 - x0);
  ydiff := ABS (y1 - y0);
  IF ydiff <= xdiff THEN
    IF x0 <= x1 THEN
      drawln (s, x0, y0, x1, y1, PlotShallow)
    ELSE
      drawln (s, x1, y1, x0, y0, PlotShallow)
    END
  ELSE
    IF y0 <= y1 THEN
      drawln (s, y0, x0, y1, x1, PlotSteep)
    ELSE
      drawln (s, y1, x1, y0, x0, PlotSteep)
    END
  END
END DrawLine;

CONST u0 = -299;
      u1 = 300;
      v0 = -20;
      v1 = 379;
CONST Kx = 4.0;
      Ky = 0.1;
VAR s : DrawingSurface;
    i : INTEGER;
    t : REAL;
    x0, y0, x1, y1 : REAL;
    x, y, u, v : REAL;
BEGIN
  InitializeDrawingSurface (s, u0, v0, u1, v1);

  (* Draw a parabola. *)
  FOR i := -101 TO 100 DO
    t := VAL (REAL, i);  x0 := Kx * t;  y0 := Ky * t * t;
    t := VAL (REAL, i + 1);  x1 := Kx * t;  y1 := Ky * t * t;
    DrawLine (s, x0, y0, x1, y1)
  END;

  (* Draw normals to that parabola. The parabola has equation y=A*x*x,
     where A=Ky/(Kx*Kx). Therefore the slope at x is dy/dx=2*A*x. The
     slope of the normal is the negative reciprocal, and so equals
     -1/(2*A*x)=-(Kx*Kx)/(2*Ky*(Kx*t))=-Kx/(2*Ky*t). *)
  FOR i := -101 TO 101 DO
    t := VAL (REAL, i);
    x := Kx * t;  y := Ky * t * t; (* (x,y) = a point on the parabola *)
    IF ABS (t) <= 0.000000001 THEN (* (u,v) = a normal vector *)
      u := 0.0;  v := 1.0
    ELSE
      u := 1.0;  v := -Kx / (2.0 * Ky * t)
    END;
    x0 := x - (1000.0 * u);  y0 := y - (1000.0 * v);
    x1 := x + (1000.0 * u);  y1 := y + (1000.0 * v);
    DrawLine (s, x0, y0, x1, y1);
  END;

  WriteTransparencyMask (s)
END Xiaolin_Wu_Task.

Here is a shell script that compiles the program, runs it, and (using Netpbm commands) makes a PNG using the outputted mask.

#!/bin/sh

# Set GM2 to wherever you have a GNU Modula-2 compiler.
GM2="/usr/x86_64-pc-linux-gnu/gcc-bin/13/gm2"

${GM2} -g -fbounds-check -fiso xiaolin_wu_line_algorithm_Modula2.mod
./a.out > alpha.pgm
ppmmake rgb:5C/06/8C 600 400 > bg.ppm
ppmmake rgb:E2/E8/68 600 400 > fg.ppm
pamcomp -alpha=alpha.pgm fg.ppm bg.ppm | pamtopng > image.png
Output:
A parabola, and lines normal to the parabola, forming a semicubic parabola as their envelope. A shade of yellow on a background shade of purple.

Nim

Simple translation of the Wikipedia algorithm.

import math
import imageman

template ipart(x: float): float = floor(x)
template fpart(x: float): float = x - ipart(x)
template rfpart(x: float): float = 1 - fpart(x)

const
  BG = ColorRGBF64([0.0, 0.0, 0.0])
  FG = ColorRGBF64([1.0, 1.0, 1.0])

func plot(img: var Image; x, y: int; c: float) =
  ## Draw a point with brigthness c.
  let d = 1 - c
  img[x, y] = ColorRGBF64([BG.r * d + FG.r * c, BG.g * d + FG.g * c, BG.b * d + FG.b * c])


func drawLine(img: var Image; x0, y0, x1, y1: float) =
  ## Draw an anti-aliased line from (x0, y0) to (x1, y1).

  var (x0, y0, x1, y1) = (x0, y0, x1, y1)
  let steep = abs(y1 - y0) > abs(x1 - x0)
  if steep:
    swap x0, y0
    swap x1, y1
  if x0 > x1:
    swap x0, x1
    swap y0, y1

  let dx = x1 - x0
  let dy = y1 - y0
  var gradient = dy / dx
  if dx == 0:
    gradient = 1

  # Handle first endpoint.
  var xend = round(x0)
  var yend = y0 + gradient * (xend - x0)
  var xgap = rfpart(x0 + 0.5)
  let xpxl1 = xend.toInt
  let ypxl1 = yend.toInt
  if steep:
    img.plot(ypxl1, xpxl1, rfpart(yend) * xgap)
    img.plot(ypxl1 + 1, xpxl1, fpart(yend) * xgap)
  else:
    img.plot(xpxl1, ypxl1, rfpart(yend) * xgap)
    img.plot(xpxl1, ypxl1 + 1, fpart(yend) * xgap)
  var 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)
  let xpxl2 = xend.toInt
  let ypxl2 = yend.toInt
  if steep:
    img.plot(ypxl2, xpxl2, rfpart(yend) * xgap)
    img.plot(ypxl2 + 1, xpxl2, fpart(yend) * xgap)
  else:
    img.plot(xpxl2, ypxl2, rfpart(yend) * xgap)
    img.plot(xpxl2, ypxl2 + 1, fpart(yend) * xgap)

  # Main loop.
  if steep:
    for x in (xpxl1 + 1)..(xpxl2 - 1):
      img.plot(intery.int, x, rfpart(intery))
      img.plot(intery.int + 1, x, fpart(intery))
      intery += gradient
  else:
    for x in (xpxl1 + 1)..(xpxl2 - 1):
      img.plot(x, intery.int, rfpart(intery))
      img.plot(x, intery.int + 1, fpart(intery))
      intery += gradient


when isMainModule:
  var img = initImage[ColorRGBF64](800, 800)
  img.fill(BG)
  for x1 in countup(100, 700, 60):
    img.drawLine(400, 700, x1.toFloat, 100)
  img.savePNG("xiaoling_wu.png", compression = 9)

ObjectIcon

The program puts up a window. In the window you can draw a line by left-mouse-button-press for one endpoint, and then another press for the other endpoint. You can draw multiple lines. When you leave (by pressing "q", for instance, or closing the window), the program stores the image as a PNG.

Rather than vary the color as such, I vary the opacity.

import
  graphics(Mouse, Window),
  io(stop),
  ipl.graphics(QuitEvents)

procedure main ()
  local width, height
  local done, w, event
  local x1, y1, x2, y2, press_is_active

  width := 640
  height := 480

  w := Window().
    set_size(width, height).
    set_bg("white").
    set_canvas("normal") | stop(&why)

  press_is_active := &no
  done := &no
  while /done do
  {
    if *w.pending() ~= 0 then
    {
      event := w.event()
      case event[1] of
      {
        QuitEvents() : done := &yes
        Mouse.LEFT_PRESS:
        {
          if /press_is_active then
          {
            x1 := event[2];  y1 := event[3]
            press_is_active := &yes
          }
          else
          {
            x2 := event[2];  y2 := event[3]
            draw_line (w, x1, y1, x2, y2)
            press_is_active := &no
          }
        }
      }
    }
  }

  w.get_pixels().to_file("xiaolin_wu_line_algorithm_OI.png")
end

procedure draw_line (w, x0, y0, x1, y1)
  local steep
  local dx, dy, gradient
  local xend, yend, xgap, intery
  local xpxl1, ypxl1
  local xpxl2, ypxl2
  local x

  x0 := real (x0)
  y0 := real (y0)
  x1 := real (x1)
  y1 := real (y1)

  # In Object Icon (as in Icon), comparisons DO NOT return boolean
  # values! They either SUCCEED or they FAIL. Thus the need for an
  # "if-then-else" here.
  steep := if abs (y1 - y0) > abs (x1 - x0) then &yes else &no

  if \steep then { x0 :=: y0;  x1 :=: y1 }
  if x0 > x1 then { x0 :=: x1;  y0 :=: y1 }
  dx := x1 - x0;  dy := y1 - y0
  gradient := if dx = 0 then 1.0 else dy / dx

  # Handle the first endpoint.
  xend := round (x0);  yend := y0 + (gradient * (xend - x0))
  xgap := rfpart (x0 + 0.5)
  xpxl1 := xend;  ypxl1 := ipart (yend)
  if \steep then
  {
    plot (w, ypxl1, xpxl1, rfpart (yend) * xgap)
    plot (w, ypxl1 + 1, xpxl1, fpart(yend) * xgap)
  }
  else
  {
    plot (w, xpxl1, ypxl1, rfpart (yend) * xgap)
    plot (w, xpxl1, ypxl1 + 1, fpart (yend) * xgap)
  }

  # The first y-intersection.
  intery := yend + gradient

  # Handle the second endpoint.
  xend := round (x1);  yend := y1 + (gradient * (xend - x1))
  xgap := fpart (x1 + 0.5)
  xpxl2 := xend;  ypxl2 := ipart (yend)
  if \steep then
  {
    plot (w, ypxl2, xpxl2, rfpart (yend) * xgap)
    plot (w, ypxl2 + 1, xpxl2, fpart (yend) * xgap)
  }
  else
  {
    plot (w, xpxl2, ypxl2, rfpart (yend) * xgap)
    plot (w, xpxl2, ypxl2 + 1, fpart (yend) * xgap)
  }

  if \steep then
    every x := xpxl1 + 1 to xpxl2 - 1 do
    {
      plot (w, ipart (intery), x, rfpart (intery))
      plot (w, ipart (intery) + 1, x, fpart (intery))
      intery := intery + gradient
    }
  else
    every x := xpxl1 + 1 to xpxl2 - 1 do
    {
      plot(w, x, ipart (intery), rfpart (intery))
      plot(w, x, ipart (intery) + 1, fpart (intery))
      intery := intery + gradient
    }

  return
end

procedure plot (w, x, y, c)
  w.set_fg ("black " || round (100.0 * c) || "%")
  w.draw_point (x, y)
  return
end

procedure ipart (x)
  local i
  i := integer (x)
  return (if i = x then i else if x < 0 then i - 1 else i)
end

procedure round (x)
  return ipart (x + 0.5)
end

procedure fpart (x)
  return x - ipart (x)
end

procedure rfpart (x)
  return 1.0 - fpart (x)
end
Output:

An example:

Some straight lines, antialiased. Black on white.

Pascal

Works with: Free Pascal
Library: SDL2

Based on Wikipwdia pseudocode with some optimizations and alpha handling.

program wu;
uses
  SDL2,
  math;

const
  FPS = 1000 div 60;
  SCALE = 6;

var
  win: PSDL_Window;
  ren: PSDL_Renderer;
  mouse_x, mouse_y: longint;
  origin: TSDL_Point;
  event: TSDL_Event;
  line_alpha: byte = 255;

procedure SDL_RenderDrawWuLine(renderer: PSDL_Renderer; x1, y1, x2, y2: longint);
var
  r, g, b, a, a_new: Uint8;
  gradient, iy: real;
  x, y: longint;
  px, py: plongint;

  procedure swap(var a, b: longint);
  var
    tmp: longint;
  begin
    tmp := a;
    a := b;
    b := tmp;
  end;

begin
  if a = 0 then
    exit;
  SDL_GetRenderDrawColor(renderer, @r, @g, @b, @a);
  if abs(y2 - y1) > abs(x2 - x1) then
    begin
      swap(x1, y1);
      swap(x2, y2);
      px := @y;
      py := @x;
    end
  else
    begin
      px := @x;
      py := @y;
    end;
  if x1 > x2 then
    begin
      swap(x1, x2);
      swap(y1, y2);
    end;
  x := x2 - x1;
  if x = 0 then
    x := 1;
  gradient := (y2 - y1) / x;
  iy := y1;
  for x := x1 to x2 do
    begin
      a_new := round(a * frac(iy));
      y := floor(iy);
      SDL_SetRenderDrawColor(renderer, r, g, b, a-a_new);
      SDL_RenderDrawPoint(renderer, px^, py^);
      inc(y);
      SDL_SetRenderDrawColor(renderer, r, g, b, a_new);
      SDL_RenderDrawPoint(renderer, px^, py^);
      iy := iy + gradient;
    end;
  SDL_SetRenderDrawColor(renderer, r, g, b, a);
end;

begin
  SDL_Init(SDL_INIT_VIDEO);
  win := SDL_CreateWindow('Xiaolin Wu''s line algorithm', SDL_WINDOWPOS_CENTERED, SDL_WINDOWPOS_CENTERED,
                        640, 480, SDL_WINDOW_RESIZABLE);
  ren := SDL_CreateRenderer(win, -1, 0);
  if ren = NIL then
    begin
      writeln(SDL_GetError);
      halt;
    end;
  SDL_SetRenderDrawBlendMode(ren, SDL_BLENDMODE_BLEND);
  SDL_RenderSetScale(ren, SCALE, SCALE);
  SDL_SetCursor(SDL_CreateSystemCursor(SDL_SYSTEM_CURSOR_CROSSHAIR));

  mouse_x := 0;
  mouse_y := 0;
  origin.x := 0;
  origin.y := 0;
  repeat
    while SDL_PollEvent(@event) = 1 do
      case event.type_ of
        SDL_KEYDOWN:
          if event.key.keysym.sym = SDLK_ESCAPE then
            halt;
        SDL_MOUSEBUTTONDOWN:
          begin
            origin.x := mouse_x;
            origin.y := mouse_y;
          end;
        SDL_MOUSEMOTION:
          with event.motion do
            begin
              mouse_x := x div SCALE;
              mouse_y := y div SCALE;
            end;
        SDL_MOUSEWHEEL:
          line_alpha := EnsureRange(line_alpha + event.wheel.y * 20, 0, 255);
        SDL_QUITEV:
          halt;
      end;

    SDL_SetRenderDrawColor(ren, 35, 35, 35, line_alpha);
    SDL_RenderDrawWuLine(ren, origin.x, origin.y, mouse_x, mouse_y);
    SDL_RenderPresent(ren);
    SDL_SetRenderDrawColor(ren, 255, 255, 255, 255);
    SDL_RenderClear(ren);
    SDL_Delay(FPS);
  until false;
end.

Perl

This is mostly a translation of the pseudo-code on Wikipedia, except that the $plot trick was inspired by the Raku example.

#!perl
use strict;
use warnings;

sub plot {
	my ($x, $y, $c) = @_;
	printf "plot %d %d %.1f\n", $x, $y, $c if $c;
}

sub ipart {
	int shift;
}

sub round {
	int( 0.5 + shift );
}

sub fpart {
	my $x = shift;
	$x - int $x;
}

sub rfpart {
	1 - fpart(shift);
}
 
sub drawLine {
	my ($x0, $y0, $x1, $y1) = @_;

	my $plot = \&plot;

	if( abs($y1 - $y0) > abs($x1 - $x0) ) {
		$plot = sub { plot( @_[1, 0, 2] ) };
		($x0, $y0, $x1, $y1) = ($y0, $x0, $y1, $x1);
	}

	if( $x0 > $x1 ) {
		($x0, $x1, $y0, $y1) = ($x1, $x0, $y1, $y0);
	}

	my $dx = $x1 - $x0;
	my $dy = $y1 - $y0;
	my $gradient = $dy / $dx;

	my @xends;
	my $intery;

	# handle the endpoints
	for my $xy ([$x0, $y0], [$x1, $y1]) {
		my ($x, $y) = @$xy;
		my $xend = round($x);
		my $yend = $y + $gradient * ($xend - $x);
		my $xgap = rfpart($x + 0.5);

		my $x_pixel = $xend;
		my $y_pixel = ipart($yend);
		push @xends, $x_pixel;

		$plot->($x_pixel, $y_pixel  , rfpart($yend) * $xgap);
		$plot->($x_pixel, $y_pixel+1,  fpart($yend) * $xgap);
		next if defined $intery;
		# first y-intersection for the main loop
		$intery = $yend + $gradient;
	}

	# main loop

	for my $x ( $xends[0] + 1 .. $xends[1] - 1 ) {
		$plot->($x, ipart ($intery),  rfpart($intery));
		$plot->($x, ipart ($intery)+1, fpart($intery));
		$intery += $gradient;
	}
}

if( $0 eq __FILE__ ) {
	drawLine( 0, 1, 10, 2 );
}
__END__
Output:
plot 0 1 0.5
plot 10 2 0.5
plot 1 1 0.9
plot 1 2 0.1
plot 2 1 0.8
plot 2 2 0.2
plot 3 1 0.7
plot 3 2 0.3
plot 4 1 0.6
plot 4 2 0.4
plot 5 1 0.5
plot 5 2 0.5
plot 6 1 0.4
plot 6 2 0.6
plot 7 1 0.3
plot 7 2 0.7
plot 8 1 0.2
plot 8 2 0.8
plot 9 1 0.1
plot 9 2 0.9

Phix

Library: Phix/pGUI
Library: Phix/online

You can run this online here, with caveats as below. Resize the window to show lines at any angle

--
-- 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()

PicoLisp

(scl 2)

(de plot (Img X Y C)
   (set (nth Img (*/ Y 1.0) (*/ X 1.0)) (- 100 C)) )

(de ipart (X)
   (* 1.0 (/ X 1.0)) )

(de iround (X)
   (ipart (+ X 0.5)) )

(de fpart (X)
   (% X 1.0) )

(de rfpart (X)
   (- 1.0 (fpart X)) )

(de xiaolin (Img X1 Y1 X2 Y2)
   (let (DX (- X2 X1)  DY (- Y2 Y1))
      (use (Grad Xend Yend Xgap Xpxl1 Ypxl1 Xpxl2 Ypxl2 Intery)
         (when (> (abs DY) (abs DX))
            (xchg 'X1 'Y1  'X2 'Y2) )
         (when (> X1 X2)
            (xchg 'X1 'X2  'Y1 'Y2) )
         (setq
            Grad (*/ DY 1.0 DX)
            Xend (iround X1)
            Yend (+ Y1 (*/ Grad (- Xend X1) 1.0))
            Xgap (rfpart (+ X1 0.5))
            Xpxl1 Xend
            Ypxl1 (ipart Yend) )
         (plot Img Xpxl1 Ypxl1 (*/ (rfpart Yend) Xgap 1.0))
         (plot Img Xpxl1 (+ 1.0 Ypxl1) (*/ (fpart Yend) Xgap 1.0))
         (setq
            Intery (+ Yend Grad)
            Xend (iround X2)
            Yend (+ Y2 (*/ Grad (- Xend X2) 1.0))
            Xgap (fpart (+ X2 0.5))
            Xpxl2 Xend
            Ypxl2 (ipart Yend) )
         (plot Img Xpxl2 Ypxl2 (*/ (rfpart Yend) Xgap 1.0))
         (plot Img Xpxl2 (+ 1.0 Ypxl2) (*/ (fpart Yend) Xgap 1.0))
         (for (X (+ Xpxl1 1.0)  (>= (- Xpxl2 1.0) X)  (+ X 1.0))
            (plot Img X (ipart Intery) (rfpart Intery))
            (plot Img X (+ 1.0 (ipart Intery)) (fpart Intery))
            (inc 'Intery Grad) ) ) ) )

(let Img (make (do 90 (link (need 120 99))))       # Create image 120 x 90
   (xiaolin Img 10.0 10.0 110.0 80.0)              # Draw lines
   (xiaolin Img 10.0 10.0 110.0 45.0)
   (xiaolin Img 10.0 80.0 110.0 45.0)
   (xiaolin Img 10.0 80.0 110.0 10.0)
   (out "img.pgm"                                  # Write to bitmap file
      (prinl "P2")
      (prinl 120 " " 90)
      (prinl 100)
      (for Y Img (apply printsp Y)) ) )

PureBasic

Macro PlotB(x, y, Color, b)
  Plot(x, y, RGB(Red(Color) * (b), Green(Color) * (b), Blue(Color) * (b)))
EndMacro

Procedure.f fracPart(x.f)
  ProcedureReturn x - Int(x)
EndProcedure

Procedure.f invFracPart(x.f)
  ProcedureReturn 1.0 - fracPart(x)
EndProcedure

Procedure drawAntiAliasedLine(x1.f, y1.f, x2.f, y2.f, color)
  Protected.f dx, dy, xend, yend, grad, yf, xgap, ix1, iy1, ix2, iy2
  Protected x
  
  dx = x2 - x1
  dy = y2 - y1
  If Abs(dx) < Abs(dy)
    Swap x1, y1
    Swap x2, y2
    Swap dx, dy
  EndIf
  
  If x2 < x1
    Swap x1, x2
    Swap y1, y2
  EndIf
  
  grad = dy / dx
  
  ;handle first endpoint
  xend = Round(x1, #pb_round_nearest)
  yend = y1 + grad * (xend - x1)
  xgap = invFracPart(x1 + 0.5)
  ix1 = xend  ;this will be used in the MAIN loop
  iy1 = Int(yend)
  PlotB(ix1, iy1, color, invFracPart(yend) * xgap)
  PlotB(ix1, iy1 + 1, color, fracPart(yend) * xgap)
  yf = yend + grad ;first y-intersection for the MAIN loop
  
  ;handle second endpoint
  xend = Round(x2, #pb_round_nearest)
  yend = y2 + grad * (xend - x2)
  xgap = fracPart(x2 + 0.5)
  ix2 = xend  ;this will be used in the MAIN loop
  iy2 = Int(yend)
  PlotB(ix2, iy2, color, invFracPart(yend) * xgap)
  PlotB(ix2, iy2 + 1, color, fracPart(yend) * xgap)
  ;MAIN loop
  For x = ix1 + 1 To ix2 - 1
    PlotB(x, Int(yf), color, invFracPart(yf))
    PlotB(x, Int(yf) + 1, color, fracPart(yf))
    yf + grad
  Next 
EndProcedure

Define w = 200, h = 200, img = 1
CreateImage(img, w, h) ;img is internal id of the image

OpenWindow(0, 0, 0, w, h,"Xiaolin Wu's line algorithm", #PB_Window_SystemMenu)

StartDrawing(ImageOutput(img))
  drawAntiAliasedLine(80,20, 130,80, RGB(255, 0, 0))
StopDrawing()

ImageGadget(0, 0, 0, w, h, ImageID(img))

Define event
Repeat
  event = WaitWindowEvent()
Until event = #PB_Event_CloseWindow

Python

"""Script demonstrating drawing of anti-aliased lines using Xiaolin Wu's line
algorithm

usage: python xiaolinwu.py [output-file]

"""
from __future__ import division
import sys

from PIL import Image


def _fpart(x):
    return x - int(x)

def _rfpart(x):
    return 1 - _fpart(x)

def putpixel(img, xy, color, alpha=1):
    """
    Paints color over the background at the point xy in img.
    Use alpha for blending. alpha=1 means a completely opaque foreground.
    """
    compose_color = lambda bg, fg: int(round(alpha * fg + (1-alpha) * bg))
    c = compose_color(img.getpixel(xy), color)
    img.putpixel(xy, c)

def draw_line(img, p1, p2, color):
    """Draws an anti-aliased line in img from p1 to p2 with the given color."""
    x1, y1 = p1
    x2, y2 = p2
    dx, dy = x2-x1, y2-y1
    steep = abs(dx) < abs(dy)
    p = lambda px, py: ((px,py), (py,px))[steep]

    if steep:
        x1, y1, x2, y2, dx, dy = y1, x1, y2, x2, dy, dx
    if x2 < x1:
        x1, x2, y1, y2 = x2, x1, y2, y1

    grad = dy/dx
    intery = y1 + _rfpart(x1) * grad
    def draw_endpoint(pt):
        x, y = pt
        xend = round(x)
        yend = y + grad * (xend - x)
        xgap = _rfpart(x + 0.5)
        px, py = int(xend), int(yend)
        putpixel(img, p(px, py), color, _rfpart(yend) * xgap)
        putpixel(img, p(px, py+1), color, _fpart(yend) * xgap)
        return px

    xstart = draw_endpoint(p(*p1)) + 1
    xend = draw_endpoint(p(*p2))

    for x in range(xstart, xend):
        y = int(intery)
        putpixel(img, p(x, y), color, _rfpart(intery))
        putpixel(img, p(x, y+1), color, _fpart(intery))
        intery += grad


if __name__ == '__main__':
    if len(sys.argv) != 2:
        print 'usage: python xiaolinwu.py [output-file]'
        sys.exit(-1)

    blue = (0, 0, 255)
    yellow = (255, 255, 0)
    img = Image.new("RGB", (500,500), blue)
    for a in range(10, 431, 60):
        draw_line(img, (10, 10), (490, a), yellow)
        draw_line(img, (10, 10), (a, 490), yellow)
    draw_line(img, (10, 10), (490, 490), yellow)
    filename = sys.argv[1]
    img.save(filename)
    print 'image saved to', filename

Racket

#lang racket
(require 2htdp/image)

(define (plot img x y c)
  (define c*255 (exact-round (* (- 1 c) 255)))
  (place-image
   (rectangle 1 1 'solid (make-color c*255 c*255 c*255 255))
   x y img))

(define ipart exact-floor) ; assume that a "round-down" is what we want when -ve
;;; `round` is built in -- but we'll use exact round (and I'm not keen on over-binding round)

(define (fpart n) (- n (exact-floor n)))
(define (rfpart n) (- 1 (fpart n)))

(define (draw-line img x0 y0 x1 y1)
  (define (draw-line-steeped img x0 y0 x1 y1 steep?)
    (define (draw-line-steeped-l-to-r img x0 y0 x1 y1 steep?)
      (define dx (- x1 x0))
      (define dy (- y1 y0))
      (define gradient (/ dy dx))
      
      (define (handle-end-point img x y)
        (define xend (exact-round x))
        (define yend (+ y (* gradient (- xend x))))
        (define xgap (rfpart (+ x 0.5)))
        (define ypxl (ipart yend))
        (define intery (+ yend gradient))
        
        (case steep?
          [(#t)
           (define img* (plot img ypxl xend (* xgap (rfpart yend))))
           (values (plot img* (+ ypxl 1) xend (* xgap (fpart yend))) xend intery)]
          [(#f)
           (define img* (plot img xend ypxl (* xgap (rfpart yend))))
           (values (plot img* xend (+ ypxl 1) (* xgap (fpart yend))) xend intery)]))
      
      (define-values (img-with-l-endpoint xpl1 intery) (handle-end-point img x0 y0))
      (define-values (img-with-r-endpoint xpl2 _) (handle-end-point img-with-l-endpoint x1 y1))
      
      (for/fold ((img img-with-l-endpoint)  (y intery))
        ((x (in-range (+ xpl1 1) xpl2)))
        (define y-i (ipart y))
        (values
         (case steep?
           [(#t)
            (define img* (plot img y-i x (rfpart y)))
            (plot img* (+ 1 y-i) x (fpart y))]
           [(#f)
            (define img* (plot img x y-i (rfpart y)))
            (plot img* x (+ 1 y-i) (fpart y))])
         (+ y gradient))))
    
    (if (> x0 x1)
        (draw-line-steeped-l-to-r img x1 y1 x0 y0 steep?)
        (draw-line-steeped-l-to-r img x0 y0 x1 y1 steep?)))
  
  (define steep? (> (abs (- y1 y0)) (abs (- x1 x0))))
  (define-values (img* _)
    (if steep?
        (draw-line-steeped img y0 x0 y1 x1 steep?)
        (draw-line-steeped img x0 y0 x1 y1 steep?)))
  img*)

(define img-1
(beside
 (scale 3 (draw-line (empty-scene 150 100) 12 12 138 88))
 (above
  (scale 1 (draw-line (empty-scene 150 100) 12 50 138 50))
  (scale 1 (draw-line (empty-scene 150 100) 75 12 75 88))
  (scale 1 (draw-line (empty-scene 150 100) 12 88 138 12)))))

(define img-2
  (beside
 (scale 3 (draw-line (empty-scene 100 150) 12 12 88 138))
 (above (scale 1 (draw-line (empty-scene 100 150) 50 12 50 138))
        (scale 1 (draw-line (empty-scene 100 150) 12 75 88 75))
        (scale 1 (draw-line (empty-scene 100 150) 88 12 12 138)))))

img-1
img-2
(save-image img-1 "images/xiaolin-wu-racket-1.png")
(save-image img-2 "images/xiaolin-wu-racket-2.png")

Output files: File:Xiaolin-wu-racket-1.png File:Xiaolin-wu-racket-2.png

Raku

(formerly Perl 6)

sub plot(\x, \y, \c) { say "plot {x} {y} {c}" }
 
sub fpart(\x) { x - floor(x) }
 
sub draw-line(@a is copy, @b is copy) {
    my Bool \steep = abs(@b[1] - @a[1]) > abs(@b[0] - @a[0]);
    my $plot = &OUTER::plot;
 
    if steep {
	$plot = -> $y, $x, $c { plot($x, $y, $c) }
	@a.=reverse;
	@b.=reverse;
    }
    if @a[0] > @b[0] { my @t = @a; @a = @b; @b = @t }

    my (\x0,\y0) = @a;
    my (\x1,\y1) = @b;
 
    my \dx = x1 - x0;
    my \dy = y1 - y0;
    my \gradient = dy / dx;
 
    # handle first endpoint
    my \x-end1 = round(x0);
    my \y-end1 = y0 + gradient * (x-end1 - x0);
    my \x-gap1 = 1 - round(x0 + 0.5);

    my \x-pxl1 = x-end1;   # this will be used in the main loop
    my \y-pxl1 = floor(y-end1);
    my \c1 = fpart(y-end1) * x-gap1;

    $plot(x-pxl1, y-pxl1    , 1 - c1) unless c1 == 1;
    $plot(x-pxl1, y-pxl1 + 1, c1    ) unless c1 == 0;
 
    # handle second endpoint
    my \x-end2 = round(x1);
    my \y-end2 = y1 + gradient * (x-end2 - x1);
    my \x-gap2 = fpart(x1 + 0.5);

    my \x-pxl2 = x-end2; # this will be used in the main loop
    my \y-pxl2 = floor(y-end2);
    my \c2 = fpart(y-end2) * x-gap2;
 
    my \intery = y-end1 + gradient;

    # main loop
    for (x-pxl1 + 1 .. x-pxl2 - 1)
	Z
	(intery, intery + gradient ... *)
    -> (\x,\y) {
	my \c = fpart(y);
	$plot(x, floor(y)    , 1 - c) unless c == 1;
	$plot(x, floor(y) + 1, c    ) unless c == 0;
    }

    $plot(x-pxl2, y-pxl2    , 1 - c2) unless c2 == 1;
    $plot(x-pxl2, y-pxl2 + 1, c2    ) unless c2 == 0;
}

draw-line [0,1], [10,2];
Output:
plot 0 1 1
plot 1 1 0.9
plot 1 2 0.1
plot 2 1 0.8
plot 2 2 0.2
plot 3 1 0.7
plot 3 2 0.3
plot 4 1 0.6
plot 4 2 0.4
plot 5 1 0.5
plot 5 2 0.5
plot 6 1 0.4
plot 6 2 0.6
plot 7 1 0.3
plot 7 2 0.7
plot 8 1 0.2
plot 8 2 0.8
plot 9 1 0.1
plot 9 2 0.9
plot 10 2 1

REXX

This REXX example uses the Xiaolin Wu line algorithm to draw a line (with output).

Apparently, there may be an error in the definition of the algorithm (which only manifests itself with negative numbers):
use of the   IPART   function should probably be   FLOOR.

[See the   talk   section on the Xiaolin Wu's line algorithm.]

http://en.wikipedia.org/wiki/Talk:Xiaolin_Wu%27s_line_algorithm


Also, it takes in account (that can easily be overlooked) of the note after the description of the algorithm:
Note:   If at the beginning of the routine   abs(dx) < abs(dy)   is true, then all plotting should be done with   x   and   y   reversed.

/*REXX program  plots/draws (ASCII)  a   line   using the   Xiaolin Wu  line algorithm. */
background= '·'                                  /*background character:  a middle-dot. */
    image.= background                           /*fill the array with middle-dots.     */
     plotC= '░▒▓█'                               /*characters used for plotting points. */
       EoE= 3000                                 /*EOE = End Of Earth,  er, ··· graph.  */
                     do j=-EoE  to +EoE          /*define the graph: lowest ──► highest.*/
                     image.j.0= '─'              /*define the graph's horizontal axis.  */
                     image.0.j= '│'              /*   "    "     "    verical      "    */
                     end   /*j*/
 image.0.0= '┼'                                  /*define the graph's axis origin (char)*/
parse arg xi yi xf yf .                          /*allow specifying the line-end points.*/
if xi=='' | xi==","  then xi= 1                  /*Not specified?  Then use the default.*/
if yi=='' | yi==","  then yi= 2                  /* "      "         "   "   "     "    */
if xf=='' | xf==","  then xf=11                  /* "      "         "   "   "     "    */
if yf=='' | yf==","  then yf=12                  /* "      "         "   "   "     "    */
minX=0;    minY=0                                /*use these as the limits for plotting.*/
maxX=0;    maxY=0                                /* "    "    "  "    "     "      "    */
call drawLine  xi, yi, xf, yf                    /*invoke subroutine and graph the line.*/
border=2                                         /*allow additional space (plot border).*/
minX=minX - border * 2;  maxX=maxX + border * 2  /*preserve screen's aspect ratio  {*2}.*/
minY=minY - border    ;  maxY=maxY + border
                              do     y=maxY  to minY  by -1;  $=      /*construct a row.*/
                                  do x=minX  to maxX;       $=$ || image.x.y;   end  /*x*/
                              say $              /*display the constructed row to term. */
                              end   /*y*/        /*graph is cropped by the MINs and MAXs*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
drawLine: parse arg x1,y1,x2,y2;      switchXY=0;        dx=x2-x1
                                                         dy=y2-y1
          if abs(dx)<abs(dy)  then parse value  x1 y1 x2 y2 dx dy  with  y1 x2 y2 x2 dy dx
          if x2<x1       then parse value  x1 x2 y1 y2   1   with   x2 x1 y2 y1   switchXY
          gradient=dy/dx
              xend=round(x1)             /*◄─────────────────1st endpoint.══════════════*/
              yend=y1 + gradient * (xend-x1);      xgap=1 - fpart(x1 + .5)
             xpx11=xend;       ypx11=floor(yend)
            intery=yend+gradient
          call plotXY  xpx11,  ypx11,    brite(1 - fpart(yend*xgap)), switchXY
          call plotXY  xpx11,  ypx11+1,  brite(    fpart(yend*xgap)), switchXY
              xend=round(x2)             /*◄─────────────────2nd endpoint.══════════════*/
              yend=y2 + gradient * (xend-x2);      xgap=    fpart(x2 + .5)
             xpx12=xend;       ypx12=floor(yend)
          call plotXY  xpx12,  ypx12  ,  brite(1 - fpart(yend*xgap)), switchXY
          call plotXY  xpx12,  ypx12+1,  brite(    fpart(yend*xgap)), switchXY

                do x=xpx11+1  to xpx12-1 /*◄═════════════════draw the line.═════════════*/
                !intery=floor(intery)
                call plotXY  x,  !intery  ,  brite(1 - fpart(intery)), switchXY
                call plotXY  x,  !intery+1,  brite(    fpart(intery)), switchXY
                intery=intery + gradient
                end   /*x*/
          return
/*──────────────────────────────────────────────────────────────────────────────────────*/
brite:    return substr(background || plotC, 1 + round( abs( arg(1) ) * length(plotC)), 1)
floor:    parse arg #; _=trunc(#);  return _   -   (#<0) * (#\=_)
fpart:    parse arg #;              return abs(# - trunc(#) )
round:    return   format(arg(1), , word(arg(2) 0, 1) )
/*──────────────────────────────────────────────────────────────────────────────────────*/
plotXY:   parse arg xx,yy,bc,switchYX;             if switchYX  then parse arg yy,xx
          image.xx.yy=bc;    minX=min(minX, xx);   maxX=max(maxX,xx)
                             minY=min(minY, yy);   maxY=max(maxY,yy);               return
output   when using the default inputs:
····│···············
····│···············
····│···············
····│··········█····
····│·········█·····
····│········█······
····│·······█·······
····│······█········
····│·····█·········
····│····█··········
····│···█···········
····│··█············
····│·█·············
····│█··············
····│···············
────┼───────────────
····│···············
····│···············

Ruby

Translation of: Tcl
def ipart(n); n.truncate; end
def fpart(n); n - ipart(n); end
def rfpart(n); 1.0 - fpart(n); end

class Pixmap
  def draw_line_antialised(p1, p2, colour)
    x1, y1 = p1.x, p1.y
    x2, y2 = p2.x, p2.y
 
    steep = (y2 - y1).abs > (x2 - x1).abs
    if steep
      x1, y1 = y1, x1
      x2, y2 = y2, x2
    end
    if x1 > x2
      x1, x2 = x2, x1
      y1, y2 = y2, y1
    end
    deltax = x2 - x1
    deltay = (y2 - y1).abs
    gradient = 1.0 * deltay / deltax
 
    # handle the first endpoint
    xend = x1.round
    yend = y1 + gradient * (xend - x1)
    xgap = rfpart(x1 + 0.5)
    xpxl1 = xend
    ypxl1 = ipart(yend)
    put_colour(xpxl1, ypxl1, colour, steep, rfpart(yend)*xgap)
    put_colour(xpxl1, ypxl1 + 1, colour, steep, fpart(yend)*xgap)
    itery = yend + gradient
 
    # handle the second endpoint
    xend = x2.round
    yend = y2 + gradient * (xend - x2)
    xgap = rfpart(x2 + 0.5)
    xpxl2 = xend
    ypxl2 = ipart(yend)
    put_colour(xpxl2, ypxl2, colour, steep, rfpart(yend)*xgap)
    put_colour(xpxl2, ypxl2 + 1, colour, steep, fpart(yend)*xgap)
 
    # in between
    (xpxl1 + 1).upto(xpxl2 - 1).each do |x|
      put_colour(x, ipart(itery), colour, steep, rfpart(itery))
      put_colour(x, ipart(itery) + 1, colour, steep, fpart(itery))
      itery = itery + gradient
    end
  end

  def put_colour(x, y, colour, steep, c)
    x, y = y, x if steep
    self[x, y] = anti_alias(colour, self[x, y], c)
  end

  def anti_alias(new, old, ratio)
    blended = new.values.zip(old.values).map {|n, o| (n*ratio + o*(1.0 - ratio)).round}
    RGBColour.new(*blended)
  end
end

bitmap = Pixmap.new(500, 500)
bitmap.fill(RGBColour::BLUE)
10.step(430, 60) do |a|
  bitmap.draw_line_antialised(Pixel[10, 10], Pixel[490,a], RGBColour::YELLOW)
  bitmap.draw_line_antialised(Pixel[10, 10], Pixel[a,490], RGBColour::YELLOW)
end
bitmap.draw_line_antialised(Pixel[10, 10], Pixel[490,490], RGBColour::YELLOW)

Scala

Uses Bitmap#Scala.

import java.awt.Color
import math.{floor => ipart, round, abs}

case class Point(x: Double, y: Double) {def swap = Point(y, x)}

def plotter(bm: RgbBitmap, c: Color)(x: Double, y: Double, v: Double) = {
  val X = round(x).toInt
  val Y = round(y).toInt
  val V = v.toFloat
  // tint the existing pixels
  val c1 = c.getRGBColorComponents(null)
  val c2 = bm.getPixel(X, Y).getRGBColorComponents(null)
  val c3 = (c1 zip c2).map{case (n, o) => n * V + o * (1 - V)}
  bm.setPixel(X, Y, new Color(c3(0), c3(1), c3(2)))
}

def drawLine(plotter: (Double,Double,Double) => _)(p1: Point, p2: Point) {
  def fpart(x: Double) = x - ipart(x)
  def rfpart(x: Double) = 1 - fpart(x)
  def avg(a: Float, b: Float) = (a + b) / 2

  val steep = abs(p2.y - p1.y) > abs(p2.x - p1.x)
  val (p3, p4) = if (steep) (p1.swap, p2.swap) else (p1, p2)
  val (a, b) = if (p3.x > p4.x) (p4, p3) else (p3, p4)
  val dx = b.x - a.x
  val dy = b.y - a.y
  val gradient = dy / dx
  var intery = 0.0

  def endpoint(xpxl: Double, yend: Double, xgap: Double) {
    val ypxl = ipart(yend)
    if (steep) {
      plotter(ypxl,   xpxl, rfpart(yend) * xgap)
      plotter(ypxl+1, xpxl,  fpart(yend) * xgap)
    } else {
      plotter(xpxl, ypxl  , rfpart(yend) * xgap)
      plotter(xpxl, ypxl+1,  fpart(yend) * xgap)
    }
  }

  // handle first endpoint
  var xpxl1 = round(a.x);
  {
    val yend = a.y + gradient * (xpxl1 - a.x)
    val xgap = rfpart(a.x + 0.5)
    endpoint(xpxl1, yend, xgap)
    intery = yend + gradient
  }

  // handle second endpoint
  val xpxl2 = round(b.x);
  {
    val yend = b.y + gradient * (xpxl2 - b.x)
    val xgap = fpart(b.x + 0.5)
    endpoint(xpxl2, yend, xgap)
  }

  // main loop
  for (x <- (xpxl1 + 1) to (xpxl2 - 1)) {
    if (steep) {
      plotter(ipart(intery)  , x, rfpart(intery))
      plotter(ipart(intery)+1, x,  fpart(intery))
    } else {
      plotter(x, ipart (intery),  rfpart(intery))
      plotter(x, ipart (intery)+1, fpart(intery))
    }
    intery = intery + gradient
  }
}

Example:

Test line drawing in various directions including vertical, horizontal, 45° and oblique (such lines are drawn multiple times to test swapped parameters).

val r = 120
val img = new RgbBitmap(r*2+1, r*2+1)
val line = drawLine(plotter(img, Color.GRAY)_)_
img.fill(Color.WHITE)
for (angle <- 0 to 360 by 30; θ = math toRadians angle; θ2 = θ + math.Pi) {
  val a = Point(r + r * math.sin(θ), r + r * math.cos(θ))
  val b = Point(r + r * math.sin(θ2), r + r * math.cos(θ2))
  line(a, b)
}
javax.imageio.ImageIO.write(img.image, "png", new java.io.File("XiaolinWuLineAlgorithm.png"))
Output:

View the PNG, available at the following URL because RosettaCode image uploads were disabled: https://lh5.googleusercontent.com/GxBAHV4nebuO1uiKboKc6nQmmtlJV47jPwVZnQHcbV7TKm0kjdKfKteclCfxmSdFJnSKvYYoB5I

Scheme

Translation of: ATS
Works with: CHICKEN Scheme version 5.3.0
Works with: Gauche Scheme version 0.9.12

This program is written for R7RS. For CHICKEN you need the r7rs and srfi-160 eggs.

;;;-------------------------------------------------------------------

(import (scheme base))
(import (scheme file))
(import (scheme inexact))
(import (scheme process-context))
(import (scheme write))

;; (srfi 160 f32) is more properly known as (scheme vector f32), but
;; is not part of R7RS-small. The following will work in both Gauche
;; and CHICKEN Schemes.
(import (srfi 160 f32))

;;;-------------------------------------------------------------------

(define-record-type <color>
  (make-color r g b)
  color?
  (r color-r)
  (g color-g)
  (b color-b))

;;; See https://yeun.github.io/open-color/
(define violet9 (make-color (/ #x5F 255.0)
                            (/ #x3D 255.0)
                            (/ #xC4 255.0)))

;;;-------------------------------------------------------------------

(define-record-type <drawing-surface>
  (drawing-surface% u0 v0 u1 v1 pixels)
  drawing-surface?
  (u0 u0%)
  (v0 v0%)
  (u1 u1%)
  (v1 v1%)
  (pixels pixels%))

(define (make-drawing-surface u0 v0 u1 v1)
  (unless (and (<= u0 u1) (<= v0 v1))
    (error "illegal drawing-surface corners"))
  (let ((width (- u1 u0 -1))
        (height (- v1 v0 -1)))
    (let ((pixels (make-f32vector (* width height) 0.0)))
      (drawing-surface% u0 v0 u1 v1 pixels))))

;;; In calls to drawing-surface-ref and drawing-surface-set! indices
;;; outside the drawing_surface are allowed. Such indices are treated
;;; as if you were trying to draw on the air.

(define (drawing-surface-ref s x y)
  (let ((u0 (u0% s))
        (v0 (v0% s))
        (u1 (u1% s))
        (v1 (v1% s)))
    (if (and (<= u0 x) (<= x u1) (<= v0 y) (<= y v1))
        (f32vector-ref (pixels% s)
                       (+ (* (- x u0) (- v1 v0 -1)) (- v1 y)))
        +nan.0)))

(define (drawing-surface-set! s x y opacity)
  (let ((u0 (u0% s))
        (v0 (v0% s))
        (u1 (u1% s))
        (v1 (v1% s)))
    (when (and (<= u0 x) (<= x u1) (<= v0 y) (<= y v1))
      (f32vector-set! (pixels% s)
                      (+ (* (- x u0) (- v1 v0 -1)) (- v1 y))
                      opacity))))

(define (write-PAM s color)
  ;; Write a Portable Arbitrary Map to the current output port, using
  ;; the given color as the foreground color and the drawing-surface
  ;; values as opacities.

  (define (float->byte v) (exact (round (* 255 v))))

  (define r (float->byte (color-r color)))
  (define g (float->byte (color-g color)))
  (define b (float->byte (color-b color)))

  (define w (- (u1% s) (u0% s) -1))
  (define h (- (v1% s) (v0% s) -1))
  (define opacities (pixels% s))

  (define (loop x y)
    (cond ((= y h) )
          ((= x w) (loop 0 (+ y 1)))
          (else
           (let ((alpha (float->byte
                         (f32vector-ref opacities (+ (* x h) y)))))
             (write-bytevector (bytevector r g b alpha))
             (loop (+ x 1) y)))))

  (display "P7") (newline)
  (display "WIDTH ") (display (- (u1% s) (u0% s) -1)) (newline)
  (display "HEIGHT ") (display (- (v1% s) (v0% s) -1)) (newline)
  (display "DEPTH 4") (newline)
  (display "MAXVAL 255") (newline)
  (display "TUPLTYPE RGB_ALPHA") (newline)
  (display "ENDHDR") (newline)
  (loop 0 0))

;;;-------------------------------------------------------------------

(define (ipart x) (exact (floor x)))
(define (iround x) (ipart (+ x 0.5)))
(define (fpart x) (- x (floor x)))
(define (rfpart x) (- 1.0 (fpart x)))

(define (plot s x y opacity)
  ;; One might prefer a more sophisticated function than mere
  ;; addition. Here, however, the function is addition.
  (let ((combined-opacity (+ opacity (drawing-surface-ref s x y))))
    (drawing-surface-set! s x y (min combined-opacity 1.0))))

(define (drawln% s x0 y0 x1 y1 steep)
  (let* ((dx (- x1 x0))
         (dy (- y1 y0))
         (gradient (if (zero? dx) 1.0 (/ dy dx)))

         ;; Handle the first endpoint.
         (xend (iround x0))
         (yend (+ y0 (* gradient (- xend x0))))
         (xgap (rfpart (+ x0 0.5)))
         (xpxl1 xend)
         (ypxl1 (ipart yend))
         (_ (if steep
                (begin
                  (plot s ypxl1 xpxl1 (* (rfpart yend) xgap))
                  (plot s (+ ypxl1 1) xpxl1 (* (fpart yend) xgap)))
                (begin
                  (plot s xpxl1 ypxl1 (* (rfpart yend) xgap))
                  (plot s xpxl1 (+ ypxl1 1) (* (fpart yend) xgap)))))

         ;; The first y-intersection.
         (intery (+ yend gradient))

         ;; Handle the second endpoint.
         (xend (iround x1))
         (yend (+ y1 (* gradient (- xend x1))))
         (xgap (fpart (+ x1 0.5)))
         (xpxl2 xend)
         (ypxl2 (ipart yend))
         (_ (if steep
                (begin
                  (plot s ypxl2 xpxl2 (* (rfpart yend) xgap))
                  (plot s (+ ypxl2 1) xpxl2 (* (fpart yend) xgap)))
                (begin
                  (plot s xpxl2 ypxl2 (* (rfpart yend) xgap))
                  (plot s xpxl2 (+ ypxl2 1) (* (fpart yend) xgap))))))

    ;; Loop over the rest of the points.
    (if steep
        (do ((x (+ xpxl1 1) (+ x 1))
             (intery intery (+ intery gradient)))
            ((= x xpxl2))
          (plot s (ipart intery) x (rfpart intery))
          (plot s (+ (ipart intery) 1) x (fpart intery)))
        (do ((x (+ xpxl1 1) (+ x 1))
             (intery intery (+ intery gradient)))
            ((= x xpxl2))
          (plot s x (ipart intery) (rfpart intery))
          (plot s x (+ (ipart intery) 1) (fpart intery))))))

(define (draw-line s x0 y0 x1 y1)
  (let ((xdiff (abs (- x1 x0)))
        (ydiff (abs (- y1 y0))))
    (if (<= ydiff xdiff)
        (if (<= x0 x1)
            ;; R7RS lets you say #false and #true, as equivalents of
            ;; #f and #t. (To support such things as #false and #true,
            ;; the "r7rs" egg for CHICKEN Scheme 5 comes with a
            ;; special reader.)
            (drawln% s x0 y0 x1 y1 #false)
            (drawln% s x1 y1 x0 y0 #false))
        (if (<= y0 y1)
            (drawln% s y0 x0 y1 x1 #true)
            (drawln% s y1 x1 y0 x0 #true)))))

;;;-------------------------------------------------------------------

(define u0 0)
(define v0 0)
(define u1 999)
(define v1 749)

(define PI (* 4.0 (atan 1.0)))
(define PI/180 (/ PI 180.0))

(define (cosdeg theta) (cos (* theta PI/180)))
(define (sindeg theta) (sin (* theta PI/180)))

(define s (make-drawing-surface u0 v0 u1 v1))

;; The values of theta are exactly representable in either binary or
;; decimal floating point, and therefore the following loop will NOT
;; do the angle zero twice. (If you might stray from exact
;; representations, you must do something different, such as increment
;; an integer.)
(let ((x0 (inexact (* (/ 380 640) u1)))
      (y0 (inexact (* (/ 130 480) v1))))
  (do ((theta 0.0 (+ theta 5.0)))
      ((<= 360.0 theta))
    (let ((cos-theta (cosdeg theta))
          (sin-theta (sindeg theta)))
      (let ((x1 (+ x0 (* cos-theta 1200.0)))
            (y1 (+ y0 (* sin-theta 1200.0))))
        (draw-line s x0 y0 x1 y1)))))

(define args (command-line))
(unless (= (length args) 2)
  (parameterize ((current-output-port (current-error-port)))
    (display (string-append "Usage: " (car args) " FILENAME"))
    (newline)
    (display (string-append "       " (car args) " -"))
    (newline) (newline)
    (display (string-append "The second form writes the PAM file"
                            " to standard output."))
    (newline)
    (exit 1)))
(if (string=? (cadr args) "-")
    (write-PAM s violet9)
    (with-output-to-file (list-ref args 1)
      (lambda () (write-PAM s violet9))))

;;;-------------------------------------------------------------------
Output:

The output of the program is a Portable Arbitrary Map, defining a transparent image of the drawn lines. Here is an example of making a complete PNG image from such a PAM (using CHICKEN Scheme):

$ csc -O5 -R r7rs -X r7rs xiaolin_wu_line_algorithm.scm
$ ./xiaolin_wu_line_algorithm image.pam
$ pamgradient lightgray lightblue lightcyan lightgray 1000 750 | pamcomp image.pam - | pamtopng > image.png

Here is what I get:

A violet starburst on a light bluish gradient background.


I thought it would be amusing to modify the code to use the R7RS-small macro system, and so I made the following second version. Note also the variable steep is gone.

;;;-------------------------------------------------------------------

(import (scheme base))
(import (scheme file))
(import (scheme inexact))
(import (scheme process-context))
(import (scheme write))

;; (srfi 160 f32) is more properly known as (scheme vector f32), but
;; is not part of R7RS-small. The following will work in both Gauche
;; and CHICKEN Schemes.
(import (srfi 160 f32))

;;;-------------------------------------------------------------------

(define-record-type <color>
  (make-color r g b)
  color?
  (r color-r)
  (g color-g)
  (b color-b))

;;; See https://yeun.github.io/open-color/
(define violet9 (make-color (/ #x5F 255.0)
                            (/ #x3D 255.0)
                            (/ #xC4 255.0)))

;;;-------------------------------------------------------------------

(define-record-type <drawing-surface>
  (drawing-surface% u0 v0 u1 v1 pixels)
  drawing-surface?
  (u0 u0%)
  (v0 v0%)
  (u1 u1%)
  (v1 v1%)
  (pixels pixels%))

(define (make-drawing-surface u0 v0 u1 v1)
  (unless (and (<= u0 u1) (<= v0 v1))
    (error "illegal drawing-surface corners"))
  (let ((width (- u1 u0 -1))
        (height (- v1 v0 -1)))
    (let ((pixels (make-f32vector (* width height) 0.0)))
      (drawing-surface% u0 v0 u1 v1 pixels))))

;;; In calls to drawing-surface-ref and drawing-surface-set! indices
;;; outside the drawing_surface are allowed. Such indices are treated
;;; as if you were trying to draw on the air.

(define (drawing-surface-ref s x y)
  (let ((u0 (u0% s))
        (v0 (v0% s))
        (u1 (u1% s))
        (v1 (v1% s)))
    (if (and (<= u0 x) (<= x u1) (<= v0 y) (<= y v1))
        (f32vector-ref (pixels% s)
                       (+ (* (- x u0) (- v1 v0 -1)) (- v1 y)))
        +nan.0)))

(define (drawing-surface-set! s x y opacity)
  (let ((u0 (u0% s))
        (v0 (v0% s))
        (u1 (u1% s))
        (v1 (v1% s)))
    (when (and (<= u0 x) (<= x u1) (<= v0 y) (<= y v1))
      (f32vector-set! (pixels% s)
                      (+ (* (- x u0) (- v1 v0 -1)) (- v1 y))
                      opacity))))

(define (write-PAM s color)
  ;; Write a Portable Arbitrary Map to the current output port, using
  ;; the given color as the foreground color and the drawing-surface
  ;; values as opacities.

  (define (float->byte v) (exact (round (* 255 v))))

  (define r (float->byte (color-r color)))
  (define g (float->byte (color-g color)))
  (define b (float->byte (color-b color)))

  (define w (- (u1% s) (u0% s) -1))
  (define h (- (v1% s) (v0% s) -1))
  (define opacities (pixels% s))

  (define (loop x y)
    (cond ((= y h) )
          ((= x w) (loop 0 (+ y 1)))
          (else
           (let ((alpha (float->byte
                         (f32vector-ref opacities (+ (* x h) y)))))
             (write-bytevector (bytevector r g b alpha))
             (loop (+ x 1) y)))))

  (display "P7") (newline)
  (display "WIDTH ") (display (- (u1% s) (u0% s) -1)) (newline)
  (display "HEIGHT ") (display (- (v1% s) (v0% s) -1)) (newline)
  (display "DEPTH 4") (newline)
  (display "MAXVAL 255") (newline)
  (display "TUPLTYPE RGB_ALPHA") (newline)
  (display "ENDHDR") (newline)
  (loop 0 0))

;;;-------------------------------------------------------------------

(define-syntax ipart
  (syntax-rules ()
    ((_ x) (exact (floor x)))))

(define-syntax iround
  (syntax-rules ()
    ((_ x) (ipart (+ x 0.5)))))

(define-syntax fpart
  (syntax-rules ()
    ((_ x) (- x (floor x)))))

(define-syntax rfpart
  (syntax-rules ()
    ((_ x) (- 1.0 (fpart x)))))

(define-syntax plot-shallow
  (syntax-rules ()
    ((_ s x y opacity)
     ;; One might prefer a more sophisticated function than mere
     ;; addition. Here, however, the function is addition.
     (let ((combined-opacity (+ opacity (drawing-surface-ref s x y))))
       (drawing-surface-set! s x y (min combined-opacity 1.0))))))

(define-syntax plot-steep
  (syntax-rules ()
    ((_ s x y opacity)
     (plot-shallow s y x opacity))))

(define-syntax drawln%
  (syntax-rules ()
    ((_ s x0 y0 x1 y1 plot)
     (let* ((dx (- x1 x0))
            (dy (- y1 y0))
            (gradient (if (zero? dx) 1.0 (/ dy dx)))

            ;; Handle the first endpoint.
            (xend (iround x0))
            (yend (+ y0 (* gradient (- xend x0))))
            (xgap (rfpart (+ x0 0.5)))
            (xpxl1 xend)
            (ypxl1 (ipart yend))
            (_ (plot s xpxl1 ypxl1 (* (rfpart yend) xgap)))
            (_ (plot s xpxl1 (+ ypxl1 1) (* (fpart yend) xgap)))

            ;; The first y-intersection.
            (intery (+ yend gradient))

            ;; Handle the second endpoint.
            (xend (iround x1))
            (yend (+ y1 (* gradient (- xend x1))))
            (xgap (fpart (+ x1 0.5)))
            (xpxl2 xend)
            (ypxl2 (ipart yend))
            (_ (plot s xpxl2 ypxl2 (* (rfpart yend) xgap)))
            (_ (plot s xpxl2 (+ ypxl2 1) (* (fpart yend) xgap))))

       ;; Loop over the rest of the points.
       (do ((x (+ xpxl1 1) (+ x 1))
            (intery intery (+ intery gradient)))
           ((= x xpxl2))
         (plot s x (ipart intery) (rfpart intery))
         (plot s x (+ (ipart intery) 1) (fpart intery)))))))

(define (draw-line s x0 y0 x1 y1)
  (let ((xdiff (abs (- x1 x0)))
        (ydiff (abs (- y1 y0))))
    (if (<= ydiff xdiff)
        (if (<= x0 x1)
            (drawln% s x0 y0 x1 y1 plot-shallow)