# fields, Tools for spatial data # Copyright 2004-2007, Institute for Mathematics Applied Geosciences # University Corporation for Atmospheric Research # Licensed under the GPL -- www.gpl.org/licenses/gpl.html "WQS.periodic" <- function (x) { if (!is.matrix(x)) x <- matrix(x, nrow = length(x), ncol = 1) D.smooth <- c(-1, 3, 3, -1) D.rough <- c(-1, 3, -3, 1) D.le <- matrix(c(D.smooth, D.rough), nrow = 2, ncol= 4, byrow = TRUE) D.re <- D.le n <- dim(x)[1] m <- dim(x)[2] tmp <- matrix(NA, nrow = n, ncol = m) tmp[1:2, ] <- D.le %*% x[c(n, 1:3), ] tmp[(n - 1):n, ] <- D.re %*% x[c((n - 2):n,1), ] stuff <- n - 4 indx <- seq(2, stuff, 2) tmp[indx + 1, ] <- D.smooth[1] * x[indx, ] + D.smooth[2] * x[indx + 1, ] + D.smooth[3] * x[indx + 2, ] + D.smooth[4] * x[indx + 3, ] tmp[indx + 2, ] <- D.rough[1] * x[indx, ] + D.rough[2] * x[indx + 1, ] + D.rough[3] * x[indx + 2, ] + D.rough[4] * x[indx + 3, ] tmp[c(seq(1, n, 2), seq(2, n, 2)), ] } "WQS.periodic.basis" <- function (N, cut.n = 8) { x <- diag(1, N) nrow <- nrow(x) ncol <- ncol(x) nn <- nrow temp <- x nn <- cut.n * 2 while (nn <= nrow) { temp[1:nn, ] <- WQSi.periodic(temp[1:nn, ]) nn <- nn * 2 } return(temp) } "WQS.periodic.T" <- function (x) { if (!is.matrix(x)) x <- matrix(x, nrow = length(x), ncol = 1) D.smooth <- c(3, -3, -1, -1) D.rough <- c(-1, 1, 3, 3) D.le <- matrix(D.smooth, nrow = 1, ncol = 4, byrow = TRUE) D.re <- matrix(D.rough, nrow = 1, ncol = 4, byrow = TRUE) n <- dim(x)[1] m <- dim(x)[2] x[c(seq(1, n, 2), seq(2, n, 2)), ] <- x tmp <- matrix(NA, nrow = n, ncol = m) tmp[n, ] <- D.le %*% x[c(n-1,n,1:2), ] tmp[1, ] <- D.re %*% x[c((n - 1),n,1,2), ] stuff <- n - 3 indx <- seq(1, stuff, 2) tmp[indx + 1, ] <- D.smooth[1] * x[indx, ] + D.smooth[2] * x[indx + 1, ] + D.smooth[3] * x[indx + 2, ] + D.smooth[4] * x[indx + 3, ] tmp[indx + 2, ] <- D.rough[1] * x[indx, ] + D.rough[2] * x[indx + 1, ] + D.rough[3] * x[indx + 2, ] + D.rough[4] * x[indx + 3, ] tmp } "WQS2d.cylinder" <- function (x, transpose = FALSE, byX=TRUE) { if( byX){ if (!transpose) t(WQS(t(WQS.periodic(x)))) else { WQS.periodic.T(t(WQS.T(t(x)))) } } else{ if (!transpose) t(WQS.periodic(t(WQS(x)))) else { WQS.T(t(WQS.periodic.T(t(x)))) } } } "WQS2di.cylinder" <- function (x, transpose = FALSE, byX=TRUE) { if( byX){ if (!transpose) t(WQSi(t(WQSi.periodic(x)))) else { WQSi.periodic.T(t(WQSi.T(t(x)))) } } else{ if (!transpose) t(WQSi.periodic(t(WQSi(x)))) else { WQSi.T(t(WQSi.periodic.T(t(x)))) } } } "WQSi.periodic" <- function (x) { if (!is.matrix(x)) x <- matrix(x, nrow = length(x), ncol = 1) D.smooth <- c(3, -3, 1, 1)/16 D.rough <- c(1, -1, 3, 3)/16 D.le <- matrix(D.smooth, nrow = 1, ncol = 4, byrow = TRUE) D.re <- matrix(D.rough, nrow = 1, ncol = 4, byrow = TRUE) n <- dim(x)[1] m <- dim(x)[2] x[c(seq(1, n, 2), seq(2, n, 2)), ] <- x tmp <- matrix(NA, nrow = n, ncol = m) tmp[n, ] <- D.le %*% x[c(n-1,n,1:2), ] tmp[1, ] <- D.re %*% x[c((n - 1),n,1,2), ] stuff <- n - 3 indx <- seq(1, stuff, 2) tmp[indx + 1, ] <- D.smooth[1] * x[indx, ] + D.smooth[2] * x[indx + 1, ] + D.smooth[3] * x[indx + 2, ] + D.smooth[4] * x[indx + 3, ] tmp[indx + 2, ] <- D.rough[1] * x[indx, ] + D.rough[2] * x[indx + 1, ] + D.rough[3] * x[indx + 2, ] + D.rough[4] * x[indx + 3, ] tmp } "WQSi.periodic.T" <- function (x) { if (!is.matrix(x)) x <- matrix(x, nrow = length(x), ncol = 1) D.smooth <- c(1, 3, 3, 1)/16 D.rough <- c(1, 3, -3, -1)/16 D.le <- matrix(c(D.smooth, D.rough), nrow = 2, ncol= 4, byrow = TRUE) D.re <- D.le n <- dim(x)[1] m <- dim(x)[2] tmp <- matrix(NA, nrow = n, ncol = m) tmp[1:2, ] <- D.le %*% x[c(n, 1:3), ] tmp[(n - 1):n, ] <- D.re %*% x[c((n - 2):n,1), ] stuff <- n - 4 indx <- seq(2, stuff, 2) tmp[indx + 1, ] <- D.smooth[1] * x[indx, ] + D.smooth[2] * x[indx + 1, ] + D.smooth[3] * x[indx + 2, ] + D.smooth[4] * x[indx + 3, ] tmp[indx + 2, ] <- D.rough[1] * x[indx, ] + D.rough[2] * x[indx + 1, ] + D.rough[3] * x[indx + 2, ] + D.rough[4] * x[indx + 3, ] tmp[c(seq(1, n, 2), seq(2, n, 2)), ] } "Wtransform.cylinder.image" <- function (x, inv = FALSE, transpose = FALSE, byX=TRUE, cut.min = 8) { # if (transpose){ inv <- !inv} n <- dim(x)[1] m <- dim(x)[2] if (n > m) { flip <- TRUE temp <- t(x) n <- dim(temp)[1] m <- dim(temp)[2] # # if the calculations are done on the transpose of the image then # switch the choice of which dimension has the periodic basis # byX <- !byX } else { flip <- FALSE temp <- x } if (n > m) stop(" number of columns of x must >= to number of\nrows") nn <- n mm <- m # test for dimensions "close" to dyadic (see help file) if( dyadic.2check( mm,nn,cut.min)==FALSE) {stop("error in length of column or row dimensions")} if (!inv) { while (nn > cut.min) { if (!transpose) { temp[1:nn, 1:mm] <- WQS2d.cylinder(temp[1:nn,1:mm], byX=byX) } else { temp[1:nn, 1:mm] <- WQS2di.cylinder(temp[1:nn, 1:mm], transpose = TRUE, byX=byX) } nn <- nn/2 mm <- mm/2 } } if (inv) { NN <- n MM <- m while (NN > cut.min) { NN <- NN/2 MM <- MM/2 } nn <- NN * 2 mm <- MM * 2 while (nn <= n) { if (!transpose) { temp[1:nn, 1:mm] <- WQS2di.cylinder(temp[1:nn, 1:mm], byX=byX) } else { temp[1:nn, 1:mm] <- WQS2d.cylinder(temp[1:nn, 1:mm], transpose = TRUE, byX=byX) } nn <- nn * 2 mm <- mm * 2 } } if (flip) { return(t(temp)) } else { return(temp) } }