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
    spheres3d(X, radius=0.1, col='black', add=TRUE)
    segments3d(X[p,], col='black', lwd=5, add=TRUE)
    
    # 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))
}

# load data from file
X = read.table("symphysis.txt",header=TRUE)

# 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='')
contour(x,x,S,add=TRUE)
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)

Homework 11. IM3 - bending energy