# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"plot.Wimage" <- function(x, cut.min, graphics.reset = TRUE, 
    common.range = FALSE, color.table = tim.colors(128), Nlevel = NULL, 
    with.lines = FALSE, omd.width = 0.2, ...) {
    m <- nrow(x)
    n <- ncol(x)
    info <- Wimage.info(m, n, cut.min)
    if (is.null(Nlevel)) {
        Nlevel <- info$Lmax
    }
    old.par <- par(no.readonly = TRUE)
    par(mar = c(0, 1, 1, 0), ...)
    # split the screen vertically into the main panel space and
    # a strip to put the legends.
    omd.bottom <- par()$cxy[2]
    figs <- rbind(c(0, 1 - omd.width, omd.bottom, 1), c(1 - omd.width, 
        1, omd.bottom, 1))
    split.screen(figs)
    # now split first screen into the panel
    ind <- split.screen(c(Nlevel + 1, 3), screen = 1)
    # draw smooth basis coefficients.
    zr.common <- range(x, na.rm = TRUE)
    S1 <- info$S[1]:info$S[2]
    S2 <- info$S[3]:info$S[4]
    M <- 1:info$L[1, 1]
    N <- 1:info$L[1, 2]
    add.boxes <- function() {
        if (with.lines) {
            xline(c(0.5, M + 0.5), col = "white", lwd = 0.5)
            yline(c(0.5, N + 0.5), col = "white", lwd = 0.5)
        }
    }
    # move to first plot
    # smooth basis functions
    screen(ind[1])
    image(M, N, x[S1, S2], zlim = zr.common, xaxt = "n", yaxt = "n", 
        xlab = "", ylab = "", col = color.table)
    add.boxes()
    # full image
    screen(ind[3])
    image(1:m, 1:n, x, xaxt = "n", yaxt = "n", col = color.table, 
        xlab = "", ylab = "")
    # save the ranges for each level to draw the legend strips later
    save.zr <- matrix(NA, Nlevel, 2)
    screen.off <- ind[3]
    KK <- 1
    for (lev in (1:Nlevel)) {
        M <- 1:info$L[lev, 1]
        N <- 1:info$L[lev, 2]
        H1 <- info$H[lev, 1]:info$H[lev, 2]
        H2 <- info$H[lev, 3]:info$H[lev, 4]
        V1 <- info$V[lev, 1]:info$V[lev, 2]
        V2 <- info$V[lev, 3]:info$V[lev, 4]
        D1 <- info$Di[lev, 1]:info$Di[lev, 2]
        D2 <- info$Di[lev, 3]:info$Di[lev, 4]
        if (common.range) {
            zr <- zr.common
        }
        else {
            zr <- range(c(x[H1, H2], x[V1, V2], x[D1, D2]), na.rm = TRUE)
            save.zr[lev, ] <- zr
        }
        screen(screen.off + KK)
        image(M, N, x[H1, H2], zlim = zr, xaxt = "n", yaxt = "n", 
            col = color.table, xlab = "", ylab = "")
        add.boxes()
        KK <- KK + 1
        screen(screen.off + KK)
        image(M, N, x[V1, V2], zlim = zr, xaxt = "n", yaxt = "n", 
            col = color.table, xlab = "", ylab = "")
        add.boxes()
        KK <- KK + 1
        screen(screen.off + KK)
        image(M, N, x[D1, D2], zlim = zr, xaxt = "n", yaxt = "n", 
            col = color.table, xlab = "", ylab = "")
        add.boxes()
        KK <- KK + 1
    }
    screen(2)
    if (common.range) {
        image.plot(zlim = zr, legend.only = TRUE, smallplot = c(0.2, 
            0.3, 0.1, 0.5), col = color.table, graphics.reset = TRUE)
    }
    else {
        ind2 <- split.screen(c(Nlevel + 1, 1), screen = 2)
        screen(ind2[1])
        image.plot(zlim = zr, legend.only = TRUE, smallplot = c(0.2, 
            0.3, 0.1, 0.9), col = color.table)
        for (lev in 1:Nlevel) {
            screen(ind2[lev + 1])
            image.plot(zlim = save.zr[lev, ], legend.only = TRUE, 
                col = color.table, smallplot = c(0.2, 0.3, 0.1, 
                  0.9), graphics.reset = TRUE)
        }
    }
    # reset margins
    if (graphics.reset) {
        close.screen(all = TRUE)
        par(cex.axis = 1)
        par(old.par)
        invisible()
    }
    else {
        return(matrix(ind, ncol = 3, byrow = TRUE))
    }
}
"W.i2s" <- function(ind, m, cut.min) {
    info <- W.info(m, cut.min = cut.min)
    # figure out to which level the index addresses.
    M <- rep(0, length(ind))
    for (kk in 1:info$Lmax) {
        M <- ifelse(ind > info$offset[kk], kk, M)
    }
    # compute off for levels taking into account
    # father wavelets and indices that are too large
    off <- rep(0, length(ind))
    M.good <- (M > 0) & (M <= info$Lmax)
    off[M.good] <- info$offset[M[M.good]]
    off[ind > m] <- NA
    list(level = M, i = ind - off, m = m, cut.min = cut.min)
}
"Wimage.i2s" <- function(ind, m, n, cut.min) {
    if (is.matrix(ind)) {
        tempI <- ind[, 1]
        tempJ <- ind[, 2]
    }
    else {
        tempJ <- ceiling(ind/m)
        tempI <- ind - (ceiling(ind/m) - 1) * m
    }
    #cat( tempI, tempJ, fill=T)
    info <- Wimage.info(m = m, n = n, cut.min = cut.min)
    tempM <- W.i2s(tempI, m = m, cut.min = info$S[2])
    tempN <- W.i2s(tempJ, m = n, cut.min = info$S[4])
    # find H V or Di.
    flavor <- rep(1, length(ind))
    flavor <- ifelse(tempM$level > tempN$level, 1, 2)
    flavor[tempM$level == tempN$level] <- 3
    level <- ifelse(tempM$level > tempN$level, tempM$level, tempN$level)
    Sind <- level == 0
    # this change makes the array indexing work below
    flavor[Sind] <- NA
    level[Sind] <- NA
    I <- tempI - info$offset.m[cbind(level, flavor)]
    J <- tempJ - info$offset.n[cbind(level, flavor)]
    #
    #set bad values to NAs
    I <- ifelse((I < 1) | (I > info$L[level, 1]), NA, I)
    #cat( level, J,info$L[level,2], fill=TRUE)
    J <- ifelse((J < 1) | (J > info$L[level, 2]), NA, J)
    # reset level and flavor for fathers
    flavor[Sind] <- 0
    level[Sind] <- 0
    I[Sind] <- tempI[Sind]
    J[Sind] <- tempJ[Sind]
    list(i = I, j = J, level = level, flavor = flavor, cut.min = cut.min, 
        m = m, n = n)
}
"Wimage.info" <- function(m = 128, n = m, cut.min = 4) {
    NN <- n
    MM <- m
    level <- 1
    while (min(c(NN, MM)) > cut.min) {
        level <- level + 1
        NN <- NN/2
        MM <- MM/2
    }
    nlevel <- level - 1
    #   print( nlevel)
    V <- Di <- H <- matrix(NA, nlevel, ncol = 4)
    N <- matrix(NA, nlevel, ncol = 2)
    S <- matrix(NA, 1, ncol = 4)
    offset.m <- matrix(NA, nlevel, ncol = 3)
    offset.n <- matrix(NA, nlevel, ncol = 3)
    n2 <- NN
    m2 <- MM
    n3 <- NN * 2
    m3 <- MM * 2
    n1 <- 1
    m1 <- 1
    level <- 1
    S <- c(1, m2, 1, n2)
    while (n3 <= n & m3 <= m) {
        #     cat( level, fill=TRUE)
        #     cat( m1,m2,m3, fill=TRUE)
        #     cat( n1,n2,n3, fill=TRUE)
        # H and V are defined so that an ordinary image plot
        # gives the expected orientationd for basis functions
        #
        H[level, ] <- c((m2 + 1), m3, n1, n2)
        V[level, ] <- c(m1, m2, (n2 + 1), n3)
        Di[level, ] <- c((m2 + 1), m3, (n2 + 1), n3)
        offset.m[level, ] <- c(H[level, 1], V[level, 1], Di[level, 
            1]) - 1
        offset.n[level, ] <- c(H[level, 3], V[level, 3], Di[level, 
            3]) - 1
        N[level, ] <- c((m3 - m2), c(n3 - n2))
        level <- level + 1
        n2 <- n3
        n3 <- n2 * 2
        m2 <- m3
        m3 <- m2 * 2
    }
    list(m = m, n = n, cut.min = cut.min, S = S, H = H, V = V, 
        Di = Di, L = N, Lmax = nlevel, offset.m = offset.m, offset.n = offset.n)
}
"Wimage.s2i" <- function(i, j, level, flavor, m, n, 
    cut.min, mat = TRUE) {
    ind <- Wimage.info(m = m, n = n, cut.min = cut.min)
    if (is.list(i)) {
        j <- i$j
        level <- i$level
        flavor <- i$flavor
        i <- i$i
    }
    if (length(i) != length(j)) {
        stop("i and j must be same length")
    }
    if (length(level) == 1) {
        level <- rep(level, length(i))
    }
    if (length(flavor) == 1) {
        flavor <- rep(flavor, length(i))
    }
    #
    # reset structure elements to NA that are out of bounds w/r to level
    Sind <- level == 0
    level[(level > ind$Lmax) | Sind] <- NA
    flavor[Sind] <- NA
    hold <- i + ind$offset.m[cbind(level, flavor)] + (j + ind$offset.n[cbind(level, 
        flavor)] - 1) * ind$m
    # handle the father wavelet cases
    hold[Sind] <- i[Sind] + (j[Sind] - 1) * ind$m
    # reset structure elements to NA that are out of bounds w/r to position
    hold[(i < 1) | (j < 1)] <- NA
    hold[i > ind$L[level, 1]] <- NA
    hold[j > ind$L[level, 2]] <- NA
    hold
    if (mat) {
        return(cbind(hold - (ceiling(hold/m) - 1) * m, ceiling(hold/m)))
    }
    else {
        return(hold)
    }
}
"W.info" <- function(m = 128, cut.min = 4) {
    MM <- m
    level <- 1
    while (MM > cut.min) {
        level <- level + 1
        MM <- MM/2
    }
    nlevel <- level - 1
    #   print( nlevel)
    N <- matrix(NA, nlevel, ncol = 1)
    H <- matrix(NA, nlevel, ncol = 2)
    S <- matrix(NA, nrow = 1, ncol = 2)
    offset <- matrix(NA, nlevel, ncol = 1)
    m2 <- MM
    m3 <- MM * 2
    m1 <- 1
    level <- 1
    S <- c(1, m2)
    off <- m2
    while (m3 <= m) {
        H[level, ] <- c((m2 + 1), m3)
        offset[level, 1] <- off
        N[level, 1] <- c((m3 - m2))
        off <- off + m3 - m2
        level <- level + 1
        m2 <- m3
        m3 <- m2 * 2
    }
    list(m = m, cut.min = cut.min, S = S, H = H, L = N, Lmax = nlevel, 
        offset = offset)
}
"W.s2i" <- function(i, level, m, cut.min) {
    # find the offsets in the array for different levels
    ind <- W.info(m, cut.min = cut.min)
    # some checks
    if (length(level) == 1) {
        level <- rep(level, length(i))
    }
    if (length(level) != length(i)) 
        stop("number of levels different from number of positions")
    # indicator for smoothed (father) basis functions
    indS <- level == 0
    # reset level for fathers to one
    level[indS] <- 1
    # reset structure elements to NA that are out of bounds w/r to level
    level[(level > ind$Lmax)] <- NA
    # add postion in level to the level offset
    hold <- ifelse(indS, i, i + ind$offset[level])
    # catch any postions that exceed the allowed limits
    hold[(i > ind$L[level]) | (i < 1)] <- NA
    hold
}