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
#-----------------------------------------------------------------------------
im <- magick::image_read("../../static/img/dithering/Portal_Companion_Cube.jpg") %>%
  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
)