複素数版・四元数版〜Least Square Conformal Mapping
曲面への共形変換マッピング?複素数・四元数のそれぞれで最小二乗法?
- 作者: ryamada
- 発売日: 2015/03/21
- メディア: Kindle版
- この商品を含むブログ (1件) を見る
--- title: "４元数を用いた曲面の共形変換マッピング" author: "ryamada" date: "Wednesday, March 18, 2015" output: html_document --- この文書は、平面の格子を曲面に張り付ける話である。 # ２次元平面と共形変換 まずは、平面に格子を作って描く。 ```{r,echo=FALSE} library(knitr) ``` ```{r} # xy二次元平面のグリッド点を作る x1 <- y1 <- seq(from=-10,to=10,length=31) X <- as.matrix(expand.grid(x1,y1)) # 点の数 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$という関数にしてみよう。 ```{r} p <- 12+1i*7;k <- 0.4 w <- X[,1] + 1i*X[,2] new.w <- (w+p)^k 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]) ``` この変形の特徴は、各所で角度が保存されていることであり、共形変換を角度保存の変換と呼ぶ。 # ３次元空間に埋め込まれた曲面 上の例では、２次元平面内で複素数を用いて格子をぐにゃりと曲げたが、二次元平面内に収まっている。 これに３次元座標も与えてやると、３次元空間内の曲面ができる。 ```{r,echo=FALSE} library(rgl) ``` ```{r setup,echo=FALSE} knit_hooks$set(rgl = hook_rgl) ``` ```{r} p <- 12+1i*7;k <- 0.4 w <- X[,1] + 1i*X[,2] new.w <- (w+p)^k new.z <- 0.4*exp(-0.1*Mod(w)^2) new.X <- cbind(Re(new.w),Im(new.w),new.z) #new.X[v.free,] <- new.X[v.free,] + rnorm(length(new.X[v.free,]),0,0.005) ``` ```{r,rgl=TRUE} plot3d(new.X,col=col) segments3d(new.X[c(t(lattice.edge)),]) ``` # 三角化 格子は座標系の変形を視覚化するのに有用だが、局面データの離散的な取扱いのためには、三角形メッシュにするのが好都合である。 Rではgeometryパッケージを使うと三角化が可能である。 ```{r} 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]) ``` ```{r,rgl=TRUE} plot3d(new.X,col=col) triangles3d(new.X[c(t(tc)),],col=2) ``` # 三角形の面積とベクトルのクロス積と面の向き 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)$で考える場合とで、クロス積の値は符号が反転する。 ```{r} 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次元ベクトルを表し、同様に符号付き面積を計算することもできる。 ４元数に慣れる意味で、その計算も行っておく。 ```{r} 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) ``` すべての三角形について符号付き面積を計算してみる。 ```{r} 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) } my.area.tri.q <- function(h){ h12 <- h[1]-h[2] h13 <- h[1]-h[3] tmp <- h12*h13 i(tmp)+j(tmp)+k(tmp) } ``` 符号付き面積の符号で三角形を塗り分けてみると、以下のようになる。 ```{r,rgl=TRUE} 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=rep(col,each=3)) ``` すべての三角形の向きをそろえたいので、万一、面積の符合があっていなかった場合には、その値をもとに、裏向きの三角形の頂点順序を変えることにする。 ```{r,rgl=TRUE} 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=rep(col,each=3)) ``` # 最小二乗法による法線ベクトル推定と共形変換推定 その１ （２ｘ面の数）ｘ（２ｘ点の数）行列を作る。 面と点とのペアごとに２ｘ２行列を取り、それを配置したもの。 この２ｘ２行列を説明する。 今、ある三角形の面積の２倍を$A$とし、 その三角形の頂点の一つをvとする。 vを除く２頂点を結ぶ辺をeとする。 ここで、辺eを、この三角形を乗せた平面上の適当な座標系で$(e_x,e_y)$と表すことにする。 このとき、この三角形と頂点とのペアのための２ｘ２行列を、$ \begin{pmatrix} e_x & -e_y \\ e_y & e_x \end{pmatrix}/\sqrt{A}$$ とする。 $(e_x,e_y)$を算出するためには、三角形を乗せた座標系を計算する必要があるから、それをする関数を作る。 ```{r} # ３頂点座標に相当する、４元数のベクトルを引数とする 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 # 正規直交基底であることの確認 ``` ```{r} # 頂点座標の四元数表現 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) } ``` 基底が得られたので、三角形ごとに、辺を基底で２次元座標表現し、それを、（２ｘ三角形の数）ｘ（２ｘ点の数）の行列に納める。 ```{r} 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) } ``` 今、このようにして作った行列Mを列について二つに分割する。 Mの列は、点に対応して定まっているが、今、点を２群に分けることとする。 点のうち、曲面の周辺の点を固定点、それ以外の点を非固定点とする。 固定点については、曲面上にはりつける二次元座標を定めて与えることとする。 非固定点については、推定の対象とする。 このときMのうち、固定点に対応する列だけを抜き出したものを$M_k$、非固定点のそれを$M_u$とする。 ここで、固定点について与える、二次元座標を、タンデムに並べたベクトルを$b_k$とし、それに対応する、非固定点のベクトル（値は未知）を$b_u$とする。 このとき、曲面に張り付けた２次元座標が共形変換的である、ということは$|(M_k,M_u) (b_k,b_u)^T|^2=0$であることが知られている。 そして、それを完全に満足するように、未知の$b_u$を得ることはできないので、推定することにするのだが、$|(M_k,M_u) (b_k,b_u)^T|^2$がなるべく小さくなるようにするという推定法をとることにすれば、それは最小二乗法である。 今、$|(M_k,M_u) (b_k,b_u)^T|^2=0$と$M_u b_u = - M_k b_k$とは同じことなので、結局、既知の$M_u,M_k,b_k$を用いて、未知の$b_u$を求める問題で、これは、$b_u = (M_u^T M_u)^{-1}M_u^T (-M_k b_k)$である。 以下では、このようにして$b_u$を求めている。そして$b_k$が固定点の二次元座標をベクトル化したものであったのと同様に$b_u$は非固定点の２次元座標の推定値であるので、それを曲面に割り付けて図示したものである。 ```{r} # v.fixは固定点の番号。予め、周辺の点のIDを登録しておいたもの # v.freeは非固定点の番号 # Mは点の数ｘ２列あるので、固定点・非固定点に対応する列番号を書き出す col.fix <- sort(c(1+(v.fix-1)*2,v.fix*2)) col.free <- sort(c(1+(v.free-1)*2,v.free*2)) # 2つの行列に分ける M.fix <- M[,col.fix] M.free <- M[,col.free] # 固定点の曲面３次元座標のうち、xy平面二次元座標を、固定座標をしている 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) ``` このようにして推定された非固定点の曲面上２次元座標と、固定点のそれとを合わせて、市松模様で描く。 三角形の色は、３頂点の座標の色の多数決で決めてあるので、色境界線は凸凹するが、三角形の３頂点座標から、三角形の任意の点の座標を線形補間すれば、より丁寧な塗り分けは可能。 ```{r,rgl=TRUE} X.c <- new.X[,1:2] X.c[v.free,] <- matrix(b.free,byrow=TRUE,ncol=2) ans <- X.c col3 <- apply(ans%/%0.5,1,sum)%%2 +1 plot3d(new.X,col=col3) plot3d(new.X) col4 <- apply(matrix(col3[tc],ncol=3),1,median) triangles3d(new.X[c(t(tc)),],col=rep(col4,each=3)) ans.1 <- ans col4.1 <- col4 ``` この方法は、曲面に直交座標系を張り付けている方法であり、直交座標系は、複素平面と見ることもできるから、推定している曲面座標は複素座標とも言える。 # 最小二乗法による法線ベクトル推定と共形変換推定 その２ その１で示した方法は、複素数の代りに４元数を用いる方法である。 こちらでは、推定にあたって、より多くの情報（具体的には法線ベクトルの方向と、その方向の長さ（曲率に相当））を推定しつつ、それを使って二次元座標化する方法である。 原理や証明はすべて飛ばして、手続きをRで実装することで、まず概要をつかむこととする。 手続きは以下の通り。 (1) ２つの行列F,Dを作る。 (2) 周辺点に制約条件を与える。この制約は、１つの４元数である。４元数の実部は、接平面における伸縮の程度を表し、虚部は法線ベクトルであって、法線方向とその方向への伸び縮み具合を表す。 (3) それらを用いて、4元数を与えなかった点に4元数を推定する。その際、最小二乗法を用いる。実部も虚部も隣接する点同士で滑らかに変化するように推定することになる。したがって、接平面での伸び縮みと法線方向とその方向の伸び縮み（曲率）が滑らかに変化するように曲面を推定する (4) 推定した４元数から、法線ベクトルと、曲面に沿った共形変換後座標を算出する。曲面に沿った２次元座標は、接面方向の伸び縮みと法線方向の伸び縮みとする。推定した情報のうち、曲面上の座標に使用しなかった情報は、法線ベクトルの方向成分だけであるが、これは、「曲面について考えるとき、変化しない要素」であるから、結果として、推定したもののうち、「曲面に沿って変化する要素のすべて」を用いて２次元座標化しものとなっている ## 行列Fの作成 行列Fは次のような行列である。 三角形を乗せる平面を張る正規直行基底に対応する四元数成分を対角成分とする行列であって、行数は三角形の数の2倍、列数は三角形の数の４倍である。 三角形の3頂点の3次元座標が得られているとき、座標を4元数化すると、次の要領で、簡単に法線ベクトルと、三角形平面を張る正規直行基底を得ることができる。 三角形の２辺に対応するベクトルの外積が法線ベクトル方向であることと、いったん法線ベクトルが決まれば、それと三角形の１辺とに直交するベクトルが同じ手続きでとれる。 ３頂点座標を３個の４元数で与え、正規直行基底を返す関数が以下の通り。 第１，２列が三角形平面を張り、第３列が法線ベクトルである。 ```{r} library(Matrix) make.F <- function(Xq,tc){ n.t <- length(tc[,1]) F <- Matrix(0,2*n.t,4*n.t) for(i in 1:n.t){ tmp.row <- (1+(i-1)*2):(i*2) tmp.col <- (1+(i-1)*4):(i*4) tmp.mat <- as.matrix(my.orthobasis(Xq[tc[i,]])) F[tmp.row,tmp.col] <- t(tmp.mat[,1:2]) } F } F <- make.F(Xq,tc) ``` ## 行列Dの作成 行列Dは次のような行列である。 行数は三角形の数の４倍、列数は点の数の４倍である。 4倍であることを除けば、 三角形の数ｘ点の数の組である。 この三角形と点とのペアごとに４ｘ４行列を定め、それを単位ブロックとして並べたものである。 ### 四元数とその４ｘ４実行列表示 ここで、四元数の４ｘ４行列表示について確認する。 四元数は１つの実部と３つの虚部からなる４要素の数であって、四要素の間で加減ができるし、乗算もできるが、積の順序を入れ替えると値が変わる。 この演算の特徴を持つのが４ｘ４行列なので、１つの四元数は１つの４ｘ４行列とみなすことができる。 つまり、行列Dは三角形の数ｘ点の数の行列であって、個々の要素が四元数なので、実行列で表すために、個々の要素を４ｘ４行列化したものになっている。 ### ４ｘ４行列の構成と並べ方 三角形と点とのペアは、点が三角形の構成要素である場合とそうでない場合とに分けられる。 構成要素でない場合には、４ｘ４行列は０行列である。 構成要素である場合には、その三角形のその点における対向辺ベクトル(その向きは、面の向きに注意して並べた頂点の順序によるものとする)と、三角形の面積によって定まる４ｘ４行列である。 対向辺ベクトルの３成分によって４ｘ４ベクトルを定め、それを三角形の面積の2倍の平方根で割ったものを、その行列とする。 ３要素からの４ｘ４行列作成は次の通り。 を虚部とする4元数(実部は0)を考え、4元数を実線形代数で計算するときに、用いられる行列となる。 具体的にはベクトル$(x1,x2,x3)$に対して $$ \begin{pmatrix} 0 & -x1 & -x2 & -x3 \\ x1 & 0 & -x3 & x2 \\ x2 & x3 & 0 & -x1 \\ x3 & -x2 & x1 & 0 \end{pmatrix} $$ となる。 以下は、4元数から、４ｘ４行列を作る関数である。 ```{r} Q.mat <- function(h){ matrix(c(h[1],h[2],h[3],h[4],-h[2],h[1],h[4],-h[3],-h[3],-h[4],h[1],h[2],-h[4],h[3],-h[2],h[1]),4,4) } h <- c(0,1,2,3) h.mat <- Q.mat(h) h.mat ``` 三角形の面積はすでに求めた通りである。 ### 行列Dの作成関数 ```{r} # new.X 頂点座標、Xq その四元数表現 # tc 向きのそろった三角形頂点IDのリスト # As 各三角形の面積の２倍 # basis.tri 各三角形を接平面とする正規直交基底 make.D <- function(Xq,tc,As){ n.t <- length(tc[,1]) n.v <- length(Xq) D <- Matrix(0,4*n.t,4*n.v) for(i in 1:n.t){ s <- tc[i,] v32 <- Xq[s[3]]-Xq[s[2]] v13 <- Xq[s[1]]-Xq[s[3]] v21 <- Xq[s[2]]-Xq[s[1]] D[(1+(i-1)*4):(i*4),(1+(s[1]-1)*4):(s[1]*4)] <- Q.mat(as.vector(v32))/sqrt(As[i]) D[(1+(i-1)*4):(i*4),(1+(s[2]-1)*4):(s[2]*4)] <- Q.mat(as.vector(v13))/sqrt(As[i]) D[(1+(i-1)*4):(i*4),(1+(s[3]-1)*4):(s[3]*4)] <- Q.mat(as.vector(v21))/sqrt(As[i]) } D } ``` ## 行列FD M=FDなる行列を考える。 この行列のサイズは、Fが(２ｘ三角形数)ｘ(４ｘ三角形数)の行列であり、 Dが(４ｘ三角形数）ｘ(４ｘ点の数)の行列であることから、 Mは(２ｘ三角形数)ｘ(４ｘ点の数)の行列である。 Fは三角形を張る直交する単位ベクトルの４元数からなる。 他方、Dは三角形の辺ベクトルの４元数表現を三角形の面積で補正したものである。 面積補正はひとまず無視して考える。 それらの積であるFDは、三角形と点のペアごとに定まる、２つずつの四元数からなる。 この２つの四元数の実部は、三角形平面上での座標値に相当する。 他方、虚部の３次元ベクトルは、法線方向のベクトルである（三角形の辺ベクトルは三角形に乗っているので、三角形上の２ベクトルのクロス積であるから、法線方向ベクトルになる）。 つまりFDは、各点について、その点が所属する三角形に関して、三角形の面積と対向辺の長さで重みづけをした上で、接平面成分の２直交方向と、平面と直交する方向の成分（を２つに分けて）、取り出したものとなっている。 この接平面成分はすでに実数なので、そのままに評価し、法線方向成分はベクトルなので、法線方向ベクトルとの内積をとることで実数評価してやることにしたい。 そのために、点に対して、実１成分、虚３成分の計４成分（一つの四元数）を対応付ける。 ここで調整するべきことの一つが、各点の法線ベクトルの向きである。 なぜなら、点は複数の三角形に帰属するが、帰属三角形の法線方向が異なるので、どの方向を向くのが最適かの調整をしなくてはならないからである。 その上で、接平面方向と法線方向との伸び縮みを曲面全体で滑らかにする必要がある。 以下では、この点ごとの四元数のうち、一部の点の４元数は与え、残りの点の四元数を推定する。 ## 制約条件 たとえばの条件として、曲面の周囲の点に次のような制約をつける。 以下の例では、曲面の周辺の法線方向が$(0,0,1)$として、考えている。 ３次元プロットで描いているところの$(x,y)$座標を用いて、$(x,0,0,y)$なる長さ4のベクトルを割り付ける。ここで、$y$から(0,0,y)$を作っている点が、法線方向が$(0,0,1)$であるという仮定に想定する。 非周辺の点については、長さ4のベクトルが未知であるとする。 ## 最小二乗法による推定 制約条件の説明で、点の数ｘ４の値を、一部は与え、残りは未知であるとした。 その全体をベクトルVとし、$M-V$ができるだけ0ベクトルに近くなるように未知の値を推定したい。 0ベクトルとの近さを最小二乗法で測ることにすれば、線形代数計算で求めることができる。 具体的には、Mを、制約する点に対応する列と、そうでない列とに分け、ベクトルVも同様に分けることにし、それぞれ$FD = (M_{unknown},M_{known})$、$V=(V_{unknown},V_{known})$とする。 このとき$|M_{unknown}V_{unknown}-M_{known}D_{known}|^2$の最小化問題が、Vの未知部分の最小二乗法的な解となる。 その解は $$ V_{u} = ((M_{u})^T M_{u})^{-1} (M_{u})^T (-M_{k}V_{k}) $$ これを算出する計算する。 ```{r} # v.fix and v.free are known and unknown vs, respectively D <- make.D(Xq,tc,As) M <- F %*% D # 固定点、非固定点の点IDより、対応する列IDを作成 col.free <- sort(c(1+4*(v.free-1),2+4*(v.free-1),3+4*(v.free-1),4*v.free)) col.fix <- sort(c(1+4*(v.fix-1),2+4*(v.fix-1),3+4*(v.fix-1),4*v.fix)) # 行列Mを固定・非固定で二分する M.u <- M[,col.free] M.k <- M[,col.fix] # 固定点の座標から４ベクトルを作り、それをタンデムにつないだベクトルを作る V.k <- c(rbind(new.X[v.fix,1],rep(0,length(v.fix)),rep(0,length(v.fix)),-new.X[v.fix,2])) # 非固定点の長さ４のベクトルの最小二乗解を出す V.u <- solve(t(M.u)%*%M.u)%*%t(M.u) %*% (-M.k %*% V.k) ``` ## 結果の解釈 推定の結果、すべての点は４つの値を持つこととなった。 ４つの値のうちの1つ目と残りの３つに分ける。 三つ組みの値を3次元ベクトルとみなし、その長さをとることで、 $ (a,\sqrt(b^2+c^2+d^2))$ という二次元座標が得られる。 この二次元座標を、3次元空間に埋め込まれた曲面上の二次元座標とすると、それが「最小二乗法に基づく共形変換座標」である。 ```{r} X.c <- new.X[,1:2] # V.uを行列化する V.u.m <- Matrix(as.vector(V.u),byrow=TRUE,ncol=4) X.c[v.free,1] <- V.u.m[,1] X.c[v.free,2] <- sqrt(apply(V.u.m[,2:4]^2,1,sum))*sign(-V.u.m[,4]) col2 <- rep(1,n.v) col2[v.fix] <- 2 ``` ```{r,rgl=TRUE} col3 <- apply(X.c%/%0.5,1,sum)%%2 +1 plot3d(new.X,col=col3) col4 <- apply(matrix(col3[tc],ncol=3),1,median) triangles3d(new.X[c(t(tc)),],col=rep(col4,each=3)) col4.2 <- col4 ``` 以上の４元数処理を関数かする。 ３次元座標を納めた列数３の行列Xと、固定点ID v.fixを与える。 固定点の２次元座標をX.oriとして与える。 ```{r} my.conformal.norm <- function(X,tc,v.fix,X.ori=X[v.fix,1:2]){ n.t <- length(tc[,1]) # 三角形の個数 As <- rep(0,n.t) for(i in 1:n.t){ As[i] <- my.area.tri(X[tc[i,],]) } ura.tri <- which(As<0) tc[ura.tri,2:3] <- tc[ura.tri,c(3,2)] As <- abs(As) Xq <- X[,1]*Hi + X[,2]*Hj + X[,3]*Hk F <- make.F(Xq,tc) D <- make.D(Xq,tc,As) M <- F %*% D col.free <- sort(c(1+4*(v.free-1),2+4*(v.free-1),3+4*(v.free-1),4*v.free)) col.fix <- sort(c(1+4*(v.fix-1),2+4*(v.fix-1),3+4*(v.fix-1),4*v.fix)) M.u <- M[,col.free] M.k <- M[,col.fix] #V.k <- c(rbind(new.X[v.fix,1],rep(0,length(v.fix)),rep(0,length(v.fix)),-new.X[v.fix,2])) V.k <- c(rbind(X.ori[,1],rep(0,length(v.fix)),rep(0,length(v.fix)),-X.ori[,2])) V.u <- solve(t(M.u)%*%M.u)%*%t(M.u) %*% (-M.k %*% V.k) X.c <- X[,1:2] # V.uを行列化する V.u.m <- Matrix(as.vector(V.u),byrow=TRUE,ncol=4) X.c[v.fix,1:2] <- X.ori X.c[v.free,1] <- V.u.m[,1] X.c[v.free,2] <- -V.u.m[,4] return(X.c) } #x1 <- y1 <- seq(from=-10,to=10,length=31) Xori <- as.matrix(expand.grid(x1,y1)) X.fix <- Xori[v.fix,] out <- my.conformal.norm(new.X,tc,v.fix) #plot(out,pch=20) #segments(out[lattice.edge[,1],1],out[lattice.edge[,1],2],out[lattice.edge[,2],1],out[lattice.edge[,2],2]) ``` ４元数を用いた市松模様。 ```{r,rgl=TRUE} col3 <- apply(out%/%0.5,1,sum)%%2 +1 plot3d(new.X,col=col3,main="quaternion") col4 <- apply(matrix(col3[tc],ncol=3),1,median) triangles3d(new.X[c(t(tc)),],col=rep(col4,each=3)) ``` 複素数での市松模様と少し違う。 ```{r,rgl=TRUE} plot3d(new.X,col=col3,main="complex") triangles3d(new.X[c(t(tc)),],col=rep(col4.1,each=3)) ``` ４元数と複素数の場合でのそれぞれの推 ```{r} plot(X.c[v.free,],ans[v.free,],pch=20,cex=0.1,xlab="quaternion",ylab="complex") ``` # 形を変えてやってみる ```{r} # xy二次元平面のグリッド点を作る x1 <- y1 <- seq(from=-10,to=10,length=41) X <- as.matrix(expand.grid(x1,y1)) # 点の数 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,col=col) ``` ```{r,rgl=TRUE} X.comp <- X[,1] + 1i*X[,2] args <- Arg(X.comp) mods <- Mod(X.comp) new.mods <- abs(mods-mean(mods)) new.X <- X new.X[,1] <- X[,1] new.X[,2] <- X[,2] z <- cos(X[,1]/10)+exp(-0.1*new.mods^2) new.X <- cbind(X,z) new.X[v.free,] <- new.X[v.free,] + rnorm(length(new.X[v.free,]),0,0.01) plot3d(new.X,col=col) ``` ```{r,rgl=TRUE} library(geometry) tc <- delaunayn(X[,1:2]) out <- my.conformal.norm(new.X,tc,v.fix) col3 <- apply(out%/%3,1,sum)%%2 +3 plot3d(new.X,col=col3,main="quaternion") col4 <- apply(matrix(col3[tc],ncol=3),1,median) triangles3d(new.X[c(t(tc)),],col=rep(col4,each=3)) ```
