Sobel edge detection on raster data in R returns lines instead of edges in r

Question

I'm using the lidR, raster, and terra packages in R to perform a canopy height model (CHM) analysis on a LiDAR dataset, and I'm attempting to use Sobel edge detection on my CHM data.

library("lidR")
library("raster")
library(rgl)
library(RANN)
library(terra)
library(sf)
library(EBImage)


las <- readLAS("lasfile")
st_crs(las) <- 25832
lasNormal <- filter_poi(las, Classification != 2)



chm <-
  rasterize_canopy(las, res = 0.5, p2r(0.3, na.fill = NULL))
gf <- focalWeight(chm, .3, "Gauss")
chm <- focal(chm, w = gf)

# Plot the smoothed CHM in grayscale
plot(chm, col = grey.colors(255))


# Set the CRS of the raster
terra::crs(chm) <- "+init=epsg:25832"
# Convert the raster to a matrix
chm_matrix <- raster::as.matrix(chm)
chm_im <- as.cimg(chm_matrix)

conv2D <- function(img, fltr) {
  # Dimensions of image and filter
  img_dim <- dim(img)
  fltr_dim <- dim(fltr)
  
  # Padding size
  pad <- floor(fltr_dim / 2)
  
  # Add zero padding to the image
  padded_img <- matrix(0, nrow = img_dim[1] + 2 * pad[1], ncol = img_dim[2] + 2 * pad[2])
  padded_img[(pad[1] + 1):(pad[1] + img_dim[1]), (pad[2] + 1):(pad[2] + img_dim[2])] <- img
  
  # Output image
  output <- matrix(0, nrow = img_dim[1], ncol = img_dim[2])
  
  # Apply the filter
  for (i in 1:img_dim[1]) {
    for (j in 1:img_dim[2]) {
      region <- padded_img[i:(i + 2 * pad[1]), j:(j + 2 * pad[2])]
      output[i, j] <- sum(region * fltr)
    }
  }
  
  return(output)
}

# Now, use this function to compute the edge detection
sh <- matrix(c(-1, 0, 1, -2, 0, 2, -1, 0, 1), nrow = 3)
sv <- t(sh)
xh <- conv2D(chm_matrix, sh)
xv <- conv2D(chm_matrix, sv)
edge <- sqrt(xh^2 + xv^2)



edge_raster <- terra::rast(edge)

# Copy geographic information from the original raster
terra::ext(edge_raster) <- terra::ext(chm)
terra::crs(edge_raster) <- terra::crs(chm)

# Plot the raster
terra::plot(edge_raster,col = grey.colors(255))

However, the result I'm getting is a raster image with lines, whereas I was expecting to see edges from the edge detection process. Can someone please help me understand why I'm getting this result and how I might be able to achieve the edge detection I'm looking for?

Answer 1

I'm not familiar with the packages you are using other than EBImage. Perhaps the Sobel object you generated is not an Image object? I've extracted your example image (less axes and labels) and applied a Sobel filter function meant to return an image. Perhaps this will do what you want.

# Assuming EBImage has been loaded

# Sobel function
  sobel <- function (x) 
  {
    sh <- matrix(c(1, 2, 1, 0, 0, 0, -1, -2, -1), nrow = 3)
    sv <- t(sh)
    xh <- filter2(x, sh, boundary = "replicate")
    xv <- filter2(x, sv, boundary = "replicate")
    return(sqrt(xh^2 + xv^2))
  }

# image
  x <- readImage(file.choose()) # extracted image saved locally
  img <- abind(x, sobel(x), along = 3)
  dev.new(width = 6, height = 1.9)
  plot(img, all = TRUE)