Error Diffusion Dithering in R

Introduction

I wanted to experiment with some dithering algorithms within R, but I also wanted complete control of the error diffusion.

After reading Tanner Helland’s articale on dithering, I decided to implement some generic dithering code in order to experiment.

Packages

• magick - for image preparation
• EBImage - easiest(?) way to get image data out of `magick`

To Do

• Serpentine error diffusion (i.e. alternating traversal direction every row)

Original Image

``````#-----------------------------------------------------------------------------
# Using imagemagick to load and convert the image to grayscale
#-----------------------------------------------------------------------------
magick::image_convert(type = 'grayscale') %>%
magick::image_scale(geometry="50%")

#-----------------------------------------------------------------------------
# Extract the matrix of values
#-----------------------------------------------------------------------------
m <- magick::as_EBImage(im)@.Data

#-----------------------------------------------------------------------------
# Show the original image
#-----------------------------------------------------------------------------
plot(EBImage::as.Image(m), interpolate = FALSE)``````

Floyd Steinberg dither

For dithering a grayscale image (values 0 to 255) down to a binary image, each pixel is considered in turn such that

• the pixel is set to 0 or 255 - depending which is closer to the current value.
• the difference between this current value and the new set value is the error.
• distribute this error among the surrounding pixels.

There are plenty of good tutorials about the mechanics of dithering, but I thought I’d try my hand at visualising what’s going on:

``````#-----------------------------------------------------------------------------
# Floyd Steinberg error diffusion matrix.
# The position of the NA indicates the position of the current pixel.
# The other values indicate the fraction of the error distributed to that
# pixel
#-----------------------------------------------------------------------------
floyd_steinberg <- matrix(
c(0, NA, 7,
3,  5, 1),
ncol=3, byrow = TRUE) / 16

dithered_matrix <- dither(m, floyd_steinberg)
plot(EBImage::as.Image(dithered_matrix), interpolate = FALSE)``````

Jarvis, Judice, and Ninke Dithering

``````#-----------------------------------------------------------------------------
# Jarvis Judice Ninke
#-----------------------------------------------------------------------------
jarvis_judice_ninke <- matrix(
c(0, 0, NA, 7, 5,
3, 5,  7, 5, 3,
1, 3,  5, 3, 1),
ncol=5, byrow = TRUE) / 48

dithered_matrix <- dither(m, jarvis_judice_ninke)
plot(EBImage::as.Image(dithered_matrix), interpolate = FALSE)``````

Appendix - the `dither()` function

``````#-----------------------------------------------------------------------------
#' Add a border of zeros around a matrix
#'
#' @param m matrix
#' @param border border width
#-----------------------------------------------------------------------------
expand_matrix <- function(m, border) {
em <- matrix(0, nrow = nrow(m) + 2*border, ncol=ncol(m) + 2*border)
em[border + seq(nrow(m)), border + seq(ncol(m))] <- m
em
}

#-----------------------------------------------------------------------------
#' Dither a numeric matrix
#'
#' @param m numeric matrix with all values in the range [0, 1]
#' @param edm the error diffusion matrix
#'
#' @return numeric matrix with dithering. all values are 0 or 1.
#-----------------------------------------------------------------------------
dither <- function(m, edm, threshold = 0.5) {

#---------------------------------------------------------------------------
# Expand the matrix so we don't have to worry about edge conditions
#---------------------------------------------------------------------------
border <- max(dim(edm))
em     <- expand_matrix(m, border = border)

#---------------------------------------------------------------------------
# Set up the error diffusion matrix and associated variables
#---------------------------------------------------------------------------
stopifnot(sum(is.na(edm)) == 1) # should be 1 NA to nominate current pixel position
edm_w           <- ncol(edm)
edm_h           <- nrow(edm)
edm_xoff        <- which(is.na(edm[1, ])) - 1
edm[is.na(edm)] <- 0

#---------------------------------------------------------------------------
# Loop over all pixels in the matrix
#---------------------------------------------------------------------------
for (row in seq(border, nrow(em) - border)) {

#-------------------------------------------------------------------------
# To which part of the matrix is the error diffusion matrix to be applied?
#-------------------------------------------------------------------------
ed_rows  <- row:(row+edm_h - 1)
ed_cols  <- (border - edm_xoff):(border - edm_xoff + edm_w - 1)

for (col in seq(border, ncol(em) - border)) {

#-----------------------------------------------------------------------
# For the current position, push it to zero or 1 and
# calculate the error
#-----------------------------------------------------------------------
val <- em[row, col]
if (val < threshold) {
error <- val
em[row, col] <- 0
} else {
error <- val - 1
em[row, col] <- 1
}

#-----------------------------------------------------------------------
# Distribute the error to the surrounding pixels
#-----------------------------------------------------------------------
em[ed_rows, ed_cols] <-  em[ed_rows, ed_cols] + edm * error

#-----------------------------------------------------------------------
# slide the diffusion matrix position one pixel to the right
#-----------------------------------------------------------------------
ed_cols <- ed_cols + 1
}
}

#---------------------------------------------------------------------------
# Return the subset of the expanded matrix which corresponds to the
# original matrix
#---------------------------------------------------------------------------
em[border + seq(nrow(m)), border + seq(ncol(m))]
}``````

Appendix: Other error diffusion matrices

``````edms <- list(
#-----------------------------------------------------------------------------
# 1-d Error Diffusion (only pushes error forwards to the next pixel)
#-----------------------------------------------------------------------------
one_d = matrix(
c(NA, 1),
ncol=2, byrow = TRUE),

#-----------------------------------------------------------------------------
# Floyd Steinberg
#-----------------------------------------------------------------------------
floyd_steinberg = matrix(
c(0, NA, 7,
3,  5, 1),
ncol=3, byrow = TRUE) / 16,

#-----------------------------------------------------------------------------
# Fake Floyd Steinberg
#-----------------------------------------------------------------------------
fake_floyd_steinberg = matrix(
c(NA, 3,
3 , 2),
ncol=2, byrow = TRUE) / 8,

#-----------------------------------------------------------------------------
# Stucki
#-----------------------------------------------------------------------------
stucki = matrix(
c(0, 0, NA, 8, 4,
2, 4, 8 , 4, 2,
1, 2, 4 , 2, 1),
ncol=5, byrow = TRUE) / 42,

#-----------------------------------------------------------------------------
# Burkes
#-----------------------------------------------------------------------------
burkes = matrix(
c(0, 0, NA, 8, 4,
2, 4,  8, 4, 2),
ncol=5, byrow = TRUE) / 32,

#-----------------------------------------------------------------------------
# Atkinson
#-----------------------------------------------------------------------------
atkinson = matrix(
c(0, NA, 1, 1,
1,  1, 1, 0,
0,  1, 0, 0),
ncol=4, byrow = TRUE) / 8,

#-----------------------------------------------------------------------------
# Jarvis Judice Ninke
#-----------------------------------------------------------------------------
jarvis_judice_ninke = matrix(
c(0, 0, NA, 7, 5,
3, 5,  7, 5, 3,
1, 3,  5, 3, 1),
ncol=5, byrow = TRUE) / 48,

#-----------------------------------------------------------------------------
# Sierra
#-----------------------------------------------------------------------------
sierra = matrix(
c(0, 0, NA, 5, 3,
2, 4,  5, 4, 2,
0, 2,  3, 2, 0),
ncol=5, byrow = TRUE) / 32,

#-----------------------------------------------------------------------------
# Two-row Sierra
#-----------------------------------------------------------------------------
two_row_sierra = matrix(
c(0, 0, NA, 4, 3,
1, 2,  3, 2, 1),
ncol=5, byrow = TRUE) / 16,

#-----------------------------------------------------------------------------
# Sierra Lite
#-----------------------------------------------------------------------------
sierra_lite = matrix(
c(0, NA, 2,
1,  1, 0),
ncol=3, byrow = TRUE) / 4
)``````