遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

pivot transformation

pivot transformation
  • こちらに生物界にある「形」を数学的に構成する方法としてpivot transformationと言うのがある
  • なかなか解らなかったのだけれど、こういうこと（らしい）
  • とても苦労したのだが、実は、求めていたものではなかったかも…（複比・卵・心臓あたりまでは求めていたものだったのだが…、特に胎児の形…。今の印象では、次元をもう１つ上げないといけない？？？要、検討）
  • 曲面がある
  • 曲面の点を次のルールで移すことで、曲面の形を変える
  • 曲面は微分可能な普通の曲面とすれば、曲面上のすべての点に接面Tが一つある
  • 今、空間にある直線hを取る
  • 曲面上の点xとその接面Tとをとる
  • Tとhとの交点をQとする
  • このQを射影変換する
  • 射影変換するにあたって、直線L上の２点X,Yが「不動点」となるような射影変換を考える
  • そのような射影変換は、X,Yともう１つh上にない点Zとが作る三角形を基本として行う
  • この三角形X,Y,Zにさらに射影変換のためのプレーヤーを追加する
  • 直線XZ上に２点L,Mを取る
  • 三角形X,Y,Zを乗せた平面上にあって、Yを通る直線vを取る
  • 射影変換は、X,Y,Zをの載せた平面上で行う
  • 直線QLとvとの交点をQ2とする
  • 直線MQ2とhとの交点をPとする
  • さらにPを通りhに垂直な面Sを取る
  • Sと直線xQとの交点をRとし、このRをpivot transformation の移った先とする…（らしい）

my.cross.line <- function(X,Y,Z,W){
	tmp1 <- (X[2]-Y[2])*(Y[1]-W[1])-(X[1]-Y[1])*(Y[2]-W[2])
	tmp2 <- (X[2]-Y[2])*(Z[1]-W[1])-(X[1]-Y[1])*(Z[2]-W[2])
	if(tmp2==0){
		tmp1 <- (X[3]-Y[3])*(Y[1]-W[1])-(X[1]-Y[1])*(Y[3]-W[3])
		tmp2 <- (X[3]-Y[3])*(Z[1]-W[1])-(X[1]-Y[1])*(Z[3]-W[3])

	}
	s <- tmp1/tmp2
	s*Z + (1-s) * W
}

my.cross.tangent <- function(A,B,X,v){
	t <- - sum((B-X)*v)/sum((A-B)*v)
	t*A+(1-t)*B

}
my.pivot.transform <- function(x,y,A=c(0,0,0),B=c(0,0,-1),C = c(-5,0,-5),r1 = 0.5,r2 =  0.3){
	L <- r1 * A + (1-r1) * C
	M <- r2 * A + (1-r2) * C
	Q1 <- my.cross.tangent(A,B,x,y)
	Q2 <- my.cross.line(B,C,L,Q1)
	P <- my.cross.line(A,B,M,Q2)
	R <- my.cross.tangent(x,Q1,P,A-B)
	return(list(L=L,M=M,Q1=Q1,Q2=Q2,P=P,R=R))
}

x <- c(3,0,-1)
y <- c(1,1,1)
A=c(0,0,0);B=c(0,0,-1);C = c(-5,0,-5)
out <- my.pivot.transform(x,y)
library(rgl)
ppts <- rbind(A,B,C,out$L,out$M,out$Q1,out$Q2,out$P,out$P,out$R,x,x+y)
ppts2 <- rbind(ppts,rep(min(ppts),3))
ppts2 <- rbind(ppts2,rep(max(ppts),3))
plot3d(ppts2,type="s", col="red", size=1)
#points3d(ppts,color="red",size=5,lit=TRUE)
lines3d(rbind(A,B),col=2)
lines3d(rbind(B,C),col=2)
lines3d(rbind(A,C),col=2)
lines3d(rbind(A,out$Q1),col=2)
lines3d(rbind(A,out$P),col=2)

lines3d(rbind(x,out$Q1),col=4)
lines3d(rbind(out$R,x),col=4)
lines3d(rbind(out$R,out$Q1),col=4)

lines3d(rbind(out$Q1,out$Q2),col=3)
lines3d(rbind(out$L,out$Q1),col=3)
lines3d(rbind(out$Q2,out$L),col=3)

lines3d(rbind(out$Q2,out$P),col=5)
lines3d(rbind(out$M,out$Q2),col=5)
lines3d(rbind(out$M,out$P),col=5)

lines3d(rbind(out$P,out$R),col=6)
lines3d(rbind(x,x+y),col=6)

rgl.planes(a=y,d= - sum(x*y),alpha=0.5)
rgl.planes(a=A-B,d= - sum(out$P*(A-B)),alpha=0.5)


##############

t <- seq(from=-6,to=6,length=200)
k1 <- 0.5
k2 <- 0.7
theta <- 4
f <- 0.05
X <- exp(k2*t)*cos(theta*t)
Y <- exp(k2*t)*sin(theta*t)
Z <- exp(k1*t)
X. <- k2 * X - theta*Y
Y. <- k2 * Y + theta*X
Z. <- k1 * Z

Z.. <-  k2 * exp(k2*t)
X.. <- - k1 * exp(k1*t) * cos(theta*t)
Y.. <- - k1 * exp(k1*t) * sin(theta*t)

XYZ.. <- cbind(X..,Y..,Z..)

sum(XYZ..[10,]*XYZ.[10,])
rots <- seq(from=0,to=1,length=100)*2*pi
#rots <- 0
XYZall <- XYZ..all <- matrix(0,0,3)
for(i in 1:length(rots)){
	Rot <- matrix(c(cos(rots[i]),sin(rots[i]),-sin(rots[i]),cos(rots[i])),2,2)
	tmpXY <- t(Rot %*% rbind(X,Y))
	tmpXY.. <- t(Rot %*% t(XYZ..[,1:2]))
	XYZall <- rbind(XYZall,cbind(tmpXY,Z))
	XYZ..all <- rbind(XYZ..all,cbind(tmpXY..,XYZ..[,3]))
}
XYZ2 <- rbind(XYZall,rep(min(XYZall)))
XYZ2 <- rbind(XYZ2,rep(max(XYZall)))


shift <- c(0,0,0)

plot3d(XYZ2)

ss <- sample(1:length(XYZall[,1]),10)
for(i in 1:length(ss)){

rgl.planes(a=XYZ..all[ss[i],],d= - sum(XYZall[ss[i],]*XYZ..all[ss[i],]),alpha=0.5)
}
rgl.planes(a=c(0,0,1),d= 0,alpha=0.5)

ret <- list()
tr.xyz <- XYZall
phi <- pi/100

RR <- matrix(c(cos(phi),sin(phi),-sin(phi),cos(phi)),2,2)

u <- 0.5
v <- 0
w <- 0.02
A <- c(0,0,0)
B <- c(0,0,-1)
C <- c(1,0,-0.5)
r1 <- 1.2
r2 <- 1.6

A[3] <- A[3]+u
B[3] <- B[3]+u
C[3] <- C[3]+u
A[1] <- A[1]+v
B[1] <- B[1]+v
C[1] <- C[1]+v

A[2] <- A[2]+w
B[2] <- B[2]+w
C[2] <- C[2]+w

for(i in 1:length(XYZall[,1])){
	tmp <- XYZall[i,]
	tmp.. <- XYZ..all[i,]
	#tmp[c(1,3)] <- RR %*% tmp[c(1,3)]
	#tmp..[c(1,3)] <- solve(RR) %*% tmp[c(1,3)]
	ret[[i]] <- my.pivot.transform(tmp+shift,tmp..,A,B,C,r1,r2)
	tr.xyz[i,] <- ret[[i]]$R
}

rng <- apply(tr.xyz,2,quantile,c(0.01,0.99))


sel <- which(tr.xyz[,1] > rng[1,1] & tr.xyz[,1] < rng[2,1] & tr.xyz[,2] > rng[1,2] & tr.xyz[,2] < rng[2,2] & tr.xyz[,3] > rng[1,3] & tr.xyz[,3] < rng[2,3])

tr.xyz2 <- tr.xyz[sel,]
tr.xyz2 <- rbind(tr.xyz[sel,],rep(min(tr.xyz[sel,]),3))
tr.xyz2 <- rbind(tr.xyz2,rep(max(tr.xyz[sel,]),3))

plot3d(tr.xyz2)
plot3d(tr.xyz[sel,])