Quando o que não é suficiente
É verdade que às vezes é important distinguir entre diferentes tipos de objetos. É um carro vindo em alta velocidade em minha direção; nesse caso, é melhor eu pular fora do caminho? Ou é um Doberman enorme (nesse caso eu provavelmente faria o mesmo)? Porém, muitas vezes na vida actual, em vez de granularidade grosseira classificaçãoo que é necessário é refinado segmentação.
Ao ampliar as imagens, não procuramos um único rótulo; em vez disso, queremos classificar cada pixel de acordo com algum critério:
Na medicina, podemos querer distinguir entre diferentes tipos de células ou identificar tumores.
Em várias ciências da terra, dados de satélite são usados para segmentar superfícies terrestres.
Para permitir o uso de planos de fundo personalizados, o software program de videoconferência deve ser capaz de distinguir o primeiro plano do plano de fundo.
A segmentação de imagens é uma forma de aprendizagem supervisionada: é necessário algum tipo de verdade básica. Aqui, ele vem em forma de máscara – uma imagem, de resolução espacial idêntica à dos dados de entrada, que designa a verdadeira classe de cada pixel. Conseqüentemente, a perda de classificação é calculada em pixels; as perdas são então somadas para produzir um agregado a ser usado na otimização.
A arquitetura “canônica” para segmentação de imagens é Rede U (em torno de 2015).
Rede U
Aqui está o protótipo da U-Internet, conforme descrito no Rönneberger et al. papel (Ronneberger, Fischer e Brox 2015).
Desta arquitetura, existem inúmeras variantes. Você pode usar diferentes tamanhos de camada, ativações, maneiras de reduzir e aumentar o tamanho e muito mais. Porém, há uma característica definidora: o formato em U, estabilizado pelas “pontes” que se cruzam horizontalmente em todos os níveis.
Resumindo, o lado esquerdo do U se assemelha às arquiteturas convolucionais usadas na classificação de imagens. Reduz sucessivamente a resolução espacial. Ao mesmo tempo, outra dimensão – a canais dimensão – é usada para construir uma hierarquia de recursos, variando de muito básicos a muito especializados.
Ao contrário da classificação, porém, a saída deve ter a mesma resolução espacial que a entrada. Assim, precisamos de aumentar novamente – isto é resolvido pelo lado direito dos EUA. Mas, como vamos chegar a um bom por pixel classificação, agora que tanta informação espacial foi perdida?
É para isso que servem as “pontes”: em cada nível, a entrada para uma camada de upsampling é um concatenação da saída da camada anterior – que passou por toda a rotina de compressão/descompressão – e alguma representação intermediária preservada da fase de downsizing. Desta forma, uma arquitetura U-Internet combina atenção aos detalhes com extração de recursos.
Segmentação de imagens cerebrais
Com a U-Internet, a aplicabilidade do domínio é tão ampla quanto a arquitetura é flexível. Aqui, queremos detectar anormalidades em exames cerebrais. O conjunto de dados, usado em Buda, Saha e Mazurowski (2019)contém imagens de ressonância magnética juntamente com imagens criadas manualmente ESTILO máscaras de segmentação de anormalidades. Está disponível em Kaggle.
Muito bem, o documento é acompanhado por um Repositório GitHub. Abaixo, acompanhamos de perto (embora não repliquemos exatamente) o código de pré-processamento e aumento de dados dos autores.
Como costuma acontecer em imagens médicas, há um notável desequilíbrio de courses nos dados. Para cada paciente, foram realizadas seções em múltiplas posições. (O número de seções por paciente varia.) A maioria das seções não apresenta lesões; as máscaras correspondentes são pretas em todos os lugares.
Aqui estão três exemplos onde as máscaras fazer indicar anormalidades:
Vamos ver se conseguimos construir uma U-Internet que gere essas máscaras para nós.
Dados
Antes de começar a digitar, aqui está um Caderno colaborativo para acompanhar convenientemente.
Nós usamos pins para obter os dados. Por favor veja esta introdução se você nunca usou esse pacote antes.
O conjunto de dados não é tão grande – inclui exames de 110 pacientes diferentes – então teremos que fazer apenas um conjunto de treinamento e validação. (Não faça isso na vida actual, pois você inevitavelmente acabará ajustando a última.)
train_dir <- "information/mri_train"
valid_dir <- "information/mri_valid"
if(dir.exists(train_dir)) unlink(train_dir, recursive = TRUE, drive = TRUE)
if(dir.exists(valid_dir)) unlink(valid_dir, recursive = TRUE, drive = TRUE)
zip::unzip(recordsdata, exdir = "information")
file.rename("information/kaggle_3m", train_dir)
# this can be a duplicate, once more containing kaggle_3m (evidently a packaging error on Kaggle)
# we simply take away it
unlink("information/lgg-mri-segmentation", recursive = TRUE)
dir.create(valid_dir)Desses 110 pacientes, mantemos 30 para validação. Mais algumas manipulações de arquivos e estamos configurados com uma bela estrutura hierárquica, com train_dir e valid_dir mantendo seus subdiretórios por paciente, respectivamente.
valid_indices <- pattern(1:size(sufferers), 30)
sufferers <- listing.dirs(train_dir, recursive = FALSE)
for (i in valid_indices) {
dir.create(file.path(valid_dir, basename(sufferers(i))))
for (f in listing.recordsdata(sufferers(i))) {
file.rename(file.path(train_dir, basename(sufferers(i)), f), file.path(valid_dir, basename(sufferers(i)), f))
}
unlink(file.path(train_dir, basename(sufferers(i))), recursive = TRUE)
}Precisamos agora de um dataset que sabe o que fazer com esses arquivos.
Conjunto de dados
Como todo torch conjunto de dados, este tem initialize() e .getitem() métodos. initialize() cria um inventário de nomes de arquivos de varredura e máscara, para ser usado por .getitem() quando ele realmente lê esses arquivos. Em contraste com o que vimos em postagens anteriores, porém, .getitem() não retorna simplesmente pares de entrada-alvo em ordem. Em vez disso, sempre que o parâmetro random_sampling é verdade, realizará amostragem ponderada, preferindo itens com lesões consideráveis. Esta opção será utilizada para o conjunto de treinamento, para compensar o desequilíbrio de courses mencionado acima.
A outra diferença entre os conjuntos de treinamento e validação é o uso de aumento de dados. As imagens/máscaras de treinamento podem ser invertidas, redimensionadas e giradas; probabilidades e valores são configuráveis.
Uma instância de brainseg_dataset encapsula toda essa funcionalidade:
brainseg_dataset <- dataset(
identify = "brainseg_dataset",
initialize = operate(img_dir,
augmentation_params = NULL,
random_sampling = FALSE) {
self$photos <- tibble(
img = grep(
listing.recordsdata(
img_dir,
full.names = TRUE,
sample = "tif",
recursive = TRUE
),
sample = 'masks',
invert = TRUE,
worth = TRUE
),
masks = grep(
listing.recordsdata(
img_dir,
full.names = TRUE,
sample = "tif",
recursive = TRUE
),
sample = 'masks',
worth = TRUE
)
)
self$slice_weights <- self$calc_slice_weights(self$photos$masks)
self$augmentation_params <- augmentation_params
self$random_sampling <- random_sampling
},
.getitem = operate(i) {
index <-
if (self$random_sampling == TRUE)
pattern(1:self$.size(), 1, prob = self$slice_weights)
else
i
img <- self$photos$img(index) %>%
image_read() %>%
transform_to_tensor()
masks <- self$photos$masks(index) %>%
image_read() %>%
transform_to_tensor() %>%
transform_rgb_to_grayscale() %>%
torch_unsqueeze(1)
img <- self$min_max_scale(img)
if (!is.null(self$augmentation_params)) {
scale_param <- self$augmentation_params(1)
c(img, masks) %<-% self$resize(img, masks, scale_param)
rot_param <- self$augmentation_params(2)
c(img, masks) %<-% self$rotate(img, masks, rot_param)
flip_param <- self$augmentation_params(3)
c(img, masks) %<-% self$flip(img, masks, flip_param)
}
listing(img = img, masks = masks)
},
.size = operate() {
nrow(self$photos)
},
calc_slice_weights = operate(masks) {
weights <- map_dbl(masks, operate(m) {
img <-
as.integer(magick::image_data(image_read(m), channels = "grey"))
sum(img / 255)
})
sum_weights <- sum(weights)
num_weights <- size(weights)
weights <- weights %>% map_dbl(operate(w) {
w <- (w + sum_weights * 0.1 / num_weights) / (sum_weights * 1.1)
})
weights
},
min_max_scale = operate(x) {
min = x$min()$merchandise()
max = x$max()$merchandise()
x$clamp_(min = min, max = max)
x$add_(-min)$div_(max - min + 1e-5)
x
},
resize = operate(img, masks, scale_param) {
img_size <- dim(img)(2)
rnd_scale <- runif(1, 1 - scale_param, 1 + scale_param)
img <- transform_resize(img, measurement = rnd_scale * img_size)
masks <- transform_resize(masks, measurement = rnd_scale * img_size)
diff <- dim(img)(2) - img_size
if (diff > 0) {
prime <- ceiling(diff / 2)
left <- ceiling(diff / 2)
img <- transform_crop(img, prime, left, img_size, img_size)
masks <- transform_crop(masks, prime, left, img_size, img_size)
} else {
img <- transform_pad(img,
padding = -c(
ceiling(diff / 2),
ground(diff / 2),
ceiling(diff / 2),
ground(diff / 2)
))
masks <- transform_pad(masks, padding = -c(
ceiling(diff / 2),
ground(diff /
2),
ceiling(diff /
2),
ground(diff /
2)
))
}
listing(img, masks)
},
rotate = operate(img, masks, rot_param) {
rnd_rot <- runif(1, 1 - rot_param, 1 + rot_param)
img <- transform_rotate(img, angle = rnd_rot)
masks <- transform_rotate(masks, angle = rnd_rot)
listing(img, masks)
},
flip = operate(img, masks, flip_param) {
rnd_flip <- runif(1)
if (rnd_flip > flip_param) {
img <- transform_hflip(img)
masks <- transform_hflip(masks)
}
listing(img, masks)
}
)Após a instanciação, vemos que temos 2.977 pares de treinamento e 952 pares de validação, respectivamente:
Como verificação de correção, vamos plotar uma imagem e uma máscara associada:
Com torché simples inspecionar o que acontece quando você altera os parâmetros relacionados ao aumento. Apenas escolhemos um par do conjunto de validação, que ainda não teve nenhum aumento aplicado, e chamamos valid_ds$ diretamente. Apenas por diversão, vamos usar parâmetros mais “extremos” aqui do que no treinamento actual. (O treinamento actual usa as configurações do repositório GitHub de Mateusz, que presumimos ter sido cuidadosamente escolhidas para um desempenho ultimate.)
img_and_mask <- valid_ds(77)
img <- img_and_mask((1))
masks <- img_and_mask((2))
imgs <- map (1:24, operate(i) {
# scale issue; train_ds actually makes use of 0.05
c(img, masks) %<-% valid_ds$resize(img, masks, 0.2)
c(img, masks) %<-% valid_ds$flip(img, masks, 0.5)
# rotation angle; train_ds actually makes use of 15
c(img, masks) %<-% valid_ds$rotate(img, masks, 90)
img %>%
transform_rgb_to_grayscale() %>%
as.array() %>%
as_tibble() %>%
rowid_to_column(var = "Y") %>%
collect(key = "X", worth = "worth", -Y) %>%
mutate(X = as.numeric(gsub("V", "", X))) %>%
ggplot(aes(X, Y, fill = worth)) +
geom_raster() +
theme_void() +
theme(legend.place = "none") +
theme(side.ratio = 1)
})
plot_grid(plotlist = imgs, nrow = 4)
Agora ainda precisamos dos carregadores de dados e nada nos impede de prosseguir para a próxima grande tarefa: construir o modelo.
batch_size <- 4
train_dl <- dataloader(train_ds, batch_size)
valid_dl <- dataloader(valid_ds, batch_size)Modelo
Nosso modelo ilustra bem o tipo de código modular que vem “naturalmente” com torch. Abordamos as coisas de cima para baixo, começando com o próprio contêiner U-Internet.
unet cuida da composição international – até onde “para baixo” vamos, diminuindo a imagem enquanto aumentamos o número de filtros, e então como vamos “para cima” novamente?
É importante ressaltar que também está na memória do sistema. Em ahead()ele monitora as saídas da camada vistas “descendo”, para serem adicionadas novamente ao “subir”.
unet <- nn_module(
"unet",
initialize = operate(channels_in = 3,
n_classes = 1,
depth = 5,
n_filters = 6) {
self$down_path <- nn_module_list()
prev_channels <- channels_in
for (i in 1:depth) {
self$down_path$append(down_block(prev_channels, 2 ^ (n_filters + i - 1)))
prev_channels <- 2 ^ (n_filters + i -1)
}
self$up_path <- nn_module_list()
for (i in ((depth - 1):1)) {
self$up_path$append(up_block(prev_channels, 2 ^ (n_filters + i - 1)))
prev_channels <- 2 ^ (n_filters + i - 1)
}
self$final = nn_conv2d(prev_channels, n_classes, kernel_size = 1)
},
ahead = operate(x) {
blocks <- listing()
for (i in 1:size(self$down_path)) {
x <- self$down_path((i))(x)
if (i != size(self$down_path)) {
blocks <- c(blocks, x)
x <- nnf_max_pool2d(x, 2)
}
}
for (i in 1:size(self$up_path)) {
x <- self$up_path((i))(x, blocks((size(blocks) - i + 1))$to(system = system))
}
torch_sigmoid(self$final(x))
}
)unet delega para dois contêineres emblem abaixo dele na hierarquia: down_block e up_block. Enquanto down_block está “apenas” ali por razões estéticas (ele delega imediatamente ao seu próprio burro de carga, conv_block), em up_block vemos as “pontes” da U-Internet em ação.
down_block <- nn_module(
"down_block",
initialize = operate(in_size, out_size) {
self$conv_block <- conv_block(in_size, out_size)
},
ahead = operate(x) {
self$conv_block(x)
}
)
up_block <- nn_module(
"up_block",
initialize = operate(in_size, out_size) {
self$up = nn_conv_transpose2d(in_size,
out_size,
kernel_size = 2,
stride = 2)
self$conv_block = conv_block(in_size, out_size)
},
ahead = operate(x, bridge) {
up <- self$up(x)
torch_cat(listing(up, bridge), 2) %>%
self$conv_block()
}
)Finalmente, um conv_block é uma estrutura sequencial contendo camadas convolucionais, ReLU e dropout.
conv_block <- nn_module(
"conv_block",
initialize = operate(in_size, out_size) {
self$conv_block <- nn_sequential(
nn_conv2d(in_size, out_size, kernel_size = 3, padding = 1),
nn_relu(),
nn_dropout(0.6),
nn_conv2d(out_size, out_size, kernel_size = 3, padding = 1),
nn_relu()
)
},
ahead = operate(x){
self$conv_block(x)
}
)Agora instancie o modelo e, possivelmente, mova-o para a GPU:
system <- torch_device(if(cuda_is_available()) "cuda" else "cpu")
mannequin <- unet(depth = 5)$to(system = system)Otimização
Treinamos nosso modelo com uma combinação de entropia cruzada e perda de dados.
Este último, embora não seja enviado com torchpode ser implementado manualmente:
calc_dice_loss <- operate(y_pred, y_true) {
easy <- 1
y_pred <- y_pred$view(-1)
y_true <- y_true$view(-1)
intersection <- (y_pred * y_true)$sum()
1 - ((2 * intersection + easy) / (y_pred$sum() + y_true$sum() + easy))
}
dice_weight <- 0.3A otimização usa gradiente descendente estocástico (SGD), juntamente com o escalonador de taxa de aprendizagem de um ciclo introduzido no contexto de classificação de imagens com tocha.
optimizer <- optim_sgd(mannequin$parameters, lr = 0.1, momentum = 0.9)
num_epochs <- 20
scheduler <- lr_one_cycle(
optimizer,
max_lr = 0.1,
steps_per_epoch = size(train_dl),
epochs = num_epochs
)Treinamento
O ciclo de treinamento segue então o esquema ordinary. Uma coisa a observar: a cada época, salvamos o modelo (usando torch_save()), para que possamos escolher mais tarde o melhor, caso o desempenho tenha piorado depois disso.
train_batch <- operate(b) {
optimizer$zero_grad()
output <- mannequin(b((1))$to(system = system))
goal <- b((2))$to(system = system)
bce_loss <- nnf_binary_cross_entropy(output, goal)
dice_loss <- calc_dice_loss(output, goal)
loss <- dice_weight * dice_loss + (1 - dice_weight) * bce_loss
loss$backward()
optimizer$step()
scheduler$step()
listing(bce_loss$merchandise(), dice_loss$merchandise(), loss$merchandise())
}
valid_batch <- operate(b) {
output <- mannequin(b((1))$to(system = system))
goal <- b((2))$to(system = system)
bce_loss <- nnf_binary_cross_entropy(output, goal)
dice_loss <- calc_dice_loss(output, goal)
loss <- dice_weight * dice_loss + (1 - dice_weight) * bce_loss
listing(bce_loss$merchandise(), dice_loss$merchandise(), loss$merchandise())
}
for (epoch in 1:num_epochs) {
mannequin$practice()
train_bce <- c()
train_dice <- c()
train_loss <- c()
coro::loop(for (b in train_dl) {
c(bce_loss, dice_loss, loss) %<-% train_batch(b)
train_bce <- c(train_bce, bce_loss)
train_dice <- c(train_dice, dice_loss)
train_loss <- c(train_loss, loss)
})
torch_save(mannequin, paste0("model_", epoch, ".pt"))
cat(sprintf("nEpoch %d, coaching: loss:%3f, bce: %3f, cube: %3fn",
epoch, imply(train_loss), imply(train_bce), imply(train_dice)))
mannequin$eval()
valid_bce <- c()
valid_dice <- c()
valid_loss <- c()
i <- 0
coro::loop(for (b in tvalid_dl) {
i <<- i + 1
c(bce_loss, dice_loss, loss) %<-% valid_batch(b)
valid_bce <- c(valid_bce, bce_loss)
valid_dice <- c(valid_dice, dice_loss)
valid_loss <- c(valid_loss, loss)
})
cat(sprintf("nEpoch %d, validation: loss:%3f, bce: %3f, cube: %3fn",
epoch, imply(valid_loss), imply(valid_bce), imply(valid_dice)))
}Epoch 1, coaching: loss:0.304232, bce: 0.148578, cube: 0.667423
Epoch 1, validation: loss:0.333961, bce: 0.127171, cube: 0.816471
Epoch 2, coaching: loss:0.194665, bce: 0.101973, cube: 0.410945
Epoch 2, validation: loss:0.341121, bce: 0.117465, cube: 0.862983
(...)
Epoch 19, coaching: loss:0.073863, bce: 0.038559, cube: 0.156236
Epoch 19, validation: loss:0.302878, bce: 0.109721, cube: 0.753577
Epoch 20, coaching: loss:0.070621, bce: 0.036578, cube: 0.150055
Epoch 20, validation: loss:0.295852, bce: 0.101750, cube: 0.748757Avaliação
Nesta execução, é o modelo ultimate que apresenta melhor desempenho no conjunto de validação. Ainda assim, gostaríamos de mostrar como carregar um modelo salvo, usando torch_load() .
Uma vez carregado, coloque o modelo em eval modo:
saved_model <- torch_load("model_20.pt")
mannequin <- saved_model
mannequin$eval()Agora, como não temos um conjunto de testes separado, já conhecemos as métricas médias fora da amostra; mas no ultimate, o que nos importa são as máscaras geradas. Vamos ver alguns, exibindo dados reais e exames de ressonância magnética para comparação.
# with out random sampling, we might primarily see lesion-free patches
eval_ds <- brainseg_dataset(valid_dir, augmentation_params = NULL, random_sampling = TRUE)
eval_dl <- dataloader(eval_ds, batch_size = 8)
batch <- eval_dl %>% dataloader_make_iter() %>% dataloader_next()
par(mfcol = c(3, 8), mar = c(0, 1, 0, 1))
for (i in 1:8) {
img <- batch((1))(i, .., drop = FALSE)
inferred_mask <- mannequin(img$to(system = system))
true_mask <- batch((2))(i, .., drop = FALSE)$to(system = system)
bce <- nnf_binary_cross_entropy(inferred_mask, true_mask)$to(system = "cpu") %>%
as.numeric()
dc <- calc_dice_loss(inferred_mask, true_mask)$to(system = "cpu") %>% as.numeric()
cat(sprintf("nSample %d, bce: %3f, cube: %3fn", i, bce, dc))
inferred_mask <- inferred_mask$to(system = "cpu") %>% as.array() %>% .(1, 1, , )
inferred_mask <- ifelse(inferred_mask > 0.5, 1, 0)
img(1, 1, ,) %>% as.array() %>% as.raster() %>% plot()
true_mask$to(system = "cpu")(1, 1, ,) %>% as.array() %>% as.raster() %>% plot()
inferred_mask %>% as.raster() %>% plot()
}Também imprimimos a entropia cruzada particular person e as perdas nos dados; relacioná-los às máscaras geradas pode fornecer informações úteis para o ajuste do modelo.
Pattern 1, bce: 0.088406, cube: 0.387786}
Pattern 2, bce: 0.026839, cube: 0.205724
Pattern 3, bce: 0.042575, cube: 0.187884
Pattern 4, bce: 0.094989, cube: 0.273895
Pattern 5, bce: 0.026839, cube: 0.205724
Pattern 6, bce: 0.020917, cube: 0.139484
Pattern 7, bce: 0.094989, cube: 0.273895
Pattern 8, bce: 2.310956, cube: 0.999824
Embora longe de serem perfeitas, a maioria dessas máscaras não é tão ruim assim – um bom resultado dado o pequeno conjunto de dados!
Conclusão
Este tem sido o nosso trabalho mais complexo torch postar até agora; no entanto, esperamos que você tenha aproveitado bem o tempo. Por um lado, entre as aplicações de aprendizagem profunda, a segmentação de imagens médicas se destaca como altamente útil para a sociedade. Em segundo lugar, arquiteturas do tipo U-Internet são empregadas em muitas outras áreas. E finalmente, vimos mais uma vez torchflexibilidade e comportamento intuitivo em ação.
Obrigado por ler!
Buda, Mateusz, Ashirbani Saha e Maciej A. Mazurowski. 2019. “Associação de subtipos genômicos de gliomas de grau inferior com características de forma extraídas automaticamente por um algoritmo de aprendizado profundo.” Computadores em Biologia e Medicina 109: 218–25. https://doi.org/https://doi.org/10.1016/j.compbiomed.2019.05.002.
Ronneberger, Olaf, Philipp Fischer e Thomas Brox. 2015. “U-Internet: Redes Convolucionais para Segmentação de Imagens Biomédicas.” CoRR abs/1505.04597. http://arxiv.org/abs/1505.04597.