R is a cool image editor #2: Dithering algorithms

Ho implementato in R alcuni algoritmi per il "dithering" delle immagini (un buon articolo introduttivo è disponibile su WikiPedia):
- Floyd-Steinberg dithering
- Bill Atkinson dithering
- Jarvis-Judice-Ninke dithering
- Sierra 2-4a dithering filter
- Stucki dithering
- Burkes dithering
- Sierra2 dithering
- Sierra3 dithering

Per ciascun algoritmo, sono presenti due funzioni. La prima consiste nel processo di convolution, mentre la seconda applica la prima funzione all'immagine caricata.
Per usare le immagini, ho usato la library rimage (per un breve commento, potete leggere il mio articolo precedente qui). Questi algoritmi sono alquanto lenti, perchè non ho implementato alcun metodo per velocizzare il convolution; convertendo il codice in linguaggio C, da usare in R, il processo diventa nettamente più rapido.
Inoltre (per questione di rapidità) queste funzioni possono essere applicate solo su immagini in scala di grigi, ma possono essere facilmente ampliate per immagini in RGB.
La prima funzione è commentata; le altre sono molto simili.

library(rimage)

y <- read.jpeg("valve.jpg")

plot(y)

---

Floyd-Steinberg dithering

	# Floyd-Steinberg dithering
	# the 2-dimensional convolution function
	FloydConvolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 1:(dim(a)[1]-1)){
	for(j in 1:(dim(a)[2]-1)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 7/16)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 3/16)
	a[i+1,j] <- a[i+1,j] + (e * 5/16)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 1/16)
	}
	}
	a
	}
	# the main function
	grey2FSdith <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 1
	dim2 <- 1
	dim1x <- 1
	dim2x <- 1
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	# this code creates a bigger matrix (image) adding 0.5 values in first/last row/col
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2, ncol=ncol(greyMatrix)+2)
	tempMatrix[2:(nrow(tempMatrix)-1),2:(ncol(tempMatrix)-1)] <- greyMatrix
	# apply the convolution function
	igrey <- FloydConvolution(tempMatrix)
	# create the image
	imagematrix(igrey[1:(nrow(igrey)-1),1:(ncol(igrey)-1)], type="grey", noclipping=TRUE)
	}

view raw Floyd-Steinberg dithering hosted with ❤ by GitHub

plot(normalize(grey2FSdith(rgb2grey(y))))

---

Bill Atkinson dithering

	# Bill Atkinson dithering
	AtkinConvolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 1/8)
	a[i,j+2] <- a[i,j+2] + (e * 1/8)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 1/8)
	a[i+1,j] <- a[i+1,j] + (e * 1/8)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 1/8)
	a[i+2,j] <- a[i+2,j] + (e * 1/8)
	}
	}
	a
	}
	grey2ATKdith <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*2, ncol=ncol(greyMatrix)+2*2)
	tempMatrix[3:(nrow(tempMatrix)-2),3:(ncol(tempMatrix)-2)] <- greyMatrix
	igrey <- AtkinConvolution(tempMatrix)
	imagematrix(igrey[3:(nrow(igrey)-2),3:(ncol(igrey)-2)], type="grey", noclipping=TRUE)
	}

view raw Bill Atkinson dithering hosted with ❤ by GitHub

plot(normalize(grey2ATKdith(rgb2grey(y))))

---

Jarvis-Judice-Ninke dithering

	# Jarvis-Judice-Ninke dithering
	JJNConvolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 7/48)
	a[i,j+2] <- a[i,j+2] + (e * 5/48)
	a[i+1,j-2] <- a[i+1,j-2] + (e * 3/48)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 5/48)
	a[i+1,j] <- a[i+1,j] + (e * 7/48)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 5/48)
	a[i+1,j+2] <- a[i+1,j+2] + (e * 3/48)
	a[i+2,j-2] <- a[i+2,j-2] + (e * 1/48)
	a[i+2,j-1] <- a[i+2,j-1] + (e * 3/48)
	a[i+2,j] <- a[i+2,j] + (e * 5/48)
	a[i+2,j+1] <- a[i+2,j+1] + (e * 3/48)
	a[i+2,j+2] <- a[i+2,j+2] + (e * 1/48)
	}
	}
	a
	}
	grey2JJNdith <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*2, ncol=ncol(greyMatrix)+2*2)
	tempMatrix[3:(nrow(tempMatrix)-2),3:(ncol(tempMatrix)-2)] <- greyMatrix
	igrey <- JJNConvolution(tempMatrix)
	imagematrix(igrey[3:(nrow(igrey)-2),3:(ncol(igrey)-2)], type="grey", noclipping=TRUE)
	}

view raw Jarvis-Judice-Ninke dithering hosted with ❤ by GitHub

plot(normalize(grey2JJNdith(rgb2grey(y))))

---

Sierra 2-4a dithering filter

	# Sierra 2-4a dithering filter
	Sierra24aConvolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 2/4)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 1/4)
	a[i+1,j] <- a[i+1,j] + (e * 1/4)
	}
	}
	a
	}
	grey2S24adith <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*1, ncol=ncol(greyMatrix)+2*1)
	tempMatrix[2:(nrow(tempMatrix)-1),2:(ncol(tempMatrix)-1)] <- greyMatrix
	igrey <- Sierra24aConvolution(tempMatrix)
	imagematrix(igrey[2:(nrow(igrey)-1),2:(ncol(igrey)-1)], type="grey", noclipping=TRUE)
	}

view raw Sierra 2-4a dithering hosted with ❤ by GitHub

plot(normalize(grey2S24adith(rgb2grey(y))))

---

Stucki dithering

	# Stucki dithering
	StuckiConvolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 8/42)
	a[i,j+2] <- a[i,j+2] + (e * 4/42)
	a[i+1,j-2] <- a[i+1,j-2] + (e * 2/42)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 4/42)
	a[i+1,j] <- a[i+1,j] + (e * 8/42)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 4/42)
	a[i+1,j+2] <- a[i+1,j+2] + (e * 2/42)
	a[i+2,j-2] <- a[i+2,j-2] + (e * 1/42)
	a[i+2,j-1] <- a[i+2,j-1] + (e * 2/42)
	a[i+2,j] <- a[i+2,j] + (e * 4/42)
	a[i+2,j+1] <- a[i+2,j+1] + (e * 2/42)
	a[i+2,j+2] <- a[i+2,j+2] + (e * 1/42)
	}
	}
	a
	}
	grey2Stucki <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*2, ncol=ncol(greyMatrix)+2*2)
	tempMatrix[3:(nrow(tempMatrix)-2),3:(ncol(tempMatrix)-2)] <- greyMatrix
	igrey <- StuckiConvolution(tempMatrix)
	imagematrix(igrey[3:(nrow(igrey)-2),3:(ncol(igrey)-2)], type="grey", noclipping=TRUE)
	}

view raw Stucki dithering hosted with ❤ by GitHub

plot(normalize(grey2Stucki(rgb2grey(y))))

---

Burkes dithering

	# Burkes dithering
	BurkesConvolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 4/32)
	a[i,j+2] <- a[i,j+2] + (e * 4/32)
	a[i+1,j-2] <- a[i+1,j-2] + (e * 2/32)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 4/32)
	a[i+1,j] <- a[i+1,j] + (e * 8/32)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 4/32)
	a[i+1,j+2] <- a[i+1,j+2] + (e * 2/32)
	}
	}
	a
	}
	grey2Burkes <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*2, ncol=ncol(greyMatrix)+2*2)
	tempMatrix[3:(nrow(tempMatrix)-2),3:(ncol(tempMatrix)-2)] <- greyMatrix
	igrey <- BurkesConvolution(tempMatrix)
	imagematrix(igrey[3:(nrow(igrey)-2),3:(ncol(igrey)-2)], type="grey", noclipping=TRUE)
	}

view raw Burkes dithering hosted with ❤ by GitHub

plot(normalize(grey2Burkes(rgb2grey(y))))

---

Sierra2 dithering

	# Sierra2 dithering
	Sierra2Convolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 4/16)
	a[i,j+2] <- a[i,j+2] + (e * 3/16)
	a[i+1,j-2] <- a[i+1,j-2] + (e * 1/16)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 2/16)
	a[i+1,j] <- a[i+1,j] + (e * 3/16)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 2/16)
	a[i+1,j+2] <- a[i+1,j+2] + (e * 1/16)
	}
	}
	a
	}
	grey2Sierra2 <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*2, ncol=ncol(greyMatrix)+2*2)
	tempMatrix[3:(nrow(tempMatrix)-2),3:(ncol(tempMatrix)-2)] <- greyMatrix
	igrey <- Sierra2Convolution(tempMatrix)
	imagematrix(igrey[3:(nrow(igrey)-2),3:(ncol(igrey)-2)], type="grey", noclipping=TRUE)
	}

view raw Sierra2 dithering hosted with ❤ by GitHub

plot(normalize(grey2Sierra2(rgb2grey(y))))

---

Sierra3 dithering

	# Sierra3 dithering
	Sierra3Convolution <- function(a){
	c <- matrix(0, nrow=dim(a)[1], ncol=dim(a)[2])
	for(i in 2:(dim(a)[1]-2)){
	for(j in 2:(dim(a)[2]-2)){
	P <- trunc(a[i,j]+0.5)
	e <- a[i,j] - P
	a[i,j] <- P
	a[i,j+1] <- a[i,j+1] + (e * 5/32)
	a[i,j+2] <- a[i,j+2] + (e * 3/32)
	a[i+1,j-2] <- a[i+1,j-2] + (e * 2/32)
	a[i+1,j-1] <- a[i+1,j-1] + (e * 4/32)
	a[i+1,j] <- a[i+1,j] + (e * 5/32)
	a[i+1,j+1] <- a[i+1,j+1] + (e * 4/32)
	a[i+1,j+2] <- a[i+1,j+2] + (e * 2/32)
	a[i+2,j-1] <- a[i+2,j-1] + (e * 2/32)
	a[i+2,j] <- a[i+2,j] + (e * 3/32)
	a[i+2,j+1] <- a[i+2,j+1] + (e * 2/32)
	}
	}
	a
	}
	grey2Sierra3 <- function(img){
	greyMatrix <- img[1:nrow(img),1:ncol(img)]
	dim1 <- 2
	dim2 <- 2
	dim1x <- 2
	dim2x <- 2
	dim1a <- dim(greyMatrix)[1]
	dim2a <- dim(greyMatrix)[2]
	tempMatrix <- matrix(0.5, nrow=nrow(greyMatrix)+2*2, ncol=ncol(greyMatrix)+2*2)
	tempMatrix[3:(nrow(tempMatrix)-2),3:(ncol(tempMatrix)-2)] <- greyMatrix
	igrey <- Sierra3Convolution(tempMatrix)
	imagematrix(igrey[3:(nrow(igrey)-2),3:(ncol(igrey)-2)], type="grey", noclipping=TRUE)
	}

view raw Sierra3 dithering hosted with ❤ by GitHub

plot(normalize(grey2Sierra3(rgb2grey(y))))

Legge di Benford, o legge della prima cifra

La distribuzione di Benford meglio nota come legge di Benford o legge della prima cifra è una distribuzione di probabilità che descrive la probabilità che un numero presente in molte raccolte di dati reali (p.es. popolazione dei comuni, quotazione delle azioni, costanti fisiche o matematiche, numero di strade esistenti nelle località) cominci con una data cifra, ad esempio "1", che secondo questa variabile casuale discreta dovrebbe essere nel 30,1% dei casi la prima cifra.
Wikipedia, retrieved 08/26/2011





Simulazione con R:

library(MASS)
benford <- function(m, n){
list <- c()

# calcola tutti i valori m^n, per n= 1, 2, ..., i, ..., n
for(i in 1:n){
list[i] <- m^i
}

# funzione per estrarre la prima cifra di un numero
bben <- function(k){
as.numeric(head(strsplit(as.character(k),'')[[1]],n=1))
}

# estrazione della prima cifra dei valori calcolati
first.digit <- sapply(list, bben)

# plot delle frequenze della prima cifra
truehist(first.digit, nbins=10, main=m)
}

par(mfrow=c(2,2))
benford(2,1000)
benford(3,640)
benford(4,500)
benford(5,440)