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

Docker で仮想環境っぽいところだったり、サーバーなどでGUI、グラフィック系のドライバーがなさそうなところでmeshlab とかそういうのを実行すると、上記のエラーがでる。
wkhtmltopdfをUbuntuで使うためのメモ | Network Project 5
を参考に、xvfb-run を経由して実行すればできるようだ。

xvfb-run meshlabserver

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

MikuHatsune2017-05-24

細胞でもなんでもいいが、視野内の物体の動きを軌跡として取得したい。
最近流行りのディープラーニングでは下調べが足りないがいろいろあるっぽい。しかし、これらの手法を使うにしても自前で開発するにしても、正解データと比較して性能評価が必要である。これをGround truth という。
共同研究者からjava ベースで動く軽快なプログラムをもらったが、javaコンパイルがうまくいかないのでR で作ってみた。
と思ったらjava が動いたのでこれはお蔵入り。動画データはこちらからパクった。

 
R でこれをやろうと思ったら、画像を表示して、対象の細胞をクリックして挫傷を取得し、1枚画像を送って、最後までやる、を繰り返すことになる。
R では画像をポチポチして座標を取得するのに、locator 関数が使える。
これで細胞の座標を取得できるが、ひだりクリックで座標を取得し、みぎクリックは空打ちになるようなので、このクリック操作により条件分岐ができる。
 
追加したい機能としては、いままでにクリックした細胞はチェックしてわかりやすくすることと、失敗したらやり直すことである。
いままでにクリックした細胞は、履歴を取っておいて画像に重ね書きすることでできる。
(小さいがバツ印がついている)

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

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

[[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)
}

勾配法による最適化

MikuHatsune2017-05-02

ある関数があって、最小化もしくは最大化したい。これらは正負を入れ替えればよいでの最小化を考える。
 
いま、適当にf(x,y)=x(x-1)(x-2)+y(y-2)(y-3)(y+7)-8 という関数があったとする。これの最小値とそのときの(x,y)の組を求めたい。
2変数なので3次元に図示するとこんな感じである。最小値はひとつ定まりそう。真っ赤なところがお盆の低いところになる。

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()


 
その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 がどのような最適化経路をたどるかを可視化しているのでわかりやすい。
 

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

MikuHatsune2017-04-05

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 では背景が濃い青だったので白抜きの文字が本当は見えている。

 
ここで、この値情報を保持したままファイルに書き出してその値が欲しいのだが、値の情報は .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 して改行コードを付け足しておく。
 
とすればなんとか取得できる。

Python でコマンドライン引数、標準入力を受け取る

結論から言えばsys.argv を使う。
http://www.yukun.info/blog/2008/07/python-command-line-arguments.html
【python】sys.argv で引数をとる - metabo346の日記
コマンドライン引数入力できるPythonプログラムをつくる - 剥いだ布団とメモランダム.old
Python を使っている。
R でもそうだが、コンパイルしないタイプの、対話型言語ではターミナルCUI (かwindows 的にはGUIとかIPython とか? よく知らないが)を使ってガチャガチャ書いて、実行ごとに演算結果を取得している。
大量の処理したいファイルがあって、いわゆるシェル的に

execute hogehoge < huga.txt

みたいにプログラムに投げたいという状況がある。このとき、python で実行できる形式に.py ファイルを書いてバッチ処理できるようにすればよい。
ここで、バッチ処理するには

python script.py

とすれば実行できるが、問題はscript.py がなんらかのファイルを受け取る必要がある場合に、上の1行にファイルパスを記述すればどうにかなるようにしたい。
 
その昔、コンピュータに全く詳しくなかったとき(今でも全然詳しくない)、コンソールにポチポチ打つことを標準入力ということを知った。
Python2で競技プログラミングする時に知っておきたいtips(入出力編) - Qiita
新約:標準入力の受け取りその1(Python3対応版) - Qiita
こうすると、input_raw (python3 ならinput)のところで、コンソールが入力受付状態になり、必要な文字を入れるとそれをpython に流すことができる。
というわけでこれを使って書いていたが、こうすると、シェルスクリプトで書いたときに、ここの部分で入力受付状態のままになって永遠に進まないようである(たぶんエンターを入力するみたいなことを実行すればいいかもしれないが面倒)。
また、この方法だと、入力受付状態のときに、tab でのパス自動補完が利かないようである。
 
というわけで、

python script.py hoge0 hoge1 hoge2 ...

というように引数をとるにはどうしたらよいか、という話で、sys.argv を使うとよいらしい。これはtab 補完も有効である。
 
上のように入力した場合、sys.argv[0] はscript.py になってしまうようなので[1:] から取ることと、ファイルパスを ~/hoge で入力したとき、os.path.expanduser をしないと /home/user/ を補完してくれないので適当につけておく。
 

import sys
import os
input_args = sys.argv[1:]
path = os.path.expanduser(...)

ボクセルデータから形を構成する

MikuHatsune2017-03-18

蛍光イメージング動画があって、そこから細胞の3次元の形を再構成したい。
手順としては、
蛍光イメージの取得:他力本願
2Dでのセグメンテーション:他力本願、大津法
セグメンテーションされた領域の抽出:openCV
セグメンテーションされた領域内で欠けている部分を埋める:openCV
z方向の何枚かを並べて3次元構成する:マーチングキューブ
という感じでやる。
 
蛍光染色した細胞はこんな感じで観察できる。

すべてのz方向では14枚の画像がある。大津法でのセグメンテーション結果は以下の通りである。
zz-007.png という連番で保存されているとする。














 
あるディレクトリに連番で画像が入っているものとして、z軸のID と時間軸のID をいい感じで処理しておく。

import numpy as np
import os
import re
import cv2
import pylab as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from skimage import measure
from skimage.draw import ellipsoid

files = os.listdir("./")
files.sort()
zid = map(lambda x: int(re.sub("-.+", "", x)), files)
zmin = min(zid)
zmax = max(zid)
tid = map(lambda x: re.sub(".+-|.png", "", x), files)

ti = 7
png = map(lambda z: "{0:02d}".format(z)+"-""{0:03d}".format(ti)+".png", range(zmin, zmax+1))
xpix, ypix, zpix = 1.24, 1.24, 3
rgb = 256 # white color

画像データをopenCVPython binding であるcv2 のimread で読み込む。imread では行列データとしてx方向のピクセル番地、y方向のピクセル番地にRBG なら3チャネル分の画素値、モノクロなら白の強度が入っている。
輪郭の抽出はfindContours ででき、contours に収められている。
contours には輪郭を構成するピクセルの番地(x,y) が入っている。輪郭の抽出の仕方には何パターンかあるが、細胞は穴があいてなく、もっとも外側の輪郭を抽出すればよい、と仮定するとRETR_EXTERNAL を指定する。また、細胞の辺縁はかなりギザギザで、近似することなく持っていてもデータを減らして近似してもよい。
 
大津法のセグメンテーションの都合で、細胞の内部領域であっても白塗りにならない部分がある。これをそのまま残しておくと次のマーチングキューブなどの3次元構成のときにノイズになりまくるので、細胞はトーラスだったり内部に空洞を持たないと仮定して、fillPoly 関数でとりあえず輪郭の内部を塗りつぶす。
contour_extract は読み込んだセグメンテーション済データの輪郭となるピクセルを抽出し、image3D 関数は輪郭とともに輪郭内部のピクセルも抽出する。引数は画像のパスで、出力はN*3 の行列である。XY の大きさの2次元行列がN*3 の行列になったので、各ピクセルのRGB 情報はセグメンテーションされて0/1 になって出力される。

def contour_extract(f):
  z = int(re.sub("-.+", "", f))
  xyz = np.array([])
  im0 = cv2.imread(f)
  im1 = cv2.cvtColor(im0, cv2.COLOR_BGR2GRAY)
  ret, thres = cv2.threshold(im1, 13, 255, cv2.THRESH_BINARY)
  contours, hierarchy = cv2.findContours(im1, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
  cv2.drawContours(im0, contours, -1, (0, rgb, 0), 1)
  a = np.sum(im0, 2)
  if z not in [zmin, zmax]:
    idx = np.where(a == rgb)
  else:
    idx = np.where(a > 0)
  xyz0 = np.c_[np.c_[idx]*np.array([xpix, ypix]), np.array([z]*len(idx[0]))*zpix]
  return xyz0

def image3D(f, fill=False):
  z = int(re.sub("-.+", "", f))
  xyz = np.array([])
  im0 = cv2.imread(f)
  im1 = cv2.cvtColor(im0, cv2.COLOR_BGR2GRAY)
  ret, thres = cv2.threshold(im1, 13, 255, cv2.THRESH_BINARY)
  contours, hierarchy = cv2.findContours(im1, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
  if fill:
    cv2.fillPoly(im0, contours, color=(rgb, rgb, rgb))
  cv2.drawContours(im0, contours, -1, (0, rgb, 0), 1)
  a = np.sum(im0, 2)
  idx = np.where(a > 0)
  xyz0 = np.c_[np.c_[idx]*np.array([xpix, ypix]), np.array([z]*len(idx[0]))*zpix]
  return xyz0


# fillPoly して読み込む
def imread_fill(filename):
  im0 = cv2.imread(filename)
  im1 = cv2.cvtColor(im0, cv2.COLOR_BGR2GRAY)
  #im2 = cv2.GaussianBlur(im1, (15, 15), 0)
  contours, hierarchy = cv2.findContours(im1, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
  cv2.fillPoly(im0, contours, color=(rgb, rgb, rgb))
  return im0

 
とりあえずプロットしてみる。

xyz = np.vstack((map(contour_extract, png)))
ori = np.vstack((map(lambda z: image3D(z, fill=True),  png)))
X, Y, Z = xyz.T
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(ori.T[0], ori.T[1], ori.T[2], color="green")
ax.scatter3D(X, Y, Z, color="blue")
ax.view_init(elev=15, azim=-50)
plt.show()

 
各zスライスのxy方向に細胞内部が存在するかしないかのデータへ変換できたので、これを3次元構成する。3次元構成にはマーチングキューブを用いる。
マーチングキューブ法はxyz 格子の各点にある値v が乗っていて、その値が閾値Iを超えるか超えないかで0 もしくは1 に変換する。この閾値が境界となって、3次元物体の表面を構成する。閾値以下(以上)の領域は等しくなるので、その境界はisosurface と言われる。
図では赤もしくは青がある閾値で隔てられ、isosurface は紫である。この紫は、赤と青のノードをつなぐエッジを適当に内分するように線が引かれる。

これを3次元xyz 格子でやると、赤と青のノードパターンがある立方体(xyz 格子がその立方体の頂点になる)の頂点で[2^8=256] パターンになる。これをすべて列挙してもいいが、回転と反転を考えると15パターンになるので、これをひたすら頑張る。

オリジナルの実装はここにあるっぽいが、R ならmisc3d パッケージ、Python ならmcubes(これとか)やskimage に入っている。今回はskimage を使ってみる。
細胞の形としては閉じてほしいから、仮にz軸の一番上(下)が細胞をぶった切っている断面だとしても、フタ(すべてが0、つまり黒い画像)をすることで閉じるようにいじってマーチングキューブを行う。
ここで、格子にはセグメンテーションの閾値で1もしくは0 になっているから、isofurface は0から1 の間の適当な数値、ここでは0.5 にしている。

 
これで3次元構成は終了。

v0 = np.dstack(map(lambda x: np.sum(imread_fill(x), axis=2), png))
v = np.zeros(tuple(np.array(np.shape(v0))+2))
idx = np.where(v0>-1)
v[ map(lambda x: x+1, idx) ] = v0[idx]

# Marchingcubes を行う
hoge = measure.marching_cubes(v, 0.5)

xyz = np.vstack((map(image3D, png)))
ori = hoge[0]
mesh = Poly3DCollection(hoge[0][hoge[1]])
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(ori.T[0], ori.T[1], ori.T[2], color="green")
#ax.scatter3D(xyz.T[0], xyz.T[1], xyz.T[2]-1, color="red")
mesh.set_edgecolor('k')
ax.add_collection3d(mesh)
plt.show()

# obj ファイルとして書き出し
with open("out.obj", "w") as w0:
  for i, j in enumerate(hoge):
    for k in j:
      line = ["v", "f"][i] + " " + " ".join(map(str, k + i)) + "\n"
      w0.write(line)

 
と思いきや、初期のマーチングキューブの実装は、簡便で高速ながら、穴(孔)があくことが知られている。というのも、15パターンの立方体の隣接の仕方次第で、孔があくパターンがあるからでる。

これを改善するために、孔を塞ぎうる追加パターンを加えることや、isosurface の面を貼るときの内分点をいじるみたいなことがあるが、既に実装してくれているこれを使ってみる。これだと物体は(たぶん必ず)閉じるらしい。
make でコンパイルしたあとで、Bin/Linux にIsoSurfaceExtraction という実行ファイルができる。これは

    • in 入力 バイナリファイル
    • out 出力 ply 形式
    • res resolution 入力のxyz 格子のそれぞれの次元数
    • dim 実際のxyz 格子の長さ この顕微鏡イメージでは1.24x1.24x3 um/pix である
    • iso isosurface の閾値

となっているので適当に設定する。

# Python でデータ加工
vdim = np.shape(v)
seq = np.zeros( np.array(vdim).prod() )
seq[ np.sum(np.c_[v.nonzero()] * [1, vdim[0], vdim[0]*vdim[1]], axis=1) ] = 1
seq.astype("uint8").tofile("np.bin")
# ターミナル
./IsoSurfaceExtraction --in np.bin --out npout.ply --res 202 202 16 --dim 1.24 1.24 3 --iso 0.5

出力は3D オブジェクトデータフォーマットのうちのひとつ、ply である。R だとRvcg, geomorph で扱えるが、geomorph のread.ply 関数だと

read.ply("npout.ply") でエラー: 
  PLY file is not ASCII format: format = binary_little_endian1.0
 追加情報:  11 件の警告がありました (警告を見るには warnings() を使って下さい) 

と言われてしまうので、Rvcg のvcgPlyRead 関数を使うとうまくいく
Python ではpymesh もしくはplyfile が使えそうだが、どちらもなぜか使えなかったのでやめた。
 
実際には、meshlab がply ファイルから各種変換、また、マーチングキューブや表面平滑化などいろいろ取りそろえていて、かつ、コマンドラインから実行できるので

meshlabserver -i npout.ply -o npout.obj

とやればply ファイルからobj ファイルへターミナルのコマンドで変換できる。
 
最終的に物体表面の離散平均曲率で色付けしてobj ファイルを読み込んでR でplot3d したらこんな感じになる。

表面が閉じた物体であるので、S^2 同相に変換できる。