Create a bidimensional regression tool - this creates the transformation function given x1, y1 and x2, y2 - these are vectors of coordinates of the initial and transformed data. Aim is to create a ‘best fit’ surface mapping \((x_1,y_1) \mapsto (x_2,y_2)\).
Code below defines all the functions needed - the key ones are transform_linear
, transform_loess
, transform_akima
and make_grid
. The first two, given a set x1
,y1
and x2
, y2
create a function that applies the best fit bidimensional transform to a new matrix of points Xin
whose first column are \(x\)-coordinates and whose second column is \(y\)-coordinates (in initial space). transform_akima
does a similar thing but the transform is an interpolation - so that the transformed space maps \((x_1,y_1)\) exactly onto \((x_2,y_2)\). make_grid
creates a grid in this two column matrix format. NA
’s are strategically inserted so that a grid is drawn if the grid is given as the argument to lines
or plot
with type='l'
.
library(akima)
ak <- function(x1,y1,x2,y2,Xin,extrap=TRUE,linear=TRUE) {
use <- apply(Xin,1,function(z) ! any(is.na(z)))
xinterp <- yinterp <- rep(NA,nrow(Xin))
xinterp[use] <- interpp(x1,y1,x2,Xin[use,1],Xin[use,2],extrap=extrap,linear=linear)$z
yinterp[use] <- interpp(x1,y1,y2,Xin[use,1],Xin[use,2],extrap=extrap,linear=linear)$z
cbind(xinterp,yinterp)
}
transform_akima <- function(x1,y1,x2,y2,extrap=TRUE,linear=TRUE) {
function(Xin){
ak(x1,y1,x2,y2,Xin,extrap,linear)
}
}
affmat <- function(x1,y1) {
N <- length(x1)
N2 <- 2*N
dMat <- matrix(0,N2,6)
top <- 1:N
btm <- top + N
dMat[top,1] <- 1
dMat[btm,2] <- 1
dMat[top,3] <- x1
dMat[btm,4] <- x1
dMat[top,5] <- y1
dMat[btm,6] <- y1
colnames(dMat) <- c("ix","iy","xx","xy","yx","yy")
as.data.frame(dMat)
}
bidim_linear <- function(x1,y1,x2,y2) {
dMat <- affmat(x1,y1)
y <- c(x2,y2)
mdl <- lm(y~ix+iy+xx+xy+yx+yy+0,data=dMat)
mdl
}
transform_linear <- function(x1,y1,x2,y2) {
mdl <- bidim_linear(x1,y1,x2,y2)
function(Xin) {
dMat <- as.data.frame(affmat(Xin[,1],Xin[,2]))
res <- predict(mdl,newdata=dMat)
top <- 1:(length(res) %/% 2)
btm <- top + length(top)
cbind(res[top],res[btm])
}
}
make_grid <- function(xmin,xmax,xstride,ymin,ymax,ystride) {
xseq <- seq(xmin,xmax,by=xstride)
yseq <- seq(ymin,ymax,by=ystride)
gr1 <- as.matrix(expand.grid(x=c(xseq,NA),y=yseq))
gr2 <- as.matrix(expand.grid(y=c(yseq,NA),x=xseq))[,2:1]
rbind(gr1,gr2)}
transform_loess <- function(x1,y1,x2,y2,...) {
mdl1 <- loess(x2~x1+y1,...)
mdl2 <- loess(y2~x1+y1,...)
function(Xin) {
dMat <- data.frame(x1=Xin[,1],y1=Xin[,2])
res <- cbind(predict(mdl1,newdata=dMat),predict(mdl2,newdata=dMat))
res
}
}
First create a test data set:
set.seed(4499864)
x1 <- runif(20)
y1 <- runif(20)
x2 <- x1 + 0.1*y1 + rnorm(20)*0.01
y2 <- y1 - 0.1*x1 + rnorm(20)*0.01
plot it …
plot(x1,y1,pch=16,asp=1)
points(x2,y2,col='red',pch=16)
Compute the linear transform
custard <- transform_linear(x1,y1,x2,y2)
apply and plot
plot(custard(make_grid(0,1,0.1,0,1,0.1)),type='l',asp=1)
Compute the loess transform
mustard <- transform_loess(x1,y1,x2,y2)
apply and plot
plot(mustard(make_grid(0,1,0.1,0,1,0.1)),type='l',asp=1)
Compute the Akima transform
bustard <- transform_akima(x1,y1,x2,y2,extrap = FALSE, linear=FALSE)
apply and plot
plot(bustard(make_grid(0,1,0.1,0,1,0.1)),type='l',asp=1)
Note that Akima ‘chops out’ grid squares not in the convex hull of the calibration points.
One alternative is to add a set of control points at the corners of the area that are unchanged by the transformation, and to re-calibrate:
cpts <- matrix(c(
0,0,
1,0,
1,1,
0,1), 4,2,byrow=TRUE)
tpts <- custard(cpts)
x1 <- c(x1,cpts[,1])
y1 <- c(y1,cpts[,2])
x2 <- c(x2,tpts[,1])
y2 <- c(y2,tpts[,2])
bustard <- transform_akima(x1,y1,x2,y2,extrap = FALSE, linear=FALSE)
plot(bustard(make_grid(0,1,0.1,0,1,0.1)),type='l',asp=1)