Exercicio #4

Bibliotecas

Para a leitura e visualização da imagem irei utilizar o pacote imager. Ademais, os pacotes dplyr e ggplot2 serão usados para manipulação de dados e visualização das áreas extraidas ao final do experimento.

library(imager)
library(ggplot2)
library(dplyr)

Leitura da Imagem

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

img <- load.image('../images/lung.jpg')

Imagem original

A imagem utilizada no experimento possui as seguintes dimensões: 1024 \(\times\) 1187.

O array 4D (Padrão CImg) que representa a imagem tem dimensão:

dim(img)

## [1] 1024 1187    1    3

Por termos uma imagem RGB, temos três canais de cores… contudo, irei extrair a imagem em escala de cinza utilizando a função imager::grayscale, pelo método LUMA, i.e, uma combinação linear das três componentes RGB.

img_gray <- img %>% grayscale %>% as.matrix

Segmentação dos Pulmões

Para segmentar as áreas dos pulmões irei utilizar kernels suavizantes e de realce. Como é a primeira vez que os utilizo, irei detalhar como os implementei:

Kernels

Gaussiano

O kernel gaussiano é baseado na seguinte equação geradora:

\[ G(s,~t) = K\times \exp\left\{-\frac{s^2 + t^2}{2\sigma^2}\right\} \]

kernel_gaussian

## function (size, sigma = 1, K = 1) 
## {
##     kern <- -(size%/%2):(size%/%2)
##     kern <- sapply(kern, function(k) exp(-(k^2)/(2 * sigma^2)))
##     kern <- outer(kern, kern, "*") * K
##     kern <- (kern)/sum(kern)
##     return(kern)
## }

A implementação gera, primeiramente, uma sequência de inteiros da divisão inteira de -size à size, por exemplo em uma janela \(3\times3\) (size = 3) a sequência seria \(\{-1, 0, 1\}\).

Após, é avaliado a parte exponencial da equação \(G(s,~t)\) em cada elemento da sequência gerada… esse resultado é possivel pela decomposição de \(G(s,~t)\) em:

\[ G(s,~t) = K\times \exp\left\{-\frac{s^2 + t^2}{2\sigma^2}\right\} = K\times\left[\exp\left\{-\frac{s^2}{2\sigma^2}\right\} \times \exp\left\{-\frac{t^2}{2\sigma^2}\right\}\right] \] Depois de calculadas as exponenciais é efetuado um produto vetorial, que é exatamente a parte entre colchetes. O resultado é então multiplicado pela constante \(K\).

Finalmente, para evitar saturação, toda a matriz gerada é padronizada pelo inverso da soma de seus coeficientes: \(\frac{G(s,~t)}{\sum_s \sum_t G(s,~t)}\).

Exemplo de kernel gaussiano com size = 3 (Janela 3x3):

kernel_gaussian(size = 3)

##            [,1]      [,2]       [,3]
## [1,] 0.07511361 0.1238414 0.07511361
## [2,] 0.12384140 0.2041800 0.12384140
## [3,] 0.07511361 0.1238414 0.07511361

Laplace e Sobel

Os kernels de laplace e sobel tem formas bem definidas, e por isso a implementação não tem “truques”, somente a forma dos kernels:

kernel_laplacian

## function (diagonals = F) 
## {
##     if (diagonals) 
##         return(-1 * matrix(c(1, 1, 1, 1, -8, 1, 1, 1, 1), 3))
##     else return(-1 * matrix(c(0, 1, 0, 1, -4, 1, 0, 1, 0), 3))
## }

kernel_sobel

## function (horizontal = T) 
## {
##     kern <- matrix(c(-1, 0, 1, -2, 0, 2, -1, 0, 1), 3)
##     if (horizontal) 
##         return(t(kern))
##     else return(-1 * kern)
## }

Exemplo de kernel laplaciano sem considerar diagonais:

kernel_laplacian(diagonals = F)

##      [,1] [,2] [,3]
## [1,]    0   -1    0
## [2,]   -1    4   -1
## [3,]    0   -1    0

Exemplo de kernel de sobel para a coordenada \(x\):

kernel_sobel(horizontal = T)

##      [,1] [,2] [,3]
## [1,]   -1    0    1
## [2,]   -2    0    2
## [3,]   -1    0    1

Função de Aplicação dos Kernels

Para aplicar os kernels ao longo das coordenadas, pixel a pixel, implementei a seguinte função:

apply_linear_kernel

## function (img, kernel) 
## {
##     pad_size <- (nrow(kernel) - 1)/2
##     img <- padding_image(img = img, ncells = pad_size)
##     img_new <- outer(X = (1 + pad_size):(nrow(img) - pad_size), 
##         Y = (1 + pad_size):(ncol(img) - pad_size), FUN = Vectorize(function(rw, 
##             cl) {
##             img_subset <- img[(rw - pad_size):(rw + pad_size), 
##                 (cl - pad_size):(cl + pad_size)]
##             result <- sum(img_subset * kernel)
##             return(result)
##         }))
##     return(minmax(img_new))
## }

A lógica da função é a seguinte:

1. Defina o tamanho da margem, baseando-se na dimensão do kernel;
2. Adicione margens de zeros à imagem;
3. Percorra pixel a pixel;
  3.1. Obtenha a submatriz centrada no pixel de referência;
  3.2. Avalie o somatório dos elementos do produto da submatriz e do kernel;
  3.3. Substitua o pixel de referência pelo resultado acima;
4. Retorne a imagem transformada.

Obs: A função se chama apply_linear_kernel por que não aplica kernels de mediana, quartil e etc., para isso criei uma outra função apply_nonlinear_kernel… mas como não usei aqui, não irei entrar no mérito.

Experimento

O melhor resultado até o momento de escrita desse trabalho foi obtido com o processamento da imagem obedecendo o seguinte fluxo:

O Fluxo começa com a obtenção do negativo da imagem, já que os pulmões estão na parte de baixa intensidade. Após, os filtros gaussiano e laplacianos são empregados para remover ruídos e evidenciar as possiveis bordas dos pulmões. Em seguida, adicionamos o filtro à imagem original e reescalonamos para o intervalo [0, 1]. Prosseguindo, utilizamos um processo de “alongamento” de intensidades no intervalo [0.18, 0.56]. Finalmente, binarizamos a imagem com um limiar de 0.18 e extraimos as regiões conectadas (filtrando para manter as maiores regiões, os dois pulmões).

O Código para obtenção do resultado final de segmentação segue:

img_gray <- 1 - img_gray

img_filter <- img_gray %>%
  apply_linear_kernel(kernel_gaussian(size = 35, sigma = 40)) %>%
  apply_linear_kernel(kernel_laplacian(diagonals = T))

img_filtered <- minmax(img_gray + img_filter)

img_bin <- img_filtered %>% 
  histogram_stretch(0.18, 0.56) %>% 
  as.cimg %>% threshold(0.18)

img_connections <- img_bin %>%
  split_connected %>% 
  purrr::keep(~ sum(.) > 100000) 

par(mfrow = c(1, 2))
temp = lapply(img_connections, plot)

img_connections_filtered <- lapply(img_connections, function(pixelset) {
  pixelset %>%
    as.data.frame %>%
    group_by(y) %>%
    filter(x == max(x) | x == min(x)) %>%
    mutate(value = 1, left = max(x) < 500) 
}) %>% Reduce(rbind, .) 

img_connections_filtered %>% 
  filter(value == 1) %>%
  ggplot(aes(x = x, y = y, color = left)) +
  geom_polygon(show.legend = F) +
  scale_y_reverse() + 
  theme_bw() +
  theme(panel.background = element_rect(fill = "black"),
        panel.border = element_blank(), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())