Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2016-06-25

MikuHatsune2016-06-25

アレのサイズの自己申告に基づくデータに見られるheaping をrstanで解析する。

rstan神がやっている人口ピラミッドのAge Heapingを階層ベイズで補正するという解析について、自己申告に基づくデータにheapingが見られるという別の話があったので、3番煎じだがrstanを使わねば!!

 

上記のブログにある分布は、ギザギザしすぎているし、データを作成するにしても難しいので竿の直径のデータを使う。

最大で5cm くらいだが、いくつかのポイントで突出している。元々は50万人の自己申告に基づくデータを収集したらしく、元データほしいと思いながらも、分布の尖っているところをlocator 関数でポチポチする。

0.017 0.000
1.062 0.111
1.140 0.590
1.238 0.093
1.365 0.119
1.414 0.198
1.521 0.119
1.863 0.198
1.970 1.707
2.088 0.303
2.215 0.224
2.517 0.320
2.586 0.608
2.723 0.355
2.908 0.425
2.996 3.190
3.152 0.555
3.279 0.634
3.455 0.555
3.504 1.218
3.670 0.651
3.846 0.555
3.963 0.643
3.992 2.195
4.139 2.056
4.197 0.198
4.324 0.276
4.461 0.154
4.529 0.381
4.656 0.119
4.774 0.398
4.852 0.093
4.940 0.041
4.989 0.241
5.116 0.041

f:id:MikuHatsune:20160624192716p:image

 

その後、線形補間して適当に水増しする。

# n にデータが入っている
l <- function(x, v1, v2){
  coef <- solve(n[c(i+1, i),], c(1, 1))
  -coef[1]/coef[2] + 1/coef[2]
}

x <- seq(0, 5, by=0.05) # x を水増し
y <- rep(0, length(x))  # x に対応する点
for(j in 2:length(x)){
  i <- tail(which(n[,1] < x[j]), 1)
  coef <- solve(n[c(i+1, i),], c(1, 1)) # ax + by = 1 を解く
  y[j] <- -coef[1]/coef[2]*x[j] + 1/coef[2]
}
plot(x, y, type="l", lwd=3, col=4, xlab="Penile diameter (cm)", ylab="人数 (万人)")
abline(v=0:5, h=0:4);abline(v=0:4+0.5, lty=3, h=0:3+0.5)

f:id:MikuHatsune:20160624192717p:image

 

さて、ここでデータについて考えよう。

局部についてのデータなので、かなりデリケートである。匿名のアンケート、データ収集といえども、羞恥心や見栄をはりたい、という気持ちはよくわかる。

さらに、局部をものさしなどの測定器具でしっかり挟んで直径を測定しない限り、ミリ単位では計測できないし、つまんだりにぎったりしたときのだいたいの大きさでcm を報告してしまうかもしれない。せめて、なんとなくわかりそうな0.5cm 刻みでは報告しそうである。

ということを考慮すると、rstan神の年齢解析から改変すると、

・下のカテゴリーから上へ過剰申告する

・2.4cm の人は、2.1とか2.2cm の人よりかは2.5cm にしてしまったほうが見栄えがいいので、よく2.5cm と報告する(これはまあ、2.1cmでもより上に報告するとかはありそうだけれども)

さらに、0.5cm 刻みだと楽なんだが、グラフでは0.5cm 刻みより外れていそうなところがあるので

・Ncm + 0.5cm + rcm のなんらかの項が存在する。

みたいなモデルを考えたい。

けどとりあえず、0.5cm 刻みで、データとしては0.05cm 刻みで報告されているデータの、とあるN.5cm より小さい5群くらいが過剰報告してくる(つまり、fromが常に下側、上のtoへ流入してくる)とする。

 

やってみたけど、全然滑らかにならなかった。もっといじりたかったけど時間がないし、こういうのは旬が大事()なので。

95% 信用区間と推定値。

f:id:MikuHatsune:20160624192718p:image

diam_idx <- which(x %in% seq(from=0, 5.5, by=0.5))
x_idx <- cut(x, diam, include.lowest=TRUE)
li <- lapply(-(5:1), function(j) data.frame(to=diam_idx, from=diam_idx + j))
diam_heap <- subset(do.call(rbind, li), from > 0)

stanmodel <- stan_model("penile_stan.stan") # rstan 神のスクリプトそのまま
dat <- list(A=length(x),
            Y=round(y*10000), J=nrow(diam_heap), From=diam_heap$from, To=diam_heap$to)
fit <- sampling(stanmodel, data=dat, chain=3, seed=123)

ex0 <- extract(fit)
qu <- apply(ex0$q, 2, quantile, c(0.025, 0.5, 0.975))*sum(y)

plot(x, y, type="n", xlab="Penile diameter (cm)", ylab="推定")
polygon(c(x, rev(x)), c(qu[1,], rev(qu[3,])), border=NA, col=grey(0.9))
lines(x, qu[2,], lty=3, lwd=3, col=4)
abline(v=0:5, h=0:4);abline(v=0:4+0.5, lty=3, h=0:3+0.5)

2016-06-24

MikuHatsune2016-06-24

ラブライブ! サンシャイン!! のAqours メンバーがラブライブのμ's の誰に相当するか考える。

ラブライブ!サンシャイン!が始まりそうである。

無印ラブライブの時はまったく注目していなかったので、アニメが始まってから(なんだこれ…)と思いながら見ていたが、3話くらいになると(ブヒィ)とか言ってた。

ラブライブは9人、サンシャインも(たぶんメンバー構造を模倣して)9人いるので、誰が誰に相当するのかを考えたい。

公式HP では、ラブライブではほのかすの隣にエリーチカがいて、主人公に対して重要そうな幼なじみポジであることりが3番手だったり、サンシャインのほうでは高海千歌が主人公っぽそうで2番手に転校生である桜内梨子、3番手に3年生で幼なじみの松浦果南がいる、という感じで、公式HP はアテになるようなならないような。

1,2,3年生が各々3人という構図、幼なじみや生徒会長属性というのは維持されているっぽいので、各々プロフィールを読むと、黒澤ダイヤが生徒会長らしい。ただし、エリーチカとは学年以外似ても似つかず、むしろ小原鞠莉のほうが設定、外見的に近い。

 

というわけで、付属情報からはなかなかよくわからないので、定量データから考える。ラブライブアイドルらしく、各メンバーのスリーサイズが設定されているので、これを使う。

てっとり早くクラスタリングをしてみよう。

綾瀬絵里と小原鞠莉は、あらゆるパラメータボインなキャラとしてすごい近いところにいるし、他には、南ことりと桜内梨子、高坂穂乃果と津島善子などが近くなっているようである。

しかし、ここで問題なのが、9*2人で近そうな人からクラスタリングされているので、対応がどうなるか、というのは考えにくい。例えば、綾瀬絵里と小原鞠莉が近いのはいいが、これらに次に近いのは松浦果南、その次は東條希である。ここをもってして松浦果南と東條希が近くて、対応すると言っていいのか。

library(gplots)
heatmap.2(as.matrix(d2), col=colorpanel(20, "blue", "grey", "red"), scale="column", Colv=FALSE, tracecol="black", hclustfun=function(x){hclust(x, "ward")})

f:id:MikuHatsune:20160624114340p:image

 

PCA をしてみる。PC1 はグラマラスさ、PC2 はギャップロリのような感じの主成分が以前と同じように取れる。例えば黒澤ルビィはPC1 が大きいが、すべてのベクトルから負方向であり、幼少ロリである。一方、綾瀬絵里と小原鞠莉はすべてがボンッキュッボンッである。

黒澤ダイヤはPC2 が負であり、身長年齢に比してスリーサイズは控えめである。一方で小泉花陽と国木田花丸は年齢身長が小さい割にグラマラスである。

というのは、特徴的なパラメータを持つキャラならばわかりやすいが、例えばPC1 もPC2 も0近辺にいる、いわば普通の特徴のキャラは対応困難である。

f:id:MikuHatsune:20160624114341p:image

 

ラブライブの9人の集合L_1に属するμ's メンバーm_i ¥in L_1, i=¥{1,2,¥dots,9¥}と、サンシャインの9人の集合L_2に属するAqours メンバーa_i ¥in L_2, i=¥{1,2,¥dots,9¥}について、m_i=a_jの一対一対応を考えたい。

と思っていたら、結局、9!の組み合わせを総当りして計算量は物理で殴ればいいので、ベタに差分をとって差が最小となる組み合わせを探す。

per <- permutations(9, 9, seq(9))
di <- mapply(function(z) sum(abs(l[per[z,],]-a)), 1:nrow(per))
l[per[which.min(di),],];a
           age height  B  W  H
東條希      17    156 90 60 82
南ことり    16    159 80 58 80
西木野真姫  15    161 78 56 83
園田海未    16    159 76 58 80
高坂穂乃果  16    157 78 58 82
星空凛      15    155 75 59 80
小泉花陽    15    156 82 60 83
絢瀬絵里    17    162 88 60 84
矢澤にこ    17    154 74 57 79
           age height  B  W  H
高海千歌    16    157 82 59 83
桜内梨子    16    160 80 58 82
松浦果南    17    162 83 58 84
黒澤ダイヤ  17    162 80 57 80
渡辺曜      16    157 82 57 81
津島善子    15    156 79 58 80
国木田花丸  15    152 83 57 83
小原鞠莉    17    163 87 60 84
黒澤ルビィ  15    154 76 56 79

主人公と思われる高海千歌は東條希に対応する。うちで9人や!とか言いそう(主人公なのに

南ことりは桜内梨子である。転校生属性に変わった。

西木野真姫は主人公の幼なじみである松浦果南である。まきちゃんは見た目はしっかり者お姉さんキャラのように見えて、にこにーのカップリング要員でしかなかったから、今後のカップリングが誰になるのか(もはやキャラ対応関係ない

園田海未は黒澤ダイヤである。スリーサイズが近いだけで属性はさらさらストレートしかかぶっていない。

高坂穂乃果は渡辺曜である。もう主人公はこの娘でいいんじゃないかな(適当

星空凛は津島善子である。津島善子はキャラデザだけみれば黒髪ロングニーソで一番のブヒ要素を持っている、と思ったら頭の中がお花畑のようで、りんちゃんと見事な対応だね()

小泉花陽は国木田花丸である。いわゆる年下ロリ巨乳枠で一致している。

綾瀬絵里と小原鞠莉はheatmap でも一致しているように、よく合っている。声優鈴木愛奈なので少し期待している。

矢澤にこは黒澤ルビィである。純粋なロリ枠だったが、年齢が3年生から1年生へとかわり、さらにロリ化している。

2016-06-11

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

COI:筆者とはなんの関係もありません。

 

読んだ。

細かいところはstan 神が書いているので、rstan ユーザーもしくは生物系解析者の立場で書いてみる。

 

まとめ

・具体例がたくさんあり、解釈の仕方も丁寧に書いてあるので、買い。

・新時代の統計の教科書には な り え な い と思う。(辛口)

・初心者向け(大嘘

 

良い点

・具体的

これに尽きると思います。

サイトからスクリプトが入手できるように、全部自分の環境で再現できる。

また、各章の解析で、RQ(Research Question)を設定しているが、例えば

RQ.2 第1群と第2群の平均値の差¥mu_{1}-¥mu_{2}の点推定。平均値の差の推定。

(ダイエットプログラムに参加した人と参加前の人では平均値に何kg の差がありますか。またその原料はどの程度の幅で確信できるでしょう。95% の確信で答えてください。)

p97 対応のある2群の差と相関の推測

RQ.4 平均値の差が基準点c より大きい確率。

(ダイエットプログラムは有料だし、参加すると辛いこともあるので、痩せても少しではベイしません。たとえばダイエットプログラムに参加した人と参加前の人の平均値¥mu_{1}-¥mu_{2} が2kg より大きい確率が70% より大きいならば参加したいです。参加すべきでしょうか、あるいは見送るべきでしょうか。)

p97 対応のある2群の差と相関の推測

という問題は、かなり、自然な思考回路であり、データ解析上、よく思うことだと思います。

古典的なt検定だと、H_0:¥mu_{1}-¥mu_{2}=0H_1:¥mu_{1}-¥mu_{2}¥not{=}0ということはわかるけれども、それ以上というかこれ以外のことは(厳密には)証明しにくいので、いいと思う。

また、分割表でも

連言命題が正しい確率

表6.5 の確率は、2つのカテゴリの比率の確率としてはそのまま解釈可能です。ただし複数の比較が同時に成り立つ確率とは異なります。

手始めに、研究上の問い「「友達」は他の誰よりも相談される比率が高い」が正しい確率を求めてみましょう。

(中略)

次に、研究上の問い「「先生」は他の誰よりも相談される比率が低い」を考えてみましょう。

(以下略)

p145 比率とクロス表の推測

というように、古典的なχ二乗検定だと、全体がだいたい同じか、何かが違うか、ということがぼんやりわかるだけで、「どことどこが比較して」とか、「組み合わせてどうか」が直接わかりやすいので、いいと思う。

 

悪い点

・「はじめての」という大嘘

初級者が読めるとは到底思えない。

 

この本は初めて統計学に入門する学生のための教科書です。

まえがき

筆者がそういうのなら、まあ、そうだと信じよう。しかし、

統計データ分析に関する予備知識はいっさい仮定していません。

うーん、まあ、本当に予備知識が0 なら難しい。分布とか、統計量とか、ベイズの定理(条件付き確率程度)とか、そういうことは少なからず必要。

数学説明には、微分積分シグマ記号・行列・ベクトル演算を使いません。

とあるように、確かに、これらの記述がない。かわりに、ベタに書き下した記述があるし、確率としてP(x|a,b)みたいな記述は頻発する。

数学の知識0 でやるのも、ぶっちゃけ言うと、不可能。

 

・R はなんだかんだで難しい

たぶん、R 初心者でいきなりrstan を始める人はいないと思うが、rstan は簡単にベイズ推定ができるとはいえ、初心者がこの本だけでいろいろやるのは難しいので、みんなググラビリティをあげよう。

 

・アンチは本当はファンじゃないのか

この本は実はアメリカ統計学会のp値はクソ記事から始まるのだが、ベイズ的な統計学をp値的な統計と対比させて本書が進むので、やはりp値的な解析の問題点や限界を知ったうえで読むべきなのである。という点では、p値はクソだと言いながらp値に理解のある大人な本なのである。

(ファンとか儲とかアンチの定義が違うと言われたらアレだが、まあ、適当)

前も言ったように、われわれの業界では、なんだかんだでp値派なので、

たぶん、お偉い人たちががんばってp値以外の他の統計量や指標、モデル、考え方を流布してくれると思うので、当面は「p値いいっすね〜(アヘ顔」って上司にいい顔をしながら、その他の手法について勉強するしかないんじゃないかな(適当

となる。ベイズ派が優勢となるまで隠れて爪を研いでいたらいいんじゃないかな(適当

なので、勉強会などで引っ張りだこになるが、p値がベイズに置き換わるのは少なくとも10年はなさそうなので、その頃にはもっと、いい教科書になりうる良書がでてるんじゃないかな(適当

2016-06-10

MikuHatsune2016-06-10

Ricci flow

読んでる。Comput Vis Image Underst. 2013 Sep 1;117(9):1107-1118.

Ricci flow という、Willmore flow とはまた違う条件で微分幾何学的物体変換を行う。

その前に下準備。

 

以前やった(離散)ガウス曲率は、1-ring下の三角形の角度と、A_{mixed}と呼ばれる傘の面積で計算していたが、これはMeyer method というらしい。

本来は、ガウス曲率は角度だけで計算できて

K_i=2¥pi-¥sum_i¥theta_i

となる。しかし、この場合は、論文の図にもあるけれども、理想的な球体のうえにランダムに点を発生させて三角形を作った時に、ガウス曲率の濃淡ができてしまう。というわけで、疎密に合わせてガウス曲率を補正したい。

というわけでA_{mixed}で割る、らしい。

 

半径r の理想的な球を作っているので、ガウス曲率は¥frac{1}{r^2}=1である。

ほとんどが1の周辺に分布している。

f:id:MikuHatsune:20160610233328p:image

 

1からかけ離れたガウス曲率をもつ頂点をプロットした。赤が高いガウス曲率、青が低い。

ズームしてみると、すごい鋭角な三角形をもつ頂点が、低いガウス曲率になっている様子。

計算はあっているようである。

f:id:MikuHatsune:20160610233329p:image

f:id:MikuHatsune:20160610233330p:image

 

理想的な球を作り、点を配置し、三角形を作る。

library(rgl)
library(onion)
my.sphere.tri.mesh <- function(n.psi=30){
  thetas <- list()
  psis <- seq(from=-pi/2,to=pi/2,length=n.psi)
  #psis <- sort(runif(n.psi, min=-pi/2, max=pi/2))
    d.psis <- psis[2]-psis[1]
    hs <- sin(psis)
    rs <- sqrt(1-hs^2)
    ls <- 2*pi*rs
    n.thetas <- floor(ls/d.psis)
    thetas[[1]] <- c(2*pi)
    for(i in 2:(n.psi-1)){
        thetas[[i]] <- seq(from=0,to=2*pi,length=n.thetas[i]+1)
        thetas[[i]] <- thetas[[i]][-(n.thetas[i]+1)]
    }
    thetas[[n.psi]] <- c(2*pi)
    sapply(thetas,length)

    bridge <- list()
    for(i in 1:(n.psi-1)){
        a <- c(thetas[[i]],2*pi)
        b <- c(thetas[[i+1]],2*pi)
        bridge[[i]] <- matrix(c(1,1),1,2)
        loop <- TRUE
        while(loop){
            n.r <- nrow(bridge[[i]])
            id.a <- bridge[[i]][n.r,1] + 1
            id.b <- bridge[[i]][n.r,2] + 1
            if(id.a > length(thetas[[i]]) & id.b > length(thetas[[i+1]])){
                if(id.a-1!=1 & id.b-1!=1){
                    bridge[[i]] <- rbind(bridge[[i]],c(1,id.b-1))
                }

                loop <- FALSE
            }else{
                if(id.a > length(thetas[[i]])){
                    tmp <- c(id.a-1,id.b)
                }else if(id.b > length(thetas[[i+1]])){
                    tmp <- c(id.a,id.b-1)
                }else{
                    if(a[id.a] < b[id.b]){
                        tmp <- c(id.a,id.b-1)
                    }else{
                        tmp <- c(id.a-1,id.b)
                    }
                }
                bridge[[i]] <- rbind(bridge[[i]],tmp)
            }
        }
    }
    xyz <- matrix(0,0,3)
    edge <- matrix(0,0,2)
    triangles <- matrix(0,0,3)
  n.triangles <- rep(0,n.psi-1)
    for(i in 1:n.psi){
        n.r <- nrow(xyz)
        if(i > 1){
            pre <- (n.r-length(thetas[[i-1]])+1):n.r
            post <- (n.r+1):(n.r+length(thetas[[i]]))
            edge <- rbind(edge,cbind(post,c(post[-1],post[1])))
            br <- bridge[[i-1]]
            new.edge <- cbind(pre[br[,1]],post[br[,2]])
            edge <- rbind(edge,new.edge)
            tmp.tri <- cbind(new.edge,rbind(new.edge[-1,],new.edge[1,]))
            tmp <- apply(tmp.tri,1,unique)
            triangles <- rbind(triangles,t(tmp))
      n.triangles[i-1] <- length(tmp[1,])
        }
        psi <- psis[i]
        theta <- thetas[[i]]
        xyz <- rbind(xyz,cbind(cos(psi) * cos(theta),cos(psi)*sin(theta),sin(psi)))

    }
    return(list(xyz=xyz,edge=edge,triangles=triangles,n.psi=n.psi,thetas=thetas,n.triangles=n.triangles))
}
knit_hooks$set(rgl = hook_rgl)
n.psi <- 25
sp.mesh <- my.sphere.tri.mesh(n.psi)
spxyz <- sp.mesh$xyz
apply(spxyz^2,1,sum)
xx <- 0.8
yy <- 0.1
large.xy1 <- which(spxyz[,1]>xx & spxyz[,2]<=yy)
large.xy2 <- which(spxyz[,1]>xx & spxyz[,2]>yy)

spxyz[large.xy1,1] <- xx + ((spxyz[large.xy1,1]-xx)*10)^2
spxyz[large.xy2,1] <- xx + ((spxyz[large.xy2,1]-xx)*10)^2
spxyz[large.xy2,2] <- yy + ((spxyz[large.xy2,1]-yy)*10)^5
plot3d(spxyz)
segments3d(spxyz[c(t(sp.mesh$edge)),])
mesh.tri <- tmesh3d(t(spxyz),t(sp.mesh$triangles),homogeneous=FALSE)

# 打点して
plot3d(sp.mesh$xyz)
# 三角形メッシュオブジェクトを作り
mesh.tri <- tmesh3d(t(sp.mesh$xyz),t(sp.mesh$triangles),homogeneous=FALSE)
# 三角形を灰色で塗る
shade3d(mesh.tri,col="gray")

四元数でごそごそしていたので、ガウス曲率の計算に四元数を入力する仕様になっている。

# cumpute Gaussian curvature of j_th vertex
# input: index of vertices of quanterion vertices object
gaussian_curvature <- function(j, vertices, faces.v){
  #for(j in seq(n_vertices)){
    plane <- t(faces.v[, apply(faces.v==j, 2, any)])
    # shuffle the order of summation of triangle areas
    p <- p1 <- 1 # start
    s <- c(j, setdiff(plane[p,], j)[1]) # 2nd vertex
    p1 <- c(p1, setdiff(which(apply(plane == tail(s, 1), 1, any)), p))
    while(p < nrow(plane) - 1){
      s <- c(s, setdiff(plane[tail(p1, 1),], s))
      p1 <- c(p1, which(sapply(apply(plane, 1, setdiff, s), length) == 1)[2])
      p <- p + 1
    }
    plane <- plane[p1,] # ordered sequentially

    circum_mat <- matrix(1, 3, 3)
    diag(circum_mat) <- -1
    theta_j <- rep(0, nrow(plane))
    circumcenter <- matrix(0, nrow(plane), 3)
    for(p in seq(nrow(plane))){
      hogemat <- t(as.matrix(vertices[plane[p,]])[-1,])
      rname <- plane[p,]
      e12 <- sweep(hogemat[!(rname %in% j),], 2, hogemat[ (rname %in% j),], "-")
      theta_j[p] <- acos( e12[1,] %*% e12[2,]/prod(sqrt(rowSums(e12^2))) )
      if(theta_j[p] < pi/2){ # non-obtuse angle compute circumcenter point
        abc <- mapply(function(k) sum(apply(hogemat[c(k%%3+1, (k+1)%%3+1),], 2, diff)^2), seq(3))
        # abc <- rowSums(rbind(apply(e12, 2, diff), e12)^2)
        abc2 <- circum_mat %*% abc * abc
        circumcenter[p,] <- colSums(diag(c(abc2)) %*% hogemat)/sum(abc2)
        # circumcenter <- (abc2[1]*hogemat[1,] + abc2[2]*hogemat[2,] + abc2[3]*hogemat[3,])/sum(abc2)
      } else { # obtuse angle
        circumcenter[p,] <- colSums(hogemat[!(rname %in% j),])/2 # center point
      }
    }
    mat <- sweep(circumcenter, 2, as.numeric(vertices[j])[-1], "-") # vector from vertex j_th
    n <- nrow(circumcenter)
    Amixed <- mapply(function(k) sqrt(prod(rowSums(mat[c(k, k%%n+1),]^2)) - (mat[k,]%*%mat[k%%n+1,])^2) / 2, seq(n))
    Gcurvature <- (2*pi - sum(theta_j))/sum(Amixed, na.rm=TRUE)
  #}
  return(Gcurvature)
}

# すべての頂点について計算
b <- mapply(function(j) gaussian_curvature(j, vertices, faces.v), seq(length(vertices)))
cut0 <- cut(b, quantile(b, c(0, 0.02, 0.98, 1)))
cols <- bluered(length(levels(cut0)))[cut0]

plot3d(mesh.tri, type="dits", axes=FALSE, xlab="", ylab="", zlab="")
shade3d(mesh.tri, col="grey")
wire3d(mesh.tri)
spheres3d(t(mesh.tri$vb), col=cols, radius=0.02)

2016-06-08

MikuHatsune2016-06-08

3次元オブジェクトをグラフにしてみる

平均曲率を用いて、3次元オブジェクトの表面を色付けした。

f:id:MikuHatsune:20160608122658p:image

凹んでいるところは赤、出っ張っているところは緑になっている。

 

いま、同じ色でひと続きになっている領域を仲間分けして、なおかつ、隣り合う領域がどうなっているかを調べたい。

ある頂点と、その頂点を含む三角形(1-ring)を選び、それに隣り合っている三角形をすべて、再度抽出する。

このとき、既に選ばれた三角形たちは、重複して抽出すると計算に時間がかかるので、適宜除く。

# 球面に色を付けた時、同じ色が付いている閉じた領域をチェックする
K0 <- NULL
for(k0 in 1:nrow(v1[[i]]$f)){
  if( !any(mapply(function(z) k0 %in% z, K0)) ){ # いままでのK に kが含まれていなかったら新規のtriangle
    print(paste("Processing", k0))
    k <- k0
    foo0 <- which(mapply(function(x) any(v1[[i]]$f[x,] %in% v1[[i]]$f[k,]), seq(nrow(v1[[i]]$f))))
    foo1 <- foo0[cols[foo0] == cols[k[1]]] # k の三角形と同じ色をもって隣り合っている三角形たち
    #plot3d(bun.v, type="n", axes=FALSE, xlab="", ylab="", zlab="")
    #shade3d(mesh.tri, col=rep(cols, each=3))
    #rgl.viewpoint(10, 30, zoom=0.7)
    m2 <- tmesh3d(t(v1[[i]]$v), t(v1[[i]]$f[foo1,]), homogeneous=FALSE)
    #wire3d(m2)
    K <- list(k, foo1)
    while( !all(K[[length(K)]] %in% K[[length(K)-1]]) ){
      k <- setdiff(K[[length(K)]], K[[length(K)-1]])
      foo0 <- which(mapply(function(x) any(v1[[i]]$f[x,] %in% v1[[i]]$f[k,]), seq(nrow(v1[[i]]$f))))
      foo1 <- foo0[cols[foo0] == cols[k[1]]]
      K <- c(K, list(foo1))
      #open3d()
      #plot3d(bun.v, type="n", axes=FALSE, xlab="", ylab="", zlab="")
      #shade3d(mesh.tri, col=rep(cols, each=3))
      m <- v1[[i]]$f[setdiff(K[[length(K)]], K[[length(K)-1]]),]
      m2 <- tmesh3d(t(v1[[i]]$v), t(v1[[i]]$f[foo1,]), homogeneous=FALSE)
      #points3d(matrix(v1[[i]]$v[m,], nc=3))
      #texts3d(matrix(v1[[i]]$v[m,], nc=3), texts=unique(c(v1[[i]]$f[k,])))
      #wire3d(m2, col=4)
    }
    K1 <- sort(unique(unlist(K)))
    K0 <- c(K0, list(K1))
  }
}

隣接する三角形が広がっていく様子。最前線しかプロットしていない。

f:id:MikuHatsune:20160608145124g:image

すべての隣り合う三角形をプロットした図。

f:id:MikuHatsune:20160608145123g:image

 

いま、ひと続きになっている領域をメッシュで色塗りしたが、このふたつの紫、水色の領域は隣り合っている。

f:id:MikuHatsune:20160608155845p:image

f:id:MikuHatsune:20160608155844p:image

隣り合う領域をadjacent matrix として、グラフ化。

背中の広い範囲の緑、首まわりの赤がおおきなノードとして表現される。赤いノードからは顔、耳2つが隣接していて、それもノードとして出ている。

ノードの大きさは、対応する領域の面積としている。

グラフ構造に落とせれば、ノードの密度や、ネットワーク解析を応用してそういう統計量もとれる。

# K0 は同じ色で塗られるtriangle のひとつづきの領域
# 隣り合う領域を探す
V0 <- mapply(function(z) unique(c(v1[[i]]$f[z,])), K0) # 含まれる頂点たち
adj_mat <- mapply(function(z2) mapply(function(z1) any(z2 %in% z1)+0, V0), V0)
g0 <- graph.adjacency(adj_mat, mode="undirected", diag=FALSE)
V(g0)$size <- mapply(function(z) sum(area0[[i]][z,]), K0) * 10
V(g0)$color <- mapply(function(z) cols[z[1]], K0)
plot(g0)

f:id:MikuHatsune:20160608155846p:image