With imager package, how to get the correlation of an image (with a kernel)?

Question

I have this pic:

I want to get the image named the "correlation" of this image with a discrete kernel, the matrix given below. I first try to use the correlate function:

kernel <- rbind(
  c( 0, 1, 1, 2, 2, 2, 1, 1, 0 )
  , c( 1,  2,  4, 5, 5, 5, 4, 2, 1 )
  , c( 1,  4,  5, 3, 0, 3, 5, 4, 1 )
  , c( 2,  5,  3, -12, -24, -12, 3, 5, 2 )
  , c( 2,  5,  0, -24, -40, -24, 0, 5, 2 )
  , c( 2,  5,  3, -12, -24, -12, 3, 5, 2 )
  , c( 1,  4,  5, 3, 0, 3, 5, 4, 1 )
  , c( 1,  2,  4, 5, 5, 5, 4, 2, 1 )
  , c( 0,  1,  1, 2, 2, 2, 1, 1, 0 )
)

library(imager)
im <- load.image("mypic.png")
img <- correlate(im, as.cimg(kernel))
plot(img)

This gives this gray image:

Now I try manually, not taking care of the borders:

RED <- red <- 255*as.matrix(R(im)[, , 1, 1])
GRE <- gre <- 255*as.matrix(G(im)[, , 1, 1])
BLU <- blu <- 255*as.matrix(B(im)[, , 1, 1])

for(r in 10:500) {
  for(c in 10:500){
    red[r, c] <- gre[r, c] <- blu[r, c] <- 0
    for(u in 1:9){
      for(v in 1:9) {
        red[r, c] <- red[r, c] + kernel[u, v] * RED[r+u-5, c+v-5]
        gre[r, c] <- gre[r, c] + kernel[u, v] * GRE[r+u-5, c+v-5]
        blu[r, c] <- blu[r, c] + kernel[u, v] * BLU[r+u-5, c+v-5]
      }
    }
  }
}

img <- as.cimg(abind::abind(red, gre, blu, along = 3L))
plot(img)

This gives a gray image as well.

The gray image is not the correct result. I can get the expected result with the magick package:

library(magick)
im <- image_read("mypic.png")
kern <- "9x9:
    0    1    1    2    2    2    1    1    0
    1    2    4    5    5    5    4    2    1
    1    4    6    3    0    3    6    4    1
    2    5    3  -12  -24  -12    3    5    2
    2    5    0  -24  -40  -24    0    5    2
    2    5    3  -12  -24  -12    3    5    2
    1    4    6    3    0    3    6    4    1
    1    2    4    5    5    5    4    2    1
    0    1    1    2    2    2    1    1    0"
image_morphology(im, kernel = kern, method = "Correlate")

I know it's the expected result because I get almost the same one with Haskell (except that the colors are inversed). So how to do with imager::correlate and manually?

Answer 1

I think the correct information is there, but the contrast is poor. You can see it better if you set normalise = TRUE

library(imager)

load.image("Klein.png") |>
  correlate(as.cimg(kernel), normalise = TRUE) |> 
  plot()

The problem seems to be that the LoG can give both positive and negative values, and we would expect the pixels where no change is occuring to be 0. Since the plot method for these images normalizes values to between 0 and 1, any negative pixels will be dark and the positive pixels are bright, with 0 values being gray.

We get something much closer to the desired result if we set all negative pixels to 0, since this will make the areas with no changes black upon normalization (effectively setting 0 as the new floor)

arr <- correlate(im, as.cimg(kernel), norm = TRUE) |> as.array()
arr[arr < 0] <- 0
as.cimg(arr) |> plot()

You can bring out more detail by raising arr to a fractional power (since the plot method by default auto-normalizes to between 0 and 1)

as.cimg(arr^0.4) |> plot()

And of course if you want an inverted image, you can do

as.cimg(-(sqrt(arr))) |> plot()

---

Edit

The main drawback of implementing this entirely in R is speed. You need to get the correlation of kernel with the 9x9 pixel square surrounding each pixel of each channel. There are almost certainly quicker ways to implement this in R, but the naive approach would be something like:

correlate_image <- function(image, kernel) {
  
  correlate_channel <- function(mat, kernel) {
    
    height_seq <- seq(nrow(kernel)) - 1
    width_seq  <- seq(ncol(kernel)) - 1
    res        <- matrix(0, ncol = ncol(mat), nrow = nrow(mat))
    c_kernel   <- c(kernel)
    
    for(i in seq(nrow(mat) - nrow(kernel) - 1)) {
      for(j in seq(ncol(mat) - ncol(kernel) - 1)) {
        submat <- c(mat[i + height_seq, j + width_seq])
        if(sd(submat) != 0) res[i, j] <- cor(submat, c_kernel)
      }
    }
    res[res < 0] <- 0
    return(res/max(res))
  }
  
  arr <- as.array(image)
  arr[,,1,1]   <- correlate_channel(arr[,,1,1], kernel)
  arr[,,1,2]   <- correlate_channel(arr[,,1,2], kernel)
  arr[,,1,3]   <- correlate_channel(arr[,,1,3], kernel)
  arr[,,1,4]   <- 1 # Set alpha channel to 1 for full opacity
  as.cimg(arr)
}

Which allows:

im <- load.image("Klein.png")

my_img <- correlate_image(im, kernel)

plot(my_img)

But takes about 5 minutes to process the image!

The Haskell version seems to use the absolute value of the convolution of the image with the kernel rather than correlation, since if we change cor(submat, c_kernel) to sum(submat * c_kernel) we get very close to the Haskell image, providing we use absolute values and clamp between 0 and 1:

my_img <- correlate_image(im, kernel)
my_img2 <- abs(my_img)
my_img2[my_img2 > 1] <- 1
plot(my_img2)