2014-11-07

円の変形

共形変換

昨日の記事で、共形変換による閉領域の単位円板対応では、外周の曲率を対応点の粗密にするとよいのでは、ということになった
さて
単位円板を単位円板に移しつつ、その外周粗密を変えるときの変換はどうするか？
外周は「ふりふり」になる様子が見える。これは細胞の変形のときの様子？

# この辺りは試行錯誤状態のコード
t <- seq(from=0,to=1,length=500) * 2 * pi
x <- cos(t)+10
y <- sin(t)+20
z <- x + 1i * y
#z <- z/max(Mod(z))
library(MCMCpack)
k <- 5
#as <- c(rdirichlet(1,rep(1,k)) * sample(c(-1,1),k,replace=TRUE))
as <- c(rdirichlet(1,rep(1,k))) 
# 複素多項式の係数を与えて共形変換する関数
my.conf.transform <- function(z,as){
	tmp <- rep(as[1],length(z))
	if(length(as)>1){
		for(i in 2:length(as)){
			tmp <- tmp + as[i] * z^(i-1)
		}
	}
	tmp
	#tmp <- tmp/max(Mod(tmp))*max(Mod(z))
}
# 単位円板の外周粗密は、外周にそって、単調増加な0-2piの列を与えることで作れる。その関数
my.f.unit.disk <- function(g){
	if(is.null(g)){
		t <- seq(from=0,to=1,length=501) * 2 * pi
		t <- t[-length(t)]
		x <- cos(t)
		y <- sin(t)
		z <- x + 1i * y
	}else{
		#g. <- sort(g)
		g. <- (g-min(g))/(max(g)-min(g))*2*pi
		g. <- g.[-length(g.)]
		x <- cos(g.)
		y <- sin(g.)
		z <- x+ 1i * y
	}
	z
}
# 外周の単調増加関数であって、連続・微分可能で、0=2piで滑らかに接続するためには、三角関数に整数倍周期を与えて、合成し、それを一様増加関数と合算する。その上で、単調増加を維持するように、三角関数部分の最大減少率を勘案すればよい
t <- seq(from=0,to=1,length=50) * 2*pi
g <- (5*sin(t) + 5*sin(t*2) )/15 + t # 5*1+5*2 = 15
plot(g)
x <- cos(t)
y <- sin(t)
z <- x + 1i * y
z <- z[-length(z)]
Z <- my.f.unit.disk(g)
Z <- Z 
plot(Z)

plot(c(z,Z+3),col=rep(1:2,each=length(z)))
# 2014/11/3の記事で、点ペアを与えて、それを十分せつめいするだけの自由度のある複素多項式を計算するのは、単なる線形代数であることを示したが、それ
my.conformal.coef <- function(zn,wn){
	n <- length(zn)
	W <- c(Re(wn),Im(wn))
	Z <- matrix(0,2*n,2*n)
	for(i in 1:n){
		tmp <- zn[i]^(0:(n-1))
		Z[i,] <- c(Re(tmp),-Im(tmp))
		Z[n+i,] <- c(Im(tmp),Re(tmp))
		Z[i,n+1] <- Z[n+i,1] <- 0
	}
	solve(Z) %*% W
	# lm(W~Z-1)

}
#zn <- z+ rnorm(length(z),0,0.1)
#wn <- Z+ rnorm(length(Z),0,0.1)
# 係数ベクトル
tmp.out <- my.conformal.coef(z,Z)
# 係数ベクトルから複素係数を作る
tmp.out.mat <- matrix(tmp.out,ncol=2)
as.est <- tmp.out.mat[,1] + 1i * tmp.out.mat[,2]
# 推定複素多項式で確かに、ＯＫであることの確認
Z. <- my.conf.transform(z,as.est)
plot(Z)
points(Z.,col=2)
# 円の内部に極座標格子を作り、その変換の様子を見る
t <- seq(from=0,to=1,length=300)*2*pi
t <- t[-1]
r <- seq(from=0,to=1,length=100)
x <- y <- c()
for(i in 1:length(r)){
	#tmp.t <- seq(from=0,to=1,length=ceiling(2*pi*r[i]*50))*2*pi
	#tmp.t <- tmp.t[-1]
	tmp.t <- t
	x <- c(x,r[i]*cos(tmp.t))
	y <- c(y,r[i]*sin(tmp.t))
}
z2 <- x + 1i * y
z2 <- c(z2,z)
plot(z2)
Z... <- my.conf.transform(z,as.est)
plot(Z...)

Z.. <- my.conf.transform(z2,as.est)
plot(Z..)
points(Z...,type="l")

par(mfcol=c(1,2))
col1 <- Mod(Z..) * length(r)
plot(Z..,type="b",cex=0.1,col=rep(1:length(t),length(r)),pch=20,asp=TRUE)
#plot(Z..,type="b",cex=0.1,col=col1,pch=20,asp=TRUE)
col2 <- (Arg(Z..)+pi)/(2*pi) * 100
plot(Z..,type="b",cex=0.1,col=rep(1:length(t),each=length(r)),pch=20,asp=TRUE)
#plot(Z..,type="b",cex=0.1,col=col2,pch=20,asp=TRUE)
par(mfcol=c(1,1))