dist-wmat-computation.R

library(raster)

# global config
INPUT_FNAME <- "/home/vlx/Work/biggis/data-julian/Auto150_georef.tif"
BAND <- 1
AGGREG_FACT <- 150
AGGREG_FUNC <- mean


# here comes the code
r <- raster(INPUT_FNAME, band = BAND)
ncol(r); nrow(r)

# aggregating pixels (downsampling)
r_aggreg <- aggregate(r, fact = AGGREG_FACT, fun = AGGREG_FUNC,
                      expand = TRUE, na.rm = TRUE)
width <- ncol(r_aggreg)
height <- nrow(r_aggreg)

print(width)
print(height)
print(width * height)
plot(r_aggreg)


# Simple convolution with circular filter.
# The circle represents neightbourhood to a particular
# distance expressed in the original CRS units.
# All the weights in the circle are the same and SUM(weights) = 1
focalw <- focalWeight(r_aggreg, d = 1, type = "circle")
plot(raster(focalw))
focalPadValue <- min(values(r_aggreg))
#focalPadValue <- mean(values(r_aggreg))
r_aggreg_conv <- focal(r_aggreg, focalw,
                       pad = TRUE, padValue = focalPadValue)
plot(r_aggreg_conv)


# Convolution with a gausian filter
# focalw <- focalWeight(r_aggreg, d = 2, type = "Gauss")
# plot(raster(focalw))
# focalPadValue <- min(values(r_aggreg))
# r_aggreg_conv <- focal(r_aggreg, focalw,
#                        pad = TRUE, padValue = focalPadValue)
# plot(r_aggreg_conv)



#r_by_col <- c(as.matrix(r_aggreg_conv)) # floating point cells
r_by_col <- as.vector(as.matrix(r_aggreg_conv), mode = "integer") # integer cells

z <- width * height # the distance matrix is a square matrix
dist <- matrix(nrow = z, ncol = z) # might be simmetrical but not necessarily :-)

# this is the computationally intensive part
for (i in 1:z) {
  for (j in 1:z) {

    row1 <- ((i - 1) %% height) + 1  # operator %% in R is modulo
    col1 <- ((i - 1) %/% height) + 1 # operator %/% in R is division without remainder

    row2 <- ((j - 1) %% height) + 1
    col2 <- ((j - 1) %/% height) + 1

    dist[i,j] <-
      abs(r_by_col[i] - r_by_col[j]) * # weight
      (row1 - row2) ^ 2 + (col1 - col2) ^ 2 # distance squared
  }
}

plot(raster(dist))

# writing the matric as CSV file without column/row names, cells delimited by space
write.table(dist, file = "dist-wmat8x8.csv", row.names = FALSE, col.names = FALSE)