1## Make RK for thin-plate splines 2mkrk.tp <- function(dm,order,mesh,weight=1) 3{ 4 ## Check inputs 5 if (!((2*order>dm)&(dm>=1))) { 6 stop("gss error: thin-plate spline undefined for the parameters") 7 } 8 if (xor(is.vector(mesh),dm==1) 9 |xor(is.matrix(mesh),dm>=2)) { 10 stop("gss error in mkrk.tp: mismatched inputs") 11 } 12 if ((min(weight)<0)|(max(weight)<=0)) { 13 stop("gss error in mkrk.tp: negative weights") 14 } 15 ## Set weights 16 if (is.vector(mesh)) N <- length(mesh) 17 else N <- dim(mesh)[1] 18 weight <- rep(weight,len=N) 19 weight <- sqrt(weight/sum(weight)) 20 ## Obtain orthonormal basis 21 phi.p <- mkphi.tp.p(dm,order) 22 nnull <- choose(dm+order-1,dm) 23 s <- NULL 24 for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env)) 25 s <- qr(weight*s) 26 if (s$rank<nnull) { 27 stop("gss error in mkrk.tp: insufficient normalizing mesh for thin-plate spline") 28 } 29 q <- qr.Q(s) 30 r <- qr.R(s) 31 ## Set Q^{T}E(|u_{i}-u_{j}|)Q 32 rk.p <- mkrk.tp.p(dm,order) 33 pep <- weight*t(weight*rk.p$fun(mesh,mesh,rk.p$env,out=TRUE)) 34 pep <- t(q)%*%pep%*%q 35 ## Create the environment 36 env <- list(dim=dm,order=order,weight=weight, 37 phi.p=phi.p,rk.p=rk.p,q=q,r=r,mesh=mesh,pep=pep) 38 ## Create the rk function 39 fun <- function(x,y,env,outer.prod=FALSE) { 40 ## Check inputs 41 if (env$dim==1) { 42 if (!(is.vector(x)&is.vector(y))) { 43 stop("gss error in rk: inputs are of wrong types") 44 } 45 nx <- length(x) 46 ny <- length(y) 47 } 48 else { 49 if (is.vector(x)) x <- t(as.matrix(x)) 50 if (env$dim!=dim(x)[2]) { 51 stop("gss error in rk: inputs are of wrong dimensions") 52 } 53 nx <- dim(x)[1] 54 if (is.vector(y)) y <- t(as.matrix(y)) 55 if (env$dim!=dim(y)[2]) { 56 stop("gss error in rk: inputs are of wrong dimensions") 57 } 58 ny <- dim(y)[1] 59 } 60 ## Return the results 61 nnull <- choose(env$dim+env$order-1,env$dim) 62 if (outer.prod) { 63 phix <- phiy <- NULL 64 for (nu in 1:nnull) { 65 phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env)) 66 phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env)) 67 } 68 phix <- backsolve(env$r,phix,transpose=TRUE) 69 phiy <- backsolve(env$r,phiy,transpose=TRUE) 70 ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE) 71 ex <- env$weight*ex 72 ex <- t(env$q)%*%ex 73 ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE) 74 ey <- env$weight*ey 75 ey <- t(env$q)%*%ey 76 env$rk.p$fun(x,y,env$rk.p$env,out=TRUE)-t(phix)%*%ey- 77 t(ex)%*%phiy+t(phix)%*%env$pep%*%phiy 78 } 79 else { 80 N <- max(nx,ny) 81 phix <- phiy <- NULL 82 for (nu in 1:nnull) { 83 phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env)) 84 phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env)) 85 } 86 phix <- backsolve(env$r,phix,transpose=TRUE) 87 phix <- matrix(phix,nnull,N) 88 phiy <- backsolve(env$r,phiy,transpose=TRUE) 89 phiy <- matrix(phiy,nnull,N) 90 ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE) 91 ex <- env$weight*ex 92 ex <- t(env$q)%*%ex 93 ex <- matrix(ex,nnull,N) 94 ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE) 95 ey <- env$weight*ey 96 ey <- t(env$q)%*%ey 97 ey <- matrix(ey,nnull,N) 98 fn1 <- function(x,n) x[1:n]%*%x[n+(1:n)] 99 fn2 <- function(x,pep,n) t(x[1:n])%*%pep%*%x[n+(1:n)] 100 env$rk.p$fun(x,y,env$rk.p$env)-apply(rbind(phix,ey),2,fn1,nnull)- 101 apply(rbind(phiy,ex),2,fn1,nnull)+ 102 apply(rbind(phix,phiy),2,fn2,env$pep,nnull) 103 } 104 } 105 ## Return the function and the environment 106 list(fun=fun,env=env) 107} 108 109## Make phi function for thin-plate splines 110mkphi.tp <- function(dm,order,mesh,weight) 111{ 112 ## Check inputs 113 if (!((2*order>dm)&(dm>=1))) { 114 stop("gss error: thin-plate spline undefined for the parameters") 115 } 116 if (xor(is.vector(mesh),dm==1) 117 |xor(is.matrix(mesh),dm>=2)) { 118 stop("gss error in mkphi.tp: mismatched inputs") 119 } 120 if ((min(weight)<0)|(max(weight)<=0)) { 121 stop("gss error in mkphi.tp: negative weights") 122 } 123 ## Set weights 124 if (is.vector(mesh)) N <- length(mesh) 125 else N <- dim(mesh)[1] 126 weight <- rep(weight,len=N) 127 weight <- sqrt(weight/sum(weight)) 128 ## Create the environment 129 phi.p <- mkphi.tp.p(dm,order) 130 nnull <- choose(dm+order-1,dm) 131 s <- NULL 132 for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env)) 133 s <- qr(weight*s) 134 if (s$rank<nnull) { 135 stop("gss error in mkphi: insufficient normalizing mesh for thin-plate spline") 136 } 137 r <- qr.R(s) 138 env <- list(dim=dm,order=order,phi.p=phi.p,r=r) 139 ## Create the phi function 140 fun <- function(x,nu,env) { 141 nnull <- choose(env$dim+env$order-1,env$dim) 142 phix <- NULL 143 for(i in 1:nnull) 144 phix <- rbind(phix,env$phi.p$fun(x,i,env$phi.p$env)) 145 t(backsolve(env$r,phix,transpose=TRUE))[,nu+1] 146 } 147 ## Return the function and the environment 148 list(fun=fun,env=env) 149} 150 151## Make pseudo RK for thin-plate splines 152mkrk.tp.p <- function(dm,order) 153{ 154 ## Check inputs 155 if (!((2*order>dm)&(dm>=1))) { 156 stop("gss error: thin-plate spline undefined for the parameters") 157 } 158 ## Create the environment 159 if (dm%%2) { 160 theta <- gamma(dm/2-order)/2^(2*order)/pi^(dm/2)/gamma(order) 161 } 162 else { 163 theta <- (-1)^(dm/2+order+1)/2^(2*order-1)/pi^(dm/2)/ 164 gamma(order)/gamma(order-dm/2+1) 165 } 166 env <- list(dim=dm,order=order,theta=theta) 167 ## Create the rk.p function 168 fun <- function(x,y,env,outer.prod=FALSE) { 169 ## Check inputs 170 if (env$dim==1) { 171 if (!(is.vector(x)&is.vector(y))) { 172 stop("gss error in rk: inputs are of wrong types") 173 } 174 } 175 else { 176 if (is.vector(x)) x <- t(as.matrix(x)) 177 if (env$dim!=dim(x)[2]) { 178 stop("gss error in rk: inputs are of wrong dimensions") 179 } 180 if (is.vector(y)) y <- t(as.matrix(y)) 181 if (env$dim!=dim(y)[2]) { 182 stop("gss error in rk: inputs are of wrong dimensions") 183 } 184 } 185 ## Return the results 186 if (outer.prod) { 187 if (env$dim==1) { 188 fn1 <- function(x,y) abs(x-y) 189 d <- outer(x,y,fn1) 190 } 191 else { 192 fn2 <- function(x,y) sqrt(sum((x-y)^2)) 193 d <- NULL 194 for (i in 1:dim(y)[1]) d <- cbind(d,apply(x,1,fn2,y[i,])) 195 } 196 } 197 else { 198 if (env$dim==1) d <- abs(x-y) 199 else { 200 N <- max(dim(x)[1],dim(y)[1]) 201 x <- t(matrix(t(x),env$dim,N)) 202 y <- t(matrix(t(y),env$dim,N)) 203 fn <- function(x) sqrt(sum(x^2)) 204 d <- apply(x-y,1,fn) 205 } 206 } 207 power <- 2*env$order-env$dim 208 switch(1+env$dim%%2, 209 env$theta*d^power*log(ifelse(d>0,d,1)), 210 env$theta*d^power) 211 } 212 ## Return the function and the environment 213 list(fun=fun,env=env) 214} 215 216## Make pseudo phi function for thin-plate splines 217mkphi.tp.p <- function(dm,order) 218{ 219 ## Check inputs 220 if (!((2*order>dm)&(dm>=1))) { 221 stop("gss error: thin-plate spline undefined for the parameters") 222 } 223 ## Create the environment 224 pol.code <- NULL 225 for (i in 0:(order^dm-1)) { 226 ind <- i; code <- NULL 227 for (j in 1:dm) { 228 code <- c(code,ind%%order) 229 ind <- ind%/%order 230 } 231 if (sum(code)<order) pol.code <- cbind(pol.code,code) 232 } 233 env <- list(dim=dm,pol.code=pol.code) 234 ## Create the phi function 235 fun <- function(x,nu,env) { 236 if (env$dim==1) x <- as.matrix(x) 237 if (env$dim!=dim(x)[2]) { 238 stop("gss error in phi: inputs are of wrong dimensions") 239 } 240 apply(t(x)^env$pol.code[,nu],2,prod) 241 } 242 ## Return the function and the environment 243 list(fun=fun,env=env) 244} 245 246## Make RK for spherical splines 247mkrk.sphere <- function(order) 248{ 249 ## Create the environment 250 env <- list(order=order) 251 ## Create the rk function 252 fun <- function(x,y,env,outer.prod=FALSE) { 253 ##% Check the inputs 254 if (is.vector(x)) x <- t(as.matrix(x)) 255 if (is.vector(y)) y <- t(as.matrix(y)) 256 if (!(is.matrix(x)&is.matrix(y))) { 257 stop("gss error in rk: inputs are of wrong types") 258 } 259 if ((dim(x)[2]!=2)|(dim(y)[2]!=2)) { 260 stop("gss error in rk: inputs are of wrong dimensions") 261 } 262 if ((max(abs(x[,1]),abs(y[,1]))>90)|(max(abs(x[,2]),abs(y[,2]))>180)) { 263 stop("gss error in rk: inputs are out of range") 264 } 265 ##% Convert to radian 266 lat.x <- x[,1]/180*pi; lon.x <- x[,2]/180*pi 267 lat.y <- y[,1]/180*pi; lon.y <- y[,2]/180*pi 268 ##% Return the result 269 rk <- function(lat.x,lon.x,lat.y,lon.y,order) { 270 z <- cos(lat.x)*cos(lat.y)*cos(lon.x-lon.y)+sin(lat.x)*sin(lat.y) 271 W <- ifelse(z<1-10^(-10),(1-z)/2,0) 272 A <- ifelse(W>0,log(1+1/sqrt(W)),0) 273 C <- ifelse(W>0,2*sqrt(W),0) 274 switch(order-1, 275 (A*4*W*(3*W-1)+6*W*(1-C)+1)/2, 276 (W*W*(A*((840*W-720)*W+72)+420*W*(1-C)+220*C-150)-4*W+3)/12, 277 (W*W*W*(A*(((27720*W-37800)*W+12600)*W-600)+ 278 (13860*(1-C)*W+14280*C-11970)*W-2772*C+1470)+ 279 15*W*W-3*W+5)/30) - 1/(2*order-1) 280 } 281 if (outer.prod) { 282 zz <- NULL 283 for (i in 1:length(lat.y)) 284 zz <- cbind(zz,rk(lat.x,lon.x,lat.y[i],lon.y[i],env$order)) 285 } 286 else zz <- rk(lat.x,lon.x,lat.y,lon.y,env$order) 287 zz 288 } 289 ## Return the function and the environment 290 list(fun=fun,env=env) 291} 292