2014-01-11

Distance Cartogram

Distance Cartogramのことが気になった
ある標本セットについて２つの距離空間情報があったときに、片方の距離空間情報で描かれている地図をもう片方の距離空間の情報で歪ませる、という話。
たとえば、個体の物理的存在位置が一つ目の距離空間で、個体の遺伝的距離がもう一つの距離空間にあるとする
歪ませた地図表現をCartogramと言うらしい(こちら)
アルゴリズムも色々と提案されているようだけれど、ちょっと自分でも考えてみたい(アルゴリズムのいい説明サイトも簡単には見つからないようだ)
片方の距離空間で地図が描けているときに、もう片方の距離空間で「基準点」間の距離がサンプリングされたら、その基準点に動かして、その基準点の間の点は「補間」する、というのが一番簡単な話だろうか
もっとも簡単な補間は１次線形補間だろうから(複雑にしたければ、ここでスムージングを考慮してスプライン関数とかにするのだろう)、それでやってみる
２次元地図だけでは、ちょっとつまらないので、任意次元を念頭にしてやってみる
基準点をとる、というのは、三角測量だから、任意次元でやる、ということはTriangulationで、それはドロネー図を描くこと
ひとたびドロネー図が描かれれば、元の空間の点は、どれか一つの三角形(三角錐…高次元化)に属するから、その三角形を見つけてやって、三角形の頂点に関して係数を出してやればよい
以下のソースは、書き散らし。ぐねぐねした絵を作るために複素関数のカラー表示をまずやって、それを適当に歪ませる。また、歪ませた空間に適当に基準点を取って、そのドロネーを出して、線型補間してある
もしかすると、はてなブログのランキングサイト「TopHatenaR」の図はarea cartogramなのかも・・・(ブログ間の似ている程度を距離として距離行列を作り、それに基づき２次元にMDSか何かをして、その上で、各ブログの面積に「被リンク」等の重みを入れて歪ませている？？？(以下がTopHatenaRの地図の部分拡大図)

歪ませる前のグネグネ図と歪ませた後のグネグネ図を並置

採用点の動きが大きすぎて三角形で対応できなくなったら(移動した後の三角形にオーバーラップや裏返しが発生したら)どうするのか、とか、の問題を無視してある。その場合は、いきなりやらずに、位相を維持するようにソロ―っと動かしたりするのかもしれないけれど、それをどうやるのかは全くわかっていない。おそらくそのあたりがCartogramアルゴリズムの本質なのでしょう
ボロノイ図の関する関連記事はこちら
複素関数のカラー表示に関する関連記事はこちら

# 複素関数の表示のための色々
# zは複素数のベクトル
# int0,int1はintensityの上下限、sat0,sat1はSaturation(彩度)の上下限
my.hsv <- function(z,int0=0.6,sat0=0.3,int1=1,sat1=1){
# 複素数の偏角
	arg <- Arg(z)
	s <- which(arg<0)
	arg[s] <- arg[s]+2*pi
# 複素数の絶対値
	r <- Mod(z)
# 絶対値が非常に大きくてもそこそこの色になるように対数変換
	s <- which(r>1)
	r[s] <- log(r[s])
# 絶対値で周期性が出るように4のmod
	r. <- 4*(r%%1)
	k <- floor(r.)
	r. <- r.-k
# 明度が上限、明度が下限、彩度が上限、彩度が下限の４パターンを
# 4のmodに対応づける
# 明度・彩度を動かすときは、複素数の絶対値で１次線形変化
	inten <- sat <- rep(0,length(r))
	s <- which(k==0)
	inten[s] <- int1
	sat[s] <- sat1-(sat1-sat0)*r.[s]
	s <- which(k==1)
	inten[s] <- int1-(int1-int0)*r.[s]
	sat[s] <- sat0
	s <- which(k==2)
	inten[s] <- int0
	sat[s] <- sat1-(sat1-sat0)*(1-r.[s])
	s <- which(k==3)
	inten[s] <- int1-(int1-int0)*(1-r.[s])
	sat[s] <- sat1

	return(cbind(arg,inten,sat))
}
my.hsv2rgb <- function(h,s,v){
# 色相の6 のmodでぐるりの情報を作る
	hi <- floor(h/(2*pi)*6)
	hi[which(hi==6)] <- 0
# 色相のぐるりの余りをfに入れ、それと明度・彩度とでp,q,tという３変数を決める
# ３変数を色相からの値を取らせる１つの原色を除いた２原色の値を定めるために使う
# 使い方は巡回させることでうまいことやる
	f <- (h/(2*pi)*6) %%1
	p <- v*(1-s)
	q <- v *(1-f*s)
	t <- v *(1-(1-f)*s)
	r <- g <- b <- rep(0,length(h))
	s <- which(hi==0)
		r[s] <- v[s];g[s] <- t[s]; b[s] = p[s];
	s <- which(hi==1)
		r[s] <- q[s];g[s] <- v[s]; b[s] = p[s];
	s <- which(hi==2)
		r[s] <- p[s];g[s] <- v[s]; b[s] = t[s];
	s <- which(hi==3)

		r[s] <- p[s];g[s] <- q[s]; b[s] = v[s];
	s <- which(hi==4)
		r[s] <- t[s];g[s] <- p[s]; b[s] = v[s];
	s <- which(hi==5)
		r[s] <- v[s];g[s] <- p[s]; b[s] = q[s];
	return(cbind(r,g,b))
}
# 三角形・三角錐…の内部かどうかの判定
judge.inside <- function(x,V){
	n <- length(x)
	K <- matrix(V[1:n,],ncol=n)
	K. <- t(K)-V[n+1,]
	tmp <- solve(K.,x-V[n+1,])
	all(tmp>=0 & sum(tmp)<=1)
}
# 内部かどうかの判定のために、頂点ごとの係数を返す
judge.inside.2 <- function(x,V){
	n <- length(x)
	K <- matrix(V[1:n,],ncol=n)
	K. <- t(K)-V[n+1,]
	tmp <- solve(K.,x-V[n+1,])
	c(tmp,1-sum(tmp))
	#all(tmp>=0 & sum(tmp)<=1)
}

# ぐねぐねなカラー絵を作る
x <- seq(from=-4,to=4,len=100)
xx <- as.matrix(expand.grid(x,x))
z <- xx[,1]+1i * xx[,2]
# 複素関数
my.f <- function(z){
	(z^2-1)*(z-2-1i)^2/(z^2+2+2*1i)
}
# 複素数に色を付ける
w <- my.f(z)
hsv <- my.hsv(w,int0=0.1,sat0=0.1,int1=1,sat1=1)
col <- my.hsv2rgb(hsv[,1],hsv[,3],hsv[,2])

plot(xx,pch=20,col=rgb(col[,1],col[,2],col[,3]))

# 適当に点を取る
n.pt <- 20
pt.id <- sample(1:length(xx[,1]),n.pt)
X <- xx[pt.id,]

# その点のドロネー図
library(geometry)
#X <- matrix(rnorm(50),ncol=2)
delaunay.X <- delaunayn(X)
plot(X)
delaunay.X.2 <- rbind(delaunay.X[,1:2],delaunay.X[,2:3],delaunay.X[,c(3,1)])
segments(X[delaunay.X.2[,1],1],X[delaunay.X.2[,1],2],X[delaunay.X.2[,2],1],X[delaunay.X.2[,2],2])

# 距離関係を歪ませる
X. <- X[,1] + 1i * X[,2]
mod.X <- Mod(X.)
arg.X <- Arg(X.)
#new.mod.X <- mod.X
new.mod.X <- mod.X^2/max(mod.X)
#new.mod.X[which(X[,1]>0)] <- new.mod.X[which(X[,1]>0)]^2
new.arg.X <- arg.X + mod.X * 0.1
X.. <- new.mod.X * (cos(new.arg.X) + 1i * sin(new.arg.X))

X. <- X
delaunay.X <- delaunayn(X)

plot(X.)
delaunay.X.2 <- rbind(delaunay.X[,1:2],delaunay.X[,2:3],delaunay.X[,c(3,1)])
segments(X.[delaunay.X.2[,1],1],X.[delaunay.X.2[,1],2],X.[delaunay.X.2[,2],1],X.[delaunay.X.2[,2],2])
X.. <- cbind(Re(X..),Im(X..))
par(mfcol=c(1,2))
plot(X.)
plot(X..)
par(mfcol=c(1,1))
par(mfcol=c(1,2))
plot(X.)
segments(X.[delaunay.X.2[,1],1],X.[delaunay.X.2[,1],2],X.[delaunay.X.2[,2],1],X.[delaunay.X.2[,2],2])

plot(X..)
segments(X..[delaunay.X.2[,1],1],X..[delaunay.X.2[,1],2],X..[delaunay.X.2[,2],1],X..[delaunay.X.2[,2],2])
par(mfcol=c(1,1))

# 標本点の真座標に基づいて補間座標を出す
new.xx <- matrix(NA,length(xx[,1]),length(xx[1,]))
# すべてのドロネー三角ごとにループする
for(i in 1:length(delaunay.X[,1])){
	tmp.v <- delaunay.X[i,]
	V <- X[tmp.v,]
	# 扱っているドロネー三角の内部の点を確認
	tmp <- t(apply(xx,1,judge.inside.2,V))
	tmp.more.0 <- apply(tmp>=0,1,all)
	#tmp.less.1 <- (apply(tmp,1,sum)) <= 1
	#insiders <- which(tmp.more.0 & tmp.less.1)
	insiders <- which(tmp.more.0)
	# 内部点の座標を線形補間
	if(length(insiders) > 0){
		new.xx[insiders,] <- tmp[insiders,] %*% X..[tmp.v,]
	}
	
	
}
xlim <- ylim <- range(c(c(xx),c(X..)))
par(mfcol=c(1,2))
plot(xx,pch=20,col=rgb(col[,1],col[,2],col[,3]),xlim=xlim,ylim=ylim)
plot(new.xx,pch=20,col=rgb(col[,1],col[,2],col[,3]),xlim=xlim,ylim=ylim)
par(mfcol=c(1,1))