下載資料

請你先到資料科學競賽網站選擇第七個任務，登入後下載4份資料，包含之前任務沒有的「supplement.zip」，裡面放著影像與標註資料。

– 與之前所有任務相同，這個任務的資料下載是需要經過申請的，請你找助教申請帳號。

讀取資料

最開始我們還是需要先將「train.csv」、「test.csv」以及「sample_submission.csv」讀進來：

library(data.table)

train_dat <- fread('train.csv', data.table = FALSE)
test_dat <- fread('test.csv', data.table = FALSE)
submit_dat <- fread('sample_submission.csv', data.table = FALSE)

如果你是使用伺服器版的Rstudio，那直接上傳上傳「supplement.zip」就能得到資料夾「image」及「label」了。或是你可以先解壓縮放到你的project的資料夾下。

檔案格式

電腦斷層影像是對病人做連續切面的描述：

F02

因此，影像是以「一組」的方式呈現的：

F01

我們可以透過這種方式把第一組「P001」中的第10張影像「U010」給讀進來，它會呈現成這樣

library(OpenImageR)
library(imager)

img_dir <- 'image/'
img_pos <- 10L

img_path <- paste0(img_dir, train_dat[img_pos, 'PID'], '/', train_dat[img_pos, 'UID'], '.png')
img <- readImage(img_path)
par(mar = rep(0, 4))
plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(img, 0, 1, 1, 0, interpolate=FALSE)

標籤格式

在三軍總醫院的標註系統中，我們在標註結束後取得的資料是以「json」的格式儲存的。

– 它的檔案格式有點複雜，但如果你想要在原來那張圖上面做堆疊，你可以透過這種方式做到：

library(jsonlite)

# Settings

label_dir <- 'label/'
img_pos <- 10L

# Read data

label_dat <- fread('label/label_list.csv', data.table = FALSE)

json_path <- paste0(label_dir, train_dat[img_pos, 'PID'], '/', train_dat[img_pos, 'UID'], '.json')
anno_list <- fromJSON(json_path)

# Show image

par(mar = rep(0, 4))
plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(img, 0, 1, 1, 0, interpolate=FALSE)

# Show annotations

for (m in 1:length(anno_list[['shapes']][['tagAnno']])) {
  
  # Processing
  
  seg_data <- anno_list[['shapes']][['segments']][[m]]
  
  x_seq <- t(seg_data[,c(1,3)]) %>% as.numeric()
  x_seq <- x_seq / anno_list[['size']][1]
  
  y_seq <- t(seg_data[,c(2,4)]) %>% as.numeric()
  y_seq <- y_seq / anno_list[['size']][2]    
  
  poly_col <- label_dat[label_dat[,'structure'] %in% anno_list[['shapes']][['tagAnno']][m], 'color']
  poly_col <- paste0(poly_col, 'A0')
  polygon(x = x_seq, y = y_seq, col = poly_col, border = poly_col)
  
}

訓練時我們必須將格式做一些處理，處理的範例如下：

library(jsonlite)
library(tiff)

# Settings

label_dir <- 'label/'
img_pos <- 10L

# Read data

json_path <- paste0(label_dir, train_dat[img_pos, 'PID'], '/', train_dat[img_pos, 'UID'], '.json')
anno_list <- fromJSON(json_path)

# Get info

current_dir <- paste0(label_dir, train_dat[img_pos,'PID'], '/', train_dat[img_pos,'UID'], '/')
if (!dir.exists(current_dir)) {dir.create(current_dir, recursive = TRUE)}

obj_names <- unique(anno_list[['shapes']][['tagAnno']])

for (m in 1:length(obj_names)) {
  
  # Plotting
  
  tiff(filename = paste0(current_dir, obj_names[m], '.tiff'), width = 256, height = 256)
  
  par(mar = rep(0, 4))
  plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
  
  # Processing
  
  used_id <- which(anno_list[['shapes']][['tagAnno']] %in% obj_names[m])
  if (length(used_id) > 1) {message(original_data[img_pos,'PID'], '/', original_data[img_pos,'UID'], ':', obj_names[m])}
  
  for (q in used_id) {
    
    seg_data <- anno_list[['shapes']][['segments']][[q]]
    
    x_seq <- t(seg_data[,c(1,3)]) %>% as.numeric()
    x_seq <- x_seq / anno_list[['size']][1]
    
    y_seq <- t(seg_data[,c(2,4)]) %>% as.numeric()
    y_seq <- y_seq / anno_list[['size']][2]    
    
    polygon(x = x_seq, y = y_seq, col = '#000000', border = '#000000')
    
  }
  
  dev.off()
  
  # Binarization
  
  seg_img <- readImage(paste0(current_dir, obj_names[m], '.tiff'))
  if (dim(seg_img)[3] %in% 3) {seg_img <- seg_img %>% rgb_2gray()}
  if (prop.table(table(seg_img))['1'] > 0.5) {seg_img <- 1 - (seg_img > 0.5)} else {seg_img <- (seg_img > 0.5) + 0}
  seg_img <- matrix(seg_img, nrow = 256, ncol = 256)
    
  writeTIFF(what = seg_img, where = paste0(current_dir, obj_names[m], '.tiff'), compression = 'none')
  
}

這個過程之前就已經做完了，讓我們試著讀取這些檔案進行套疊：

library(magrittr)

# Settings

label_dir <- 'label/'
img_dir <- 'image/'

img_pos <- 10L

# Read data

label_dat <- fread('label/label_list.csv', data.table = FALSE)

# Show image

img_path <- paste0(img_dir, train_dat[img_pos, 'PID'], '/', train_dat[img_pos, 'UID'], '.png')
img <- readImage(img_path)

par(mar = rep(0, 4))
plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(img, 0, 1, 1, 0, interpolate = FALSE)

# Read annotations

current_dir <- paste0(label_dir, train_dat[img_pos, 'PID'], '/', train_dat[img_pos, 'UID'], '/')
anno_files <- list.files(current_dir) %>% gsub('\\.tiff', '', .)

# Show annotations

for (m in 1:length(anno_files)) {
  
  anno_img <- readImage(paste0(current_dir, anno_files[m], '.tiff'))
  anno_img <- (anno_img > 0.5)
  
  poly_col <- label_dat[label_dat[,'structure'] %in% anno_files[m], 'color']
  poly_col <- paste0(poly_col, 'A0')  
  
  anno_raster <- c('#00000000', poly_col)[anno_img + 1L] %>% array(., dim = dim(anno_img))
  rasterImage(anno_raster, 0, 1, 1, 0, interpolate = FALSE)
  
}

標籤編碼

有了每張圖的標籤之後，我們要對資料及標籤進行資料重整，編碼方式如下：

library(data.table)
library(OpenImageR)
library(magrittr)

# Settings

img_dir <- 'image/'
label_dir <- 'label/'
process_dir <- 'processed/'

# Read data

label_dat <- fread('label/label_list.csv', data.table = FALSE)

# Patient based processing

if (!dir.exists(process_dir)) {dir.create(process_dir)}
PIDs <- unique(train_dat[,'PID'])

tab_list <- list()

pb <- txtProgressBar(max = length(PIDs), style = 3)

for (i in 1:length(PIDs)) {
  
  current_path <- paste0(process_dir, PIDs[i], '.RData')
  sub_train_dat <- train_dat[train_dat[,'PID'] %in% PIDs[i],]
  
  x_array <- array(0, dim = c(256, 256, nrow(sub_train_dat), 3))
  y_array <- array(0L, dim = c(256, 256, nrow(sub_train_dat), nrow(label_dat)))
  
  for (j in 1:nrow(sub_train_dat)) {
    
    img_path <- paste0(img_dir, PIDs[i], '/', sub_train_dat[j, 'UID'], '.png')
    x_array[,,j,] <- readImage(img_path)
    
    current_label_dir <- paste0(label_dir, PIDs[i], '/', sub_train_dat[j, 'UID'], '/')
    label_files <- list.files(current_label_dir) %>% gsub('.tiff', '', .)
    
    if (length(label_files) > 0) {
      
      for (k in 1:length(label_files)) {
        
        label_path <- paste0(label_dir, PIDs[i], '/', sub_train_dat[j, 'UID'], '/', label_files[k], '.tiff')
        label_img <- readImage(label_path)
        pos_id <- which(label_dat[,'structure'] %in% label_files[k])
        if (dim(label_img)[3] %in% 3) {label_img <- label_img %>% rgb_2gray()}
        y_array[,,j,pos_id] <- (label_img > 0.5) + 0L
        
      }
      
    }
    
  }
  
  tab_list[[i]] <- apply(y_array, 3:4, sum)
  save(x_array, y_array, file = current_path)
  setTxtProgressBar(pb, i)
  
}

close(pb)

那這樣我們就能產生執行器了，需要注意的是我們需要用undersampling來解決類別不平衡問題，所以需要一個mask來做遮蔽：

library(mxnet)
library(abind)

# Settings

fold_pos <- 50L

img_dir <- 'image/'
process_dir <- 'processed/'
PIDs <- unique(train_dat[,'PID'])

neg_pixel <- c(461, 465, 238, 241, 51, 54, 29, 89, 133)
neg_pixel <- neg_pixel * fold_pos

# Read data

my_iterator_core <- function (batch_size, random_crop = 224) {
  
  batch <- 0
  batch_per_epoch <- nrow(train_dat) / batch_size
  
  reset <- function() {batch <<- 0}
  
  iter.next <- function() {
    batch <<- batch+1
    if (batch > batch_per_epoch) {return(FALSE)} else {return(TRUE)}
  }
  
  value <- function() {
    batch_PID <- sample(PIDs, 1, replace = TRUE)
    load(paste0(process_dir, batch_PID, '.RData'))
    batch_idxID <- sample(1:dim(x_array)[3], batch_size, replace = TRUE)
    batch_x_list <- list()
    batch_y_list <- list()
    batch_mask_list <- list()
    for (j in 1:batch_size) {
      random_row <- sample(0:(dim(x_array)[1] - random_crop), 1)
      random_col <- sample(0:(dim(x_array)[2] - random_crop), 1)
      batch_x_list[[j]] <- x_array[random_row+1:random_crop,random_col+1:random_crop,batch_idxID[j],]
      batch_y_list[[j]] <- y_array[random_row+1:random_crop,random_col+1:random_crop,batch_idxID[j],]
      sample_mask_list <- list()
      for (k in 1:length(neg_pixel)) {
        sub_mask_array <- array(0L, dim = c(random_crop, random_crop))
        neg_mask_pos <- sample(1:length(sub_mask_array), neg_pixel[k])
        sub_mask_array[neg_mask_pos] <- 1L
        pos_mask_pos <- which(y_array[random_row+1:random_crop,random_col+1:random_crop,batch_idxID[j],k] == 1)
        sub_mask_array[pos_mask_pos] <- fold_pos
        mirror_mask_pos <- which(y_array[random_row+1:random_crop,(dim(y_array)[2] + 1) - (random_col+1:random_crop),batch_idxID[j],k] == 1)
        sub_mask_array[mirror_mask_pos] <- fold_pos
        sample_mask_list[[k]] <- sub_mask_array
      }
      batch_mask_list[[j]] <- abind(sample_mask_list, along = 3)
    }
    batch_x <- abind(batch_x_list, along = 4)
    batch_y <- abind(batch_y_list, along = 4)
    batch_mask <- abind(batch_mask_list, along = 4)
    data <- mx.nd.array(batch_x)
    label <- mx.nd.array(batch_y)
    mask <- mx.nd.array(batch_mask)
    return(list(data = data, label = label, mask = mask))
  }
  
  return(list(reset = reset, iter.next = iter.next, value = value, batch_size = batch_size, batch = batch))
}

my_iterator_func <- setRefClass("Custom_Iter",
                                fields = c("iter", "batch_size", "random_crop"),
                                contains = "Rcpp_MXArrayDataIter",
                                methods = list(
                                  initialize = function (iter, batch_size = 8, random_crop = 224){
                                    .self$iter <- my_iterator_core(batch_size = batch_size, random_crop = random_crop)
                                    .self
                                  },
                                  value = function () {
                                    .self$iter$value()
                                  },
                                  iter.next = function () {
                                    .self$iter$iter.next()
                                  },
                                  reset = function () {
                                    .self$iter$reset()
                                  },
                                  finalize = function () {
                                  }
                                )
)

my_iter <- my_iterator_func(iter = NULL, batch_size = 16, random_crop = 224)

讓我們來看看他的成果：

label_dat <- fread('label/label_list.csv', data.table = FALSE)

my_iter$reset()
my_iter$iter.next()

## [1] TRUE

batch_data <- my_iter$value()

par(mar = rep(0, 4), mfrow = c(1, 1))

plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(as.raster(as.array(batch_data$data)[,,,1]), 0, 1, 1, 0, interpolate=FALSE)
for (k in 1:nrow(label_dat)) {
  current_col <- paste0(label_dat[k,'color'], 'A0')
  current_img <- as.array(batch_data$label)[,,k,1]
  if (!table(current_img)['1'] %in% NA) {
    current_img <- c('#00000000', current_col)[current_img + 1]
    current_img <- array(current_img, dim = c(224, 224))
    rasterImage(current_img, 0, 1, 1, 0, interpolate=FALSE)
  }
}

我們可以看看label與mask的結果，做以比較：

par(mar = rep(0, 4), mfcol = c(3, 6))

for (k in 1:nrow(label_dat)) {
  plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
  rasterImage(as.raster(as.array(batch_data$data)[,,,1]), 0, 1, 1, 0, interpolate=FALSE)
  current_col <- paste0(label_dat[k,'color'], 'A0')
  current_img <- as.array(batch_data$label)[,,k,1]
  current_img[current_img > 1] <- 1
  if (!table(current_img)['1'] %in% NA) {
    current_img <- c('#00000000', current_col)[current_img + 1]
    current_img <- array(current_img, dim = c(224, 224))
    rasterImage(current_img, 0, 1, 1, 0, interpolate = FALSE)
  }
}

for (k in 1:nrow(label_dat)) {
  plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
  rasterImage(as.raster(as.array(batch_data$data)[,,,1]), 0, 1, 1, 0, interpolate=FALSE)
  current_col <- paste0(label_dat[k,'color'], 'A0')
  current_img <- as.array(batch_data$mask)[,,k,1]
  current_img[current_img > 1] <- 1
  if (!table(current_img)['1'] %in% NA) {
    current_img <- c('#00000000', current_col)[current_img + 1]
    current_img <- array(current_img, dim = c(224, 224))
    rasterImage(current_img, 0, 1, 1, 0, interpolate = FALSE)
  }
}

模型訓練

接著就能進行模型訓練了，我們以最基礎版的U-net來看看：

library(mxnet)

data <- mx.symbol.Variable('data')
bn_data <- mx.symbol.BatchNorm(data = data, fix.gamma = TRUE, name = 'bn_data')

down_conv1 <- mx.symbol.Convolution(data = bn_data, kernel = c(7, 7), pad = c(3, 3), num_filter = 32, no.bias = TRUE, name = 'down_conv1')
down_convbn1 <- mx.symbol.BatchNorm(data = down_conv1, fix.gamma = FALSE, name = 'down_convbn1')
down_convrelu1 <- mx.symbol.Activation(data = down_convbn1, act_type = "relu", name = 'down_convrelu1')
down_convpool1 <- mx.symbol.Pooling(data = down_convrelu1, pool_type = "avg", kernel = c(2, 2), stride = c(2, 2), name = 'down_convpool1')

down_conv2 <- mx.symbol.Convolution(data = down_convpool1, kernel = c(5, 5), pad = c(2, 2), num_filter = 64, no.bias = TRUE, name = 'down_conv2')
down_convbn2 <- mx.symbol.BatchNorm(data = down_conv2, fix.gamma = FALSE, name = 'down_convbn2')
down_convrelu2 <- mx.symbol.Activation(data = down_convbn2, act_type = "relu", name = 'down_convrelu2')
down_convpool2 <- mx.symbol.Pooling(data = down_convrelu2, pool_type = "avg", kernel = c(2, 2), stride = c(2, 2), name = 'down_convpool2')

down_conv3 <- mx.symbol.Convolution(data = down_convpool2, kernel = c(3, 3), pad = c(1, 1), num_filter = 128, no.bias = TRUE, name = 'down_conv3')
down_convbn3 <- mx.symbol.BatchNorm(data = down_conv3, fix.gamma = FALSE, name = 'down_convbn3')
down_convrelu3 <- mx.symbol.Activation(data = down_convbn3, act_type = "relu", name = 'down_convrelu3')
down_convpool3 <- mx.symbol.Pooling(data = down_convrelu3, pool_type = "avg", kernel = c(2, 2), stride = c(2, 2), name = 'down_convpool3')

down_conv4 <- mx.symbol.Convolution(data = down_convpool3, kernel = c(3, 3), pad = c(1, 1), num_filter = 256, no.bias = TRUE, name = 'down_conv4')
down_convbn4 <- mx.symbol.BatchNorm(data = down_conv4, fix.gamma = FALSE, name = 'down_convbn4')
down_convrelu4 <- mx.symbol.Activation(data = down_convbn4, act_type = "relu", name = 'down_convrelu4')
down_convpool4 <- mx.symbol.Pooling(data = down_convrelu4, pool_type = "avg", kernel = c(2, 2), stride = c(2, 2), name = 'down_convpool4')

down_conv5 <- mx.symbol.Convolution(data = down_convpool4, kernel = c(3, 3), pad = c(1, 1), num_filter = 512, no.bias = TRUE, name = 'down_conv5')
down_convbn5 <- mx.symbol.BatchNorm(data = down_conv5, fix.gamma = FALSE, name = 'down_convbn5')
down_convrelu5 <- mx.symbol.Activation(data = down_convbn5, act_type = "relu", name = 'down_convrelu5')
down_convpool5 <- mx.symbol.Pooling(data = down_convrelu5, pool_type = "avg", kernel = c(2, 2), stride = c(2, 2), name = 'down_convpool5')

down_conv6 <- mx.symbol.Convolution(data = down_convpool5, kernel = c(3, 3), pad = c(1, 1), num_filter = 512, no.bias = TRUE, name = 'down_conv6')
down_convbn6 <- mx.symbol.BatchNorm(data = down_conv6, fix.gamma = FALSE, name = 'down_convbn6')
down_convrelu6 <- mx.symbol.Activation(data = down_convbn6, act_type = "relu", name = 'down_convrelu6')

up_deconv1 <- mx.symbol.Deconvolution(data = down_convrelu6, kernel = c(2, 2), stride = c(2, 2), num_filter = 512, name = 'up_deconv1')
up_deconvbn1 <- mx.symbol.BatchNorm(data = up_deconv1, fix.gamma = FALSE, name = 'up_deconvbn1')
up_deconvrelu1 <- mx.symbol.Activation(data = up_deconvbn1, act_type = "relu", name = 'up_deconvrelu1')
up_conv1 <- mx.symbol.Convolution(data = up_deconvrelu1, kernel = c(3, 3), pad = c(1, 1), num_filter = 512, no.bias = TRUE, name = 'up_conv1')
up_convbn1 <- mx.symbol.BatchNorm(data = up_conv1, fix.gamma = FALSE, name = 'up_convbn1')
up_convrelu1 <- mx.symbol.Activation(data = up_convbn1, act_type = "relu", name = 'up_convrelu1')

up_concat2 <- mx.symbol.concat(data = list(up_convrelu1, down_convrelu5), num.args = 2, dim = 1, name = 'up_concat2')
up_deconv2 <- mx.symbol.Deconvolution(data = up_concat2, kernel = c(2, 2), stride = c(2, 2), num_filter = 256, name = 'up_deconv2')
up_deconvbn2 <- mx.symbol.BatchNorm(data = up_deconv2, fix.gamma = FALSE, name = 'up_deconvbn2')
up_deconvrelu2 <- mx.symbol.Activation(data = up_deconvbn2, act_type = "relu", name = 'up_deconvrelu2')
up_conv2 <- mx.symbol.Convolution(data = up_deconvrelu2, kernel = c(3, 3), pad = c(1, 1), num_filter = 256, no.bias = TRUE, name = 'up_conv2')
up_convbn2 <- mx.symbol.BatchNorm(data = up_conv2, fix.gamma = FALSE, name = 'up_convbn2')
up_convrelu2 <- mx.symbol.Activation(data = up_convbn2, act_type = "relu", name = 'up_convrelu2')

up_concat3 <- mx.symbol.concat(data = list(up_convrelu2, down_convrelu4), num.args = 2, dim = 1, name = 'up_concat3')
up_deconv3 <- mx.symbol.Deconvolution(data = up_concat3, kernel = c(2, 2), stride = c(2, 2), num_filter = 256, name = 'up_deconv3')
up_deconvbn3 <- mx.symbol.BatchNorm(data = up_deconv3, fix.gamma = FALSE, name = 'up_deconvbn3')
up_deconvrelu3 <- mx.symbol.Activation(data = up_deconvbn3, act_type = "relu", name = 'up_deconvrelu3')
up_conv3 <- mx.symbol.Convolution(data = up_deconvrelu3, kernel = c(3, 3), pad = c(1, 1), num_filter = 256, no.bias = TRUE, name = 'up_conv3')
up_convbn3 <- mx.symbol.BatchNorm(data = up_conv3, fix.gamma = FALSE, name = 'up_convbn3')
up_convrelu3 <- mx.symbol.Activation(data = up_convbn3, act_type = "relu", name = 'up_convrelu3')

up_concat4 <- mx.symbol.concat(data = list(up_convrelu3, down_convrelu3), num.args = 2, dim = 1, name = 'up_concat4')
up_deconv4 <- mx.symbol.Deconvolution(data = up_concat4, kernel = c(2, 2), stride = c(2, 2), num_filter = 256, name = 'up_deconv4')
up_deconvbn4 <- mx.symbol.BatchNorm(data = up_deconv4, fix.gamma = FALSE, name = 'up_deconvbn4')
up_deconvrelu4 <- mx.symbol.Activation(data = up_deconvbn4, act_type = "relu", name = 'up_deconvrelu4')
up_conv4 <- mx.symbol.Convolution(data = up_deconvrelu4, kernel = c(3, 3), pad = c(1, 1), num_filter = 288, no.bias = TRUE, name = 'up_conv4')
up_convbn4 <- mx.symbol.BatchNorm(data = up_conv4, fix.gamma = FALSE, name = 'up_convbn4')
up_convrelu4 <- mx.symbol.Activation(data = up_convbn4, act_type = "relu", name = 'up_convrelu4')

up_deconv5 <- mx.symbol.Deconvolution(data = up_convrelu4, kernel = c(2, 2), stride = c(2, 2), num_filter = 288, num_group = 9, name = 'up_deconv5')
up_deconvbn5 <- mx.symbol.BatchNorm(data = up_deconv5, fix.gamma = FALSE, name = 'up_deconvbn5')
up_deconvrelu5 <- mx.symbol.Activation(data = up_deconvbn5, act_type = "relu", name = 'up_deconvrelu5')
up_conv5 <- mx.symbol.Convolution(data = up_deconvrelu5, kernel = c(3, 3), pad = c(1, 1), num_filter = 288, num_group = 9, no.bias = TRUE, name = 'up_conv5')
up_convbn5 <- mx.symbol.BatchNorm(data = up_conv5, fix.gamma = FALSE, name = 'up_convbn5')
up_convrelu5 <- mx.symbol.Activation(data = up_convbn5, act_type = "relu", name = 'up_convrelu5')

linear_pred <- mx.symbol.Convolution(data = up_convrelu5, kernel = c(3, 3), pad = c(1, 1), num_filter = 9, name = 'linear_pred')
logistic_pred <- mx.symbol.Activation(data = linear_pred, act.type = 'sigmoid', name = 'logistic_pred')

# CE loss

label <- mx.symbol.Variable(name = 'label')
mask <- mx.symbol.Variable(name = 'mask')

eps <- 1e-8
ce_loss_pos <- mx.symbol.broadcast_mul(mx.symbol.log(logistic_pred + eps), label)
ce_loss_neg <- mx.symbol.broadcast_mul(mx.symbol.log(1 - logistic_pred + eps), 1 - label)
ce_loss_sum <- mx.symbol.broadcast_add(ce_loss_pos, ce_loss_neg)
ce_loss_weight <- mx.symbol.broadcast_mul(ce_loss_sum, mask)
sum_mask <- mx.symbol.sum(mask, axis = 0:3)
sum_ce_loss <- mx.symbol.sum(ce_loss_weight, axis = 0:3)
mean_ce_loss <- mx.symbol.broadcast_div(sum_ce_loss, sum_mask)
ce_loss <- mx.symbol.MakeLoss(0 - mean_ce_loss, name = 'ce_loss')

指定優化器後就能直接訓練：

my_optimizer <- mx.opt.create(name = "sgd", learning.rate = 0.005, momentum = 0.9, wd = 1e-4)

my.eval.metric.loss <- mx.metric.custom(
  name = "ce-loss", 
  function(real, pred) {
    return(as.array(pred))
  }
)

mx.set.seed(0)

model <- mx.model.FeedForward.create(symbol = ce_loss, X = my_iter, optimizer = my_optimizer,
                                     eval.metric = my.eval.metric.loss,
                                     input.names = c('data'), output.names = c('label', 'mask'),
                                     array.batch.size = 16, ctx = mx.gpu(), num.round = 100)

model$symbol <- logistic_pred

輸出預測

讓我們試著做一張圖片的輸出看看：

# Generate image array

test_x_array <- array(0, dim = c(256, 256, 3, 1))
img_path <- paste0(img_dir, test_dat[4, 'PID'], '/', test_dat[4, 'UID'], '.png')
test_x_array[,,,1] <- readImage(img_path)

# Predicting

test_y_array <- predict(model, test_x_array, ctx = mx.gpu())

# Show image

par(mar = rep(0, 4))
plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(test_x_array[,,,1], 0, 1, 1, 0, interpolate=FALSE)

# Show annotation

for (k in 1:nrow(label_dat)) {
  current_col <- paste0(label_dat[k,'color'], 'A0')
  current_img <- (test_y_array[,,k,1] > 0.5) + 0L
  if (!table(current_img)['1'] %in% NA) {
    current_img <- c('#00000000', current_col)[current_img + 1]
    current_img <- array(current_img, dim = c(256, 256))
    rasterImage(current_img, 0, 1, 1, 0, interpolate=FALSE)
  }
}

在產生答案前，我們需要先知道我們的label格式，這是病人「P012」的「U004」：

head(submit_dat)

##    PID  UID row col   label
## 1 P012 U004 149 128 medulla
## 2 P012 U004 150 128 medulla
## 3 P012 U004 151 128 medulla
## 4 P012 U004 147 129 medulla
## 5 P012 U004 148 129 medulla
## 6 P012 U004 149 129 medulla

讓我們把它圖像化表示一下，當然這個標籤我們有做了一些隨機的處理，所以位置是錯的：

library(OpenImageR)
library(imager)

# Settings

img_dir <- 'image/'

# Read data

label_dat <- fread('label/label_list.csv', data.table = FALSE)
img_path <- paste0(img_dir, 'P012/U004.png')
img <- readImage(img_path)

# Show image

par(mar = rep(0, 4))
plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(img, 0, 1, 1, 0, interpolate=FALSE)

# Show annotation

obj_names <- unique(submit_dat[,'label'])

for (k in 1:length(obj_names)) {
  current_col <- paste0(label_dat[label_dat[,'structure'] %in% obj_names[k],'color'], 'A0')
  current_img <- matrix('#00000000', nrow = 256, ncol = 256)
  pos_matrix <- submit_dat[submit_dat[,'label'] %in% obj_names[k], c('row', 'col')] %>% as.matrix()
  current_img[pos_matrix] <- current_col
  rasterImage(current_img, 0, 1, 1, 0, interpolate=FALSE)
}

讓我們對所有預測資料進行預測，並將結果儲存成這樣的格式：

# Generate image array

test_x_array <- array(0, dim = c(256, 256, 3, nrow(test_dat)))

for (i in 1:nrow(test_dat)) {
  img_path <- paste0(img_dir, test_dat[i, 'PID'], '/', test_dat[i, 'UID'], '.png')
  test_x_array[,,,i] <- readImage(img_path)
}

# Predicting

test_y_array <- predict(model, test_x_array, ctx = mx.gpu())

# Processing

cut_p <- 0.5

submit_dat_list <- list()

for (i in 1:nrow(test_dat)) {
  
  for (j in 1:nrow(label_dat)) {
    
    sub_y_array <- test_y_array[,,j,i]
    pos_info <- which(sub_y_array > cut_p, arr.ind = TRUE)
    
    if (nrow(pos_info) > 0) {
      
      sub_submit_dat <- data.frame(PID = test_dat[i,'PID'],
                                   UID = test_dat[i,'UID'],
                                   row = pos_info[,'row'],
                                   col = pos_info[,'col'],
                                   label = label_dat[j,'structure'],
                                   stringsAsFactors = FALSE)
      
      submit_dat_list[[length(submit_dat_list) + 1]] <- sub_submit_dat
      
    }
    
  }
  
}

my_submit_dat <- do.call('rbind', submit_dat_list)
write.csv(my_submit_dat, file = 'my_submission.csv', na = '', row.names = FALSE, quote = FALSE)

我們大致上能取得0.3151的成績，這並不算非常理想，確實有許多進步空間。

3D卷積

要改進這個圖像分割網路還有非常多能改善的空間，如：

更複雜的資料擴增(翻轉)
適當的切點
使用「驗證組」對訓練過程做出指引
交叉驗證
網路結構的修正，更深的網路等
多階段的學習率修正策略

這些你應該都有機會能自己進行修正，但我們這裡介紹一個非常重要的想法，那就是對於「電腦斷層」這種3維影像而言，我們其實可以設計3D卷積神經網路進行運算。
我們需要修正一下疊代器，讓他輸出3D圖像，注意他的維度：

library(mxnet)
library(abind)

# Settings

fold_pos <- 50L

img_dir <- 'image/'
process_dir <- 'processed/'
PIDs <- unique(train_dat[,'PID'])

neg_pixel <- c(461, 465, 238, 241, 51, 54, 29, 89, 133)
neg_pixel <- neg_pixel * fold_pos

# Read data

my_iterator_core <- function (batch_size, random_crop = 224, z_crop = 20) {
  
  batch <- 0
  batch_per_epoch <- ceiling(length(PIDs) / batch_size * 30 / z_crop)
  
  reset <- function() {batch <<- 0}
  
  iter.next <- function() {
    batch <<- batch+1
    if (batch > batch_per_epoch) {return(FALSE)} else {return(TRUE)}
  }
  
  value <- function() {
    batch_PID <- sample(PIDs, batch_size, replace = TRUE)
    batch_x_list <- list()
    batch_y_list <- list()
    batch_mask_list <- list()
    for (j in 1:batch_size) {
      load(paste0(process_dir, batch_PID[j], '.RData'))
      random_row <- sample(0:(dim(x_array)[1] - random_crop), 1)
      random_col <- sample(0:(dim(x_array)[2] - random_crop), 1)
      random_z <- sample(0:(dim(x_array)[3] - z_crop), 1)
      batch_x_list[[j]] <- x_array[random_row+1:random_crop,random_col+1:random_crop,random_z+1:z_crop,]
      batch_y_list[[j]] <- y_array[random_row+1:random_crop,random_col+1:random_crop,random_z+1:z_crop,]
      sample_mask_list <- list()
      for (k in 1:length(neg_pixel)) {
        sub_mask_array <- array(0L, dim = c(random_crop, random_crop, z_crop))
        neg_mask_pos <- sample(1:length(sub_mask_array), neg_pixel[k] * z_crop)
        sub_mask_array[neg_mask_pos] <- 1L
        pos_mask_pos <- which(y_array[random_row+1:random_crop,random_col+1:random_crop,random_z + 1:z_crop,k] == 1)
        sub_mask_array[pos_mask_pos] <- fold_pos
        mirror_mask_pos <- which(y_array[random_row+1:random_crop, (dim(y_array)[2] + 1) - (random_col+1:random_crop),random_z + 1:z_crop,k] == 1)
        sub_mask_array[mirror_mask_pos] <- fold_pos
        sample_mask_list[[k]] <- sub_mask_array
      }
      batch_mask_list[[j]] <- abind(sample_mask_list, along = 4)
    }
    batch_x <- abind(batch_x_list, along = 5)
    batch_y <- abind(batch_y_list, along = 5)
    batch_mask <- abind(batch_mask_list, along = 5)
    data <- mx.nd.array(batch_x)
    label <- mx.nd.array(batch_y)
    mask <- mx.nd.array(batch_mask)
    return(list(data = data, label = label, mask = mask))
  }
  
  return(list(reset = reset, iter.next = iter.next, value = value, batch_size = batch_size, batch = batch))
}

my_iterator_func <- setRefClass("Custom_Iter",
                                fields = c("iter", "batch_size", "random_crop", "z_crop"),
                                contains = "Rcpp_MXArrayDataIter",
                                methods = list(
                                  initialize = function (iter, batch_size = 8, random_crop = 224, z_crop = 20){
                                    .self$iter <- my_iterator_core(batch_size = batch_size, random_crop = random_crop, z_crop = z_crop)
                                    .self
                                  },
                                  value = function () {
                                    .self$iter$value()
                                  },
                                  iter.next = function () {
                                    .self$iter$iter.next()
                                  },
                                  reset = function () {
                                    .self$iter$reset()
                                  },
                                  finalize = function () {
                                  }
                                )
)

my_iter <- my_iterator_func(iter = NULL, batch_size = 1, random_crop = 224, z_crop = 20)

接著，讓我們指定模型結構。3D卷積其實並不難實現！

library(mxnet)

data <- mx.symbol.Variable('data')
bn_data <- mx.symbol.BatchNorm(data = data, fix.gamma = TRUE, name = 'bn_data')

down_conv1 <- mx.symbol.Convolution(data = bn_data, kernel = c(7, 7, 1), pad = c(3, 3, 0), num_filter = 32, no.bias = TRUE, name = 'down_conv1')
down_convbn1 <- mx.symbol.BatchNorm(data = down_conv1, fix.gamma = FALSE, name = 'down_convbn1')
down_convrelu1 <- mx.symbol.Activation(data = down_convbn1, act_type = "relu", name = 'down_convrelu1')
down_convpool1 <- mx.symbol.Pooling(data = down_convrelu1, pool_type = "avg", kernel = c(2, 2, 1), stride = c(2, 2, 1), name = 'down_convpool1')

down_conv2 <- mx.symbol.Convolution(data = down_convpool1, kernel = c(5, 5, 1), pad = c(2, 2, 0), num_filter = 64, no.bias = TRUE, name = 'down_conv2')
down_convbn2 <- mx.symbol.BatchNorm(data = down_conv2, fix.gamma = FALSE, name = 'down_convbn2')
down_convrelu2 <- mx.symbol.Activation(data = down_convbn2, act_type = "relu", name = 'down_convrelu2')
down_convpool2 <- mx.symbol.Pooling(data = down_convrelu2, pool_type = "avg", kernel = c(2, 2, 1), stride = c(2, 2, 1), name = 'down_convpool2')

down_conv3 <- mx.symbol.Convolution(data = down_convpool2, kernel = c(3, 3, 1), pad = c(1, 1, 0), num_filter = 128, no.bias = TRUE, name = 'down_conv3')
down_convbn3 <- mx.symbol.BatchNorm(data = down_conv3, fix.gamma = FALSE, name = 'down_convbn3')
down_convrelu3 <- mx.symbol.Activation(data = down_convbn3, act_type = "relu", name = 'down_convrelu3')
down_convpool3 <- mx.symbol.Pooling(data = down_convrelu3, pool_type = "avg", kernel = c(2, 2, 1), stride = c(2, 2, 1), name = 'down_convpool3')

down_conv4 <- mx.symbol.Convolution(data = down_convpool3, kernel = c(3, 3, 1), pad = c(1, 1, 0), num_filter = 256, no.bias = TRUE, name = 'down_conv4')
down_convbn4 <- mx.symbol.BatchNorm(data = down_conv4, fix.gamma = FALSE, name = 'down_convbn4')
down_convrelu4 <- mx.symbol.Activation(data = down_convbn4, act_type = "relu", name = 'down_convrelu4')
down_convpool4 <- mx.symbol.Pooling(data = down_convrelu4, pool_type = "avg", kernel = c(2, 2, 1), stride = c(2, 2, 1), name = 'down_convpool4')

down_conv5 <- mx.symbol.Convolution(data = down_convpool4, kernel = c(3, 3, 3), pad = c(1, 1, 1), num_filter = 512, no.bias = TRUE, name = 'down_conv5')
down_convbn5 <- mx.symbol.BatchNorm(data = down_conv5, fix.gamma = FALSE, name = 'down_convbn5')
down_convrelu5 <- mx.symbol.Activation(data = down_convbn5, act_type = "relu", name = 'down_convrelu5')
down_convpool5 <- mx.symbol.Pooling(data = down_convrelu5, pool_type = "avg", kernel = c(2, 2, 1), stride = c(2, 2, 1), name = 'down_convpool5')

down_conv6 <- mx.symbol.Convolution(data = down_convpool5, kernel = c(3, 3, 3), pad = c(1, 1, 1), num_filter = 512, no.bias = TRUE, name = 'down_conv6')
down_convbn6 <- mx.symbol.BatchNorm(data = down_conv6, fix.gamma = FALSE, name = 'down_convbn6')
down_convrelu6 <- mx.symbol.Activation(data = down_convbn6, act_type = "relu", name = 'down_convrelu6')

up_deconv1 <- mx.symbol.Deconvolution(data = down_convrelu6, kernel = c(2, 2, 1), stride = c(2, 2, 1), num_filter = 512, name = 'up_deconv1')
up_deconvbn1 <- mx.symbol.BatchNorm(data = up_deconv1, fix.gamma = FALSE, name = 'up_deconvbn1')
up_deconvrelu1 <- mx.symbol.Activation(data = up_deconvbn1, act_type = "relu", name = 'up_deconvrelu1')
up_conv1 <- mx.symbol.Convolution(data = up_deconvrelu1, kernel = c(3, 3, 3), pad = c(1, 1, 1), num_filter = 512, no.bias = TRUE, name = 'up_conv1')
up_convbn1 <- mx.symbol.BatchNorm(data = up_conv1, fix.gamma = FALSE, name = 'up_convbn1')
up_convrelu1 <- mx.symbol.Activation(data = up_convbn1, act_type = "relu", name = 'up_convrelu1')

up_concat2 <- mx.symbol.concat(data = list(up_convrelu1, down_convrelu5), num.args = 2, dim = 1, name = 'up_concat2')
up_deconv2 <- mx.symbol.Deconvolution(data = up_concat2, kernel = c(2, 2, 1), stride = c(2, 2, 1), num_filter = 256, name = 'up_deconv2')
up_deconvbn2 <- mx.symbol.BatchNorm(data = up_deconv2, fix.gamma = FALSE, name = 'up_deconvbn2')
up_deconvrelu2 <- mx.symbol.Activation(data = up_deconvbn2, act_type = "relu", name = 'up_deconvrelu2')
up_conv2 <- mx.symbol.Convolution(data = up_deconvrelu2, kernel = c(3, 3, 3), pad = c(1, 1, 1), num_filter = 256, no.bias = TRUE, name = 'up_conv2')
up_convbn2 <- mx.symbol.BatchNorm(data = up_conv2, fix.gamma = FALSE, name = 'up_convbn2')
up_convrelu2 <- mx.symbol.Activation(data = up_convbn2, act_type = "relu", name = 'up_convrelu2')

up_concat3 <- mx.symbol.concat(data = list(up_convrelu2, down_convrelu4), num.args = 2, dim = 1, name = 'up_concat3')
up_deconv3 <- mx.symbol.Deconvolution(data = up_concat3, kernel = c(2, 2, 1), stride = c(2, 2, 1), num_filter = 256, name = 'up_deconv3')
up_deconvbn3 <- mx.symbol.BatchNorm(data = up_deconv3, fix.gamma = FALSE, name = 'up_deconvbn3')
up_deconvrelu3 <- mx.symbol.Activation(data = up_deconvbn3, act_type = "relu", name = 'up_deconvrelu3')
up_conv3 <- mx.symbol.Convolution(data = up_deconvrelu3, kernel = c(3, 3, 1), pad = c(1, 1, 0), num_filter = 256, no.bias = TRUE, name = 'up_conv3')
up_convbn3 <- mx.symbol.BatchNorm(data = up_conv3, fix.gamma = FALSE, name = 'up_convbn3')
up_convrelu3 <- mx.symbol.Activation(data = up_convbn3, act_type = "relu", name = 'up_convrelu3')

up_concat4 <- mx.symbol.concat(data = list(up_convrelu3, down_convrelu3), num.args = 2, dim = 1, name = 'up_concat4')
up_deconv4 <- mx.symbol.Deconvolution(data = up_concat4, kernel = c(2, 2, 1), stride = c(2, 2, 1), num_filter = 256, name = 'up_deconv4')
up_deconvbn4 <- mx.symbol.BatchNorm(data = up_deconv4, fix.gamma = FALSE, name = 'up_deconvbn4')
up_deconvrelu4 <- mx.symbol.Activation(data = up_deconvbn4, act_type = "relu", name = 'up_deconvrelu4')
up_conv4 <- mx.symbol.Convolution(data = up_deconvrelu4, kernel = c(3, 3, 1), pad = c(1, 1, 0), num_filter = 288, no.bias = TRUE, name = 'up_conv4')
up_convbn4 <- mx.symbol.BatchNorm(data = up_conv4, fix.gamma = FALSE, name = 'up_convbn4')
up_convrelu4 <- mx.symbol.Activation(data = up_convbn4, act_type = "relu", name = 'up_convrelu4')

up_deconv5 <- mx.symbol.Deconvolution(data = up_convrelu4, kernel = c(2, 2, 1), stride = c(2, 2, 1), num_filter = 288, num_group = 9, name = 'up_deconv5')
up_deconvbn5 <- mx.symbol.BatchNorm(data = up_deconv5, fix.gamma = FALSE, name = 'up_deconvbn5')
up_deconvrelu5 <- mx.symbol.Activation(data = up_deconvbn5, act_type = "relu", name = 'up_deconvrelu5')
up_conv5 <- mx.symbol.Convolution(data = up_deconvrelu5, kernel = c(3, 3, 1), pad = c(1, 1, 0), num_filter = 288, num_group = 9, no.bias = TRUE, name = 'up_conv5')
up_convbn5 <- mx.symbol.BatchNorm(data = up_conv5, fix.gamma = FALSE, name = 'up_convbn5')
up_convrelu5 <- mx.symbol.Activation(data = up_convbn5, act_type = "relu", name = 'up_convrelu5')

linear_pred <- mx.symbol.Convolution(data = up_convrelu5, kernel = c(3, 3, 1), pad = c(1, 1, 0), num_filter = 9, name = 'linear_pred')
logistic_pred <- mx.symbol.Activation(data = linear_pred, act.type = 'sigmoid', name = 'logistic_pred')

# CE loss

label <- mx.symbol.Variable(name = 'label')
mask <- mx.symbol.Variable(name = 'mask')

eps <- 1e-8
ce_loss_pos <- mx.symbol.broadcast_mul(mx.symbol.log(logistic_pred + eps), label)
ce_loss_neg <- mx.symbol.broadcast_mul(mx.symbol.log(1 - logistic_pred + eps), 1 - label)
ce_loss_sum <- mx.symbol.broadcast_add(ce_loss_pos, ce_loss_neg)
ce_loss_weight <- mx.symbol.broadcast_mul(ce_loss_sum, mask)
sum_mask <- mx.symbol.sum(mask, axis = 0:4)
sum_ce_loss <- mx.symbol.sum(ce_loss_weight, axis = 0:4)
mean_ce_loss <- mx.symbol.broadcast_div(sum_ce_loss, sum_mask)
ce_loss <- mx.symbol.MakeLoss(0 - mean_ce_loss, name = 'ce_loss')

指定優化器後就能直接訓練：

my_optimizer <- mx.opt.create(name = "sgd", learning.rate = 0.005, momentum = 0.9, wd = 1e-4)

my.eval.metric.loss <- mx.metric.custom(
  name = "ce-loss", 
  function(real, pred) {
    return(as.array(pred))
  }
)

mx.set.seed(0)

model <- mx.model.FeedForward.create(symbol = ce_loss, X = my_iter, optimizer = my_optimizer,
                                     eval.metric = my.eval.metric.loss,
                                     input.names = c('data'), output.names = c('label', 'mask'),
                                     array.batch.size = 1, ctx = mx.gpu(), num.round = 100)

model$symbol <- logistic_pred

訓練完成後，讓我們試著做一張圖片的輸出看看：

# Generate image array

test_PIDs <- unique(test_dat[,'PID'])

use_Patient <- 1

sub_test_dat <- test_dat[test_dat[,'PID'] %in% test_PIDs[use_Patient],]
test_x_array <- array(0, dim = c(256, 256, nrow(sub_test_dat), 3, 1))

for (j in 1:nrow(sub_test_dat)) {

  img_path <- paste0(img_dir, test_PIDs[use_Patient], '/', sub_test_dat[j, 'UID'], '.png')
  test_x_array[,,j,,] <- readImage(img_path)

}

# Predicting

test_y_array <- predict(model, test_x_array, ctx = mx.gpu())

# Show image

use_cut <- 4

par(mar = rep(0, 4), mfrow = c(1, 1))
plot(NA, xlim = c(0.04, 0.96), ylim = c(0.96, 0.04), axes = FALSE)
rasterImage(test_x_array[,,use_cut,,1], 0, 1, 1, 0, interpolate=FALSE)

# Show annotation

for (k in 1:nrow(label_dat)) {
  current_col <- paste0(label_dat[k,'color'], 'A0')
  current_img <- (test_y_array[,,use_cut,k,1] > 0.5) + 0L
  if (!table(current_img)['1'] %in% NA) {
    current_img <- c('#00000000', current_col)[current_img + 1]
    current_img <- array(current_img, dim = c(256, 256))
    rasterImage(current_img, 0, 1, 1, 0, interpolate=FALSE)
  }
}

讓我們對所有預測資料進行預測，並將結果儲存成剛剛的格式：

# Get predictions

cut_p <- 0.5
test_PIDs <- unique(test_dat[,'PID'])

submit_dat_list <- list()

for (i in 1:length(test_PIDs)) {

  # Generate image array

  sub_test_dat <- test_dat[test_dat[,'PID'] %in% test_PIDs[i],]
  test_x_array <- array(0, dim = c(256, 256, nrow(sub_test_dat), 3, 1))

  for (j in 1:nrow(sub_test_dat)) {

    img_path <- paste0(img_dir, test_PIDs[i], '/', sub_test_dat[j,'UID'], '.png')
    test_x_array[,,j,,] <- readImage(img_path)

  }

  # Predicting

  test_y_array <- predict(model, test_x_array, ctx = mx.gpu())

  # Processing

  for (j in 1:nrow(sub_test_dat)) {

    for (k in 1:nrow(label_dat)) {

      sub_y_array <- test_y_array[,,j,k,1]
      pos_info <- which(sub_y_array > cut_p, arr.ind = TRUE)

      if (nrow(pos_info) > 0) {

        sub_submit_dat <- data.frame(PID = sub_test_dat[1,'PID'],
                                     UID = sub_test_dat[j,'UID'],
                                     row = pos_info[,'row'],
                                     col = pos_info[,'col'],
                                     label = label_dat[k,'structure'],
                                     stringsAsFactors = FALSE)

        submit_dat_list[[length(submit_dat_list) + 1]] <- sub_submit_dat

      }

    }

  }

}

# Write out

my_submit_dat <- do.call('rbind', submit_dat_list)
write.csv(my_submit_dat, file = 'my_submission.csv', na = '', row.names = FALSE, quote = FALSE)

我們大致上能取得0.3580的成績，相較於2D卷積網路有很大的提升，你可以基於3D卷積網路進一步的進行優化！

Lesson 7 電腦斷層影像物件分割挑戰

下載資料

讀取資料

檔案格式

標籤格式

標籤編碼

模型訓練

輸出預測

3D卷積