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

蛍光イメージング動画があって、そこから細胞の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

画像データをopenCV のPython 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 格子の各点にある値 が乗っていて、その値が閾値を超えるか超えないかで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 したらこんな感じになる。

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