外接中心

• こちらにも書いた通り、三角化した後、その外接中心を取ることは便利
• 外接中心をさらっと線形代数処理で出すには、こちらのpydecのそれがよい
• Barycentric coordinateで出してそのあと通常座標にすればよい
• それをRで書き換えると：

my.circ.center.bary <- function(pts){
	d <- length(pts[1,])
	n.pt <- length(pts[,1])
	A <- 2* pts %*% t(pts)
	A <- cbind(A,rep(1,n.pt))
	A <- rbind(A,c(rep(1,n.pt),0))
	b <- c(apply(pts^2,1,sum),1)
	solve(A,b)[-(n.pt+1)]
}

pts <- matrix(c(0,0,4,0),byrow=TRUE,ncol=2)
my.circ.center.bary(pts)

my.circ.center2 <- function(pts){
	tmp <- my.circ.center.bary(pts)
	apply(pts * tmp,2,sum) 
}

my.circ.center2(pts)
pts <- rbind(v1,v2,v3)
my.circum.center(v1,v2,v3) # この関数は後述
kkk <- my.circ.center2(pts)

library(geometry)
library(onion)
# ついでに３次元でhttps://ja.wikipedia.org/wiki/%E5%A4%96%E6%8E%A5%E5%86%86　にあるような出し方をする
my.circum.center <- function(v1,v2,v3){
	e1 <- v2-v1
	e2 <- v3-v2
	e3 <- v1-v3
	S <- sqrt(sum(extprod3d(e1,e2)^2))/2
	
	A <- sum(e1^2)
	B <- sum(e2^2)
	C <- sum(e3^2)
	
	(A*(B+C-A)*v3 + B*(C+A-B)*v1 + C*(A+B-C)*v2)/(16*S^2)
}

v1 <- runif(3)
v2 <- runif(3)
v3 <- runif(3)
circm <- my.circum.center(v1,v2,v3)
# 三角形頂点が外接中心から等距離にあることを確認
sum((circm-v1)^2)
sum((circm-v2)^2)
sum((circm-v3)^2)

# ついでに、それを使って三角メッシュのLaplacianのつくりが、三角形のcotangent経由で出すか、外接中心を使ったDEC経由で出すかの比較をしておく
my.check.cot <- function(v1,v2,v3,v4){
	alpha <- acos(sum((v1-v2)*(v3-v2))/sqrt(sum((v1-v2)^2))/sqrt(sum((v3-v2)^2)))
	beta <- acos(sum((v1-v4)*(v3-v4))/sqrt(sum((v1-v4)^2))/sqrt(sum((v3-v4)^2)))
	ret1 <- (1/tan(alpha) + 1/tan(beta))/2
	print(1/tan(alpha))
	print(1/tan(beta))
	circm1 <- my.circum.center(v1,v2,v3)
	circm2 <- my.circum.center(v1,v3,v4)
	
	L1 <- sqrt(sum((circm1-(v3+v1)/2)^2))*sign(1/tan(alpha))
	L2 <- sqrt(sum((circm2-(v3+v1)/2)^2))*sign(1/tan(beta))
	ret2 <- (L1 + L2)/sqrt(sum((v3-v1)^2))
	#ret2 <- sqrt(sum((circm1-circm2)^2))
	return(list(ret1,ret2))
}

n.iter <- 1000
diffs <- rep(NA,n.iter)
for(i in 1:n.iter){
v1 <- runif(3)
v2 <- runif(3)
v3 <- runif(3)
v4 <- runif(3)

out <- my.check.cot(v1,v2,v3,v4)
diffs[i] <- out[[1]]-out[[2]]
}