lunedì 29 agosto 2011

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)
}


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)
}


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)
}


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)
}


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)
}


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)
}


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)
}


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)
}


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



venerdì 26 agosto 2011

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)