# Shape analysis¶

Series of homeworks from the course "Shape Analysis" teched by Prof. PaedDr. RNDr. Stanislav Katina, Ph.D. at Masaryk University - Science faculty.

In [3]:
%load_ext rmagic


# Homework 1. Tetragon rotation¶

In [4]:
%%R
seal = function(X){
return (cbind(rbind(X, X[1,])))
}

cplot = function(...){
# clean plot
plot(..., asp=1, axes=FALSE, xlab='', ylab='')
}

rotate = function(X, a){
T = matrix(c(
cos(a), sin(a),
-sin(a), cos(a)),
byrow=TRUE, nrow=2)
return(X %*% T)
}

# tetragon
X = rbind(c(-1,0), c(0,1), c(1,0), c(0,-1))

# plot tetragon
plot(seal(X), type='b', asp=1, pch=4, axes=FALSE, ylab='', xlab='')
lines(seal(rotate(X, 1)), type='b', asp=1, pch=5, lty=2)
lines(seal(rotate(X, -1)), type='b', asp=1, pch=6, lty=3)
legend('bottomleft', c('original', 'clockwise', 'counter-clockwise'), lty=1:3)


# Homework 2. Procrustes superimposition of two tetragons¶

In [5]:
%%R

rotate = function(X, a){
# rotation matrix
T = matrix(c(
cos(a), sin(a),
-sin(a), cos(a)),
byrow=TRUE, nrow=2)
return(X %*% T)
}

translate = function(X, t){
# translation by given vector
X = cbind(X, rep(1, dim(X)[1]))
T = matrix(c(
1,0,t[1],
0,1,t[2],
0,0,1
), byrow=TRUE, nrow=3
)
T = t(T)
return((X %*% T)[,c(1,2)])
}

scale = function(X, k){
# scale by constant or vector
if(typeof(k) == "double"){
k = rep(k, dim(X)[1])
}

T = matrix(c(
k[1],0,
0,k[2]
), byrow=TRUE, nrow=2
)
return(X %*% T)
}

cs = function(X){
# centroid size
X = opt_translation(X)
return(sqrt(sum(diag(X %*% t(X)))))
}

centroid = function(X){
# centroids
return(colMeans(X))
}

opt_translation = function(X){
# optimal translation
return(translate(X, -centroid(X)))
}

opt_scaling = function(X, Y){
# optimal scaling from X to Y
k = rep(cs(Y) / cs(X), 2)
return(scale(X, k))
}

opt_rotation = function(X, Y){
# optimal rotation from X to Y
s = svd(t(X) %*% Y)
G = s$u %*% s$v
return(X %*% G)
}

In [6]:
%%R -w 800
fplot = function(X1, X2, sub=''){
plot(seal(X1), type='b', asp=1, axes=FALSE, ylab='', xlab='', pch=1:4, sub=sub)
lines(seal(X2), type='b', lty=2, asp=1, axes=FALSE, ylab='', xlab='', pch=1:4)
}

X1 = rbind(c(-1,0), c(0,1), c(1,0), c(0,-1))
X2 = rbind(c(-1,0), c(0,0.7), c(1,0), c(0,-0.9))

# transform X2
X2 = translate(scale(rotate(X2, 1), 0.5), c(0.2,0.2))

par(mfrow=c(2,2))
# original coordinates
fplot(X1, X2, sub='original coordinates')

# center coordinates
fplot(X1, opt_translation(X2), sub='centered coordinates')

# centered and scaled coordinates
fplot(X1, opt_scaling(opt_translation(X2), X1), sub='center and scaled coordinates')

# centered, scaled and rotated coordinates
fplot(X1, opt_rotation(opt_scaling(opt_translation(X2), X1), X1), sub='centered, scaled and rotated coordinates')


# Homework 3. Procrustes superimposition of n tetragons¶

In [7]:
%%R
library(MASS)

apply3d = function(X, fun){
# apply function on matrices in third dimension
for(i in 1:dim(X)[3]){
X[,,i] = fun(X[,,i])
}
return(X)
}

procrust = function(M, threshold=0.00001){
# center and scale all matrices
M = apply3d(M, function(x) opt_translation(x))
M = apply3d(M, function(x) scale(x, 1./cs(x)))

# Procrust cycle
Xp = M[,,1]

while(TRUE){
last_Xp = Xp

# rotate to reference matrix
M = apply3d(M, function(x) opt_rotation(x, Xp))

# new reference matrix as a mean
Xp = apply(M, c(1,2), mean)

# break
cat('Procrust error:', sum(abs(last_Xp - Xp)), '\n')
if(sum(abs(last_Xp - Xp)) < threshold){
break
}
}
return(M)
}

gen_similar = function(Y, m=10){
# generate similar configuration matrices
M = array(0, dim=c(dim(Y)[1], dim(Y)[2], m))

for(k in 1:dim(Y)[1]){
M[k,,] = t(mvrnorm(n=m, mu=Y[k,], Sigma=diag(1,2) * 0.1^2))
}
return(M)
}

flatten_matrix = function(X){
return (cbind(as.vector(X[,1,]), as.vector(X[,2,])))
}

# original tetragons
X1 = rbind(c(-1,0), c(0,1), c(1,0), c(0,-1))
X2 = rbind(c(-1,0), c(0,2), c(1,0), c(0,-1))

M = gen_similar(X1, m=1000)
N = gen_similar(X2, m=1000)


### Original coordinates¶

In [8]:
%%R -w 800
fplot = function(M, ...){
plot(flatten_matrix(M), type = 'p', col = 'gray', pch=16, asp=1, axes=FALSE, ylab='', xlab='', ...)
# mean coordinates
X = apply(M, c(1,2), mean)
lines(seal(X), type='b', pch=1:4)
abline(v=0, col='red')
}

par(mfrow=c(1,2))
fplot(M, sub='original coordinates', xlim = c(-1.5, 1.5), ylim = c(-1.5, 2.5))
fplot(N, sub='original coordinates', xlim = c(-1.5, 1.5), ylim = c(-1.5, 2.5))


### Procrustes coordinates¶

In [9]:
%%R -w 800

PM = procrust(M)
PN = procrust(N)

par(mfrow=c(1,2))
fplot(PM, sub='Procrustes coordinates', ylim = c(-0.6, 0.7))
fplot(PN, sub='Procrustes coordinates', ylim = c(-0.6, 0.7))

Procrust error: 0.1011805
Procrust error: 0.0002462531
Procrust error: 6.217558e-07
Procrust error: 0.1635371
Procrust error: 0.0002536229
Procrust error: 3.983855e-07


### Rotated Procrustes coordinates¶

In [10]:
%%R -w 800
align_with_x = function(M, point1, point2){
# rotate shapes so that mean landmarks 2 and 4 lies on x-axis
v = rowMeans(M[point1,,]) - rowMeans(M[point2,,])
angle = atan(v[1] / v[2])
return(apply3d(M, function(x) rotate(x, angle)))
}

par(mfrow=c(1,2))
fplot(align_with_x(PM, 2, 4), sub='Rotated Procrustes coordinates', ylim = c(-0.6, 0.7))
fplot(align_with_x(PN, 2, 4), sub='Rotated Procrustes coordinates', ylim = c(-0.6, 0.7))


### Residual deviation of Procrustes coordinates¶

Let $SD(x)$ denote standard deviation in $x$ direction, $SD(y)$ sd in $y$ direction and $d$ distance from centroid. Then $k = \frac{\frac{1}{2} (SD(x) + SD(y))}{\sqrt{1 - d^2}}$ is supposed to be constant for all landmarks.

In [11]:
%%R
procrust_res = function(M){
sdx = apply(M[,1,],1,sd)
sdy = apply(M[,2,],1,sd)
cx = rowMeans(M[,1,])
cy = rowMeans(M[,2,])
c = apply(M, 2, mean)
d = sqrt((cx - c[1])^2 + (cy - c[2])^2)
sqrt_d = sqrt(1-d^2)
df = data.frame(d=d, sdx=sdx, sdy=sdy, '1-sqrt_d'=sqrt_d, k=0.5 * (sdx + sdy) / sqrt_d)
colnames(df) = c('d', 'SD(x)', 'SD(y)', '1-sqrt(d)', 'k')
return(df)
}

procrust_res(PN)

          d      SD(x)      SD(y) 1-sqrt(d)          k
1 0.3942879 0.02949443 0.03048001 0.9189870 0.03263073
2 0.6718438 0.02080161 0.02060487 0.7406929 0.02795118
3 0.3952984 0.02945890 0.02913935 0.9185527 0.03189705
4 0.4806373 0.02711922 0.02805762 0.8769195 0.03146060


# Homework 4. Shearing¶

In [12]:
%%R -w 1000
shear = function(X, ax=0, ay=0){
T = matrix(c(
1, tan(ay),
tan(ax), 1),
byrow=TRUE, nrow=2)
return(X %*% T)
}

# tetragon
X = rbind(c(0,0), c(0,1), c(1,1), c(1,0))

# shearing
par(mfrow=c(1,4))
plot(seal(X), type='b', asp=1, axes=FALSE, xlim=c(0,1.5),ylim=c(0,1.5), ylab='', xlab='', sub='original coordinates')
plot(seal(shear(X, ax=0.5)), type='b', asp=1, axes=FALSE, xlim=c(0,1.5), ylim=c(0,1.5), ylab='', xlab='', sub='horizontal shearing')
plot(seal(shear(X, ay=0.5)), type='b', asp=1, axes=FALSE, xlim=c(0,1.5), ylim=c(0,1.5), ylab='', xlab='', sub='vertical shearing')
plot(seal(shear(X, ax=0.5, ay=0.5)), type='b', asp=1, axes=FALSE, xlim=c(0,1.5), ylim=c(0,1.5), ylab='', xlab='', sub='horizontal and vertical shearing')


## Cube shearing¶

In [13]:
%%R
library(rgl)

shear3d_x = function(X, ax){
T = diag(3)
T[2,1] = tan(ax)
return(X %*% T)
}

shear3d_xy = function(X, ax, ay){
T = diag(3)
T[2,1] = tan(ax)
T[1,2] = tan(ay)
return(X %*% T)
}

shear3d_xyz = function(X, ax, ay, az){
T = diag(3)
T[3,1] = tan(ax)
T[3,2] = tan(ay)
T[2,3] = tan(az)
return(X %*% T)
}

save_3d_image = function(X, p, filename){
open3d()
rgl.viewpoint(theta=0,phi=0,fov=30)
bg3d('white')
par3d(windowRect=c(100,100,600,600))

# draw

# save as a file
snapshot3d(filename)
}

# cube
X1 = c(-1,-1,-1)
X2 = c(-1,-1,1)
X3 = c(-1,1,-1)
X4 = c(-1,1,1)
X5 = c(1,-1,-1)
X6 = c(1,-1,1)
X7 = c(1,1,-1)
X8 = c(1,1,1)
X = rbind(X1,X2,X3,X4,X5,X6,X7,X8)

# connect vertices
p = c()
for(i in 1:7){
for(j in (i+1):8){
if(sum(abs(X[i,] - X[j,])) == 2){
p = c(p, c(i,j))
}
}
}

save_3d_image(X, p, 'original.png')
save_3d_image(shear3d_x(X, 0.5), p, 'x.png')
save_3d_image(shear3d_xy(X, 0.5, 0.5), p, 'xy.png')
save_3d_image(shear3d_xyz(X, 0.5, 0.5, 0.5), p, 'xyz.png')

In [14]:
from IPython.display import Image


### Original¶

In [15]:
Image(filename='original.png', width=300, height=300)

Out[15]:

### Shearing in $(x_1, x_2)$ plane by $x_1$¶

In [16]:
Image(filename='x.png', width=300, height=300)

Out[16]:

### Shearing in $(x_1, x_2)$ plane¶

In [17]:
Image(filename='xy.png', width=300, height=300)

Out[17]:

### Shearing by all axes¶

In [18]:
Image(filename='xyz.png', width=300, height=300)

Out[18]:

# Homework 5. Tetragon reflection¶

In [19]:
%%R -w 1000
reflection = function(X, ax=1){
if(ax==2){
T = matrix(c(
1, 0,
0, -1),
byrow=TRUE, nrow=2)
} else {
T = matrix(c(
-1, 0,
0, 1),
byrow=TRUE, nrow=2)
}
return(X %*% T)
}

X = rbind(c(-1,0), c(-1,1), c(1,1.2), c(1,0))
par(mfrow=c(1,2))
plot(seal(X), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='original coordinates', pch=1:4)
abline(v=0, col='gray')
plot(seal(reflection(X, ax=1)), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='reflected coordinates', pch=1:4)
abline(v=0, col='gray')


## Reflection and relabeling of tetragon¶

In [20]:
%%R -w 800
Q = matrix(0,4,4)
Q[1,4] = 1
Q[2,3] = 1
Q[3,2] = 1
Q[4,1] = 1

X = rbind(c(-1,0), c(-1,1), c(1,1.2), c(1,0))
XR = reflection(X, ax=1)
print(Q)
par(mfrow=c(1,2))
plot(seal(X), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='original coordinates', pch=1:4)
abline(v=0, col='gray')
plot(seal(Q %*% XR), type='b', asp=1, axes=FALSE, ylab='', xlab='', sub='reflected and relabeled coordinates', pch=1:4)
abline(v=0, col='gray')

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    1
[2,]    0    0    1    0
[3,]    0    1    0    0
[4,]    1    0    0    0


# Homework 6. IM1¶

### 1. Is IM1 identical with IM1a?¶

After writing out model IM1 we get equivalent system of equations \begin{align} y &\sim Sw + 1_k c + xa \\ 0 &\sim 1_k w \\ 0 &\sim x^T a \end{align} model IM1a can be breaken down into the same system. The difference between IM1 a IM1a is only in interchange of rows/columns and they are equivalent to each other.

### 2. Why does the matrix $X_p$ has that form in IM1?¶

See first question. Matrix $X_p$ looks like that to fulfill conditions $0 = 1_k w$ a $0 = x^T a$.

### 3. How do you calculate $\hat{\beta}$ a $\hat{\beta^*}$?¶

Assuming full rank we get

\begin{align} \hat{\beta} &= (X_p^T X_p)^{-1} X_p^T y_p \\ \hat{\beta^*} &= (X_p^{*T} X_p^*)^{-1} X_p^{*T} y_p \end{align}

### 4. What function would you use to calculate the estimate?¶

model = lm(yp ~ Xp)
beta = summary(model)$coefficients  ### 5. Describe the relationship of$\hat{y}$to$y$?¶ $$\hat{y} = X_p\hat{\beta} = X_p (X_p^T X_p)^{-1} X_p^T y = Hy,$$ where$H$is called hat matrix or projective matrix In [21]: %%R -w 800 S_matrix = function(x,y){ # distane matrix k = length(x) S = 1/12 * abs(x %*% t(rep(1,k)) - rep(1,k) %*% t(x))^3 } IM1 = function(x, y){ k = length(x) # distance matrix S = S_matrix(x,y) # design matrix Xp = matrix(0,k+2,k+2) Xp[3:(k+2),3:(k+2)] = S Xp[k+1,1:k] = rep(1,k) Xp[k+2,1:k] = x Xp[1:k,k+1] = rep(1,k) Xp[1:k,k+2] = x # yp yp = c(y, c(0,0)) # linear model model = lm(yp ~ Xp - 1) # estiamtion of y yh = predict(model)[1:k] return(yh) } # sinus function n = 101 x = seq(0,1,length=n) y = sin(2*pi*x^3)^3 + 0.1 * rnorm(n) # IM1 model par(mfrow=c(1,2)) yh = IM1(x, y) plot(x, y, sub='IM1 model') lines(x, yh, col='blue') # spline method sp = spline(x,y,method='natural') plot(x, y, sub='spline(x,y,method=natural)') lines(sp$x, sp$y, col='red')  ### Radial basis function of spline IM1¶ In [22]: %%R library(rgl) save_3d = function(S, x, filename){ slim = range(S) ncols = 50 colorlut = terrain.colors(ncols + 1) # height color lookup table foo = trunc(ncols * (S - slim[1]) / (slim[2] - slim[1])) + 1 col = colorlut[foo] # assign colors to heights for each point open3d() rgl.viewpoint(theta=0,phi=30,fov=60) par3d(windowRect=c(100,100,600,600)) rgl.surface(x,x,S, color=col, back='lines') bg3d("white") decorate3d() snapshot3d(filename) } S = 10 * S_matrix(x,y) save_3d(S, x, 'basal.png')  In [23]: from IPython.display import Image Image(filename='basal.png')  Out[23]: In [24]: %%R S = S_matrix(x,y) image(x, x, S, asp=1, col=terrain.colors(30), axes=FALSE, xlab='', ylab='') contour(x,x,S,add=TRUE) abline(a=0,b=1, col='red', lw=2)  # Homework 7. PRM1¶ In [25]: %%R -w 800 -h 800 sqrt_matrix = function(a){ # hack positive definiteness by replacing negative values with zeros a.eig <- eigen(a) eig_values = pmax(a.eig$values, 0)
a.sqrt <- a.eig$vectors %*% diag(sqrt(eig_values)) %*% solve(a.eig$vectors)
return(a.sqrt)
}

PRM1 = function(x, y, lam, method=1, return_model=FALSE){
# method == 0 -> estimation by hat matrix
# method == 1 -> estimation with lm function
k = length(x)

# matrix Sp
S = S_matrix(x,x)
Sp = matrix(0,k+2,k+2)
Sp[3:(k+2),3:(k+2)] = S

# matrix X
R = sqrt_matrix(Sp)

# matrix Xdm
Xdm = cbind(rep(1,k), x, S)

if(method == 0){
# matrix H
H = Xdm %*% solve(t(Xdm) %*% Xdm + lam * Sp) %*% t(Xdm)

# estimate y
yh = H %*% y
return(yh)

} else if(method == 1){
# matrix Xp
Xp = rbind(Xdm, sqrt(lam) * R)

# vector yp
yp = c(y, rep(0, k+2))

# linear model
model = lm(yp ~ Xp - 1)
if(return_model){
return(model)
}else{
yh = predict(model)[1:k]
return(yh)
}
}
}

polyfit = function(x, y, d=12){
return(predict(lm(y ~ poly(x, d, raw=TRUE))))
}

lam = 4.774251e-6

par(mfrow=c(2,2))
plot(x, y, sub='smooth.spline')
sp = smooth.spline(x,y,all.knots = TRUE, cv=TRUE)
lines(sp$x, sp$y, col='blue')

plot(x, y, sub='PRM1')
lines(x, PRM1(x,y,lam=lam,method=0), col='blue')

plot(x, y, sub='12 degree polynomial')
lines(x, polyfit(x,y,d=12), col='blue')

plot(x, y, sub='PRM1a')
lines(x, PRM1(x,y,lam=lam,method=1), col='blue')

In [26]:
%%R -w 1000 -h 700
library(psych)

MSS = function(x,y,lam,fun){
# return MSS error for estimation PRM1 (assuming you know the original function)
yh = PRM1(x,y,lam)
# real function values
yr = fun(x)
return(mean((yh - yr)^2))
}

OCV = function(x, y, lam){
k = length(x)
model = PRM1(x, y, lam, return_model=TRUE)
hat = influence(model)$hat[1:k] yh = predict(model)[1:k] return(mean((yh - y)^2 / (1 - hat)^2)) } GCV = function(x, y, lam){ k = length(x) model = PRM1(x, y, lam, return_model=TRUE) hat = influence(model)$hat[1:k]

yh = predict(model)[1:k]
return(k*(sum((yh - y)^2) / tr(diag(k) - hat)^2))
}

plot_cv = function(title, x, y, objfun, yfun,...){
fun = function(lam) objfun(x,y,lam,...)

# find minimum lambda
f = function(lam) objfun(x,y,exp(lam), ...)
opt = optimize(f=f, lower=-20, upper=0)
opt_lam = exp(opt$minimum) # plot surrounding lambdas lambdas = seq(opt_lam /5, opt_lam * 5, length.out=30) error = lapply(lambdas, fun) # plot lambdas plot(lambdas, error, type='l', ylab=paste(title, '(lambda)'), xlab='lambda', sub=paste('lambda(MSS)=', sprintf("%.7f",opt_lam))) points(opt_lam, fun(opt_lam), pch=19, lwd=5) # plot spline plot(x, y, sub=paste('penalized regression spline', title)) lines(x, PRM1(x,y,opt_lam,method=1), col='blue') lines(x, yfun(x), col='red') legend('topleft', c('E[y]', 'odhad y'), col=c('red', 'blue'), lty=1) } # MSS fun = function(x) sin(2*pi*x^3)^3 par(mfcol=c(2,3)) plot_cv('MSS', x, y, MSS, fun, fun) plot_cv('OCV', x, y, OCV, fun) plot_cv('GCV', x, y, GCV, fun)  # Homework 8. IM1 (Curve resampling)¶ In [27]: %%R -w 800 cumchord = function(X){ # cumulative chord distance return(cumsum(sqrt(apply((X-rbind(X[1,],X[-(nrow(X)),]))^2,1,sum)))) } resample = function(X, method, eps=1e-8){ x = X[,1] y = X[,2] i = 0 while(TRUE){ i = i + 1 # cumulative chordal distance CH <- cumchord(cbind(y,x)) # if distance is small enough, stop dCH = diff(CH) error = sum(abs(dCH - mean(dCH))) if(error < eps) break # new splines if(method == 'spline'){ x = spline(CH, x, method="natural", n=50)$y
y = spline(CH, y, method="natural", n=50)$y } else if(method == 'approx'){ x = spline(CH, x, n=50)$y
y = spline(CH, y, n=50)$y } } cat(i, 'iterations\n') return(cbind(x,y)) } interpolation = function(X, n=50){ CH <- cumchord(cbind(X[,2],X[,1])) x = spline(CH, X[,1], n=n)$y
y = spline(CH, X[,2], n=n)\$y
return(cbind(x,y))
}

# draw resampled curve
par(mfrow=c(1,3))
plot(interpolation(X), type='l', asp=1, axes=FALSE, xlab='', ylab='', sub='interpolation')
points(X[,1], X[,2])
plot(resample(X, method='spline'), type='o', asp=1, axes=FALSE, xlab='', ylab='', sub='spline', col='blue')
plot(resample(X, method='approx'), type='o', asp=1, axes=FALSE, xlab='', ylab='', sub='approx', col='red')

7 iterations
7 iterations


# Homework 9. IM3¶

In [28]:
%%R
# phi distance function
phi = function(x){
d = sqrt(sum(x^2))
if(d == 0){
return(0)
}else{
return(d^2 * log(d))
}
}

S_matrix3 = function(X,Y){
# return distance matrix between rows of X a Y
m = dim(X)[1]
n = dim(Y)[1]
S = matrix(, nrow = m, ncol = n)
for(i in 1:m){
for(j in 1:n){
S[i,j] = phi(X[i,] - Y[j,])
}
}
return(S)
}

S = S_matrix3(X,X)
save_3d(S, seq(0,2000,length=dim(S)[1]), 'basal3d.png')


### Radial basis function of IM3¶

In [29]:
from IPython.display import Image
Image(filename='basal3d.png')

Out[29]:
In [30]:
%%R
S = S_matrix3(X,X)
x = seq(0,2000,length=dim(S)[1])
image(x, x, S, asp=1, col=terrain.colors(30), axes=FALSE, xlab='', ylab='')
abline(a=0,b=1, col='red', lw=2)


# Homework 10. IM3 - TPS net¶

In [31]:
%%R -w 1000
net = function(X, kx=20, margin=0.2){
x1 = min(X[,1])
x2 = max(X[,1])
y1 = min(X[,2])
y2 = max(X[,2])

# add 1/5 margin to all directions
x_range = x2 - x1
y_range = y2 - y1
x1 = x1 - margin * x_range
x2 = x2 + margin * x_range
y1 = y1 - margin * y_range
y2 = y2 + margin * y_range

x_grid = seq(x1,x2,length=kx)

xmin = min(x_grid)
xmax = max(x_grid)
hy = (xmax - xmin) / (kx - 1)

y_grid = seq(y1,y2,by=hy)
ky = length(y_grid)

Xcp <- as.matrix(expand.grid(x_grid,y_grid))
return(Xcp)
}

draw_net = function(X, Xcp, kx=20, net_type='lines', type='b', ...){
# show the shape
X = cbind(rbind(X, X[1,]))
plot(X, asp=1, type=type, xlim=range(Xcp[,1]), ylim=range(Xcp[,2]), axes=FALSE, xlab='', ylab='', pch=19, ...)

# show net
if(net_type == 'lines'){
ky = dim(Xcp)[1] / kx
for(i in 1:ky){
lines(Xcp[1:kx + (i-1)*kx,], lwd=0.3)
}
for(i in 1:kx){
lines(Xcp[(1:ky)*kx-i+1,], lwd=0.3)
}
} else if(net_type == 'points'){
points(Xcp, pch=19, cex=0.2)
}
}

# shapes
X1 = rbind(c(-0.5,-0.5),c(-0.5,0.5),c(0.5,0.5),c(0.5,-0.5))
X2 = rbind(c(-0.5,-0.65),c(-0.5,0.65),c(0.5,0.65),c(0.5,-0.65))

# TPS net
par(mfrow=c(1,4))
draw_net(X1, net(X1, kx=20), kx=20, net_type='points')
draw_net(X2, net(X2, kx=20), kx=20, net_type='points')
draw_net(X1, net(X1, kx=20), kx=20, net_type='lines')
draw_net(X2, net(X2, kx=20), kx=20, net_type='lines')

In [32]:
%%R -w 800
shear = function(X, ax=0, ay=0){
T = matrix(c(
1, tan(ay),
tan(ax), 1),
byrow=TRUE, nrow=2)
return(X %*% T)
}

X = rbind(c(-0.5,-0.5), c(-0.5,0.5), c(0.5,0.5), c(0.5, -0.5))
XN = net(X)
par(mfrow=c(1,2))
draw_net(shear(X, ax=0.5), shear(XN, ax=0.5))
draw_net(shear(X, ay=0.5), shear(XN, ay=0.5))

In [33]:
%%R -w 1000
extrapolation = function(X,Y,l){
return(Y + l * (Y - X))
}

IM3_Xp = function(X,Y){
k = dim(X)[1]
S = S_matrix3(X,X)
Xp = rbind(cbind(S, rep(1,k), X),
c(rep(1,k), rep(0,3)),
cbind(t(X), matrix(0,2,3)))
return(Xp)
}

IM3 = function(X,Y,Xcp){
k = dim(X)[1]

if(is.vector(Y)){
# IM2 model
Yp = c(Y,rep(0,3))
} else {
# IM3 model
Yp = rbind(Y,matrix(0,3,2))
}

Xp = IM3_Xp(X,Y)

# linear model a its coefficients
model = lm(Yp ~ Xp - 1)
beta = coef(model)

# new transformed matrix Xcp
Scp = S_matrix3(Xcp,X)

Xpcp = rbind(cbind(Scp, rep(1,k), Xcp),
c(rep(1,k), rep(0,3)),
cbind(t(X), matrix(0,2,3)))

Ypcp = Xpcp %*% beta

return(Ypcp[1:dim(Xcp)[1],])
}

plot_extrapolation = function(X,Y,l, only_net=FALSE, kx=20, ...){
Z = extrapolation(X,Y,l)

# shapes
if(!only_net){
cplot(seal(X), type='b', lty=2, sub=paste('l =', l))
lines(seal(Z), type='b')

# arrows
for(i in 1:dim(Z)[1]){
arrows(X[i,1], X[i,2], Z[i,1], Z[i,2], code=2, length=0.1, col='blue')
}
}

# net
Zcp = IM3(X, Z, net(X, kx=kx))
draw_net(Z, Zcp, sub=paste('l =', l), kx=kx, ...)
}

X1 = rbind(c(0,-0.5), c(-0.5,0), c(0,0.5), c(0.5,0))
X2 = rbind(c(0,-0.51), c(-0.51,0), c(0,0.4), c(0.51,0))

par(mfrow=c(1,2))
plot_extrapolation(X1, X2, 0)
plot_extrapolation(X1, X2, 1)

In [34]:
%%R -w 800
par(mfrow=c(1,3))
plot_extrapolation(X1, X2, 1, only_net=TRUE)
plot_extrapolation(X1, X2, 4, only_net=TRUE)
plot_extrapolation(X1, X2, 8, only_net=TRUE)