Bibliotecas

Para a leitura e visualização da imagem irei utilizar o pacote imager. Para a criação dos histogramas irei utilizar a biblioteca de visualização de dados ggplot2. Ademais, os pacotes dplyr e tidyr serão usados para manipulação de dados no final do experimento, juntamente da biblioteca de tabelas kableExtra.

library(imager)
library(ggplot2)
library(dplyr)
library(tidyr)
library(kableExtra)

Leitura da Imagem

Vamos ler a imagem através da função imager::load.image:

flower <- load.image('../images/flower.jpg')

plot(flower)

Imagem original

A imagem utilizada no experimento possui as seguintes dimensões: 264 \(\times\) 398. Um aspecto interessante do objeto criado pela função imager::load.image é que a forma como a imagem é representada segue o padrão CImg, i.e, um array 4D que consiste das seguintes informações:

\[ \begin{bmatrix}x~~(\text{largura}) & y~~(\text{comprimento}) & z~~(\text{tempo}) & c~~(\text{canal})\end{bmatrix} \]

o array 4D que representa a imagem tem dimensão:

dim(flower)

## [1] 264 398   1   3

Por se tratar de uma imagem estática, temos apenas 1 frame (\(z\)). Ademais, por padrão, o esquema de cores do objeto criado é o RGB, dessa forma temos 3 canais (\(c\)) de cores.

Extração da Componente RED

Como o esquema de cores da figura é o RGB, então para extrair o canal RED da imagem precisamos da primeira componente da quarta dimensão (\(c\)):

flower_red <- flower[,,,1]

Binarização

A função de binarização da imagem é bem simples, basicamente criamos uma matriz lógica, da mesma dimensão da imagem, que contem TRUE se a intensidade do pixel for menor ou igual aoo valor de treshold (limiar) e FALSE caso contrário. Após a definição dessa matriz lógica, efetuei uma conversão para uma matriz numérica, simplesmente subtraindo o valor 1 da matriz lógica, já que o R interpreta TRUE ou FALSE também como 1 e 0. Dessa forma, a função retorna uma matriz de mesma dimensão que contém 0 se a intensidade do pixel for menor ou igual ao limiar, e 1 caso contrário.

binarize <- function(image, threshold) {
  img_ix <- image <= threshold 
  return(1 - img_ix)
}

Para binarizar a imagem da componente RED extraida anteriormente, precisamos definir o threshold pela qual a função discrimina a intensidade dos pixels.

Vamos avaliar o comportamento da intensidade dos pixels através de seu respectivo histograma:

flower_red %>% as.vector() %>%
  qplot(xlab = 'Intensidade do Pixel', ylab = 'Frequência', 
        fill = I("grey"), col = I("black"))  +
  scale_y_continuous(labels = function(x) {
    format(x, big.mark = ".", decimal.mark = ",", scientific = FALSE)}) +
  scale_x_continuous(breaks = seq(0, 1, 0.1)) + theme_bw()

Histograma da intensidade dos pixels no canal RED

Aqui podemos notar um comportamento conhecido como zero inflated, i.e, um acumulo significativo de zeros na freuquência dos dados. Esse padrão se deve ao fato de a imagem conter em sua maioria pixeis pretos (intensidade zero).

Uma forma de amenizar esse comportamento na visualização é plotar apenas as intensidades positivas (intensidade > 0):

flower_red %>% as.vector() %>% .[. > 0] %>%
  qplot(xlab = 'Intensidade do Pixel', ylab = 'Frequência', 
        fill = I("grey"), col = I("black"))  +
  scale_y_continuous(labels = function(x) {
    format(x, big.mark = ".", decimal.mark = ",", scientific = FALSE)}) + 
  scale_x_continuous(breaks = seq(0, 1, 0.1)) +
  geom_vline(xintercept = 0.12, color = 'red', linetype = 'dashed') + theme_bw()

Histograma da intensidade positiva dos pixels no canal RED
Ponto crítico \(\approx\) 0,12 em destaque

Perceba que um ponto crítico no histograma está situado na região de intensidade \(\approx\) 0.12, sugerindo que esse ponto possa ser uma boa escolha para o limiar da função de binarização.

Para fins de simplicidade utilizei o valor de treshold \(=\) 0.12 na função de binarização:

flower_red_bin <- binarize(flower_red, 0.12)

plot(flower)
plot(flower_red_bin %>% as.cimg)

Imagem original (à esquerda) e imagem binarizada (à direita)

Aparentemente a escolha do valor do treshold foi satisfatória, já que pela comparação das figuras não é possivel notar uma perda “abrupta” de informação da área alvo.

Adjacência

Nessa seção vamos implementar quatro métodos de adjacência: N4, ND, N8 e M.

Para os métodos N4 e ND, as funções tem etapas bastante parecidas:

1. definir a adjacência ao ponto inicial dado;

2. definir uma matriz de adjacência da forma: [x, y, intensidade];

3. obter o subconjunto de intensidade = 1;

4. retornar a matriz resultante.

adjacency_N4 <- function(img, locX, locY) {
  N4 <- matrix(c(c(locX-1, locY, img[locX-1, locY]), c(locX, locY-1, img[locX, locY-1]),
                 c(locX, locY+1, img[locX, locY+1]), c(locX+1, locY), img[locX+1, locY]), 
               ncol = 3, byrow = T)
  N4 <- N4[N4[,3] == 1,] 
  
  if(length(N4) == 0) return(NULL)
  else if(is.null(dim(N4))) return(t(N4))
  else return(N4)
}
adjacency_ND <- function(img, locX, locY) {
  ND <- matrix(c(c(locX-1, locY-1, img[locX-1, locY-1]), c(locX+1, locY-1, img[locX+1, locY-1]),
                 c(locX-1, locY+1, img[locX-1, locY+1]), c(locX+1, locY+1), img[locX+1, locY+1]), 
               ncol = 3, byrow = T)
  ND <- ND[ND[,3] == 1,] 
  
  if(length(ND) == 0) return(NULL)
  else if(is.null(dim(ND))) return(t(ND))
  else return(ND)
}

Para o método N8, utilizamos a definição N8 = N4 \(\cup\) ND.

adjacency_N8 <- function(img, locX, locY) {
  N8 <- rbind(adjacency_N4(img, locX, locY), adjacency_ND(img, locX, locY))
  
  return(N8)
}

Finalmente, o método de adjacência M é bastante parecido com o método N8, porém temos uma pequena alteração:

1. definir a adjacência ao ponto inicial dado (N4 e ND);

2. obter o subconjunto de ND que não tenha pontos com coordenadas \(x\) ou \(y\) já presentes em N4;

3. definir uma matriz de adjacência da forma: [x, y, intensidade];

4. obter o subconjunto de intensidade = 1;

5. retornar a matriz resultante.

adjacency_M <- function(img, locX, locY) {
  N4 <- adjacency_N4(img, locX, locY)
  ND <- adjacency_ND(img, locX, locY)
  
  if(is.null(N4) & is.null(ND)) return(NULL)
  if(is.null(N4)) return(ND)
  if(nrow(N4) > 2 | is.null(ND)) return(N4)
  
  ND_M <- !(ND[, 1] %in% N4[, 1] | ND[, 2] %in% N4[, 2])
  M <- rbind(N4, ND[ND_M,])
  
  return(M)
}

Para evitar problemas de busca fora dos limites da imagem pela função de adjacência, defini uma função de expansão da imagem, i.e, adição de margens com intensidade 0 à imagem.

augment_img <- function(img, cells) {
  img_zeros <- matrix(0, nrow = nrow(img) + 2*cells, ncol = ncol(img) + 2*cells)
  img_zeros[1:nrow(img) + cells, 1:ncol(img) + cells] <- img
  return(img_zeros)
}

A função de adjcência é então formada por um laço do tipo while que percorre cada ponto adjacente ao ponto de referência, armazenando suas respectivas coordenadas [\(x\), \(y\)] na matriz adj_mx. O laço é interrompido caso não haja mais pontos adjacentes na matriz adj_mx a serem explorados, i.e, não há mais pontos de conexão disponíveis.

adjacency <- function(img, locX, locY, type, aug_cells = 2) {
  img <- augment_img(img, aug_cells)
  
  adj_fn <- list('N4' = adjacency_N4, 'ND' = adjacency_ND, 'N8' = adjacency_N8, 'M' = adjacency_M)
  adj_mx <- array(c(locX, locY), dim = c(1, 2), dimnames = list(NULL, c('x','y')))
  adj_fn <- adj_fn[[type]]
  
  i <- 1
  while(i <= nrow(adj_mx)) {
    adj_arr <- adj_fn(img, adj_mx[i, 1], adj_mx[i, 2])[, -3]
    adj_mx = unique(rbind(adj_mx, adj_fn(img, adj_mx[i, 1], adj_mx[i, 2])[, -3]))
    i = i + 1
  }
  
  return(adj_mx - aug_cells)
}

Para o ponto inicial, irei utilizar a coordenada 100, 130 como referência para o início da busca do algoritmo.

plot(flower_red_bin %>% as.cimg)
points(x = 100, y = 130, col = 'red')

Imagem binarizada
Coordenada [100, 130] em destaque

Por fim, aplicamos as funções de adjacência na imagem binarizada e obtemos a matriz de coordenadas conectadas ao ponto [100, 130]:

adj_N4 <- adjacency(flower_red_bin, locX = 100, locY = 130, type = 'N4')
adj_ND <- adjacency(flower_red_bin, locX = 100, locY = 130, type = 'ND')
adj_N8 <- adjacency(flower_red_bin, locX = 100, locY = 130, type = 'N8')
adj_M <- adjacency(flower_red_bin, locX = 100, locY = 130, type = 'M')

par(mfrow = c(1, 4))
flower_red_bin_N4 <- adj_N4 %>%
  as_tibble() %>%
  mutate(value = 1) %>%
  complete(x = 1:264, y = 1:398, fill = list(value = 0)) %>%
  as.cimg(dims = c(264, 398, 1, 1))
flower_red_bin_ND <- adj_ND %>%
  as_tibble() %>%
  mutate(value = 1) %>%
  complete(x = 1:264, y = 1:398, fill = list(value = 0)) %>%
  as.cimg(dims = c(264, 398, 1, 1))
flower_red_bin_N8 <- adj_N8 %>%
  as_tibble() %>%
  mutate(value = 1) %>%
  complete(x = 1:264, y = 1:398, fill = list(value = 0)) %>%
  as.cimg(dims = c(264, 398, 1, 1))
flower_red_bin_M <- adj_M %>%
  as_tibble() %>%
  mutate(value = 1) %>%
  complete(x = 1:264, y = 1:398, fill = list(value = 0)) %>%
  as.cimg(dims = c(264, 398, 1, 1))

plot(flower_red_bin_N4)
plot(flower_red_bin_ND)
plot(flower_red_bin_N8)
plot(flower_red_bin_M)

Área conectada ao ponto [100, 130] pelos métodos de adjacência
N4 (figura 1), ND (figura 2), N8 (figura 3) e M (figura 4)

Dimensões da Imagem

Vamos tentar estimar as dimensões da imagem através de três medidas de distância: euclidiana, manhattan (cityblock) e chebyshev (chessboard). Contudo, precisamos obter as coordenadas extremas de cada imagem obtida pela aplicação dos métodos de adjacência:

minh_N4 <- adj_N4[which.min(adj_N4[,1]),]
maxh_N4 <- adj_N4[which.max(adj_N4[,1]),]
minw_N4 <- adj_N4[which.min(adj_N4[,2]),]
maxw_N4 <- adj_N4[which.max(adj_N4[,2]),]
minh_ND <- adj_ND[which.min(adj_ND[,1]),]
maxh_ND <- adj_ND[which.max(adj_ND[,1]),]
minw_ND <- adj_ND[which.min(adj_ND[,2]),]
maxw_ND <- adj_ND[which.max(adj_ND[,2]),]
minh_N8 <- adj_N8[which.min(adj_N8[,1]),]
maxh_N8 <- adj_N8[which.max(adj_N8[,1]),]
minw_N8 <- adj_N8[which.min(adj_N8[,2]),]
maxw_N8 <- adj_N8[which.max(adj_N8[,2]),]
minh_M <- adj_M[which.min(adj_M[,1]),]
maxh_M <- adj_M[which.max(adj_M[,1]),]
minw_M <- adj_M[which.min(adj_M[,2]),]
maxw_M <- adj_M[which.max(adj_M[,2]),]

Vamos verificar os pontos extremos de cada imagem:

par(mfrow = c(1, 4))

plot(flower_red_bin_N4)
points(x = minh_N4[1], y = minh_N4[2], col = 'red')
points(x = maxh_N4[1], y = maxh_N4[2], col = 'red')
points(x = minw_N4[1], y = minw_N4[2], col = 'red')
points(x = maxw_N4[1], y = maxw_N4[2], col = 'red')

plot(flower_red_bin_ND)
points(x = minh_ND[1], y = minh_ND[2], col = 'red')
points(x = maxh_ND[1], y = maxh_ND[2], col = 'red')
points(x = minw_ND[1], y = minw_ND[2], col = 'red')
points(x = maxw_ND[1], y = maxw_ND[2], col = 'red')

plot(flower_red_bin_N8)
points(x = minh_N8[1], y = minh_N8[2], col = 'red')
points(x = maxh_N8[1], y = maxh_N8[2], col = 'red')
points(x = minw_N8[1], y = minw_N8[2], col = 'red')
points(x = maxw_N8[1], y = maxw_N8[2], col = 'red')

plot(flower_red_bin_M)
points(x = minh_M[1], y = minh_M[2], col = 'red')
points(x = maxh_M[1], y = maxh_M[2], col = 'red')
points(x = minw_M[1], y = minw_M[2], col = 'red')
points(x = maxw_M[1], y = maxw_M[2], col = 'red')

Pontos extremos nas imagens obtidas pelos métodos de adjacência
N4 (figura 1), ND (figura 2), N8 (figura 3) e M (figura 4)

Apesar das similaridades entre os pontos extremos da imagem, vamos estimar as dimensões da imagem nas quatro alternativas acima.

Vamos definir as funções de distâncias:

euclidean <- function(x, y) {
  return(sqrt((x[1] - y[1])^2 + (x[2] - y[2])^2) %>% round)
}
cityblock <- function(x, y) {
  return(abs(x[1] - y[1]) + abs(x[2] - y[2]) %>% round)
}
chessboard <- function(x, y) {
  return(max(abs(x[1] - y[1]), abs(x[2] - y[2])) %>% round)
}

…sem mistérios até aqui.

A estimativa das dimensões da imagem é dada pela distância, dado um método de distância, de dois pontos extremos (mínimo e máximo). O código abaixo cria uma matriz de dados com as dimensões estimadas dos quatro métodos de adjacência utilizando as três funções de distância.

dimensions <- data.frame(
  distance = c('euclidean', 'cityblock', 'chessboard'),
  N4 = c(paste0('[', euclidean(minw_N4, maxw_N4), ', ', euclidean(minh_N4, maxh_N4), ']'),
         paste0('[', cityblock(minw_N4, maxw_N4), ', ', cityblock(minh_N4, maxh_N4), ']'),
         paste0('[', chessboard(minw_N4, maxw_N4), ', ', chessboard(minh_N4, maxh_N4), ']')),
  ND = c(paste0('[', euclidean(minw_ND, maxw_ND), ', ', euclidean(minh_ND, maxh_ND), ']'),
         paste0('[', cityblock(minw_ND, maxw_ND), ', ', cityblock(minh_ND, maxh_ND), ']'),
         paste0('[', chessboard(minw_ND, maxw_ND), ', ', chessboard(minh_ND, maxh_ND), ']')),
  N8 = c(paste0('[', euclidean(minw_N8, maxw_N8), ', ', euclidean(minh_N8, maxh_N8), ']'),
         paste0('[', cityblock(minw_N8, maxw_N8), ', ', cityblock(minh_N8, maxh_N8), ']'),
         paste0('[', chessboard(minw_N8, maxw_N8), ', ', chessboard(minh_N8, maxh_N8), ']')),
  M = c(paste0('[', euclidean(minw_M, maxw_M), ', ', euclidean(minh_M, maxh_M), ']'),
         paste0('[', cityblock(minw_M, maxw_M), ', ', cityblock(minh_M, maxh_M), ']'),
         paste0('[', chessboard(minw_M, maxw_M), ', ', chessboard(minh_M, maxh_M), ']'))
  
)

Finalmente, as dimensões estimadas da imagem são mostradas abaixo:

dimensions %>%
  tibble::column_to_rownames('distance') %>%
  kbl(row.names = T, align = 'cccc', 
      caption = 'Dimensões estimadas da região alvo na imagem') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Bom, é isso…

Até a proxima 👋.