Search code examples
rfractals

Box-counting in R


I am trying to use box-counting in R to calculate the fractal dimension of an object represented by a binary matrix (a 512 x 512 matrix entirely populated by 1s and 0s). Naturally, I found this question, but the question is not similar to mine and the answer is lackluster. I'm not myself familiar with the mathematical details of box-counting, but am familiar with its utility in the calculation of fractal dimension. That being said, much of the documentation behind some functions in R claiming to calculate the box-counting dimension are explained in ways that are a bit over my head (from what I've read online, I'm not the only one). I'm curious as to whether or not anyone on this site has had experience with these programs, and if they could share their wisdom?

Edit

The data I am trying to analyze is in the form of a binary matrix, with the only nonzero entries being 1's. An example of such a matrix is as follows:

A = matrix(c(0,0,0,0,0,0,1,0,0,1,1,0,0,1,1,1),byrow=T,nrow=4,ncol=4)

I'm only inquiring about square matrices. Although this matrix is small, the matrices I am dealing with are much larger (512x512 as stated earlier). The ideal function would take as an input a matrix, as one above, and output the box-counting dimension of the figure contained in said matrix.


Solution

  • This code should work:

    # create a dataset to calculate de Hausdorff-Besicovitch dimension
    mat <- matrix(runif(512*512),nrow = 512,ncol = 512)
    
    mat[mat<=0.5] <- 0
    mat[mat>0.5] <- 1
    
    cant <- sum(mat)
    
    fragment <- rep(2,10)**(0:9)
    Table <- data.frame(Delta = rep(512,10)/(fragment ), N = fragment**2)
    Table$LogDelta <- log(Table$Delta)
    
    for(i in 2:10){
      delta_aux <- Table$Delta[i]
    
      for(j in 1:fragment [i]){
        row_id <- ((j-1)*delta_aux+1):(j*delta_aux)
        for(k in 1:fragment [i]){
          col_id <- ((k-1)*delta_aux+1):(k*delta_aux)
          if(sum(mat[row_id,col_id]) == 0){
            Table$N[i] <- Table$N[i] - 1
          }
        }
      }
    }
    
    Table$LogN <- log(Table$N)
    lm_dim <- lm(Table$LogN ~ Table$LogDelta)
    
    plot(Table$LogN ~ Table$LogDelta)
    abline(lm_dim)
    
    print('The box-counting dimension is:')
    print(-lm_dim$coefficients[2])
    
    # without the borders
    Table <- Table[2:nrow(Table),]
    lm_dim <- lm(Table$LogN ~ Table$LogDelta)
    
    plot(Table$LogN ~ Table$LogDelta)
    abline(lm_dim)
    
    print('The box-counting dimension is:')
    print(-lm_dim$coefficients[2])
    

    Note that, this code is done only for matrixes of 512X512, for different dimensions you will have to change it. Also, you may want to do the linear fitting without the two borders because that will give you an error in the dimensions.