Xiaolin Wu's line algorithm: Difference between revisions
m
→{{header|Wren}}: Changed to Wren S/H
m (→{{header|Wren}}: Changed to Wren S/H) |
|||
| (17 intermediate revisions by 3 users not shown) | |||
Line 994:
var x : int
in
for (x := succ xpxl1; x <>
begin
plot (s, ipart intery, x, rfpart intery);
Line 1,005:
var x : int
in
for (x := succ xpxl1; x <>
begin
plot (s, x, ipart intery, rfpart intery);
Line 1,240:
LINE X%*2, Y%*2, X%*2, Y%*2
ENDPROC</syntaxhighlight>
=={{header|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.
<syntaxhighlight lang="c">// 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);
}
</syntaxhighlight>
=={{header|C}}==
Line 1,547 ⟶ 1,703:
}
</syntaxhighlight>
=={{header|Common Lisp}}==
{{trans|Scheme}}
The program outputs a transparency mask in Portable Gray Map format. It draws lines normal to a tractrix. They trace out a catenary.
<syntaxhighlight lang="lisp">
;;; 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)
;;;-------------------------------------------------------------------
</syntaxhighlight>
Here is a script to run the program and produce a PNG image.
<syntaxhighlight lang="sh">
#!/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
</syntaxhighlight>
{{out}}
[[File:Xiaolin wu line algorithm CommonLisp.png|thumb|none|alt=A catenary as the evolute of a tractrix. Reddish gradient on bluish gradient.]]
=={{header|D}}==
Line 1,660 ⟶ 1,999:
im2.savePPM6("xiaolin_lines2.ppm");
}</syntaxhighlight>
=={{header|Fortran}}==
{{trans|Common Lisp}}
The program outputs a Portable Gray Map representing a transparency mask. The mask is full of straight lines, drawn by a solution of this Rosetta Code task. Some of the lines draw a piecewise approximation of an ellipse, and others draw the normals of the ellipse. The normals form an envelope that is the evolute of the ellipse.
<syntaxhighlight lang="fortran">
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
</syntaxhighlight>
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.
<syntaxhighlight lang="shell">
#!/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
</syntaxhighlight>
{{out}}
[[File:Xiaolin wu line algorithm Fortran.png|thumb|none|alt=An ellipse and its normals (forming an evolute as their envelope), all of it in goofy colors.]]
=={{header|FreeBASIC}}==
Line 2,718 ⟶ 3,352:
image = ConstantImage[Black, {100, 100}];
Fold[DrawLine[#1, {20, 10}, #2] &, image, AngleVector[{20, 10}, {75, #}] & /@ Subdivide[0, Pi/2, 10]]</syntaxhighlight>
=={{header|MATLAB}}==
{{trans|Julia}}
<syntaxhighlight lang="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
</syntaxhighlight>
=={{header|Modula-2}}==
{{trans|Fortran}}
The program outputs a transparency map in Portable Gray Map format. It draws lines normal to a parabola. The envelope formed is a semicubic parabola.
<syntaxhighlight lang="Modula2">
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.
</syntaxhighlight>
Here is a shell script that compiles the program, runs it, and (using Netpbm commands) makes a PNG using the outputted mask.
<syntaxhighlight lang="sh">
#!/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
</syntaxhighlight>
{{out}}
[[File:Xiaolin wu line algorithm Modula2.png|thumb|none|alt=A parabola, and lines normal to the parabola, forming a semicubic parabola as their envelope. A shade of yellow on a background shade of purple.]]
=={{header|Nim}}==
Line 4,158 ⟶ 5,144:
(error "illegal drawing-surface corners"))
(let ((width (- u1 u0 -1))
(let ((pixels (make-f32vector (* width height) 0.0)))
(drawing-surface% u0 v0 u1 v1 pixels))))
Line 4,167 ⟶ 5,153:
(define (drawing-surface-ref s x y)
(
(
(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)
(
(
(when (and (<= u0 x) (<= x u1) (<= v0 y) (<= y
(f32vector-set! (pixels% s)
(+ (* (- x u0) (- v1 v0 -1)) (- v1 y))
opacity))))
(define (write-PAM s color)
Line 4,266 ⟶ 5,258:
(do ((x (+ xpxl1 1) (+ x 1))
(intery intery (+ intery gradient)))
((= x
(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
(plot s x (ipart intery) (rfpart intery))
(plot s x (+ (ipart intery) 1) (fpart intery))))))
Line 4,348 ⟶ 5,340:
Here is what I get:
[[File:Xiaolin wu line algorithm SCM.2023.04.28.14.08.55.png|thumb|none|alt=A violet starburst on a light bluish gradient background.]]
I thought it would be amusing to modify the code to use the R7RS-small macro system, and so I made the following second version. Note also the variable <code>steep</code> is gone.
<syntaxhighlight lang="scheme">
;;;-------------------------------------------------------------------
(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)
(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)))))
;;;-------------------------------------------------------------------
(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))))
;;;-------------------------------------------------------------------
</syntaxhighlight>
=={{header|Sidef}}==
Line 4,624 ⟶ 5,848:
{{trans|Kotlin}}
{{libheader|DOME}}
<syntaxhighlight lang="
import "dome" for Window
import "math" for Math
| |||