Intro

Having scanned a space (a forest) in voxels, coding as 0 non empty cells, The procedure computes the percolation in each of the 3 dimensions, identifying percolating clusters (in that dimension) and their smallest section (orthogonal to that direction). Moreover, total number of cells is computed for each percolatin (in any direction) cluster. The procedure is tested using a random matrix as input.

Set up

library(tidyverse)
## -- Attaching packages ------------- tidyverse 1.2.1 --
## v ggplot2 3.2.1     v purrr   0.3.2
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## -- Conflicts ---------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(mmand)

Functions

# --- subset 'a', taking slice 'i', in dimensions 's' ---
slice <- function(a, s, i){
  a[which(slice.index(a, c(s)) == i, arr.ind = TRUE)]
}

# --- labels of percolating clusters (along dimension 's') ---
slice_intersect <- function(a, s){
  t1 <- table(slice(a, s, 1))
  t2 <- table(slice(a, s, dim(a)[s]))
  intersect(unlist(dimnames(t1)), unlist(dimnames(t2)))
}

# --- smallest section of cluster 'cl' orthogonal to dimension 's' ---
min_cl_sect <- function(c, s, cl){
  min_sect <- NA
  for(i in 1:dim(c)[s]){
    sect <- table(slice(c, s, i))[cl]
    if(is.na(min_sect)) min_sect <- sect
    if(min_sect > sect) min_sect <- sect
  }
  return(min_sect)
}

# --- percolation statistics ---
#     a list with tibbles:
#        min_sect_t: percolating clusters, by dimension, with min_sect
#        cl_ncell_t: number of cells for each perc. cl
percolation_sts <- function(space, k = NA){
  if(is.na(k)){
    # kernel - von Newmann 3d connections
    k <- array(c(
      0, 0, 0, 
      0, 1, 0, 
      0, 0, 0,
      
      0, 1, 0, 
      1, 1, 1, 
      0, 1, 0,
      
      0, 0, 0, 
      0, 1, 0, 
      0, 0, 0),
      dim = c(3, 3, 3))
  }
  
  clusters <- components(x,k)
  # labels k-adiacent cells with the cluster id number
  # considering zero-cells as 'non-adiacent'
  
  # for each dimension: 
  #  find pecolating clusters, their volume and min section
  
  min_sect_t <- tibble()
  # dim_num, pc_lbl, pc_min_sect
  for(d in 1:length(dim(clusters))){
    pcs <- slice_intersect(clusters, d)
    if(!is_empty(pcs)) {
      for(pc in pcs){
        ms <- min_cl_sect(clusters, d, pc)
        min_sect_t <- bind_rows(min_sect_t,
          tibble(dim_num = d, pc_lbl = pc, pc_min_sect = ms))
      }
    }
  }
  
  cl_ncell_t <- tibble()
  # pc_lbl, pc_ncells
  if(nrow(min_sect_t) > 0) {
    cl_ncell_t <- tibble( 
           pc_lbl = unique(min_sect_t$pc_lbl),
           pc_ncells = table(clusters)[pc_lbl]
           )
  }
  return(list(min_sect_t = min_sect_t, cl_ncell_t))
}

test application

p <- .7
x <- array(sample(0:1, 1000, replace=T,prob=c(1-p, p)), dim=c(10, 10, 10))
percolation_sts(x)
## $min_sect_t
## # A tibble: 3 x 3
##   dim_num pc_lbl pc_min_sect
##     <int> <chr>        <int>
## 1       1 1               63
## 2       2 1               61
## 3       3 1               60
## 
## [[2]]
## # A tibble: 1 x 2
##   pc_lbl pc_ncells
##   <chr>      <int>
## 1 1            685
x %>% 
  reshape2::melt() %>% 
  filter(Var3 %in% c(1, 10)) %>% 
  ggplot() +
  facet_wrap(vars(Var3)) +
  geom_raster(aes(Var1, Var2, fill = value))