I have a function that takes a raster(image). My main goal is “For each pixel find out how far is it from the edge(low pixel value) in all direction. We count the number of step it takes to reach the edge, add it to values calculated in all the direction and take an average. For that pixel the value now becomes this average value”. I am able to achieve this with the code below”

# Function to calculate slope breaks
calculateSlopeBreaks <- function(rasterObj) {
  rows <- nrow(rasterObj)
  cols <- ncol(rasterObj)
  
  # Create a raster object to store the window sizes
  windowRaster <- raster(nrows = rows, ncols = cols)
  
  # Define the eight cardinal directions
  directions <- list(c(0,1), c(1,1), c(1,0), c(1,-1), c(0,-1), c(-1,-1), c(-1,0), c(-1,1))
  
  # Loop through each pixel in the raster object
  for (r in 1:rows) {
    for (c in 1:cols) {
      
      # Get the value of the current central pixel
      centralPixelValue <- rasterObj[r,c]
      slopeBreaks <- numeric(length(directions))
      
      # Loop through the cardinal directions
      for (d in 1:length(directions)) {
        direction <- directions[[d]]
        r_offset <- direction[1]
        c_offset <- direction[2]
        
        # Move in the current direction until a lower value is found
        for (step in 1:max(rows, cols)) {
          r_step <- r + step * r_offset
          c_step <- c + step * c_offset
          
          # Check bounds of the raster
          if (r_step < 1 || r_step > rows || c_step < 1 || c_step > cols) break
          
          # Get the value of the neighboring pixel
          pixelValue <- rasterObj[r_step, c_step]
          
          # Check if a break in slope is found
          if (pixelValue < centralPixelValue) {
            slopeBreaks[d] <- step
            break
          }
        }
      }
      
      # Calculate the mean value of breaks or assign zero if none found
      windowSize <- if (all(slopeBreaks == 0)) 0 else mean(slopeBreaks[slopeBreaks > 0])
      windowRaster[r, c] <- windowSize
    }
  }
  
  return(windowRaster)
}


# Apply the slope breaks calculation
windowSizes <- calculateSlopeBreaks(chm) 

This is slow for big images. As i have a multiple cores i want to run this process in parallel so that it is faster.

My intial idea is that i will run the calculation for each direction in parallel. And this is what i tired

process_direction <- function(rasterMatrix, direction, rows, cols) {
  results <- matrix(0, nrow = rows, ncol = cols)
  r_offset <- direction[1]
  c_offset <- direction[2]
  
  for (r in 1:rows) {
    for (c in 1:cols) {
      
      minPixelValue <- rasterMatrix[r, c]
      minStep <- 0
      allGreater <- TRUE
      isTreeEdgePixel <- FALSE
      
      r_step <- r + r_offset
      c_step <- c + c_offset
      
      if (r_step >= 1 && r_step <= rows && c_step >= 1 && c_step <= cols) {
        pixelValue <- rasterMatrix[r_step, c_step]
        
        if (pixelValue == 0) {
          isTreeEdgePixel <- TRUE
          break
        } else if (pixelValue < minPixelValue) {
          allGreater <- FALSE
          minPixelValue <- pixelValue
          minStep <- 1
        }
      }
      print(minStep)
      results[r, c] <- minStep
    }
  }
  
  return(list("results" = results, "allGreater" = allGreater, "isTreeEdgePixel" = isTreeEdgePixel))
}

rasterToMatrix <- function(rasterObj) {
  return(as.matrix(rasterObj))
}

calculateSlopeBreaksParallel <- function(rasterObj) {
  rasterMatrix <- rasterToMatrix(rasterObj)
  rows <- nrow(rasterObj)
  cols <- ncol(rasterObj)
  directions <- list(c(0,1), c(1,1), c(1,0), c(1,-1), c(0,-1), c(-1,-1), c(-1,0), c(-1,1))
  no_cores <- 16
  
  cl <- makeCluster(no_cores)
  clusterExport(cl, list("rasterMatrix", "rows", "cols"), envir = environment())
  
  all_results <- parLapply(cl, directions, process_direction, rasterMatrix=rasterMatrix, rows=rows, cols=cols)
  print(all_results)
  stopCluster(cl)
  
  windowRaster <- raster(nrows = rows, ncols = cols)
  
  for (r in 1:rows) {
    for (c in 1:cols) {
      slopeBreaks <- sapply(all_results, function(x) x$results[r, c])
      print(slopeBreaks)
      isTreeEdgePixel <- any(sapply(all_results, function(x) x$isTreeEdgePixel))
      allGreater <- all(sapply(all_results, function(x) x$allGreater))
      
      if (isTreeEdgePixel) {
        windowSize <- 0
      } else if (allGreater) {
        windowSize <- 0
      } else {
        windowSize <- mean(slopeBreaks) * 0.5
      }
      
      windowRaster[r, c] <- windowSize
    }
  }
  
  return(windowRaster)
}

windowSizes <- calculateSlopeBreaksParallel(chm)

This is just giving me a blank image. I tried some debugging and figure out that the return from the process_direction function is always zero.