2010-07-31

凸包断面

R 分割表凸包多次元球

原点を中心とする多次元単位正球Qがあるとする
Qの表面上に複数の点Tを置く
Tの各点を通り、Sに接するような接面がTの個数あるが、これによって囲まれ、原点を含むような立体は、凸包である。この凸包Sを考える
相互に直交する２つの単位ベクトルA,Bがあるとき、原点を始点とする２ベクトルA,Bが作る平面Pを考える
凸包SのPでの断面は、２次元凸包である。これを、S(P)とする
S(P)を描いてみることにする
手順
- Tの１つTiとA,Bとの内積,があるとき、Tiの面Pへの射影はT'= A +Bである。
- Tiが作る接面とPとの交点は、T'を法線ベクトルとして、原点からT'方向に距離1/|T'|の点を通る直線である。この点をTi''とする
- 今、Tの２つの点とその接面と、それがPと作る交線を考え、２本の交線の交点座標を考える
- 今、P上の座標をベクトルAを第１軸,Bを第２軸とし、Aを基準とした角座標を考える
  - Ti''の角をti,長さをri=1/|T'| とし、Tj''も同様に表すことにする
  - 交線上の点Xの角座標をtx, rxとすれば、rx cos(tx-ti) = ri, rx cos(tj-tx)=rj なる関係があることを利用して、交点の各座標は求められる
- 実際、角0の方向から、１方向に1周するようにして、順次、交線が作る平面上の凸包S(P)を構成する交線がどれになるかを確認することができる。
- 以下はそのようにして描いたものである。

#Simplex2
#最短辺を選ぶ
#ついで、現在の法線ベクトルから反時計周りにpiの法線ベクトルを持つ辺を対象とし
#交点と、現在の法線ベクトルとのなす余弦が最大になるような辺を次の辺として採用する
#最短辺に戻ったら終わり

Simplex2<-function(Ts,a,r){
 A<-Ts%*%t(a)
 R<-Ts%*%r
 L<-sqrt(A^2+R^2)
 Cs<-A/L
 Ss<-R/L
 Thetas<-ThetaAngle(Cs,Ss)
 #Thetas[which(Ss<0)]<-2*pi -Thetas[which(Ss<0)]

 TPointsX<-1/L * Cs
 TPointsY<-1/L * Ss
 xlim<-ylim<-c(-10,10)
 plot(TPointsX,TPointsY,xlim=xlim,ylim=ylim)

 
 #selected<-c(which(A==max(A)))
 selected<-which(L==max(L))[1]
 StartTheta<-Thetas[selected[1]]
 Theta2<-(Thetas-StartTheta)%%(2*pi)
 #Theta2[selected[1]]<-2*pi
 current<-selected[1]
 #VsX<-c(1/max(A))
 #VsY<-c(0)
 VsX<-c()
 VsY<-c()
 loop<-TRUE
 #first<-TRUE
 while(loop){
  #if(!first){
   #Theta2[selected[1]]<-2*pi
  #}
  #if(first){
   #first<-FALSE
  #}
  currentTheta2<-Theta2[current]
  #print(current)
  plusPiTheta2<-currentTheta2+pi
  EndTheta2<-min(2*pi,plusPiTheta2)
  #print(EndTheta2)
  if(EndTheta2==(2*pi)){ # 最初の法線からpi以上回ったら、最初の法線も考慮の対象に入れるため
   Theta2[selected[1]]<-(currentTheta2+EndTheta2)/2
  }
  Mores<-which(Theta2>currentTheta2 & Theta2<EndTheta2)


  tmpL<-sqrt((1/L[Mores]*Cs[current]-1/L[current]*Cs[Mores])^2 + 
   (1/L[Mores]*Ss[current]-1/L[current]*Ss[Mores])^2 )
   
  tmpSs<-(-1/L[current]*Cs[Mores]+1/L[Mores]*Cs[current]) / tmpL
  tmpCs<-(1/L[current]*Ss[Mores]-1/L[Mores]*Ss[current]) / tmpL
  
  if(length(tmpCs)>0){
   ttcc<-tmpCs * Cs[current] + tmpSs * Ss[current]
   tmpid<-which(ttcc==max(ttcc))[1]
   #tmpid<-which(tmpThetas==min(tmpThetas))
   tmpselected<-Mores[tmpid]
   selected<-c(selected,tmpselected)
   current<-tmpselected

  }else{
   loop<-FALSE
  }
  if(current==selected[1]){
   loop<-FALSE
  }
  #print(selected)
  #plot(1/L[selected]*Cs[selected],1/L[selected]*Ss[selected])
  #for(i in 1:length(selected)){
   #abline(1/(L[selected[i]]*Ss[selected[i]]),-Cs[selected[i]]/Ss[selected[i]],col="red")
  #}
 }
 #print(Thetas[selected])

 crossX<-crossY<-rep(0,length(selected)-1)
 for(i in 1:(length(selected)-1)){
  x1<-1/L[selected[i]]*Cs[selected[i]]
  y1<-1/L[selected[i]]*Ss[selected[i]]
  #if(i<length(selected)){
   x2<-1/L[selected[i+1]]*Cs[selected[i+1]]
   y2<-1/L[selected[i+1]]*Ss[selected[i+1]]
  #}else{
   #x2<-1/L[selected[1]]*Cs[selected[1]]
   #y2<-1/L[selected[1]]*Ss[selected[1]]
  #}
  crossX[i]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2)
  crossY[i]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2)
 }

 
 list(L=L,Cs=Cs,Ss=Ss,selected=selected,VsX=crossX,VsY=crossY)
}

df<-10


Nt<-50

# 法線ベクトルの集合
Ts<-RandomSphere(n=Nt,df=df)
Ts<-rbind(Ts,-Ts)

# 真実のモデルの方向の単位ベクトル
a<-RandomSphere(n=1,df=df)

# 輪切りにする面の数
Ns<-1000
# その単位ベクトル
rs<-RandomSphere(n=Ns,df=df)
ip<-rs %*% t(a)
rs<-t(t(rs)-t(a)%*%t(ip))

rs<-rs / sqrt(apply(rs^2,1,sum))

plot(as.data.frame(rs))
# 統計量の大きさとモデルの強さとの比として計算する数(Nr*2+1)
Nr<-8

Ratios<-10^((-Nr):Nr)

# 統計量とモデルの強さの比あたり、計算する、角度の数
Na<-12
Angles<-seq(from=0,to=1,by=1/(Na-1))*2*pi


for(i in 1:length(rs[,1])){

simplexOut<-Simplex2(Ts,a,rs[i,])


#simplexOut<-Simplex(Ts,a,r)

#simplexOut
xlim<-ylim<-c(min(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss),
max(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss))

xlim<-ylim-c(min(simplexOut$VsX,simplexOut$VsY),max(simplexOut$VsX,simplexOut$VsY)) 

xlim<-ylim<-c(-10,10)
plot(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss,xlim=xlim,ylim=ylim,pch=20)
for(i in 1:length(simplexOut$selected)){
 abline(1/(simplexOut$L[simplexOut$selected[i]]*simplexOut$Ss[simplexOut$selected[i]]),-simplexOut$Cs[simplexOut$selected[i]]/simplexOut$Ss[simplexOut$selected[i]],col="red")
}

par(new=TRUE)
plot(simplexOut$VsX,simplexOut$VsY,xlim=xlim,ylim=ylim,pch=19,col="blue")

}


df<-5


Nt<-10

# 法線ベクトルの集合
Ts<-RandomSphere(n=Nt,df=df)
Ts<-rbind(Ts,-Ts)

# 真実のモデルの方向の単位ベクトル
a<-RandomSphere(n=1,df=df)

# 輪切りにする面の数
Ns<-1000
# その単位ベクトル
rs<-RandomSphere(n=Ns,df=df)
ip<-rs %*% t(a)
rs<-t(t(rs)-t(a)%*%t(ip))

rs<-rs / sqrt(apply(rs^2,1,sum))

plot(as.data.frame(rs))
# 統計量の大きさとモデルの強さとの比として計算する数(Nr*2+1)
Nr<-8

Ratios<-10^((-Nr):Nr)

# 統計量とモデルの強さの比あたり、計算する、角度の数
Na<-12
Angles<-seq(from=0,to=1,by=1/(Na-1))*2*pi

for(i in 1:Ns){
 #単位凸包を作る
 simplexOut<-Simplex(Ts,a,rs[i,])
 xlim<-ylim<-c(-10,10)
 plot(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss,xlim=xlim,ylim=ylim,pch=20)
 for(i in 1:length(simplexOut$selected)){
  abline(1/(simplexOut$L[simplexOut$selected[i]]*simplexOut$Ss[simplexOut$selected[i]]),-simplexOut$Cs[simplexOut$selected[i]]/simplexOut$Ss[simplexOut$selected[i]],col="red")
 }

 par(new=TRUE)
 plot(simplexOut$VsX,simplexOut$VsY,xlim=xlim,ylim=ylim,pch=19,col="blue")

 #相対的な位置を変えたとき(Ratiosの値)の関数として表したい
 #アウトプットは、２交点の座標なのだが・・・虚数も含めて-b \pm sqrt(D) /2a のDとかを使って表す

}