Aim

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)\).

The code:

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
  }
}

Demo

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)