270 lines
7.5 KiB
Fortran
270 lines
7.5 KiB
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
|