Juliaでバーチャル地球儀を作ってみた

# はじめに

　Python並みに書きやすく、C++並みに速いという噂の[Julia](https://julialang.org/)に興味を持ち、試しに重た目のプログラムを書いてみました。
　やりたいことは、地球をグリグリ回して表示する地図アプリです。

# インストールしたもの

+ [Julia 1.4.2] (https://julialang.org/)
+ anaconda3にJulia用の環境を追加
+ Jupyter Notebookを追加
+ Julia用の各種パッケージ(Blink, ImageMagick, ImageView, Images, Interact, Plots, WebIO, IJulia(※Jupyter用))

普段Pythonを使っているので、anacondaで環境を分けました。
JuliaからPythonのmatplotlibなどを使うためにPython連携があると便利なようです。

# 作ったアプリ
見た目はこんな感じです。スライダーでグリグリ地球を回します。ズームもできます。
**(2020/6/28追記)とりあえず、どんなアプリか試してみたい人向けに、Windows10用のバイナリ(exe)を用意しました。**（JuliaではなくPython版です）
[GitHub:VRMapEarthApp.exe](https://github.com/okimebarun/VRMapEarth.jl/releases)

![VRMapEarth_AppImage2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597022/b286f1d4-5a41-5af0-3f68-dc1339391b65.jpeg)


#### コンセプト
開発コンセプトはこんな感じです。

+ [正距円筒図法(Equirectangular projection)](https://ja.wikipedia.org/wiki/%E6%AD%A3%E8%B7%9D%E5%86%86%E7%AD%92%E5%9B%B3%E6%B3%95)で描かれた地図を3Dの真球にマッピングする。
+ 地軸の角度と地軸を中心とした回転でグリグリ回す。
+ １点透視図法により、ズームイン・ズームアウトに対応。
+ 簡易レイトレーシング。陰影も付ける。
+ 背景も付ける。
+ Juliaの関数型言語っぽい特徴を使って、関数の合成をしまくる。

#### 幾何学的な部分
下図のような関係にある、視点、画面座標系、球面座標系、背景座標系の
視点からの直線（視線）との交点を計算しながら画面座標系にマッピングしていきます。

![レイトレーシング原理.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597022/5fed87f2-e2ac-d842-7093-185488ec313e.jpeg)


特に、重要な計算式は、視線と球面の交点(x,y,z)から地図上の緯度・経度を求める部分です。
以下の記号を使います。

+ $polerv$ : 角度$\phi$で回転された地軸
+ $polerh$ : polervに直交するベクトル
+ $q=(x,y,z)$ : 球面と視線の交点
+ $r$ : 球面の半径
+ $lat$ : 経度
+ $long$ : 緯度
+ $\phi$ : 地軸の回転
+ $\theta$ : 地軸"周り"の回転

計算順序は以下のとおりです。

(1) polerv, polerhを設定

```math
polerv = [0, cos(\phi), sin(\phi)] \\
polerh = [0, sin(\phi), -cos(\phi)] \\
```

(2) 緯度(lat)を求める

```math
c1  = 1/r * q \cdot polerv \\
lat = arccos(c1) \\
```

(3) ベクトル$q$の地軸polervとの直交成分の単位ベクトルを求める。

```math
b0 = 1/r * q - c1 * polrev \\
b = b0/||b0||
```

(4) 経度(long)を求める。

```math
cv = b \times polerh \\
c2 = cv \cdot polerv \\
c3 = b \cdot polerh \\
long = arccos(c3) sign(c2) + \theta
```

$cv$は$b$と$polerh$との外積です。$arccos$では$[0,\pi]$しか求められないので、外積の符号で補完しています。

# ソースコード
GitHubにアップロードしてあります。

[ okimebarun/VRMapEarth.jl ](https://github.com/okimebarun/VRMapEarth.jl)

特にこだわりの部分を解説します。

計算パラメータ用の構造体です。
Dictと付いていますが、辞書だと遅そうなので構造体に変えています。

```Julia
mutable struct ParamDict
    p::Float32  # z-axis of eye point
    r::Float32  # sphere radius
    v0::Int     # mid-point of v
    w0::Int     # mid-point of w
    m0::Int     # mid-point of m
    n0::Int     # mid-point of n
    z1::Float32 # z-axis of img-plane
    z2::Float32 # z-axis of background-plane (dont use now)
    msh::Int    # map size for longtitude
    msv::Int    # map size for latitude
end;
```

視線と球面との交点を求める関数です。パラメータ表示された直線$[at,bt,-t+p]$と球面の式$x^2+y^2+z^2=r^2$から導けます。

```Julia
function spxyz(d::ParamDict, a::Float32, b::Float32)
    ab1 = a^2+b^2+1
    r1 = d.p^2 -(d.p^2 - d.r^2)*ab1
    if r1 < 1e-10
        (false, 0.0, 0.0, 0.0)
    else
        t = (d.p - sqrt(r1)) / ab1
        x = a*t
        y = b*t
        z = -t + d.p
        (true, x, y, z)
    end
end
```

視線と球面との交点の(x,y,z)から緯度経度(lat,long)および、球中心から交点へのベクトルと視線ベクトルのなす角度のcosを、反射係数用に求めています。

```Julia
function polerax(d::ParamDict, x::Float32, y::Float32, z::Float32,
        long0::Int, lat0::Int)
    
    phi = pi*(Float32(lat0)/180.0+1.0)
    ax = 0
    ay = cos(phi)
    az = sin(phi)
    polerv = [ax, ay,  az] # earth's axis
    polerh = [ 0, az, -ay] # orthogonal to polerv
    
    # get latitude
    r = d.r
    c1 = [x,y,z]'*polerv/r
    lat = acos(c1*0.9999)

    # get longtitude
    # inline function
    cross3d(a,b)=[a[2]b[3]-a[3]b[2], a[3]b[1]-a[1]b[3], a[1]b[2]-a[2]b[1]]
    b  = [x,y,z]/r - polerv*c1 # projection on polerv
    b  = b/sqrt(b'*b) # unitilize
    cv = cross3d(b, polerh) # outer product
    c2 = cv'*polerv
    c3 = b' *polerh
    long = acos(c3*0.9999) * sign(c2)+ pi*Float32(long0)/180.0   
    
    # convert to 360 / 180 degree
    long2 = Float32(360.0*(long/(2pi) + 2.0) % 360.0) # Mod on float
    lat2  = Float32(180.0*( lat/(pi)  + 2.0) % 180.0)

    # reflection coefficient
    coef = [x, y, z]'*[0.0, 0.0, 1.0]/d.r
    (long2, lat2, coef)
end
```

画面座標から、以下の結果を計算する関数です。

+ 球面と視線が交わる場合、地図上の位置と反射係数
+ 球面と視線が交わらない場合、背景上の位置

関数を繋げて素直に計算できるのがいいですね。

```Julia
function convToMap(d::ParamDict, w::Int, v::Int, long0::Int, lat0::Int)
    (a,b) = imcoef(d, v, w)
    (bon, x, y, z) = spxyz(d, a, b)
    if bon
        (long, lat, coef) = polerax(d, x,y,z, long0, lat0)
        (mx, my) = mapxy(d, long, lat)
        (true, my, mx, coef)
    else
        (m, n) = backax(d, v, w)
        (false, n,  m, 0.0)
    end
end
```

画面座標系のピクセルpから地球画像もしくは背景画像の対応ピクセルを返す関数です。
RGB{N0f8}型を使うと大変簡単に書けます。

```Julia
    # convert to image
    function imconv1(p)
        # shade effect
        rgb0 = RGB{N0f8}(0.1,0.1,0.2)
        # select pixel from two images depending on the first flag
        p[1] ? RGB{N0f8}( img_earth[p[2],p[3]]*p[4] + rgb0*(1.0-p[4]) ) : img_galaxy[p[2],p[3]]
    end
```

スライダーの情報と画面座標系からカラーイメージを返す関数です。
ドット(.)を使うと、２次元配列計算がfor文なしで書ける！のがいいですね。

```Julia
    # visualize the result
    function view1(long::Int, lat::Int, z1::Float32, d1=d1, ws=ws, vs=vs)
        d1.z1 = z1
        imap1a(w,v)=convToMap(d1, w, v, long, lat)
        # Julia like mapping!
        m   = imap1a.( ws, vs' )
        img = imconv1.(m)
        colorview(RGB, img)
    end
```


# Juliaを使ってみた感想

個人的にJuliaを初めて使ってみた感想です。

#### 良かった点

+ 数学の記述っぽく素直にプログラミングできる。
+ 関数がメタっているので再利用しやすい。C++のtemplateみたいな感じか。
+ atan(Inf)などが普通に計算できる！ゼロ割気にしないでいい。
+ ドット（例：func.(x)）が凄く使いやすい。
+ ライブラリがそれなりに充実している。

#### 悪かった点

+ using SomePackageするたびにコンパイルが始まる…。遅いです。キャッシュして(-_-;)。
+ Interact.jlはJupyter NotebookではOKだったのに、Jupyter Labで使えなかった。（そのうち改善？）
+ 関数によっては型が対応してないので、いちいちキャストが必要。
+ 謎の丸め誤差が生じる場合がある。Float64使えば良かったのかな？

#### 興味が湧いた点

+ 高速化。分散計算やGPUなど、色々。
+ シミュレーション。数式が書きやすいので色々試しやすそう。
+ UIアプリ化。Interact.jlはシンプルで使いやすい。

# 参考リンク
下記のページを参考にさせて頂きました。

[Juliaの画像処理ライブラリを使ってわかったこと](https://logmi.jp/tech/articles/320153)
[Interact.jl を用いて Julia で使う対話機能を持つ可視化ソフトを作ろう](https://qiita.com/SatoshiTerasaki/items/ac48b2a2cb33cd476446)
[実例で学ぶJuliaプログラミング言語入門](http://itref.fc2web.com/lang/julia/)
[1から始める Juliaプログラミング](https://www.coronasha.co.jp/np/isbn/9784339029055/)
[Equirectangular projection](https://en.wikipedia.org/wiki/Equirectangular_projection)
[高品質なフリー画像素材](https://pixabay.com/ja/)