# Box counting procedure to determine fractal dimension
# of an image. Currently restricted to square images
# with side equal to a power of 2
#
# Rajarshi Guha <rajarshi@prsidency.com>
# September 2004

require(rimage)

i2m <- function(img) {
    m <- matrix(0, nrow=attr(img,'dim')[1], ncol=attr(img,'dim')[2])
    idx <- which(img > 0.5)
    m[idx] <- 1
    m
}


bc <- function(img, thvalue=0.5, verbose=TRUE) {

    if (attr(img,'dim')[1] != attr(img,'dim')[2]) {
        stop("Currently it only works for square images")
    }
    if (round(log2(attr(img,'dim')[1]), digits=1) != log2(attr(img,'dim')[1])) {
        stop("Currently dimensions of the image must be a power of 2")
    }
        
    bv <- c()
    powers <- 1:9
    boxsize <- 2^powers
    imgsize <- attr(img,'dim')[1]
   
    # we want the contour
    sb <- normalize(sobel(img))
    tsb <- thresholding(sb, th=thvalue, mode='fixed')
    tsb <- i2m(img)


    print('Box Size   Count   Total Box')
    for (b in boxsize) {
        bcount <- 0
        bstart <- seq(1,imgsize, by=b)
        totalbox <- 0
        for (x in bstart) {
            for (y in bstart) {
                totalbox <- totalbox+1
                xend <- x + b - 1
                yend <- y + b - 1

                # since we applied sobel filter the image
                # is inverted. So original 1's become 0's and so on
                #if (length(which( tsb[ x:xend, y:yend ] < thvalue)) > 0) {
                #    bcount <- bcount + 1
                #}
                if (any( tsb[ x:xend, y:yend ] < thvalue)) {
                    bcount <- bcount + 1
                }
            }
        }
        if (verbose == TRUE) {
            print(paste(b,bcount,totalbox))
        }
        bv <- c(bv, bcount)
    }
    plot(log(1/boxsize), log(bv))
    df <- data.frame(x=log(1/boxsize), y=log(bv))
    model <- lm(y ~  x, df)
    abline(model)
    print(paste('fractal dimension =',model$coeff[2]))
    list(fdim=model$coeff[2],model=model)
}