遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

Least Square Conformal Mapping

Least Square Conformal Mapping 共形変換 曲面


  • ３次元空間に埋め込まれた曲面を三角形分割し、曲面の周囲に２次元座標を与え、曲面の非周辺点に共形変換に近くなるように２次元座標を与える方法。「共形変換に近い」かどうかの基準は最小二乗法を採用し、その解を求めるにあたって線形代数的に一意に最適化する
  • ２次元平面と共形変換
    • まずは、平面に格子を作って描く。
# xy二次元平面のグリッド点を作る
x1 <- y1 <- seq(from=-10,to=10,length=31)
X <- as.matrix(expand.grid(x1,y1))
#X <- X %*% matrix(c(cos(pi/6),sin(pi/6),-sin(pi/6),cos(pi/6)),2,2)
# 点の数
n.v <- length(X[,1])
# 格子のエッジリストを作る
distX <- as.matrix(dist(X,method="manhattan"))
distX <- abs(distX-abs(x1[1]-x1[2])) < 0.01
lattice.edge <- which(distX,arr.ind=TRUE)
# 正方形周囲の点と内部の点とを分ける
v.fix <- which(apply((X-10)*(X+10),1,prod)==0)
v.free <- which(apply((X-10)*(X+10),1,prod)!=0)

col <- rep(1,n.v)
col[v.fix] <- 2
plot(X,pch=20,col=col)
segments(X[lattice.edge[,1],1],X[lattice.edge[,1],2],X[lattice.edge[,2],1],X[lattice.edge[,2],2])
  • 今、この平面をぐにゃりと曲げる。曲げるにあたって、複素関数を使う。
  • 二次元平面の点$(x,y)$を$w=x + i y$と見立て、$f(w)$と言う複素関数の実部と虚部とを変形後の$(x,y)座標とすることにする。

$f(w)=(w+p)^k$という関数にしてみよう。

p <- 12+1i*7;k <- 0.4
w <- X[,1] + 1i*X[,2]
new.w <- (w+p)^k
#new.w <- w
new.X <- cbind(Re(new.w),Im(new.w))
plot(new.X,pch=20,col=col)
segments(new.X[lattice.edge[,1],1],new.X[lattice.edge[,1],2],new.X[lattice.edge[,2],1],new.X[lattice.edge[,2],2])
  • この変形の特徴は、各所で角度が保存されていることであり、共形変換を角度保存の変換と呼ぶ。
  • ３次元空間に埋め込まれた曲面
  • 上の例では、２次元平面内で複素数を用いて格子をぐにゃりと曲げたが、二次元平面内に収まっている。これに３次元座標も与えてやると、３次元空間内の曲面ができる。
p <- 12+1i*7;k <- 0.4
w <- X[,1] + 1i*X[,2]
new.w <- (w+p)^k
#new.w <- w
new.z <- 0.4*exp(-0.02*Mod(w)^2)
new.z <- 0.4*exp(-0.1*(X[,1]^2+X[,2]^2))
new.X <- cbind(Re(new.w),Im(new.w),new.z)
plot3d(new.X,col=col)
segments3d(new.X[c(t(lattice.edge)),])
  • 三角化
  • 格子は座標系の変形を視覚化するのに有用だが、局面データの離散的な取扱いのためには、三角形メッシュにするのが好都合である。
  • Rではgeometryパッケージを使うと三角化が可能である。
library(geometry)
tc <- delaunayn(X[,1:2])
n.t <- length(tc[,1]) # 三角形の個数
plot(new.X,col=col,pch=20)
segments(new.X[tc[,1],1],new.X[tc[,1],2],new.X[tc[,2],1],new.X[tc[,2],2])
segments(new.X[tc[,2],1],new.X[tc[,2],2],new.X[tc[,3],1],new.X[tc[,3],2])
segments(new.X[tc[,3],1],new.X[tc[,3],2],new.X[tc[,1],1],new.X[tc[,1],2])
plot3d(new.X,col=col)
triangles3d(new.X[c(t(tc)),],col=rep(2,n.t))
  • 三角形の面積とベクトルのクロス積と面の向き
  • 3次元空間に三角形があり、３頂点座標が$(x_1,y_1,z_1),(x_2,y_2,z_2),(x_3,y_3,z_3)$であるとする。
  • 三角形の頂点の並べ方を$(v1,v2,v3)$とみるとき、頂点v1を基点として、$v2-v1$ベクトルと$v3-v1$ベクトルのクロス積を計算すると、その絶対値は三角形の面積の2倍(２つのベクトルが張る平行四辺形の面積)となり、その符号は2ベクトルの外積ベクトルが2ベクトルと右手系か左手系の関係になるかを示す。したがって、以下に示すように、$(v1,v2,v3)$の順で考える場合と、$(v1,v3,v2)$で考える場合とで、クロス積の値は符号が反転する。
v1 <- c(0,0,0);v2 <- c(2,0,0);v3 <- c(1,2,3)
v21 <- v2 - v1
v31 <- v3 - v2
# (v1,v2,v3)順
tmp <- extprod3d(v21,v31)
A1 <- sum(tmp)
# (v1,v3,v2)順
tmp <- extprod3d(v31,v21)
A2 <- sum(tmp)
print(A1)
print(A2)
  • ４元数の虚部を用いて3次元ベクトルを表し、同様に符号付き面積を計算することもできる。
  • ４元数に慣れる意味で、その計算も行っておく。
library(onion)
H <- c(Hi,Hj,Hk) # ４元数虚部基底
v1q <- sum(v1*H);v2q <- sum(v2*H);v3q <- sum(v3*H)
v21q <- v2q - v1q
v31q <- v3q - v2q
# (v1,v2,v3)順
tmp <- v21q * v31q
A1q <- i(tmp) + j(tmp) + k(tmp)

# (v1,v3,v2)順
tmp <- v31q * v21q
A2q <- i(tmp) + j(tmp) + k(tmp)

print(A1q)
print(A2q)
  • すべての三角形について符号付き面積を計算してみる。
my.area.tri <- function(vs){
  v1 <- vs[1,];v2 <- vs[2,];v3 <- vs[3,]
  v21 <- v2 - v1
  v31 <- v3 - v1
  tmp <- sum(extprod3d(v21,v31))
  sum(tmp)
}
  • 符号付き面積の符号で三角形を塗り分けてみると、以下のように、まだらになる。
As <- rep(0,n.t)
for(i in 1:n.t){
  As[i] <- my.area.tri(new.X[tc[i,],])
}

plot3d(new.X)
col <- rep(2.5,n.t) + sign(As)*0.5
triangles3d(new.X[c(t(tc)),],col=col)
  • すべての三角形の向きをそろえたいので、この符号付き面積の値をもとに、裏向きの三角形の頂点順序を変えることにする。変えたあとで、面積を計算しなおそう。
ura.tri <- which(As<0)
tc[ura.tri,2:3] <- tc[ura.tri,c(3,2)]
As <- rep(0,n.t)
for(i in 1:n.t){
  As[i] <- my.area.tri(new.X[tc[i,],])
}

plot3d(new.X)
col <- rep(2.5,n.t) + sign(As)*0.5
triangles3d(new.X[c(t(tc)),],col=col)
  • 最小二乗法による法線ベクトル推定と共形変換推定 その１
  • 複素数行列を作る
  • サイズは（２ｘ面の数）ｘ（２ｘ点の数）とすることで、実行列として作る
  • 面と点とのペアごとに２ｘ２複素行列を取り、それを配置したもの。
my.orthobasis <- function(q){
  q21 <- q[2]-q[1]
  q31 <- q[3]-q[1]
  
  v1 <- q21
  v3 <- Im(v1*q31)
  v2 <- Im(v3*v1)
  c(v1/Mod(v1),v2/Mod(v2),v3/Mod(v3))
}
v1 <- c(0,0,0);v2 <- c(2,0,0);v3 <- c(1,2,3)
H <- c(Hi,Hj,Hk) # ４元数虚部基底
v1q <- sum(v1*H);v2q <- sum(v2*H);v3q <- sum(v3*H)
q <- c(v1q,v2q,v3q)
q.basis <- my.orthobasis(q)
q.basis
t(q.basis) %*% q.basis # 正規直交基底であることの確認

# 頂点座標の四元数表現
Xq <- new.X[,1]*Hi + new.X[,2]*Hj + new.X[,3]*Hk
basis.tri <- list()
for(i in 1:n.t){
  s <- tc[i,]
  q <- Xq[s]
  basis.tri[[i]] <- my.orthobasis(q)
}
  • 疎行列を使う
library(Matrix)
M <- Matrix(0,2*n.t,2*n.v)
for(i in 1:n.t){
  s <- tc[i,]
  b <- basis.tri[[i]]
  A <- As[i]
  v32 <- Xq[s[3]]-Xq[s[2]]
  v13 <- Xq[s[1]]-Xq[s[3]]
  v21 <- Xq[s[2]]-Xq[s[1]]
  
  tmp1 <- Re(v32*b[1])
  tmp2 <- Re(v32*b[2])
  M[(1+(i-1)*2):(i*2),(1+(s[1]-1)*2):(s[1]*2)] <- c(tmp1,tmp2,-tmp2,tmp1)/sqrt(A)
  tmp1 <- Re(v13*b[1])
  tmp2 <- Re(v13*b[2])
  M[(1+(i-1)*2):(i*2),(1+(s[2]-1)*2):(s[2]*2)] <- c(tmp1,tmp2,-tmp2,tmp1)/sqrt(A)
  tmp1 <- Re(v21*b[1])
  tmp2 <- Re(v21*b[2])
  M[(1+(i-1)*2):(i*2),(1+(s[3]-1)*2):(s[3]*2)] <- c(tmp1,tmp2,-tmp2,tmp1)/sqrt(A)
}
col.fix <- sort(c(1+(v.fix-1)*2,v.fix*2))
col.free <- sort(c(1+(v.free-1)*2,v.free*2))
M.fix <- M[,col.fix]
M.free <- M[,col.free]
b.fix <- Matrix(c(t(new.X[v.fix,1:2])),ncol=1)
b.free <- solve(t(M.free)%*%M.free)%*%t(M.free)%*%(-M.fix%*%b.fix)
  • ２次元座標を作る
X.c <- new.X[,1:2]
X.c[v.free,] <- matrix(b.free,byrow=TRUE,ncol=2)
ans <- X.c
col2 <- rep(1,n.v)
col2[v.fix] <- 2
plot(ans,pch=20,col=col2)
segments(ans[lattice.edge[,1],1],ans[lattice.edge[,1],2],ans[lattice.edge[,2],1],ans[lattice.edge[,2],2])

col3 <- apply(ans%/%0.5,1,sum)%%2 +1
plot3d(new.X,col=col3,size=10,pch=15)