2015-03-14

外接円とその高次元一般化〜DECのために

三角形の外接円は３辺の垂直二等分線の交点で、もちろん、その中心から３点への距離は等しい
これを三角形の一般化である単体に一般化したものが、circumcenter
このcircumcenterはDEC(離散微分幾何)において、便利
- たとえば、三角形分割・単体分割を細かくしていくときに、そのcircumcenterに新たな点を発生させ、分割を細かくすると「良い感じ」の細分が得られる。circumcentric subdivisionと呼ぶ
- また、circumcenterを置くことは、双対を考えるときに、よい。たとえば、三角形を考えているときに、３頂点・３辺・１三角形と、その双対としての３多角形・３線分・１点が取れる。今、三角形にはその面積が「値」として与えられているときに、双対としての３多角形の面積は、三角形の面積の分割となっており、その計算には、「底辺」ｘ「高さ」が必要となるが、このときに、三角形の外接円中心を取っておくと、それぞれのアイテムの量(線分ならば長さ)を考えて、その積を考えることで、ユークリッド空間の測度計算にマッチしているから。法線ベクトルもすぐに計算できる。
- これにより、単体を構成するサブ単体のすべてに、「その全体の体積」のうちの、「この部分が支配体積 support volume」が決まる
さて、そのcircumcenterであるが、高次元単体のそれは、低次元単体のcircumcenterから、再帰的に計算すればよい
まず、２頂点、１辺である、1単体のcircumcenterは辺の中点
ついで、k-1単体のcircumcenterがわかっているときにそれらをファセットとするk単体のcircumcenterは：
- k+1個のcircumcenterから、k+1個のファセットの法線方向に伸びたところにある。これは、すべてのk+1の延長線の交点になるが、直線の交点は２直線が決まればよいので、全部のファセットで実施する必要はない
やってみよう

my.circumcenter <- function(X.,rot=FALSE){
# 頂点座標に完全一致のものがあったりすると、逆行列計算でエラーが出るので、ランダムな回転行列 R を使ってそのような位置からずれるようにするべく、
	if(rot){
		library(GPArotation)
		d <- ncol(X.)
		R <- Random.Start(d)
		X <- t(R %*% t(X.))
	}else{
		X <- X.
	}
	n.r <- nrow(X)
	n.c <- ncol(X)
# ２点の場合
	if(n.r==2){
		ret <- apply(X,2,mean)
# それ以外
	}else{
		tmp <- matrix(0,n.r,n.c)
		tmpN <- matrix(0,n.r,n.c)
		#for(i in 1:n.r){# ２ファセットの外接円中心だけで事足りる
		for(i in 1:2){
# 再帰的にファセット外接円中心を求める
			tmp.mat <- matrix(X[-i,],ncol=n.c)
			tmp[i,] <- my.circumcenter(tmp.mat,rot=rot)
# 考慮しているファセットの１点を基準点とする
			v <- tmp.mat[1,]
# ファセット外の点を基準点からの相対位置ベクトルとする
			q <- X[i,] - v
# ファセット内の基準点以外の点の相対位置ベクトルを列ベクトルとする行列
			U <- matrix(t(matrix(tmp.mat[-1,],ncol=n.c)) - v,nrow=n.c)
# ファセット外の点からの垂線の足をU %*% p + vとすると、以下が成り立つ
			p <- solve(t(U) %*% U) %*% t(U) %*% q
# これが垂線の方向ベクトル
			tmpN[i,] <- q - (U%*%p)
		}
# ２直線は、それぞれのcircumcenter + 係数ｘ垂線ベクトルなので、２直線の係数tを求める
		t <- solve(matrix(c(tmpN[1,1],tmpN[1,2],-tmpN[2,1],-tmpN[2,2]),2,2)) %*% ((-1)*c(tmp[1,1:2]-tmp[2,1:2]))
		ret <- tmp[1,] + t[1] * tmpN[1,]
	}

	if(rot){# 逆回転しておく
		return(t(R) %*% ret)
	}else{
		return(ret)
	}
}

使ってみる
- ２次元空間。三角形。原点を中心とした単位円周上に点があるときは、circumcenterは原点だから、それを利用して確かにそうなるかどうかを検算

k <- 2 # 次元
X <- matrix(rnorm(k*(k+1)),ncol=k)
X <- X/sqrt(apply(X^2,1,sum))
# 中心をずらす
X <- t(t(X) + c(1,3))
X
my.circumcenter(X) # 確かにc(1,3)が中心

次元を上げて、同じことをする

# 次元を上げて
k <- 7
X <- matrix(rnorm(k*(k+1)),ncol=k)
X <- X/sqrt(apply(X^2,1,sum))
X. <- t(t(X) + 1:k)
my.circumcenter(X.)

回転をする、というオプションで回すかそうでないか

# 回転しないとエラーになるような例
X <- matrix(c(0,0,4,0,3,2),byrow=TRUE,ncol=2)
X <- matrix(c(-runif(3),1,0,1,0,1,0,0,1,1),byrow=TRUE,ncol=3)
my.circumcenter(X)
my.circumcenter(X,rot=TRUE)

> my.circumcenter(X)
 以下にエラー solve.default(matrix(c(tmpN[1, 1], tmpN[1, 2], -tmpN[2, 1], -tmpN[2,  : 
  Lapack routine dgesv: システムは正確に特異です: U[2,2] = 0 
> my.circumcenter(X,rot=TRUE)
            [,1]
[1,] -0.07023866
[2,] -0.07023866
[3,]  0.50000000

球面上乱点から表面三角形メッシュを作り、その上にcircumcenterもプロットしてみる

library(rgl)
library(geometry)
k <- 3
Z <- matrix(rnorm(k*100),ncol=k)
Z <- Z/sqrt(apply(Z^2,1,sum))
tZ <- delaunayn(Z)
sZ = t(surf.tri(Z, delaunayn(Z)))
rgl.triangles(Z[sZ,1], Z[sZ,2], Z[sZ,3],col="gray")
points3d(Z,col=2,size=10)

sZ

X.circ.ctr <- matrix(0,length(sZ[1,]),3)
sZ2 <- sZ
sZ2 <- matrix(0,3,0)

for(i in 1:length(X.circ.ctr[,1])){
	X.circ.ctr[i,] <- my.circumcenter(Z[sZ[,i],])
}

Z.2 <- rbind(Z,X.circ.ctr)
nn <- length(Z[1,])
for(i in 1:length(sZ[1,])){
	tmp <- cbind(c(nn+i,sZ[1,i],sZ[2,i]),c(nn+i,sZ[1,i],sZ[3,i]),c(nn+i,sZ[2,i],sZ[3,i]))
	sZ2 <- cbind(sZ2,tmp)
}
rgl.triangles(Z.2[sZ2,1], Z.2[sZ2,2], Z.2[sZ2,3],col="gray")
points3d(Z,col=2,size=10)

points3d(Z.2,col=3,size=10)