Implementing the MATLAB filter2 function in R

Question

I am currently implementing the filter2 MATLAB function in R, which is a method for 2D convolution. I have made for the 2D convolution work, but how the 'valid' option in filter2 works is not quite clear to me.

The MATLAB function is described here: http://se.mathworks.com/help/matlab/ref/filter2.html

My implementation:

filter2D <- function(img, window) {
  # Algoritm for 2D Convolution
  filter_center_index_y <- median(1:dim(window)[1])
  filter_max_index_y <- dim(window)[1]
  filter_center_index_x <- median(1:dim(window)[2])
  filter_max_index_x <- dim(window)[2]

  # For each position in the picture, 2D convolution is done by 
  # calculating a score for all overlapping values within the two matrices
  x_min <- 1
  x_max <- dim(img)[2]
  y_min <- 1
  y_max <- dim(img)[1]

  df <- NULL
  for (x_val in c(x_min:x_max)){
    for (y_val in c(y_min:y_max)){
      # Distanced from cell
      img_dist_left <- x_val-1
      img_dist_right <- x_max-x_val
      img_dist_up <- y_val-1
      img_dist_down <- y_max-y_val

      # Overlapping filter cells
      filter_x_start <- filter_center_index_x-img_dist_left
      if (filter_x_start < 1) {
        filter_x_start <- 1
      }
      filter_x_end <- filter_center_index_x+img_dist_right
      if (filter_x_end > filter_max_index_x) {
        filter_x_end <- filter_max_index_x
      }
      filter_y_start <- filter_center_index_y-img_dist_up
      if (filter_y_start < 1) {
        filter_y_start <- 1
      }
      filter_y_end <- filter_center_index_y+img_dist_down
      if (filter_y_end > filter_max_index_y) {
        filter_y_end <- filter_max_index_y
      }

      # Part of filter that overlaps
      filter_overlap_matrix <- filter[filter_y_start:filter_y_end, filter_x_start:filter_x_end]

      # Overlapped image cells
      image_x_start <- x_val-filter_center_index_x+1
      if (image_x_start < 1) {
        image_x_start <- 1
      }
      image_x_end <- x_val+filter_max_index_x-filter_center_index_x
      if (image_x_end > x_max) {
        image_x_end <- x_max
      }
      image_y_start <- y_val-filter_center_index_y+1
      if (image_y_start < 1) {
        image_y_start <- 1
      }
      image_y_end <- y_val+filter_max_index_y-filter_center_index_y
      if (image_y_end > y_max) {
        image_y_end <- y_max
      }

      # Part of image that is overlapped
      image_overlap_matrix <- img[image_y_start:image_y_end, image_x_start:image_x_end]

      # Calculating the cell value
      cell_value <- sum(filter_overlap_matrix*image_overlap_matrix)
      df = rbind(df,data.frame(x_val,y_val, cell_value))
    }
  }

  # Axis labels
  x_axis <- c(x_min:x_max)
  y_axis <- c(y_min:y_max)

  # Populating matrix
  filter_matrix <- matrix(df[,3], nrow = x_max, ncol = y_max, dimnames = list(x_axis, y_axis))

  return(filter_matrix)
}

Running the method:

> image
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    7    8    9   10   11   12
[3,]   13   14   15   16   17   18
[4,]   19   20   21   22   23   24
[5,]   25   26   27   28   29   30
[6,]   31   32   33   34   35   36

> filter
     [,1] [,2] [,3]
[1,]    1    2    1
[2,]    0    0    0
[3,]   -1   -2   -1

> filter2D(image, filter)
    1   2   3   4   5   6
1 -22 -32 -36 -40 -44 -35
2 -36 -48 -48 -48 -48 -36
3 -36 -48 -48 -48 -48 -36
4 -36 -48 -48 -48 -48 -36
5 -36 -48 -48 -48 -48 -36
6  76 104 108 112 116  89 

This is the same output that filter2(image, filter) produces in matlab, however, when the option 'valid' is added the following output is generated:

-48  -48  -48  -48                                                                                                                                                                                                 
-48  -48  -48  -48                                                                                                                                                                                                 
-48  -48  -48  -48                                                                                                                                                                                                 
-48  -48  -48  -48  

It is not entirely obvious how filter2 with the 'valid' option generates this. Is it just using the center values? Or is it doing something more sophisticated?

Answer 1

Before I start, your code is actually doing 2D correlation. 2D convolution requires that you perform a 180 degree rotation on the kernel before doing the weighted sum. Correlation and convolution are in fact the same operation if the kernel is symmetric (i.e. the transpose of the kernel is equal to itself). I just wanted to make that clear before I start. Also, the documentation for filter2 does state that correlation is being performed.

---

The 'valid' option in MATLAB simply means that it should return only the outputs where the kernel fully overlaps the 2D signal when performing filtering. Because you have a 3 x 3 kernel, this means that at location (2,2) in the 2D signal for example, the kernel does not go outside of the signal boundaries. Therefore, what is returned is the extent of the filtered 2D signal where the kernel was fully inside the original 2D signal. To give you an example, if you placed the kernel at location (1,1), some of the kernel would go out of bounds. Handling out of bounds conditions when filtering can be done in many ways which may affect results and the interpretation of those results when it comes down to it. Therefore the 'valid' option is desired as you are using true information that forms the final result. You aren't interpolating or performing any estimations for data that goes beyond the borders of the 2D signal.

Simply put, you return a reduced matrix which removes the border elements. The filter being odd shaped makes this easy. You simply remove the first and last floor(M/2) rows and the first and last floor(N/2) columns where M x N is the size of your kernel. Therefore, because your kernel is 3 x 3, this means that we need to remove 1 row from the top and 1 row from the bottom, as well as 1 column from the left and 1 column from the right. This results in the -48 within a 4 x 4 grid as you see from the output of MATLAB.

Therefore, if you want to use the 'valid' option in your code, simply make sure that you remove the border elements in your result. You can do this right at the end before your return the matrix:

# Place your code here...
# ...
# ...
# Now we're at the end of your code

# Populating matrix
filter_matrix <- matrix(df[,3], nrow = x_max, ncol = y_max, dimnames = list(x_axis, y_axis))

# New - Determine rows and columns of matrix as well as the filter kernel
nrow_window <- nrow(window)
ncol_window <- ncol(window)
nrows <- nrow(filter_matrix)
ncols <- ncol(filter_matrix)

# New - Figure out where to cut off
row_cutoff <- floor(nrow_window/2)
col_cutoff <- floor(ncol_window/2)

# New - Remove out borders
filter_matrix <- filter_matrix[((1+row_cutoff):(nrows-row_cutoff)), ((1+col_cutoff):(ncols-col_cutoff))]

# Finally return matrix
return(filter_matrix)    

Example Run

Using your data:

> image <- t(matrix(c(1:36), nrow=6, ncol=6))
> image
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    7    8    9   10   11   12
[3,]   13   14   15   16   17   18
[4,]   19   20   21   22   23   24
[5,]   25   26   27   28   29   30
[6,]   31   32   33   34   35   36
> filter <- matrix(c(1,0,-1,2,0,-2,1,0,-1), nrow=3, ncol=3)
> filter
     [,1] [,2] [,3]
[1,]    1    2    1
[2,]    0    0    0
[3,]   -1   -2   -1

I ran the function and now I get:

> filter2D(image,filter)
    2   3   4   5
2 -48 -48 -48 -48
3 -48 -48 -48 -48
4 -48 -48 -48 -48
5 -48 -48 -48 -48

I think it may be important to leave the horizontal and vertical labels to be the way they are so you can explicitly see that not all of the signal is being returned, which is what the code is currently doing.... that's up to you though. I'll leave that to you to decide.