R中的盒子计数
Box-counting in R
我正在尝试使用 R 中的盒计数来计算由二进制矩阵(完全由 1 和 0 填充的 512 x 512 矩阵)表示的对象的分形维数。自然而然地找到了this题,但是题目和我的不一样,答案也很乏味。我自己并不熟悉盒子计数的数学细节,但熟悉它在分形维数计算中的实用性。话虽这么说,但 R 中一些声称计算盒计数维度的函数背后的许多文档的解释方式让我有点难以理解(根据我在网上阅读的内容,我不是唯一的)。我很好奇这个网站上是否有人有过使用这些程序的经验,他们是否可以分享他们的智慧?
编辑
我要分析的数据是二进制矩阵的形式,唯一的非零项是 1。这样一个矩阵的例子如下:
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)
我问的只是方阵。虽然这个矩阵很小,但我处理的矩阵要大得多(如前所述,512x512)。理想的函数会将矩阵作为输入,如上所示,并输出所述矩阵中包含的图形的盒计数维度。
此代码应该有效:
# 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])
请注意,此代码仅适用于 512X512 的矩阵,对于不同的维度,您将不得不更改它。此外,您可能希望在没有两个边界的情况下进行线性拟合,因为这会给您带来尺寸误差。
我正在尝试使用 R 中的盒计数来计算由二进制矩阵(完全由 1 和 0 填充的 512 x 512 矩阵)表示的对象的分形维数。自然而然地找到了this题,但是题目和我的不一样,答案也很乏味。我自己并不熟悉盒子计数的数学细节,但熟悉它在分形维数计算中的实用性。话虽这么说,但 R 中一些声称计算盒计数维度的函数背后的许多文档的解释方式让我有点难以理解(根据我在网上阅读的内容,我不是唯一的)。我很好奇这个网站上是否有人有过使用这些程序的经验,他们是否可以分享他们的智慧?
编辑
我要分析的数据是二进制矩阵的形式,唯一的非零项是 1。这样一个矩阵的例子如下:
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)
我问的只是方阵。虽然这个矩阵很小,但我处理的矩阵要大得多(如前所述,512x512)。理想的函数会将矩阵作为输入,如上所示,并输出所述矩阵中包含的图形的盒计数维度。
此代码应该有效:
# 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])
请注意,此代码仅适用于 512X512 的矩阵,对于不同的维度,您将不得不更改它。此外,您可能希望在没有两个边界的情况下进行线性拟合,因为这会给您带来尺寸误差。