第一節：物體分割任務(3)

我們先把檔案讀進來，這是Training set的圖像：

library(imager)
library(abind)
library(jpeg)

train_img_list <- list()
train_files <- list.files('ISBI/train-volume', pattern = '.jpg', full.names = TRUE)

for (i in 1:length(train_files)) {
  train_img_list[[i]] <- readJPEG(train_files[i])
}

train.x <- abind(train_img_list, along = 3)
dim(train.x) <- c(512, 512, 1, 30)

par(mar=rep(0,4), mfcol = c(2, 5))
for (i in 1:10) {
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  img = as.raster(train.x[,,,i])
  rasterImage(img, 0, 0, 1, 1, interpolate=FALSE)
}

第一節：物體分割任務(4)

我們先把檔案讀進來，這是Training set的標註資訊，要特別注意的是由於jpg檔案應該是有損壓縮，所以讀回來的值會因為壓縮的問題導致必須重新二值化：

train_label_list <- list()
train_files <- list.files('ISBI/train-labels', pattern = '.jpg', full.names = TRUE)

for (i in 1:length(train_files)) {
  train_label_list[[i]] <- readJPEG(train_files[i])
}

train.y <- abind(train_label_list, along = 3)
dim(train.y) <- c(512, 512, 1, 30)

train.y[train.y > 0.5] <- 1
train.y[train.y <= 0.5] <- 0

par(mar=rep(0,4), mfcol = c(2, 5))
for (i in 1:10) {
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  img = as.raster(train.y[,,,i])
  rasterImage(img, 0, 0, 1, 1, interpolate=FALSE)
}

第一節：物體分割任務(5)

接著就是編寫Iterator：

library(mxnet)

my_iterator_core = function(batch_size) {
  
  batch = 0
  batch_per_epoch = dim(train.y)[4]/batch_size
  
  reset = function() {batch <<- 0}
  
  iter.next = function() {
    batch <<- batch+1
    if (batch > batch_per_epoch) {return(FALSE)} else {return(TRUE)}
  }
  
  value = function() {
    idx <- 1:batch_size + (batch - 1) * batch_size
    idx[idx > dim(train.y)[4]] <- sample(1:dim(train.y)[4], sum(idx > dim(train.y)[4]))
    data <- mx.nd.array(array(train.x[,,,idx], dim = c(dim(train.x)[1:3], batch_size)))
    label <- mx.nd.array(array(train.y[,,,idx], dim = c(dim(train.y)[1:3], batch_size)))
    return(list(data = data, label = label))
  }
  
  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"),
                                contains = "Rcpp_MXArrayDataIter",
                                methods = list(
                                  initialize = function(iter, batch_size = 100){
                                    .self$iter <- my_iterator_core(batch_size = batch_size)
                                    .self
                                  },
                                  value = function(){
                                    .self$iter$value()
                                  },
                                  iter.next = function(){
                                    .self$iter$iter.next()
                                  },
                                  reset = function(){
                                    .self$iter$reset()
                                  },
                                  finalize=function(){
                                  }
                                )
)

嘗試使用一下Iterator：

my_iter <- my_iterator_func(iter = NULL, batch_size = 2)
my_iter$reset()
my_iter$iter.next()

## [1] TRUE

batch_data <- my_iter$value()

par(mar=rep(0,4), mfcol = c(2, 2))
for (i in 1:2) {
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(as.array(batch_data$data)[,,,i]), 0, 0, 1, 1, interpolate=FALSE)
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(as.array(batch_data$label)[,,,i]), 0, 0, 1, 1, interpolate=FALSE)
}

第一節：物體分割任務(6)

接著編寫Model architecture(這個網路有點深，我們需要Batch Normalization的輔助)：

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

conv1 <- mx.symbol.Convolution(data = bn_data, kernel = c(3, 3), pad = c(1, 1), num_filter = 8, no.bias = TRUE, name = 'conv1')
bn1 <- mx.symbol.BatchNorm(data = conv1, fix.gamma = FALSE, name = 'bn1')
relu1 <- mx.symbol.Activation(data = bn1, act_type = "relu", name = 'relu1')
pool1 <- mx.symbol.Pooling(data = relu1, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool1')

conv2 <- mx.symbol.Convolution(data = pool1, kernel = c(3, 3), pad = c(1, 1), num_filter = 16, no.bias = TRUE, name = 'conv2')
bn2 <- mx.symbol.BatchNorm(data = conv2, fix.gamma = FALSE, name = 'bn2')
relu2 <- mx.symbol.Activation(data = bn2, act_type = "relu", name = 'relu2')
pool2 <- mx.symbol.Pooling(data = relu2, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool2')

conv3 <- mx.symbol.Convolution(data = pool2, kernel = c(3, 3), pad = c(1, 1), num_filter = 32, no.bias = TRUE, name = 'conv3')
bn3 <- mx.symbol.BatchNorm(data = conv3, fix.gamma = FALSE, name = 'bn3')
relu3 <- mx.symbol.Activation(data = bn3, act_type = "relu", name = 'relu3')
pool3 <- mx.symbol.Pooling(data = relu3, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool3')

conv4 <- mx.symbol.Convolution(data = pool3, kernel = c(3, 3), pad = c(1, 1), num_filter = 64, no.bias = TRUE, name = 'conv4')
bn4 <- mx.symbol.BatchNorm(data = conv4, fix.gamma = FALSE, name = 'bn4')
relu4 <- mx.symbol.Activation(data = bn4, act_type = "relu", name = 'relu4')
pool4 <- mx.symbol.Pooling(data = relu4, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool4')

conv5 <- mx.symbol.Convolution(data = pool4, kernel = c(3, 3), pad = c(1, 1), num_filter = 128, no.bias = TRUE, name = 'conv5')
bn5 <- mx.symbol.BatchNorm(data = conv5, fix.gamma = FALSE, name = 'bn5')
relu5 <- mx.symbol.Activation(data = bn5, act_type = "relu", name = 'relu5')

deconv6 <- mx.symbol.Deconvolution(data = relu5, kernel = c(2, 2), stride = c(2, 2), num_filter = 64, name = 'deconv6')
bn6 <- mx.symbol.BatchNorm(data = deconv6, fix.gamma = FALSE, name = 'bn6')
relu6 <- mx.symbol.Activation(data = bn6, act_type = "relu", name = 'relu6')

deconv7 <- mx.symbol.Deconvolution(data = relu6, kernel = c(2, 2), stride = c(2, 2), num_filter = 32, name = 'deconv7')
bn7 <- mx.symbol.BatchNorm(data = deconv7, fix.gamma = FALSE, name = 'bn7')
relu7 <- mx.symbol.Activation(data = bn7, act_type = "relu", name = 'relu7')

deconv8 <- mx.symbol.Deconvolution(data = relu7, kernel = c(2, 2), stride = c(2, 2), num_filter = 16, name = 'deconv8')
bn8 <- mx.symbol.BatchNorm(data = deconv8, fix.gamma = FALSE, name = 'bn8')
relu8 <- mx.symbol.Activation(data = bn8, act_type = "relu", name = 'relu8')

deconv9 <- mx.symbol.Deconvolution(data = relu8, kernel = c(2, 2), stride = c(2, 2), num_filter = 8, name = 'deconv9')
bn9 <- mx.symbol.BatchNorm(data = deconv9, fix.gamma = FALSE, name = 'bn9')
relu9 <- mx.symbol.Activation(data = bn9, act_type = "relu", name = 'relu9')

linear_pred <- mx.symbol.Convolution(data = relu9, kernel = c(1, 1), num_filter = 1, 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')

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_mean <- 0 - mx.symbol.mean(ce_loss_pos + ce_loss_neg, axis = 0:3)
ce_loss <- mx.symbol.MakeLoss(ce_loss_mean, name = 'ce_loss')

定義Optimizer：

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

第一節：物體分割任務(7)

開始訓練：

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

mx.set.seed(0)

model_1 <- mx.model.FeedForward.create(symbol = ce_loss, X = my_iter, optimizer = my_optimizer,
                                       eval.metric = my.eval.metric.loss,
                                       array.batch.size = 2, ctx = mx.gpu(), num.round = 20)

把Testing set的Data讀進來，讓我們看看結果：

test_img_list <- list()
test_files <- list.files('ISBI/test-volume', pattern = '.jpg', full.names = TRUE)

for (i in 1:length(test_files)) {
  test_img_list[[i]] <- readJPEG(test_files[i])
}

test.x <- abind(test_img_list, along = 3)
dim(test.x) <- c(512, 512, 1, 30)

看看模型的預測結果：

model_1$symbol <- logistic_pred
pred_y.1 <- predict(model_1, test.x)
pred_y.1[pred_y.1 > 0.5] <- 1
pred_y.1[pred_y.1 <= 0.5] <- 0

par(mar = rep(0, 4), mfcol = c(3, 4))
for (i in 1:4) {
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(test.x[,,,i]), 0, 0, 1, 1, interpolate = FALSE)
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  pred_img <- pred_y.1[,,,i]
  pred_img[pred_img == 0] <- '#0000FF80'
  pred_img[pred_img == 1] <- '#FFFFFF00'
  rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(test.x[,,,i]), 0, 0, 1, 1, interpolate = FALSE)
  rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)
}

你覺得效果如何?

練習1：U-Net

物體分割用自編碼器的結構已經做的非常好了，但我們趕快查一查其他人是怎樣做這件事情的。我們很快的查到了在這個任務上最出名的是Olaf Ronneberger、Philipp Fischer以及Thomas Brox於2015年提出的研究：U-Net: Convolutional Networks for Biomedical Image Segmentation。
下面是Paper中U-Net的結構，他的重點在於除了設計直筒式的網路結構之外，更在網路結構中增加了相對應圖像大小的直通通道：

F05

這是我們第一次看著Paper實現他們所提出的演算法，試著做做看吧！記住，重要的是概念的運用！

練習1答案(1)

讓我們根據Paper提出的概念重新定義Model architecture：

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

conv1 <- mx.symbol.Convolution(data = bn_data, kernel = c(3, 3), pad = c(1, 1), num_filter = 8, no.bias = TRUE, name = 'conv1')
bn1 <- mx.symbol.BatchNorm(data = conv1, fix.gamma = FALSE, name = 'bn1')
relu1 <- mx.symbol.Activation(data = bn1, act_type = "relu", name = 'relu1')
pool1 <- mx.symbol.Pooling(data = relu1, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool1')

conv2 <- mx.symbol.Convolution(data = pool1, kernel = c(3, 3), pad = c(1, 1), num_filter = 16, no.bias = TRUE, name = 'conv2')
bn2 <- mx.symbol.BatchNorm(data = conv2, fix.gamma = FALSE, name = 'bn2')
relu2 <- mx.symbol.Activation(data = bn2, act_type = "relu", name = 'relu2')
pool2 <- mx.symbol.Pooling(data = relu2, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool2')

conv3 <- mx.symbol.Convolution(data = pool2, kernel = c(3, 3), pad = c(1, 1), num_filter = 32, no.bias = TRUE, name = 'conv3')
bn3 <- mx.symbol.BatchNorm(data = conv3, fix.gamma = FALSE, name = 'bn3')
relu3 <- mx.symbol.Activation(data = bn3, act_type = "relu", name = 'relu3')
pool3 <- mx.symbol.Pooling(data = relu3, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool3')

conv4 <- mx.symbol.Convolution(data = pool3, kernel = c(3, 3), pad = c(1, 1), num_filter = 64, no.bias = TRUE, name = 'conv4')
bn4 <- mx.symbol.BatchNorm(data = conv4, fix.gamma = FALSE, name = 'bn4')
relu4 <- mx.symbol.Activation(data = bn4, act_type = "relu", name = 'relu4')
pool4 <- mx.symbol.Pooling(data = relu4, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool4')

conv5 <- mx.symbol.Convolution(data = pool4, kernel = c(3, 3), pad = c(1, 1), num_filter = 128, no.bias = TRUE, name = 'conv5')
bn5 <- mx.symbol.BatchNorm(data = conv5, fix.gamma = FALSE, name = 'bn5')
relu5 <- mx.symbol.Activation(data = bn5, act_type = "relu", name = 'relu5')

deconv6 <- mx.symbol.Deconvolution(data = relu5, kernel = c(2, 2), stride = c(2, 2), num_filter = 64, name = 'deconv6')
bn6 <- mx.symbol.BatchNorm(data = deconv6, fix.gamma = FALSE, name = 'bn6')
relu6 <- mx.symbol.Activation(data = bn6, act_type = "relu", name = 'relu6')

concat7 <- mx.symbol.concat(data = list(relu6, relu4), num.args = 2, dim = 1, name = 'concat7')
deconv7 <- mx.symbol.Deconvolution(data = concat7, kernel = c(2, 2), stride = c(2, 2), num_filter = 32, name = 'deconv7')
bn7 <- mx.symbol.BatchNorm(data = deconv7, fix.gamma = FALSE, name = 'bn7')
relu7 <- mx.symbol.Activation(data = bn7, act_type = "relu", name = 'relu7')

concat8 <- mx.symbol.concat(data = list(relu7, relu3), num.args = 2, dim = 1, name = 'concat8')
deconv8 <- mx.symbol.Deconvolution(data = concat8, kernel = c(2, 2), stride = c(2, 2), num_filter = 16, name = 'deconv8')
bn8 <- mx.symbol.BatchNorm(data = deconv8, fix.gamma = FALSE, name = 'bn8')
relu8 <- mx.symbol.Activation(data = bn8, act_type = "relu", name = 'relu8')

concat9 <- mx.symbol.concat(data = list(relu8, relu2), num.args = 2, dim = 1, name = 'concat9')
deconv9 <- mx.symbol.Deconvolution(data = concat9, kernel = c(2, 2), stride = c(2, 2), num_filter = 8, name = 'deconv9')
bn9 <- mx.symbol.BatchNorm(data = deconv9, fix.gamma = FALSE, name = 'bn9')
relu9 <- mx.symbol.Activation(data = bn9, act_type = "relu", name = 'relu9')

concat10 <- mx.symbol.concat(data = list(relu9, relu1), num.args = 2, dim = 1, name = 'concat10')
linear_pred <- mx.symbol.Convolution(data = concat10, kernel = c(1, 1), num_filter = 1, 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')

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_mean <- 0 - mx.symbol.mean(ce_loss_pos + ce_loss_neg, axis = 0:3)
ce_loss <- mx.symbol.MakeLoss(ce_loss_mean, name = 'ce_loss')

練習1答案(2)

開始訓練：

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

mx.set.seed(0)

model_2 <- mx.model.FeedForward.create(symbol = ce_loss, X = my_iter, optimizer = my_optimizer,
                                       eval.metric = my.eval.metric.loss,
                                       array.batch.size = 2, ctx = mx.gpu(), num.round = 20)

使用model_2做預測：

model_2$symbol <- logistic_pred
pred_y.2 <- predict(model_2, test.x)
pred_y.2[pred_y.2 > 0.5] <- 1
pred_y.2[pred_y.2 <= 0.5] <- 0

讓我們比較一下兩個模型的預測結果(左邊是自編碼器的結構、右邊是U-Net)：

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

for (i in 1:2) {
  
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(test.x[,,,i]), 0, 0, 1, 1, interpolate = FALSE)
  
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  pred_img <- pred_y.1[,,,i]
  pred_img[pred_img == 0] <- '#0000FF80'
  pred_img[pred_img == 1] <- '#FFFFFF00'
  rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)
  
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(test.x[,,,i]), 0, 0, 1, 1, interpolate = FALSE)
  rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)
  
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(test.x[,,,i]), 0, 0, 1, 1, interpolate = FALSE)
  
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  pred_img <- pred_y.2[,,,i]
  pred_img[pred_img == 0] <- '#0000FF80'
  pred_img[pred_img == 1] <- '#FFFFFF00'
  rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)
  
  plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
  rasterImage(as.raster(test.x[,,,i]), 0, 0, 1, 1, interpolate = FALSE)
  rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)
  
}

看Paper的效果如何?多試試不同圖片的預測結果吧！

練習1引申(1)

物體分割最有名的研究其實是Jonathan Long、Evan Shelhamer以及Trevor Darrell於2015年提出的研究：Fully Convolutional Networks for Semantic Segmentation

F06

基本上這篇研究的Model architecture其實與自編碼器的結構大同小異，但融合了一些U-Net的特點，比較特別的是論文中比較了幾種不同的上採樣結構：

F07

最終，研究者認為上採樣3次的結構最理想：

F08

讓我們實驗看看結果如何！

練習1引申(2)

這邊會用到一個新的函數「mx.symbol.UpSampling」，讓我們有限度的修改我們的結構以增加可比較性：

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

conv1 <- mx.symbol.Convolution(data = bn_data, kernel = c(3, 3), pad = c(1, 1), num_filter = 8, no.bias = TRUE, name = 'conv1')
bn1 <- mx.symbol.BatchNorm(data = conv1, fix.gamma = FALSE, name = 'bn1')
relu1 <- mx.symbol.Activation(data = bn1, act_type = "relu", name = 'relu1')
pool1 <- mx.symbol.Pooling(data = relu1, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool1')

conv2 <- mx.symbol.Convolution(data = pool1, kernel = c(3, 3), pad = c(1, 1), num_filter = 16, no.bias = TRUE, name = 'conv2')
bn2 <- mx.symbol.BatchNorm(data = conv2, fix.gamma = FALSE, name = 'bn2')
relu2 <- mx.symbol.Activation(data = bn2, act_type = "relu", name = 'relu2')
pool2 <- mx.symbol.Pooling(data = relu2, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool2')

conv3 <- mx.symbol.Convolution(data = pool2, kernel = c(3, 3), pad = c(1, 1), num_filter = 32, no.bias = TRUE, name = 'conv3')
bn3 <- mx.symbol.BatchNorm(data = conv3, fix.gamma = FALSE, name = 'bn3')
relu3 <- mx.symbol.Activation(data = bn3, act_type = "relu", name = 'relu3')
pool3 <- mx.symbol.Pooling(data = relu3, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool3')

conv4 <- mx.symbol.Convolution(data = pool3, kernel = c(3, 3), pad = c(1, 1), num_filter = 64, no.bias = TRUE, name = 'conv4')
bn4 <- mx.symbol.BatchNorm(data = conv4, fix.gamma = FALSE, name = 'bn4')
relu4 <- mx.symbol.Activation(data = bn4, act_type = "relu", name = 'relu4')
up4 <- mx.symbol.UpSampling(data = relu4, num_args = 1, scale = 2, sample_type = 'bilinear', num_filter = 64, name = 'up4')
pool4 <- mx.symbol.Pooling(data = relu4, pool_type = "max", kernel = c(2, 2), stride = c(2, 2), name = 'pool4')

conv5 <- mx.symbol.Convolution(data = pool4, kernel = c(3, 3), pad = c(1, 1), num_filter = 128, no.bias = TRUE, name = 'conv5')
bn5 <- mx.symbol.BatchNorm(data = conv5, fix.gamma = FALSE, name = 'bn5')
relu5 <- mx.symbol.Activation(data = bn5, act_type = "relu", name = 'relu5')
up5 <- mx.symbol.UpSampling(data = relu5, num_args = 1, scale = 4, sample_type = 'bilinear', num_filter = 128, name = 'up5')

concat6 <- mx.symbol.concat(data = list(relu3, up4, up5), num.args = 3, dim = 1, name = 'concat6')
deconv6 <- mx.symbol.Deconvolution(data = concat6, kernel = c(2, 2), stride = c(2, 2), num_filter = 16, name = 'deconv6')
bn6 <- mx.symbol.BatchNorm(data = deconv6, fix.gamma = FALSE, name = 'bn6')
relu6 <- mx.symbol.Activation(data = bn6, act_type = "relu", name = 'relu6')

deconv7 <- mx.symbol.Deconvolution(data = relu6, kernel = c(2, 2), stride = c(2, 2), num_filter = 8, name = 'deconv7')
bn7 <- mx.symbol.BatchNorm(data = deconv7, fix.gamma = FALSE, name = 'bn7')
relu7 <- mx.symbol.Activation(data = bn7, act_type = "relu", name = 'relu7')

linear_pred <- mx.symbol.Convolution(data = relu7, kernel = c(1, 1), num_filter = 1, 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')

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_mean <- 0 - mx.symbol.mean(ce_loss_pos + ce_loss_neg, axis = 0:3)
ce_loss <- mx.symbol.MakeLoss(ce_loss_mean, name = 'ce_loss')

練習1引申(3)

實驗看看結果如何吧：

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

mx.set.seed(0)

model_3 <- mx.model.FeedForward.create(symbol = ce_loss, X = my_iter, optimizer = my_optimizer,
                                       eval.metric = my.eval.metric.loss,
                                       array.batch.size = 2, ctx = mx.gpu(), num.round = 20)

使用model_3做預測：

model_3$symbol <- logistic_pred
pred_y.3 <- predict(model_3, test.x)
pred_y.3[pred_y.3 > 0.5] <- 1
pred_y.3[pred_y.3 <= 0.5] <- 0

讓我們比較一下三個模型的預測結果(左邊是自編碼器的結構、中邊是U-Net、右邊是FCN-4s)：

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

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
rasterImage(as.raster(test.x[,,,1]), 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
pred_img <- pred_y.1[,,,1]
pred_img[pred_img == 0] <- '#0000FF80'
pred_img[pred_img == 1] <- '#FFFFFF00'
rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
rasterImage(as.raster(test.x[,,,1]), 0, 0, 1, 1, interpolate = FALSE)
rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
rasterImage(as.raster(test.x[,,,1]), 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
pred_img <- pred_y.2[,,,1]
pred_img[pred_img == 0] <- '#0000FF80'
pred_img[pred_img == 1] <- '#FFFFFF00'
rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
rasterImage(as.raster(test.x[,,,1]), 0, 0, 1, 1, interpolate = FALSE)
rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
rasterImage(as.raster(test.x[,,,1]), 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
pred_img <- pred_y.3[,,,1]
pred_img[pred_img == 0] <- '#0000FF80'
pred_img[pred_img == 1] <- '#FFFFFF00'
rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)

plot(NA, xlim = 0:1, ylim = 0:1, xaxt = "n", yaxt = "n", bty = "n")
rasterImage(as.raster(test.x[,,,1]), 0, 0, 1, 1, interpolate = FALSE)
rasterImage(pred_img, 0, 0, 1, 1, interpolate = FALSE)

看來Paper也不見得說的都對，有時候還是要多做實驗比較一下！

第二節：物件識別模型的前期演進(1)

物體分割與物件識別兩個任務乍看之下好像物體分割更為困難，但實際上物體分割因為很容易可以想出預測目標的格式，因此物件識別其實難上非常多！

$F09$

– 給各位一個思考的時間，你能想像一下要怎樣，你要注意的是預測函數需要「固定」輸出數目，因此如何預測「不固定」數量的預測框其實就是一個大難題！

第二節：物件識別模型的前期演進(2)

目前為止我們對於圖像的分類任務應該還算是得心應手，一個比較簡單的想法就是透過滑動視窗(Sliding Windows)的概念搜索物體在哪。

– 舉例來說我們可以像下圖這樣對每一個候選框進行分類任務，看看這個框是否有包含我們希望預測的物件：

F10

在重複使用不同長寬比/大小的候選框預測後，可能會有許多重複的框都包含我們想要的物體，我們可以保留機率最高的並且移除與其重複過多的眶，如下：

但這個方法具有幾個明顯的問題，首先我們要搜索的框的數目非常多(框的長寬比及大小組合理論上是無限多的)，而將每一個框內的圖都使用卷積神經網路做預測顯然要浪費非常多的時間，因此這個方法原則上是不可行的。

第二節：物件識別模型的前期演進(3)

一個比較開創性的工作是Ross Girshick、Jeff Donahue、Trevor Darrell以及Jitendra Malik在2013年底所提出的Region CNN(R-CNN)，計算細節在論文Rich feature hierarchies for accurate object detection and semantic segmentation內有詳述，這應該是第一個使用卷積神經網路做物件識別任務的模型。

F12

剛剛滑動視窗搜索策略的邏輯本質上沒有錯，主要的問題其實是在於需要搜索的候選框數目是無限多的，因此解決的方式也很直覺：既然深度神經網路有能力做出任意預測，那何不除了預測存在物體的機率外，同時預測候選框的實際長寬及位置！

– 模型的概念如下圖所示：

透過這種方式，至少我們能夠減少候選框的數目到一個還算合理的數字(雖然還是非常多)，從而實現物件識別模型。

第二節：物件識別模型的前期演進(4)

我們這裡稍微說明一下RCNN的實現細節，他並不是真的使用了滑動視窗的搜索策略，他的候選框決定的方式是透過Selective Search，這是基於Felzenszwal於2004年發表的論文：Efficient Graph-Based Image Segmentation

– 他的概念大概是透過簡單的色階變換決定圖像內有多少可分割的物體，而RCNN就是透過這種方式決定候選框有哪些：

F14

你可以想像一下整個R-CNN的訓練過程有多麻煩，首先你要對你手上的每一張圖片做出候選框的預測，並且比較這些候選框是否能代表真實框，最後再把所有樣本拿去訓練一個CNN做分類，再把所有正向樣本拿去做迴歸任務預測框的偏移量，整個過程說多麻煩就有多麻煩。

– 更麻煩的是，在做預測的時候同樣要把整個步驟重新跑過一遍，可見他還有非常大的改進空間！

第二節：物件識別模型的前期演進(5)

Ross Girshick在2015年又發表了Fast R-CNN來改進前面R-CNN的缺點。

– 整個R-CNN最大的問題在於為什麼要做卷積神經網路的運算這麼多次，難道過程不能整合?

Fast R-CNN就是根據這個想法做出改進，他把候選框的圈選放到了一個卷積神經網路之後，然而使用一個叫ROI Pooling的方法整合候選框內的資訊，最後再使用另一個卷積神經網路做出候選框的類別及邊框預測：

F15

第二節：物件識別模型的前期演進(6)

Fast R-CNN跟R-CNN的主要不同在於下面3點：

用來提取特徵的卷積神經網絡是作用在整個圖片上，而不是各個候選框上，而且這個卷積網絡也參與訓練過程。
候選框的搜索是在卷積神經網絡的輸出上，而不是原始圖片上。
在R-CNN裡，我們將形狀各異的提議區域resize後使用同樣的形狀來進行特徵提取。而在Fast R-CNN中整個過程為了整合進神經網路的推理之中，引入了興趣區域池化層(Region of Interest Pooling，RoI Pooling)來resize每個候選區域。

F16

– 這裡要特別說明一下Fast R-CNN的候選框提取方法，由於相較於R-CNN速度加快了不少，Fast R-CNN預先定義了一系列非常多的候選框(或稱作錨框，anchor box)。

– 與R-CNN最大的不同是對於R-CNN而言，每張圖片的候選框都不一樣，而對於Fast R-CNN而言每張圖片的候選框都一樣。

這些改進除了讓推理速度加快許多之外，更重要的是整個過程可以被一個卷積神經網路給包含在內，從而簡化訓練過程。

第二節：物件識別模型的前期演進(7)

如果說Fast R-CNN還有值得改進的地方，那就是選擇候選框的方式還可以再更精簡一點！想不到甚麼更好的方法選框?那就把這個過程也交給卷積神經網絡吧！

– 在Fast R-CNN出來後不久，Ross Girshick又與他的同事Kaiming He(何愷明)、Shaoqing Ren以及Jian Sun發表了Faster R-CNN：Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks

讓我們看看從R-CNN、Fast R-CNN到Faster R-CNN一系列的演進：

F17

第二節：物件識別模型的前期演進(8)

Faster R-CNN的主要創新點在於Region Proposal Networks，這個網路主要是透過卷積網路進行全局掃描，並且預先定義好k個anchor boxes，對其進行預測，預測輸出為anchor boxes的偏移量：

F18

你可以把Faster R-CNN的Region Proposal Networks想像成是這樣的，那就是Fast R-CNN中是對每個預先定義好的anchor boxes都進行後續的分類以及偏移預測，由於任務相當複雜所以後面的網路是很大的，從而導致速度仍然很慢。

– 而在Faster R-CNN中，他把這個過程拆成兩個步驟，首先先使用Region Proposal Networks做相對簡單的任務(判斷anchor box是否存在以及預測偏移量)，之後再用後面的網路進行後續的分類以及偏移預測。

Faster R-CNN透過一系列的改進後，是目前為止在物件識別任務中精準度最高的模型。然而他的運算過程仍然相當繁瑣複雜(這是兩階段模型的通病)，從而導致檢測速度仍然過慢無法實際應用，因此之後的研究試圖把物件識別任務整合至同個卷積神經路中，在確保準確度接近Faster R-CNN的同時提升數倍的推理速度。

第三節：一階段物件識別模型介紹(1)

在物件識別模型上另一個開創性的工作當屬Joseph Redmon、Santosh Divvala、Ross Girshick以及Ali Farhadi在2016年發表的YOLO：You Only Look Once: Unified, Real-Time Object Detection
這是一個突破性的邏輯轉換，在此之前的所有物件識別模型網路，總是希望能夠先找出框大概的位置在哪，接著把該範圍內的圖像做ROI Pooling，之後再做出預測。

– 而YOLO的邏輯在於他不再試圖把框大概的位置找出來，而是直接找尋「物件中心」，再透過「物件中心」的資訊直接預測出框的大小及位置。

模型的邏輯在於，卷積神經網路最後的特徵圖一定能夠代表回原圖相對應的位置(這點我們在第9課Class Activation Mapping的實驗中已經說明了)，因此我們可以把原始圖像想像成這個樣子：

$F19$

那麼對於每一個Gird，他將直接預測出2個候選框，並且直接預測出他的機率：

F20

這邊有個問題給各位思考一下，根據你對卷積神經網路的了解，物件中心的Gird是否真的「有可能」預測一個遠大於他大小的物件框?他所擁有的資訊足夠嗎?

第三節：一階段物件識別模型介紹(2)

所以以剛剛那張圖片為例，負責狗(橘色)、腳踏車(綠色)、汽車(紅色)的Gird就分別如下：

那知道這是否是某個物件的中心之外，更重要的是還需要知道：

框的中心位置(x座標與y座標)
框的長度及寬度
該框是屬於哪一個類別的物件

所以對於每一個Gird，假設網路需要他每預測1個框，那他就必須提供[5(與框有關的資訊)+C(物件類別數)]個數值。

F22

怎麼對每一個Gird都做出預測呢?使用1×1的卷積核就可以囉！

第三節：一階段物件識別模型介紹(3)

YOLO相較於R-CNN系列的模型捨棄了Region Proposal的過程，而僅僅用了非常有限的資訊(中心的資訊)預測目標框的資訊以及物件類別，這當然不可避免的嚴重危害了準確度。因此後面的模型就試圖解決資訊不足的問題！
一個基於YOLO概念所提出的模型：SSD(Single Shot MultiBox Detector)，他試圖利用不同尺度的特徵圖做出類似於YOLO事情，成功的把準確度提升至與Faster R-CNN比肩的程度：

F23

由於不同大小的物件可以由不同尺度的Grid所負責，透過這種方式那一個Gird所擁有的資訊就較為足夠了，做出的預測也較為精準。

– 除了多尺度的預測之外，還有一個與YOLO不同的地方在於他不是憑空預測框的資訊，而是透過預先定義好一系列anchor boxes，再去預測目標框與anchor boxex的偏移量。

第三節：一階段物件識別模型介紹(4)

SSD的概念已經不錯了，如果還說有什麼問題那就是負責小物件的Gird因為為於網路的較淺層的地方，特徵抽取量可能不夠，而為了解決這個問題Facebook團隊提出了：特徵金字塔網路(Feature Pyramid Networks，FPN)。

– 該論文為Feature Pyramid Networks for Object Detection，看看作者群你會發現怎麼又是Ross Girshick以及Kaiming He

我們這裡直接介紹特徵金字塔網路的續作：RetinaNet，他們主要的差別是訓練過程的部分，整體結構是相當類似的：

F24

這個突破其實並不大，其實就是把反卷積器與Fully Convolutional Networks運用到這上面，但他也確實提升了大量的準確度：

F25

家庭作業：利用訓練好的YOLO模型實現物體檢測

物件識別模型由於正處於研究熱點，因此進展非常迅速，有任何好的idea馬上就會被運用(抄襲)並且被後起之秀超越，我們這裡使用一個近期集大成的物件識別模型：YOLO v3

– 這個模型是用COCO Detection Challenge中的資料集訓練而成的，抓取物件的標的共有80項，你可以到這裡下載完整的語法及模型，我們可以透過裡面的語法進行預測，多嘗試一下不同的圖片感受他的威力吧！

這週的回家作業就是透過任何管道了解YOLO v3是怎樣預測出邊框的，你必須一步一步參考裡面的語法了解它是如何做到的！

– 需要注意的是，這並不是原版的YOLO v3模型，這是一個較小的模型，準確度有略為降低，重要的是你是否有辦法從語法中、論文內、網路資源中找到這個最先進的模型是怎樣執行物件識別任務的！

家庭作業提示

這份語法非常的長，你除了該多看看paper之外(YOLO v3的PAPER)，你可以依據下列提示完成任務：

由於這個模型的預測式包含了3個部分(分別是下採用8倍、16倍以及32倍)，因此一般的predict函數沒有辦法使用，所以這個語法包含了一個「my_predict」函數，你必須優先了解它的輸出格式為何！
接著考慮到我們需要把輸出格式重新編碼成物件框的格式，我們需要讀取anchor_boxs (yolo v3).RData，這裡面記錄了9個anchor box的長寬資訊。有了這些anchor box之後，我們是利用下面這個方式進行解碼：

F26

接著注意Decode_fun之內，除了第一步解碼的部分，第二步是要把多餘的物件框移除(與高機率框重複太多的)，這裡我們移除的依據是使用IoU的大小作為依據，而IoU的定義如下，計算方法交給IoU_function：

F27

把預測輸出還原回框後，下一步就只是顯示圖片了，函數「Show_img」負責顯示圖片！

需要注意的是，你必須先學會使用Model，至於怎樣訓練這個Model那又是另一回事了。很明顯的針對這個任務我們要設計新的損失函數，所以有關訓練過程的部分可以通通跳過，先理解他的預測方式更為重要！

結語

今天的課程帶各位同學了解到物體分割與物件識別模型結構，其中物件識別模型我們主要講述了他的發展過程以及推理邏輯。

– 至於怎麼訓練一個物件識別模型，這還會遇到非常多的難題，並且非常考驗你的程式能力！

你可能會非常驚嘆於深度學習模型的潛力，還記得第一節課我們提到的目標嗎?我們希望能夠發展一個預測函數從而建立「任意$x$」與「任意$y$」的映射過程。這在最開始的時候是非常難想像的，但課程至此你應該會開始覺得這真的是有可能實現的了！
我們下一節課再帶領大家一步一步實現一個簡單任務：皮卡丘識別任務！由於下一節課的過程真的非常困難，我建議大家先下載這個Github上的範例：MxNetR-YOLO，試著先把整個皮卡丘識別的過程完整的跑過一輪，並了解其中每一行語法的意義：

F26

深度學習理論與實務

前言

第一節：物體分割任務(1)

第一節：物體分割任務(2)