お試しメモ

• 空間にグラフがあるとする
• 空間に点があるとする
• 点からグラフへの最短距離を描く
• こちらで作った、「凸包」への垂線の足を利用する
• ループを回しているのでちょっと重いかも

• 垂線の足がどの辺に乗っているかで色分けすると、空間が塗り分けて見えやすそう

# Space
d<-2 # dimension
# Graph
n.node<-4 # number of nodes
# 適当に0,1の対称行列(グラフの連結行列)を作る
m<-matrix(runif(n.node^2),n.node,n.node)
m<-m+t(m)
mean.m<-mean(m)
m<-(m>mean.m)
diag(m)<-0
m
# coodinates of nodes
x<-matrix(runif(n.node*d),ncol=d)
x



# coordinates of points in the space
n.pt<-1000
a<-matrix(runif(n.pt*d),ncol=d)
# 垂線の足を計算する
# グラフへの垂線の足は複数である可能性がある
# 下の例では、グラフのノードが何度も選ばれている例しかないが
# 厳密には、複数の辺にたまたま同長の垂線の足が取れることがある
#feet<-t(apply(a,1,calc.foot,x=x,m=m))
feet<-apply(a,1,calc.foot,x=x,m=m)

# グラフを赤で、垂線を緑で描く
plot(x,col=2,xlim=c(0,1),ylim=c(0,1))
	m2<-m
	m2[lower.tri(m2)] <- 0
	edges<-which(m2==1,arr.ind=TRUE)
for(i in 1:length(edges[,1])){
	segments(x[edges[i,1],1],x[edges[i,1],2],x[edges[i,2],1],x[edges[i,2],2],col=2)
}
for(i in 1:n.pt){
	#for(j in 1:length(feet[[i]][,1])){
	for(j in 1:length(feet[[i]]$edge)){
		#segments(a[i,1],a[i,2],feet[[i]][i,1],feet[i,2],col=3)
		# sum(m2)がエッジの本数なので、それごとに色を変える
		segments(a[i,1],a[i,2],feet[[i]]$coods[j,1],feet[[i]]$coods[j,2],col=rainbow(sum(m2))[feet[[i]]$edge[j]])
		#print(feet[[i]]$edge[j])
		#print(rainbow(100)[feet[[i]]$edge[j]])
	}
	
}

###############
# 使用関数
# ベクトルのノルム
nrm<-function(A){
	sqrt(sum(A^2))
}
# AからXsが作る「hyperplane」への垂線の足を計算する
# ２点が作る「hyperplane」は直線
# 垂線の足がXsが作る凸包(２点の場合は線分)の上にあるかどうかを判定する
MinDistance2<-function(A,Xs,bothside=TRUE){
 df<-length(A)
 N<-length(Xs[,1])
 Xs1<-t(t(Xs)-A)
 C<-apply(Xs1,2,sum)/N
 Xs2<-t(t(Xs1)-C)
# \sum_{i=1}^N a_i v_i と表したときに \sum_{i=1}^N a_i=1という制約がある
# v_iは与えられている凸包の頂点座標(ただし、重心を原点にとりなおしている
# 垂線の足をbとすると
# Ma=bという連立方程式を解いてaを求めるたい
# ただし凸包はN次元ではなくてN-1次元なので、N次元のうち１次元分は取り除
いて
# その代わり係数の和が１という制約等式を入れれば
# あとは行列を使って連立方程式を解くだけなので、Rに任せる

 Xs3<-Xs2 %*% t(Xs1)
 Xs4<-rbind(Xs3[1:(N-1),],rep(1,N))
 coefs<-solve(Xs4,c(rep(0,(N-1)),1))
 Apr<-t(Xs1) %*% coefs
 L<-nrm(Apr)
 Apr<-Apr+A
 inside<-FALSE
 if(length(coefs[coefs<0])==0 || (bothside & length(coefs[coefs>0])==0) ){# すべての要素が0以上の条件
  inside<-TRUE
 }
 list(A=A,Xs=Xs,distance=L,Apr=Apr,coefs=coefs,inside=inside)

}

# function to get foot of points to the graph
# 点から線分への足を計算するには、足が線分上にあるかその外側にあるかの判定が必要
calc.foot<-function(X,x,m){
	difference<-t(x)-X
	distance<-sqrt(apply(difference^2,2,sum))
	m2<-m
	m2[lower.tri(m2)] <- 0
	edges<-which(m2==1,arr.ind=TRUE)
	min.dist.per.edge<-rep(0,length(edges[,1]))
	t.per.edge<-rep(0,length(edges[,1]))
	closest<-matrix(0,length(edges[,1]),length(X))
	for(i in 1:length(edges[,1])){
		tmp.out<-MinDistance2(X,x[edges[i,],])
		if(tmp.out$inside){
			t.per.edge[i]<-tmp.out$coefs[1]
			min.dist.per.edge[i]<-tmp.out$distance
			closest[i,]<-tmp.out$Apr
		}else{
			if(distance[edges[i,1]]<distance[edges[i,2]]){
				t.per.edge[i]<-1
				min.dist.per.edge[i]<-distance[edges[i,1]]
				closest[i,]<-x[edges[i,1],]
			}else{
				min.dist.per.edge[i]<-distance[edges[i,2]]
				closest[i,]<-x[edges[i,2],]
			}
		}
	}
	min.pt<-which(min.dist.per.edge==min(min.dist.per.edge))
	if(length(min.pt)==1){
		return(list(coods=matrix(closest[min.pt,],nrow=1),edge=min.pt))
	}else{
		return(list(coods=closest[min.pt,],edge=min.pt))
	}
}