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
- See Bresenham's line algorithm for aliased lines.
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
/* 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
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).
AutoHotkey
#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
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:
D
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
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:
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
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.
- Output:
An example GIF:
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
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
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
// 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
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
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:
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:
Pascal
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
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
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
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:
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)