Visualizing what CNNs learn

In this notebook¹, we’ll take a look at how how CNNs interpret images from a visual perspective. There are several approaches one could take but we’ll look at the following three:

• Visualizing intermediate CNN outputs (“intermediate activations”). This is useful to understand how successive CNN layers transform their input, and to get a first idea of the meaning of individual CNN filters.
• Visualizing CNN filters. This is useful to understand precisely what visual pattern or concept each filter in a CNN is receptive to.
• Visualizing superpixels and heatmaps of class activation in an image. This is useful to understand which part of an image were identified as belonging to a given class, and thus allows to localize objects in images.

For the first method, activation visualization, we’ll use the small CNN that we trained from scratch in the cat vs. dog module. For the other two methods, we will use the VGG16 model that we introduced in the CNN Transfer learning notebook.

Required packages

# Initialize package
library(keras)
library(tensorflow)
library(grid)         # for producing section 1 & 2 images
library(gridExtra)    # for producing section 1 & 2 images
library(magick)       # for producing section 3 images
library(viridis)      # for producing section 3 images
library(lime)         # for producing section 3 images

# We disable eager execution to extract the gradients later on. This isn't
# required as we could also get gradients with tf$GradientTape. See this
# discussion for details https://github.com/rstudio/keras/issues/945
tf$compat$v1$disable_eager_execution()

Visualizing intermediate activations

Visualizing intermediate activations consists of displaying the feature maps that are output by various convolution and pooling layers in a network. This gives a view into how an input is decomposed unto the different filters learned by the network. These feature maps we want to visualize have 3 dimensions: width, height, and depth (aka channels). Each channel encodes relatively independent features, so the proper way to visualize these feature maps is by independently plotting the contents of every channel, as a 2D image. Let’s start by loading the model that we saved in the cat vs. dog module.

file_local <- here::here("materials", "04-computer-vision-CNNs", "cats_and_dogs_small_2.h5")
model <- load_model_hdf5(file_local)
model

Model
Model: "sequential"
__________________________________________________________________________________________________
Layer (type)                                Output Shape                           Param #        
==================================================================================================
conv2d (Conv2D)                             (None, 148, 148, 32)                   896            
__________________________________________________________________________________________________
max_pooling2d (MaxPooling2D)                (None, 74, 74, 32)                     0              
__________________________________________________________________________________________________
conv2d_1 (Conv2D)                           (None, 72, 72, 64)                     18496          
__________________________________________________________________________________________________
max_pooling2d_1 (MaxPooling2D)              (None, 36, 36, 64)                     0              
__________________________________________________________________________________________________
conv2d_2 (Conv2D)                           (None, 34, 34, 128)                    73856          
__________________________________________________________________________________________________
max_pooling2d_2 (MaxPooling2D)              (None, 17, 17, 128)                    0              
__________________________________________________________________________________________________
conv2d_3 (Conv2D)                           (None, 15, 15, 128)                    147584         
__________________________________________________________________________________________________
max_pooling2d_3 (MaxPooling2D)              (None, 7, 7, 128)                      0              
__________________________________________________________________________________________________
flatten (Flatten)                           (None, 6272)                           0              
__________________________________________________________________________________________________
dropout (Dropout)                           (None, 6272)                           0              
__________________________________________________________________________________________________
dense (Dense)                               (None, 512)                            3211776        
__________________________________________________________________________________________________
dense_1 (Dense)                             (None, 1)                              513            
==================================================================================================
Total params: 3,453,121
Trainable params: 3,453,121
Non-trainable params: 0
__________________________________________________________________________________________________

This will be the input image we will use – a picture of a dog, not part of images that the network was trained on. This image is particularly interesting because it also has a picture of a human but if you were to predict the class of this image our model correctly predicts a dog.

# Preprocesses the image into a 4D tensor
img_path <- here::here("materials", "data", "dogs-vs-cats", "test", "dogs", "dog.1508.jpg")
img <- image_load(img_path, target_size = c(150, 150))
img_tensor <- image_to_array(img)
img_tensor <- array_reshape(img_tensor, c(1, 150, 150, 3))
img_tensor <- img_tensor / 255

dim(img_tensor)

[1]   1 150 150   3

plot(as.raster(img_tensor[1,,,]))

In order to extract the feature maps you want to look at, you’ll create a Keras model that takes batches of images as input, and outputs the activations of all convolution and pooling layers. To do this, we will use the keras_model() function (recall keras_model() was introduced in the collaborative filtering notebook.

layer_outputs <- lapply(model$layers[1:8], function(layer) layer$output)
activation_model <- keras_model(inputs = model$input, outputs = layer_outputs)

When fed an image input, this model returns the activation output values for each of the 8 hidden layers in the convolution section of our model.

activations <- activation_model %>% predict(img_tensor)

For instance, this is the activation of the first convolution layer for our dog image. Note how the dimension align to the dimensions of the first layer in our CNN model.

first_layer_activation <- activations[[1]]
dim(first_layer_activation)

[1]   1 148 148  32

It’s a 148 x 148 feature map with 32 channels. Let’s visualize some of these channels.

plot_channel <- function(channel) {
  rotate <- function(x) t(apply(x, 2, rev))
  image(rotate(channel), axes = FALSE, asp = 1,
        col = terrain.colors(12))
}

The 20th channel seems to capture the outline of the dog and also the person’s face.

plot_channel(first_layer_activation[1,,,20])

The 28th channel seems to outline some of the other features in the image ( i.e. window panels).

plot_channel(first_layer_activation[1,,,28])

We can go ahead and create images for each channel in each of the 8 layers of our convolution part of our model. The following will create a new subdirectory, plot every channel in each of our 8 activation maps, stack the results in one big image tensor, with channels stacked side-by-side plot, and save them in the new directory.

dir.create("dog_activations")
image_size <- 58
images_per_row <- 16
for (i in 1:8) {
  
  layer_activation <- activations[[i]]
  layer_name <- model$layers[[i]]$name
 
  n_features <- dim(layer_activation)[[4]]
  n_cols <- n_features %/% images_per_row
 
  png(paste0("dog_activations/", i, "_", layer_name, ".png"), 
      width = image_size * images_per_row, 
      height = image_size * n_cols)
  op <- par(mfrow = c(n_cols, images_per_row), mai = rep_len(0.02, 4))
  
  for (col in 0:(n_cols - 1)) {
    for (row in 0:(images_per_row - 1)) {
      channel_image <- layer_activation[1,,,(col*images_per_row) + row + 1]
      plot_channel(channel_image)
    }
  }
  
  par(op)
  dev.off()
}

You should now have plots visualizing the first 8 feature channels:

1. 1_conv2d.png
2. 2_max_pooling2d.png
3. 3_conv2d_1.png
4. 4_max_pooling2d_1.png
5. 5_conv2d_2.png
6. 6_max_pooling2d_2.png
7. 7_conv2d_3.png
8. 8_max_pooling2d_3.png

The following plots the feature maps in the first layer and then the final 8th layer. Note how the final layer is very abstract to human eyes…but it has a lot of mathematical information interpretable by our CNN model. The activations of layers higher-up carry less and less information about the specific input being seen, and more and more information about the target (in our case, the class of the image: cat or dog).

knitr::include_graphics("dog_activations/1_conv2d.png")

knitr::include_graphics("dog_activations/8_max_pooling2d_3.png")

Visualizing CNN filters:

Another easy thing to do to inspect the filters learned by convnets is to display the visual pattern that each filter is meant to respond to. This can be done with gradient ascent in input space: applying gradient descent to the value of the input image of a convnet so as to maximize the response of a specific filter, starting from a blank input image. The resulting input image would be one that the chosen filter is maximally responsive to.

The process is simple: we will build a loss function that maximizes the value of a given filter in a given convolution layer, then we will use stochastic gradient descent to adjust the values of the input image so as to maximize this activation value. For instance, here’s a loss for the activation of filter 0 in the layer “block3_conv1” of the VGG16 network, pre-trained on ImageNet:

model <- application_vgg16(
  weights = "imagenet",
  include_top = FALSE
  )

layer_name <- "block3_conv1"
filter_index <- 1
layer_output <- get_layer(model, layer_name)$output
loss <- k_mean(layer_output[,,,filter_index])

To implement gradient descent, we will need the gradient of this loss with respect to the model’s input. To do this, we will use the k_gradients() Keras backend function. In TF 2.0 a change was made where k_gradients() can only be ran when eager execution is turned off, which is why we ran tf$compat$v1$disable_eager_execution() at the top of the script. See https://github.com/rstudio/keras/issues/945 for more details.

# The call to k_gradients returns an R list of tensors (of size 1 in this case). 
# Hence, you keep only the first element—which is a tensor.
grads <- k_gradients(loss, model$input)[[1]]

A non-obvious trick to use for the gradient descent process to go smoothly is to normalize the gradient tensor, by dividing it by its L2 norm (the square root of the average of the square of the values in the tensor). This ensures that the magnitude of the updates done to the input image is always within a same range.

# Add 1e-5 before dividing to avoid accidentally dividing by 0
grads <- grads / (k_sqrt(k_mean(k_square(grads))) + 1e-5)

Now you need a way to compute the value of the loss tensor and the gradient tensor, given an input image. You can define a Keras backend function to do this: iterate() is a function that takes a tensor (as a list of tensors of size 1) and returns a list of two tensors: the loss value and the gradient value.

iterate <- k_function(list(model$input), list(loss, grads))
c(loss_value, grads_value) %<-% iterate(list(array(0, dim = c(1, 150, 150, 3))))

At this point we can define an R loop to do stochastic gradient descent:

# Starts from a gray image with some noise
input_img_data <- array(runif(150 * 150 * 3), dim = c(1, 150, 150, 3)) * 20 + 128

# Runs gradient ascent for 40 step
step <- 1
for (i in 1:40) { 
  # Computes the loss value and gradient value
  c(loss_value, grads_value) %<-% iterate(list(input_img_data))
  
  # Adjusts the input image in the direction that maximizes the loss
  input_img_data <- input_img_data + (grads_value * step)  
}

The resulting image tensor is a floating-point tensor of shape (1, 150, 150, 3), with values that may not be integers within [0, 255]. Hence you need to post-process this tensor to turn it into a displayable image. You do so with the following straightforward utility function.

deprocess_image <- function(x) {
  dms <- dim(x)

    # Normalizes the tensor: centers on 0., ensures that std is 0.1
  x <- x - mean(x)
  x <- x / (sd(x) + 1e-5)
  x <- x * 0.1
  
  # Clips to [0, 1]
  x <- x + 0.5
  x <- pmax(0, pmin(x, 1))
  
  # Returns with the original image dimensions
  array(x, dim = dms)
}

Now you have all the pieces. Let’s put them together into an R function that takes as input a layer name and a filter index, and returns a valid image tensor representing the pattern that maximizes the activation of the specified filter.

generate_pattern <- function(layer_name, filter_index, size = 150) {
  
  # Builds a loss function that maximizes the activation of the nth filter of the layer under consideration
  layer_output <- model$get_layer(layer_name)$output
  loss <- k_mean(layer_output[,,,filter_index])
  
  # Computes the gradient of the input picture with regard to this loss
  grads <- k_gradients(loss, model$input)[[1]]
  
  # Normalization trick: normalizes the gradient
  grads <- grads / (k_sqrt(k_mean(k_square(grads))) + 1e-5)
  
  # Returns the loss and grads given the input picture
  iterate <- k_function(list(model$input), list(loss, grads))
  
  # Starts from a gray image with some noise
  input_img_data <-
  array(runif(size * size * 3), dim = c(1, size, size, 3)) * 20 + 128

  # Runs gradient ascent for 40 steps
  step <- 1
  for (i in 1:40) {
    c(loss_value, grads_value) %<-% iterate(list(input_img_data))
    input_img_data <- input_img_data + (grads_value * step)
    }
  
  img <- input_img_data[1,,,]
  deprocess_image(img)
  }

Now we can start visualizing every single filter in every layer. For simplicity, we will only look at the first 64 filters in each layer, and will only look at the first layer of the following convolution blocks: block1_conv1, block2_conv1, block3_conv1, block4_conv1. We will arrange the outputs on a 8x8 grid of filter patterns.

layers <- c("block1_conv1", "block2_conv1", "block3_conv1", "block4_conv1")
dir.create("vgg_filters")

for (layer_name in layers) {
size <- 140
  png(paste0("vgg_filters/", layer_name, ".png"),
      width = 8 * size, height = 8 * size)
  grobs <- list()
  for (i in 0:7) {
    for (j in 0:7) {
      pattern <- generate_pattern(layer_name, i + (j*8) + 1, size = size)
      grob <- rasterGrob(pattern,
                         width = unit(0.9, "npc"),
                   height = unit(0.9, "npc"))
      grobs[[length(grobs) + 1]] <- grob
    }
    }
  grid.arrange(grobs = grobs, ncol = 8)
  dev.off()
}

The previous code results in a new directory /vgg_filters that contains four images. Below is the image for the first CNN block filter and the second image is the fourth block. These filter visualizations tell you a lot about how CNN layers see the world: each layer in a CNN learns a collection of filters such that their inputs can be expressed as a combination of the filters. You can see that earlier filters (i.e. block1_conv1) encode more simple directional edges and colors whereas later filters (i.e. block4_conv1) start to resemble textures found in natural images (i.e. eyes, feathers, balls).

knitr::include_graphics(c("vgg_filters/block1_conv1.png", "vgg_filters/block4_conv1.png"))

Visualizing heatmaps of class activation

This last technique demonstrates how to highlight the parts of an image that is likely to be the most important in how the model classified the image. For example, in our dog image, which contains a human, we can see that our model classifies the image, with 0.89 probability, as a dog (in this model the target is encoded as 0-cat and 1-dog). So one question may arise, what about this image causes our model to predict this? What pixels in the image are most relevant?

# predicting dog picture
predict(model, img_tensor)

          [,1]
[1,] 0.8906708

Manual class activation map

A class-activation heatmap is a 2D grid of scores associated with a specific output class, computed for every location in any input image, indicating how important each location is with respect to the class under consideration. For instance, given an image fed into a cat-versus-dog CNN, CAM visualization allows you to generate a heatmap for the class “cat,” indicating how cat-like different parts of the image are, and also a heatmap for the class “dog,” indicating how dog-like parts of the image are.

The specific implementation you’ll use is the one described in “Grad-CAM: Visual Explanations from Deep Networks via Gradient-based Localization.”footnote:[Ramprasaath R. Selvaraju et al., Cornell University Library, March 21, 2017, https://arxiv.org/abs/1610.02391.]. It consists of taking the output feature map of a convolution layer, given an input image, and weighing every channel in that feature map by the gradient of the class with respect to the channel. Intuitively, one way to understand this trick is that you’re weighting a spatial map of “how intensely the input image activates different channels” by “how important each channel is with regard to the class,” resulting in a spatial map of “how intensely the input image activates the class.”

Let’s demonstrate with our dog image.

# This is the prediction vector for our dog
dog_output <- model$output[, 1]

# The is the output feature map of the `conv2d_3` layer, the last
# convolutional layer in our Cats vs Dogs model that we imported earlier
last_conv_layer <- model %>% get_layer("conv2d_3")

# This is the gradient of the "dog" class with regard to the output
# feature map of `conv2d_3`
grads <- k_gradients(dog_output, last_conv_layer$output)[[1]]

# This is a vector of shape (128,), where each entry is the mean
# intensity of the gradient over a specific feature map channel
pooled_grads <- k_mean(grads, axis = c(1, 2, 3))

# This function allows us to access the values of the quantities we just defined:
# `pooled_grads` and the output feature map of `conv2d_3`,
# given a sample image
iterate <- k_function(list(model$input),
                      list(pooled_grads, last_conv_layer$output[1,,,]))

# These are the values of these two quantities, as arrays,
# given our sample image of two elephants
c(pooled_grads_value, conv_layer_output_value) %<-% iterate(list(img_tensor))

# We multiply each channel in the feature map array by
# "how important this channel is" with regard to the dog class
for (i in 1:128) {
  conv_layer_output_value[,,i] <- 
    conv_layer_output_value[,,i] * pooled_grads_value[[i]] 
}

# The channel-wise mean of the resulting feature map
# is our heatmap of class activation
heatmap <- apply(conv_layer_output_value, c(1,2), mean)

We can now create the heatmap and save it to disk for later. Note that the heatmap width and height need to be the same as our preprocessed image, which we cropped to 150x150.

heatmap <- pmax(heatmap, 0) 
heatmap <- heatmap / max(heatmap)
write_heatmap <- function(heatmap, filename, width = 150, height = 150,
                          bg = "white", col = terrain.colors(12)) {
  png(filename, width = width, height = height, bg = bg)
  op = par(mar = c(0,0,0,0))
  on.exit({par(op); dev.off()}, add = TRUE)
  rotate <- function(x) t(apply(x, 2, rev))
  image(rotate(heatmap), axes = FALSE, asp = 1, col = col)
}
write_heatmap(heatmap, "dog_heatmap.png") 

knitr::include_graphics("dog_heatmap.png")

Now we can use the magick package to generate an image that superimposes the original image with the heatmap we just obtained. We can see that the heatmap is focused on the dogs face and, more specifically, on the left eye and upper part of the snout of the dog.

library(magick) 
library(viridis) 
# Read the original elephant image and it's geometry
image <- image_read(img_path)
info <- image_info(image) 
geometry <- sprintf("%dx%d!", info$width, info$height) 
# Create a blended / transparent version of the heatmap image
pal <- col2rgb(viridis(20), alpha = TRUE) 
alpha <- floor(seq(0, 255, length = ncol(pal))) 
pal_col <- rgb(t(pal), alpha = alpha, maxColorValue = 255)
write_heatmap(heatmap, "dog_overlay.png", 
              width = 14, height = 14, bg = NA, col = pal_col) 
# Overlay the heatmap
image_read("dog_overlay.png") %>% 
  image_resize(geometry, filter = "quadratic") %>% 
  image_composite(image, operator = "blend", compose_args = "20") %>%
  plot() 

Superpixels with LIME

The LIME package offers a similar capability referred to as superpixels. Superpixels is the process of segmenting an image. We can use this concept to help identify parts of an image that explain our model’s prediction.

plot_superpixels(img_path)

First, we’ll create a function that preprocesses our image.

image_prep <- function(x) {
  arrays <- lapply(x, function(path) {
  img <- image_load(path, target_size = c(150, 150))
  x <- image_to_array(img)
  x <- array_reshape(x, c(1, dim(x)))
  x <- x / 255
  })
  do.call(abind::abind, c(arrays, list(along = 1)))
}

Next, we’ll create a lime explainer object that takes our image path, our cats_and_dogs_small_2.h5 model object, and the image preprocessing function. Also, since keras does not supply class labels in the model object we pass a vector that relates the label to the 0-1 predictions.

model_labels <- c("0" = "cat", "1" = "dog")
explainer <- lime(img_path, as_classifier(model, model_labels), image_prep)

We can now run the explain() function which implements the methodology in https://arxiv.org/abs/1602.04938. See https://lime.data-imaginist.com/ for more details. In this example, we are searching for the 35 superpixels that help to explain our model’s prediction. The current settings result in an explanation fit (\(R^2\)) of 0.7. You could tweak the parameter settings for explain() to try increase the explanation fit. We see that the face of the dog is identified as one of the superpixels that helps to explain the correct prediction.

explanation <- explain(
  img_path, 
  explainer,
  n_labels = 2,
  n_features = 35,
  n_superpixels = 35,
  weight = 10,
  background = "white"
  )

plot_image_explanation(explanation)

An alternative view that puts more focus on the relevant superpixels, but removes the context can be seen by using display = ‘block’:

plot_image_explanation(explanation, display = 'block')

---

1. Note that much of the content in this module comes from Chapter 5 of Deep Learning with R and also from this notebook by the author.↩