Hatena::ブログ(Diary)

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

Prima Project

2017-10-31

Docker でmeshlabserver を使おうと思ったらcannot connect to X server とエラーが出る

Docker で仮想環境っぽいところだったり、サーバーなどでGUI、グラフィック系のドライバーがなさそうなところでmeshlab とかそういうのを実行すると、上記のエラーがでる。

wkhtmltopdfをUbuntuで使うためのメモ | Network Project 5

を参考に、xvfb-run を経由して実行すればできるようだ。

xvfb-run meshlabserver

2017-08-14

2017-05-24

MikuHatsune2017-05-24

手動で細胞の動きを追跡したい

細胞でもなんでもいいが、視野内の物体の動きを軌跡として取得したい。

最近流行りのディープラーニングでは下調べが足りないがいろいろあるっぽい。しかし、これらの手法を使うにしても自前で開発するにしても、正解データと比較して性能評価が必要である。これをGround truth という。

共同研究者からjava ベースで動く軽快なプログラムをもらったが、java のコンパイルがうまくいかないのでR で作ってみた。

と思ったらjava が動いたのでこれはお蔵入り。動画データはこちらからパクった。

f:id:MikuHatsune:20170524192744g:image

 

R でこれをやろうと思ったら、画像を表示して、対象の細胞をクリックして挫傷を取得し、1枚画像を送って、最後までやる、を繰り返すことになる。

R では画像をポチポチして座標を取得するのに、locator 関数が使える。

これで細胞の座標を取得できるが、ひだりクリックで座標を取得し、みぎクリックは空打ちになるようなので、このクリック操作により条件分岐ができる。

 

追加したい機能としては、いままでにクリックした細胞はチェックしてわかりやすくすることと、失敗したらやり直すことである。

いままでにクリックした細胞は、履歴を取っておいて画像に重ね書きすることでできる。

(小さいがバツ印がついている)

f:id:MikuHatsune:20170524193023p:image

失敗したらやり直すのはちょっとめんどうだったので、画像をすべてのタイムフレームでクリックしたあとに気に入らなかったら全部捨てる仕様になっている。

f:id:MikuHatsune:20170524192849p:image

 

出力は細胞ごとにリストで、タイムフレームごとでのピクセル単位の(重心)座標が返ってくる。集中力が続く限り連続してやれる。

f:id:MikuHatsune:20170524193147p:image

[[1]]
      [,1]     [,2]       [,3]
 [1,]    1 64.97317   6.472511
 [2,]    2 63.64829  19.304606
 [3,]    3 59.67362  25.343239
 [4,]    4 64.02683  38.552749
 [5,]    5 65.54098  53.083210
 [6,]    6 61.18778  68.934621
 [7,]    7 61.18778  81.011887
 [8,]    8 53.04918  92.145617
 [9,]    9 53.23845 101.203566
[10,]   10 53.80626 108.940565
[11,]   11 53.80626 118.375929
[12,]   12       NA         NA
[13,]   13       NA         NA
[14,]   14       NA         NA

[[2]]
      [,1]     [,2]      [,3]
 [1,]    1 65.35171  6.283804
 [2,]    2 64.02683 18.738484
 [3,]    3 59.10581 24.965825
 [4,]    4 46.61401 34.401189
 [5,]    5 42.82861 49.497771
 [6,]    6 44.15350 61.386330
 [7,]    7 43.01788 77.237741
 [8,]    8       NA        NA
 [9,]    9       NA        NA
[10,]   10       NA        NA
[11,]   11       NA        NA
[12,]   12       NA        NA
[13,]   13       NA        NA
[14,]   14       NA        NA

この出力を再利用すれば、途中から作業を再開することができる。ただし、履歴が増えすぎると重い。

 

あとはバッチ処理にしてしまおうと思ったが、バッチ処理ではplot したときにプロットGUI が出てこないので諦めている。

Ground truth が必要といったが、この作業をひたすらやるのも諦めている。

###
#
# Rscript --slave --vanilla trackerR.R args0 args1 args2
#
###

library(EBImage)
library(png)
library(jpeg)
args <- commandArgs(TRUE)
if(length(args) < 2){
  cat(
"command line args are
args[0]: image sequence folder. format is currently .png
args[1]: output filename. space separated txt. ID,t,x,y
args[2]: optional. filepath of temporaly result of tracking. you can restart from this file.\n"
      )
  do <- FALSE
} else if (length(args) > 3){
  cat(
"args is less than 4.\n"
)
  do <- FALSE
} else {
  do <- TRUE
}

input_dir <- args[1]
out_name  <- args[2]
if(length(args) == 3){
  tmp_res <- read.table(args[3], header=TRUE)
  final <- split(tmp_res[,-1], tmp_res$ID)
} else {
  final <- list()
}

txt0 <- "Do you continue next tracking ?\nYes: Left click\nNo: Right click"
txt1 <- "Do you accept your tracking ?\nYes: Left click\nNo: Right click"

imgs <- list.files(input_dir, full.names=TRUE)
ans <- TRUE
accept <- FALSE
if(do){
  while(ans){
    while(!accept){
      res <- matrix(0, length(imgs), 3)
      for(i in seq(imgs)){
        f <- readImage(imgs[i])
        xleft <- 1
        ybottom <- nrow(f)
        xright <- ncol(f)
        ytop <- 1
        par(mar=rep(0, 4))
        image(t(0), xlim=c(xleft, xright), ylim=c(ybottom, ytop), col="white")
        rasterImage(f, xleft, ybottom, xright, ytop)
        legend("topright", legend=paste0("t = ", i), text.col="yellow", cex=1.5)
        if(length(final) > 0){
          for(j in seq(final)){
            points(final[[j]][i, 2], final[[j]][i, 3], col="yellow", pch=4, cex=1, font=2)
          }
        }
        xy <- unlist(locator(1))
        res[i, ] <- switch(is.numeric(xy)+1, c(i, NA, NA), c(i, xy))
      }
      text((xleft+xright)/2, (ybottom+ytop)/2, txt1, col="yellow", pos=3, cex=2, font=2)
      accept <- is.list(locator(1))
    }
    final <- c(final, list(res))
    image(t(0), xlim=c(xleft, xright), ylim=c(ybottom, ytop), col="white")
    rasterImage(f, xleft, ybottom, xright, ytop)
    text((xleft+xright)/2, (ybottom+ytop)/2, txt0, col="yellow", pos=3, cex=2, font=2)
    ans <- is.list(locator(1))
    accept <- FALSE
  }
  output <- cbind(rep(seq(final), each=length(imgs)), do.call(rbind, final))
  colnames(output) <- c("ID", "t", "x", "y")
  write.table(output, out_name, quote=FALSE, row.names=FALSE)
}

2017-05-02

MikuHatsune2017-05-02

勾配法による最適化

ある関数があって、最小化もしくは最大化したい。これらは正負を入れ替えればよいでの最小化を考える。

 

いま、適当にf(x,y)=x(x-1)(x-2)+y(y-2)(y-3)(y+7)-8 という関数があったとする。これの最小値とそのときの(x,y)の組を求めたい。

2変数なので3次元に図示するとこんな感じである。最小値はひとつ定まりそう。真っ赤なところがお盆の低いところになる。

f:id:MikuHatsune:20170502171100p:image

x <- y <- seq(-10, 10, length=300)
f <- function(x, y) x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8
xy <- expand.grid(x, y)
z <- f(xy[,1], xy[,2])
i <- 100
cols <- rainbow(i)[cut(z, quantile(z, seq(0, 1, length=i-1)))]
plot3d(xy[,1], xy[,2], z, col=cols, xlab="x", ylab="y", zlab="z")

ちなみにPython の3D プロットは機嫌が悪いようで使えなかった。イラッ☆

離散的に解いた最小値は

xy[which.min(z),]; min(z)
           Var1      Var2
22864 -5.785953 -4.916388
[1] -1245.717

 

その1:解析的に解く

関数が決まっていて、微分が出来て、微分=0 が解析的に解けるならばそれで解いたほうが精確であるし、それが最小値であることが保証される。

高校数学までで言えば、f(x)x の関数のとき、x で微分して導関数が0となるときのx が得られて、増減表から凹になる部分を探す、というのがセオリー。

大学数学でいえば、パラメータがいくつかあるときに、各パラメータで微分(偏微分)して偏微分たちの連立方程式が0 となるパラメータセットを考えればよい。

今回のf(x,y) は微分できそうな感じで用意したので、微分しよう。微分にはsympy を用いる。

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
x, y = symbols("x y")
f = x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8

parms = [x, y]                           # 変数
deriv = map(lambda z: diff(f, z), parms) # 各変数での偏微分
sol = solve(deriv, parms)                # 偏微分したものを連立方程式で解く

z = map(lambda z: f.subs({x:z[0], y:z[1]}).evalf(), sol)
z = np.array(z, dtype=np.complex)

xy = np.array(sol, dtype=np.complex)

sol はsympy で解いた偏微分方程式のオブジェクトで、ところどころに虚数I が入り乱れている。sympy のオブジェクトにはsubs 関数で、辞書形式{x:x0, y:y0} と入れれば、代入したことになる。このとき、sympy 自体がまだ関数であればdoit() をすることで計算され、計算した値が欲しければevalf() で数値を得られる。

sympy 形式でI で表示されている場合、numpy に戻すにはdtype=np.complex を指定する。

すると、

z
array([   11.74548233 +6.77626358e-21j,   -10.78588513 +1.33638236e-51j,
        -565.66379583 +3.34095589e-49j,     4.81250937 +6.77626358e-21j,
         -17.71885808 -8.45930031e-49j,  -572.59676879 -5.13170824e-49j,
        -668.34781483 +6.77626358e-21j,  -690.87918229 +4.37999317e-47j,
       -1245.75709299 +4.41326909e-47j])
xy
array([[ 0.43596395 +2.71050543e-20j,  0.83707381 +1.05879118e-22j],
       [ 0.43596395 +2.71050543e-20j,  2.56096376 -5.29395592e-23j],
       [ 0.43596395 +2.71050543e-20j, -4.89803757 +3.38813179e-21j],
       [ 1.58881707 +5.29395592e-23j,  0.83707381 +1.05879118e-22j],
       [ 1.58881707 +5.29395592e-23j,  2.56096376 -5.29395592e-23j],
       [ 1.58881707 +5.29395592e-23j, -4.89803757 +3.38813179e-21j],
       [-5.77478101 +2.64697796e-23j,  0.83707381 +1.05879118e-22j],
       [-5.77478101 +2.64697796e-23j,  2.56096376 -5.29395592e-23j],
       [-5.77478101 +2.64697796e-23j, -4.89803757 +3.38813179e-21j]])

となり、z をみると最小値-1245 というのがあって、それに対応するxy は(-5.775, -4.898) であり、悪くない。

 

その2:勾配法(急速降下法)

関数の式は立式できても、微分がなかなか解けないという状況はある。そのとき、「変数を少し動かして、もともと最小化したかった関数がより小さくなるようにする」というのを計算機的にゴリ押しする。

アルゴリズム的に書けば、i 番目のパラメータをi+1 番目にいい感じに更新するときの具合が、関数から定まる勾配によって支配され、

w^{(i+1)} ¥leftarrow w^{(i)}-¥eta¥nabla E(w^{(i)})

と書ける。ここkで、¥nabla が(偏微分から求まる)勾配、w がパラメータセット、¥eta が勾配に応じてどれくらい動かすかの学習率である。

学習率はハイパーパラメータなので、事前にいい感じに設定してやる。アルゴリズムが進化すればこのハイパーパラメータすらも勝手に決まっていく。

¥eta=0.005 とし、(0,0) から始めると、10回くらいの繰り返しで最小値に行っている。この手法は遅い法で、しかも局所解に落ち込むこともある。

解析的に微分できない場合は、微分の定義に戻って、f(x,y) から微小距離¥Delta だけ動かして、¥Delta ¥to ¥infty の極限をとるのを模倣して、¥Delta=10^{-8} くらいで近似すればいい。この方法は、解析的に微分できる場合でも、2通りで試して正解だったらちゃんとできている、というのがわかるので、やっておくとよいらしい。

ここでは、f(x,y)¥to f(x+¥Delta, y+¥Delta) の微分ではなく、演算誤差の精度の都合上、本当に見たい点(x,y) を挟んだf(x-¥Delta,y-¥Delta)¥to f(x+¥Delta, y+¥Delta) の差分を見ている。なので割り算するときも2¥Delta である。

# 最小化する関数
def J(x, y):
  return x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8

# 勾配
# 偏微分で与えられるならそれでもよい
def gradient_descent(x, y, h=10e-8):
  dx = (J(x+h, y) - J(x-h, y)) / (2*h)
  dy = (J(x, y+h) - J(x, y-h)) / (2*h)
  return np.array([dx, dy])

xy = [0, 0]
np.array(map(lambda z: z.subs({x:xy[0], y:xy[1]}).evalf(), deriv))

def GD(init, eta=0.005, epoch=100):
  res = np.zeros((epoch, len(init)))
  res[0] = init
  cost = np.zeros(epoch)
  cost[0] = J(init[0], init[1])
  for i in range(1, epoch):
    init -= eta * gradient_descent(init[0], init[1])
    res[i] = init
    cost[i] = J(init[0], init[1])
  return res, cost

s1 = GD([0, 0])

plt.subplot(121)
plt.plot(s1[0][:,0], s1[0][:,1])
plt.xlim(-10, 1)
plt.ylim(-10, 1)
plt.title("x and y coordinate")
plt.subplot(122)
plt.plot(s1[1])
plt.title("Value to be minimized")
plt.show()

f:id:MikuHatsune:20170502171136p:image

 

その3:共役勾配法

共役なベクトル列がうまく選べると、急速降下法はその1点の勾配の向きまかせに動いていたが、ちゃんと最適化方向に進行方向が決まるらしい。しかもどれくらい行くべきかの学習率も勝手に決まる。実装力の低いPython novice はライブラリを使う。

scipy にoptimize という一連の最小化最適化ライブラリがあるので、いい感じに使う。書き方はここからパクる。

最小化したい関数J は、パラメータセット¥theta の関数にしておくのが簡単である。他に引数が欲しければ、*args と書いておいて関数内で各変数にばらす(後述)。

# conjugate gradient
# 急速降下法と引数の書き方を少し変える
def J(theta):
  x, y = theta
  return x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8

def gradient_descent(theta):
  h = 10e-8
  grad = np.zeros(len(theta))
  di = np.diag(np.repeat(h, len(theta)))
  for i in range(len(theta)):
    grad[i] = (J(theta+di[i]) - J(theta-di[i]))/(2*h)
  return grad

init_theta = [10, 10]
res = optimize.fmin_cg(J, init_theta, fprime=gradient_descent, full_output=True)
Optimization terminated successfully.
         Current function value: -1245.757093
         Iterations: 9
         Function evaluations: 29
         Gradient evaluations: 29

この程度の簡単さなら実行速度の恩恵は実感できないが、高速らしい。結果のres には最終的な最小値と、その時のxy の組が入っている。

disp=False にすればコンソールに出力しない、full_output=True にすればすべての出力結果を取得できる。

R ではRcgmin パッケージでできて

library(Rcgmin)
fr <- function(xy){
  x <- xy[1]
  y <- xy[2]
  z <- x*(x-2)*(x-1)*(x+8) + y*(y-2)*(y-3)*(y+7) -8
  return(z)
  }

gr <- function(xy){
  h <- 10e-8
  grad <- rep(0, length(xy))
  for(i in seq(xy)){
    h1 <- replace(rep(0, length(xy)), i, h)
    grad[i] <- (fr(xy+h1) - fr(xy-h1)) / (2*h)
  }
  return(grad)
}

Rcgmin(fn=fr, gr=gr, par=c(10, 10))
$par
[1] -5.774781 -4.898038

$value
[1] -1245.757

$counts
[1] 875 129

$convergence
[1] 1

$message
[1] "Too many function evaluations (> 866) "

 

ようやくここまで来たので、球面調和関数の最適な回転を見つけたい。

球面調和関数の係数は、もとの物体を回転させなくても、係数を回転させればよい。この回転は

c_l^m=¥sum_{m^{’}=-l}^l c_l^{m^{’}}D_{m^{’}m}^l

で得られる。これはsympy ではWigner D 行列としてある。

回転させるパラメータ空間は(¥alpha,¥beta,¥gamma) であり、これらは¥[0,2¥pi¥] の間で連続に変化するはずである。なので、最小値(と最大値)がひとつ定まることは保証されている、はず。ただし局所解であることは多々あると思う。

ここで、最適化関数J には、パラメータ¥theta=(¥alpha,¥beta,¥gamma) を入れ、L^2 ノルムの計算に必要な、球面調和関数係数の各項を*args から入力することにしている。

from sympy import *
from sympy import pi
from sympy.physics.quantum.spin import Rotation
from sympy.functions import re
from scipy import optimize
import numpy as np
import numpy.linalg as LA
import math
from scipy.special import sph_harm as _sph_harm
import pyshtools as p
import matplotlib.pyplot as plt

l, m, n, a, b, g = symbols("l m n a b g") # mp is n
WigD = Rotation.D(l, m, n, a, b, g)
v1 = [10, 1, 2, 0, 1, 5, 3, 4, 0]
r = flatten([[sum(transpose(Matrix(v1[(k**2):((k+1)**2)])) * Matrix([WigD.subs({l:k,m:j,n:n0}) for n0 in range(-k, k+1)])) for j in range(-k, k+1)] for k in range(Lmax+1)])
r = map(lambda z: z.subs({a:abg[0],b:abg[1],g:abg[2]}).doit().evalf(), r)
r = np.array(r, dtype=np.complex) # 完成品

でもこれ、shtools.rotate と結果が合わない…

diff で微分もできるけれども、計算がものすごく遅い。

というわけで、数値微分によって勾配を求めたものとする。

def tri(clm_coef):
  vec = clm_coef.coeffs
  return np.concatenate([np.fliplr(vec[1])[:,:-1], vec[0]], axis=1)

def invariantL2vec(vec):
  return np.array([np.linalg.norm(vec[range(k**2, (k+1)**2)]) for k in range(Lmax+1)])

def clm2vec(clm, Lmax=None):
  if Lmax == None:
    Lmax = clm.lmax
  return np.concatenate([tri(clm)[k][(Lmax-k):(Lmax+k+1)] for k in range(Lmax+1)])

lmax = 2
x = p.SHCoeffs.from_zeros(lmax)
ls = [0, 1, 1, 1, 2, 2, 2, 2, 2]
ms = [0, -1, 0, 1, -2, -1, 0, 1, 2]
value = [10, 1, 2, 0, 1, 5, 3, 4, 0]

for (l, m, v) in zip(ls, ms, value):
  x.set_coeffs(v, l, m)

PI = pi.evalf()
a, b, g = PI/2, PI/3, PI/4
x0 = x.rotate(a, b, g, degrees=False)

def J(theta, *args):
  clm0, clm1 = args
  #return LA.norm(clm2vec(clm0) - clm2vec(clm1.rotate(theta[0], theta[1], theta[2], degrees=False)))
  return LA.norm((clm0 - clm1.rotate(theta[0], theta[1], theta[2], degrees=False)).coeffs)

# partial derivatives
def gradient_descent(theta, *args):
  clm1, clm2 = args
  h = 10e-8
  grad = np.zeros(3)
  di = np.diag(np.repeat(h, 3))
  for i in range(3):
    grad[i] = (J(theta+di[i], clm1, clm2) - J(theta-di[i], clm1, clm2))/(2*h)
  return grad

init_theta = np.array([1, 1, 1])
optimize.fmin_cg(J, init_theta, fprime=gradient_descent, args=(x0, x))

Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 2.529874
         Iterations: 2
         Function evaluations: 17
         Gradient evaluations: 5
array([ 4.49933254, -6.11515311,  3.59954578])

初期値が(1,1,1) だと最小値が0にならなかったようである。いま、球面調和関数の係数x(¥frac{¥pi}{2},¥frac{¥pi}{3},¥frac{¥pi}{4}) だけ回転させたものがx_0 であり、最小値は0となって一致しないとおかしい。

というわけで変な局所解にハマったようである。

本当ならばこういうときに探索する領域をいい感じに飛ばす手法があるっぽいが、適当に初期値を振ることにすると

optimize.fmin_cg(J, np.random.rand(3)*2*pi.evalf(), fprime=gradient_descent, args=(x0, x))
Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 0.000000
         Iterations: 37
         Function evaluations: 94
         Gradient evaluations: 81
array([1.57079627211901, 1.04719757847306, 0.785398109529357], dtype=object)

となって、確かに(¥frac{¥pi}{2},¥frac{¥pi}{3},¥frac{¥pi}{4}) が推定されている。

ちなみに、x_0 を回してx に一致させるには、(¥alpha,¥beta,¥gamma) で回した時には(-¥gamma, -¥beta,-¥alpha) で逆変換すると戻る。

 

勾配法と最適化についてはここらあたり。

パターン認識と機械学習 上

パターン認識と機械学習 上

これは確率的勾配法やモーメント法、AdaGrad がどのような最適化経路をたどるかを可視化しているのでわかりやすい。

 

2017-04-05

MikuHatsune2017-04-05

MeshLab を使って曲率を取得するのをバッチ処理して自動化する

3次元物体の曲率を以前やっていたのだが、自前のスクリプトに重大なミスが発覚したので、プログラミング力の高い人のプログラムから借用することを考える。

離散的な3D オブジェクトから、平均曲率とガウス曲率を計算したいが、意外とない。

と思ったらMeshLab にあるようである。

.obj を読みこんで、Filters > Color creation and processing > Discrete curvature と選択すると、mean かgaussian か選べと出るので、どちらか選んでapply すると、こんな感じになる。

色は Filter > Quality measure and computations > Quality mapper applier から、カラースケールを適当にいじる。

Render > Show vert quality histogram を選ぶと、今回は頂点について曲率を計算しているので、頂点の値に応じたヒストグラムが出てくる。これはGUI では背景が濃い青だったので白抜きの文字が本当は見えている。

f:id:MikuHatsune:20170405201643p:image

 

ここで、この値情報を保持したままファイルに書き出してその値が欲しいのだが、値の情報は .ply 形式しか保持してくれないようである。なので、ファイルの保存は .ply にする。

ここで、.ply にはbinary 形式とascii 形式があるので、export のときにadditional parameters のbinary encoding のチェックをはずさないと、その後のデータ処理で困る。

 

というのがGUI での一連の操作である。ここで、この処理をしないといけないファイルが複数あるとき、for loop に入れたい。

この時に使えるのがmeshlabserver コマンドである。

まず、このGUI での操作を、Filters > Show current filter script に行って、マクロ化したスクリプト .mlx ファイルを手にいれる。ここでは、Filter タブでの操作(のみ、たぶん)が記録される。

こうして、GUI 操作がマクロ化されたので、

meshlabserver -i input.obj -o output.ply -s procedure.mlx -om vq > /dev/null 2>&1

とする。入力はなんでもよく、出力に .ply を指定すれば勝手に .ply 形式で出力してくれる。-s でやりたい操作.mlx を指定して、今回は頂点の曲率の値が欲しいので、-vq で vertices quality 値を .ply ファイルに書き出しておく。

ここで、meshlabserver の標準出力がどうしても邪魔なので、/dev/null 2>&1 であの世へ送る。

 

この場合に問題になるのが、-o output.ply がどうしてもbinary 形式の .ply ファイルしか書き出さないようで、けっこうなmeshlab ユーザーが困っているようである。というわけで、binary 形式で出力した .ply をascii 形式の .ply に変換したいのだが、どうもいいのがないようなのでR で頑張る。

R のRvcg パッケージに .ply ファイルをbinary 形式でも読める関数があるので、これで .ply ファイルを読み込み、目的のquality 値だけ取り出して出力できるようにバッチ処理にしてみる。

# quality.R
library(Rvcg)
args <- commandArgs(TRUE)
ply <- vcgPlyRead(args)
cat(sprintf("%f", ply$quality))
cat("\n")

これをバッチ処理すると、コンソールにquality 値がずらずらと出力されるので、適当にファイルに保存する。

Rscript --slave --vanilla quality.R $ply.ply > `printf %04d $ply`.tmp

R でバッチ処理したものを出力すると、

[1] 1

というのが残るので、cat して改行コードを付け足しておく。

 

とすればなんとか取得できる。

f:id:MikuHatsune:20170405204033p:image