Hatena::ブログ(Diary)

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

Prima Project

2016-05-23

MikuHatsune2016-05-23

球状の物体をメルカトル図法で2次元にする

物体をいじって球にしている。いま、球になった物体があるとして、これを2次元平面に展開して眺めたい。

地球儀が世界地図になるメルカトル図法みたいな感じ。

 

もとのうさぎは、各頂点の(離散)ガウス曲率にしたがって色付けしてある。正なら赤、0近くは灰色、負は緑である。

f:id:MikuHatsune:20160523180119p:image

足の裏は平らなのでほぼ0、背中も0。凸凹に応じて赤か緑である。しっぽはでっぱっているので赤、しっぽの付け根は緑である。

f:id:MikuHatsune:20160523180120p:image

 

球に位相変換したものである。正球に近づくと、半径がRとすれば、ガウス曲率はすべての頂点で¥frac{1}{R^2}に近づく。

色をつけるとほとんど1色になるが、もとの頂点と辺で作られる三角形表面がどこに移動しているかを示すため、最初のうさぎの形に応じて色付けしたガウス曲率をそのまま維持しておくとする。

すると、球になったときに、もともと足の裏だったところとか、しっぽだったところの名残がわかりやすい。

f:id:MikuHatsune:20160523180122p:image

f:id:MikuHatsune:20160523180121p:image

 

メルカトル図法への展開は、sphereplot もしくはSphericalCubatureパッケージを使う。使ってみた感触としてはsphereplot のほうが手軽な感じがする。

library(sphereplot)
xyz <- c(1, 1, 0) # xyz 座標
s <- car2sph(x=xyz[1], y=xyz[2], z=xyz[3], deg=FALSE)

library(SphericalCubature)
rect2polar(xyz)

 

座標データをshpereplot を使って緯度・経度にする。半径はデフォルトでは地球の半径になっているので、球に変換したオブジェクトの半径にする。

半径は体積から計算できる。

s <- car2sph(x=xyz[1], y=xyz[2], z=xyz[3], deg=FALSE)
merc2d <- mercator(s[,1:2], r=ra)

 

メルカトル図法をやると、3次元である球を2次元展開しているから、地球儀でいうところの北極南極が間延びする。

しかも、どうしても連続な地面が分断してしまう。

この状態で、もとの三角形を再構成すると、こんな感じのクッソ長い三角形ができる。

f:id:MikuHatsune:20160523180124p:image

 

間延びするのは北極南極で、三角形の辺がはんぱなく伸びるから、あまりにも長い辺を持つ三角形は、塗るのをやめよう、という感じにしておくと、まあ、どうにかなる。

すべての三角形の辺の長さを計算しておく。いくらか重複があるけど気にしない。

# 距離の分布
di <- c(mapply(function(x) c(dist(merc2d[ new_v$f[x,], ])), seq(nrow(new_v$f))))
dq <- quantile(di, 0.99) # 99% の三角形は塗る

plot(merc2d, pch=16, cex=0.03)
for(j in 1:nrow(new_v$f)){
  tmp <- merc2d[ new_v$f[j, c(1:3,1)],]
  if( all(c(dist(tmp[1:3,])) < dq) ){
    polygon(tmp, col=cols[j], border=NA)
  }
}

library(MASS)
x0 <- merc2d[,1]; y0 <- merc2d[, 2]
kd <- kde2d(x0, y0, c(bandwidth.nrd(x0), bandwidth.nrd(y0)), n=1000)
contour(kd, xlab="latitude",ylab="longitude", add=TRUE, col="blue", lwd=4)

 

三角形の面積でも判定しようと思った。平面上の3点(x_1, y_1), (x_2, y_2), (x_3, y_3)で囲まれる三角形の面積A

¥begin{align}A&=|(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)|/2¥¥&=¥begin{vmatrix}x_2-x_1 & x_3-x_1 ¥¥ y_2-y_1 & y_3-y_1¥end{vmatrix}/2¥end{align}

で表されるが、辺の距離ほど分別能がよくなかったので不採用。

# 三角形の面積の分布
ai <- mapply(function(x){abs(det(merc2d[ new_v$f[x,][1:2],] - merc2d[ new_v$f[x,][3],]))/2}, seq(nrow(new_v$f)))

f:id:MikuHatsune:20160523180123p:image

カーネル密度かぶせてみたりとか。

f:id:MikuHatsune:20160523213940p:image

2016-05-19

update したらログインできなくなるという謎仕様

windows 10 よりもさらに凶悪な仕様である。

 

そもそものモチベーションは、rstan を使おうと思ったらまったくインストールできない。

以前はcran に登録されていなかったので、devtools からrstan サイトを指定してうんぬんかんぬんだったが、現在ではinstall.package から普通に入りますよと

 

入りますよと

 

 

!?!?

 

入らない。

ggplot2 がいるとかplyr がいるとかRcpp がいるとか、そういうのを逐一最新ではないバージョンを入れたりしていたのだが、やっぱり入らない。

で、原因はR3.0.2 だったからなのだが、R のメジャーアップデートというのが一筋縄では行かない。少なくともR 側からはできない(はず)。

 

R はそもそも

sudo apt-get install r-base

でいれていたが、これだと、ソフトウェアセンターに登録してあるもののうち、現行のubunutu が認識している最新のものしかダウンロードインストールしてこない(たぶん)。

というわけで新しいR をインストールできるようにする

公開鍵はやらなくてもできた。

sudo gedit /etc/apt/sources.list

Geditが起動したらテキストの最終行に下記の一行を追加する。

deb http://cran.ism.ac.jp/bin/linux/ubuntu trusty/

青文字の部分はUbuntuのバージョン名で、対応は下記のとおり。

12.04:precise

14.04:trusty

15.04:vivid

15.10:wily

16.04:xenial

 

ちなみにこれは、日本のcran ミラーといえばtsukuba だったのに、いつのまにか統数研になっているから、ミラーが変わればまた変わるんだろう(適当)。

sudo apt-get update

これでcran において最新のR がインストールできる。

 

R3.3.0 になれば、rstan がinstall.package からインストールできた。

 

が、ここで単純には終わらず、PC を再起動したら、ubuntuログイン画面で何故かログインできない。

パスワードを正しくいれても、セッションの開始に失敗しました、と出る。

こちらを参考にしてやってみる。

端末にログイン、というのは、そもそもGUI ログインできていないので、Ctrl * Alt + F1CUI ログインする。

ここでは、username を入力して、そのあと、パスワードを入力すればよい。

その後、

sudo apt-get install ubuntu-desktop

とすればよい。再インストールはやらなくてもできた。

しかし、ネットワークが繋がっていない場合は詰むんじゃないかこれ。

 

というわけで、不用意なupdate はやめよう。

2016-05-13

Eigen で行列式

体積と表面積の計算をRではなくcpp のままやれば速いということで、Eigen の協力をあおぐわけだが、Eigen の行列式は使うのが少し手間が必要で、Eigen/LUもinclude する必要がある。

matrix.determinant()

で求められる。

2016-05-10

物体の体積と表面積

物体の体積と表面積を計算したい。

三角メッシュ化されていれば、表面の三角形を足せば表面積はどうにかなる。

問題は体積。

 

The GNU Triangulated Surface Library(GTS)というもので計算できるっぽい。

sudo apt-get install libgts-dev

 

Python でも使えるようになっている。しかし、pip でインストールできなかったのでソースからインストールする。

 

例によって説明など存在しない。

 

GTS が使うのはgts ファイルという、こんな感じのファイル。

一辺10の正四面体の頂点、辺、三角形の数と

座標

エッジグラフ

エッジで構成される三角形

をひたすら並べている。

8 18 12
 0  0  0
10  0  0
10 10  0
 0 10  0
 0  0 10
10  0 10
10 10 10
 0 10 10
1 2
2 3
3 4
4 1
5 6
6 7
7 8
8 5
1 5
2 6
3 7
4 8
1 6
2 7
3 8
4 5
1 3
5 7
13 1 10
13 5 9
14 2 11
14 6 10
15 3 12
15 7 11
16 4 9
16 8 12
17 2 1
17 4 3
18 5 6
18 7 8

 

これをpython gts で体積、表面積の計算をする。

# Python 2.7
import gts
f = "cube.gts"
g = open(f)
s = gts.read(g) # gts 読み込み

s.volume() # 体積
s.area()   # 表面積

 

問題はgts ファイルというやつで、obj ファイルからの変換がめんどくさい。

gtsの変換プログラム作ったよというメーリスも動かない。

 

Rvcg パッケージには、vcgGetEdge(mesh3d obj) という関数があって、三角形エッジをグラフにしてくれるが、これを使って上のcubeをいじってみると

library(Rvcg)
e <- vcgGetEdge(mesh.tri)
edge <- mapply(function(x) which(sapply(apply(t(mesh.tri$it), 1, setdiff, e[x, 1:2]), length) == 1), seq(nrow(e)))


r <- vector("list", max(edge))
for(i in seq(ncol(edge))){
  r[edge[,i]] <- lapply(r[edge[,i]], c, i)
}
r <- do.call(rbind.data.frame, r)
colnames(r) <- NULL

こんなファイルになる。番号が少し入れ替わる。

8 18 12
 0  0  0
10  0  0
10 10  0
 0 10  0
 0  0 10
10  0 10
10 10 10
 0 10 10
1 2
1 3
1 4
1 5
1 6
2 3
2 6
2 7
3 4
3 7
3 8
4 5
4 8
5 6
5 7
5 8
6 7
7 8
1 5 7
4 5 14
6 8 10
7 8 17
9 11 13
10 11 18
3 4 12
12 13 16
1 2 6
2 3 9
14 15 17
15 16 18

 

これを使ってpython gts をやってみると、area() はできるが volume() で

RuntimeError: Surface is not orientable

と言われてエラー。

 

それでは世の中に転がっているものから obj → gts 変換をしようとして、Hira 3D viewerというものでobj ファイル読み込みからのgts 出力をしてみるのだが、gts で頂点が重複して出力されるのでどうにもならん。

だがしかし、GUI で体積は計算できる。

 

がんばって変換スクリプト書くしかなさそう。

 

と思っていたら計算できるっぽい。これとかこれとか。

4点(0,0,0), (x_1, y_1, z_1), (x_2, y_2, z_2), (x_3, y_3, z_3)で作られる四面体をひたすら足せばいいっぽくて、この物体の体積は

¥begin{vmatrix}x_1& y_1& z_1 ¥¥x_2& y_2& z_2¥¥x_3& y_3& z_3 ¥end{vmatrix}/6

と、行列式で求められる。

sum(mapply(function(x) det(bun.v[bun.f[x,],])/6, seq(nrow(bun.f))))

2016-05-09

四元数の演算

四元数内積が欲しかったのでメモ。

四元数qは2つの虚数単位i, j, (k=ij), と4つの実数a, b, c, dを用いて

q=a+ib+jc+kd

と書ける。これを実数部分(スカラー部)と虚部(ベクトル部)に分けて

q=(x,¥bf{u})

とも書く。

いま、ふたつの四元数q_1=(x,¥bf{u}), q_2=(y,¥bf{v})について、積は

q_1*q_2=(xy-¥bf{u} ¥cdot ¥bf{v},x¥bf{v}+y¥bf{u}+¥bf{u}¥times¥bf{v})

となる。

 

計算して確かめる。三角形の面積も計算してみる。

set.seed(1234)
xyz <- matrix(runif(6), 3)
q12 <- as.quaternion(rbind(0, xyz))

# xyz での内積
xyz[,1] %*% xyz[,2]

# 四元数での内積
prod(Re(q12)) - Re(q12[1] * q12[2])

# xyz での三角形の面積
sqrt(sum(mapply(function(r) det(cbind(rbind(t(xyz[c(r, r%%3+1),]), 0), 1))^2 ,seq(3))))/2
 
# 四元数での三角形の面積
Mod(Im(q12[1] * q12[2]))/2

 

cpp ではEigenライブラリを使うと爆速だった。簡単な使い方はこちらにまとまっている。