search

statgenetJimu / ShapeMove

ShapeMoveAnalysis Project
1 stars 2 forks source link

3Dのオプティカルフローについて #13

Open shigetosenoo opened 8 years ago

shigetosenoo commented 8 years ago

阪大の瀬尾です．

3Dのオプティカルフローについて調べており， Amat, F., Myers, E. & Keller, P. Fast and robust optical flow for time-lapse microscopy using super-voxels. Bioinformatics 29, 373–380 (2013). の論文が，自作のソフトウェアとITKの実装(Demon法)を比較していますので，これを再現＋αしてみようと思います．

この論文は，3Dでオプティカルフローを求める際に，SLIC(Simple linear iterative clustering)法 https://www.youtube.com/watch?v=TGaNGktTlhQ を利用したSupervoxelを考えてこれを計算に使うことで計算量を削減しているようです． この著者らの自作のソフトウェアはITKのコードを利用していないかもしれません．

shigetosenoo commented 8 years ago

私の環境はubuntu14.04でgcc version 4.8.4を利用． ソフトウェアは https://www.janelia.org/open-science/optical-flow-time-lapse-microscopy-using-super-voxels からダウンロード．

展開して出来たOpticalFlowの中で mkdir build cd build cmake .. make と実行．大量にエラーが出る．

コンパイラのバージョンによって？ヘッダ等のicludeの追加設定が必要だったり不要だったりするもよう． SLICinterface.cpp にstring.h， SLIC.cppにstring.hとlimits.h， opticalflow.cppにcstdlib， multithreadLBFGS.cppにstdio.h,string.hをinclude．

Makefileの中で"libtools -static -o"とするとエラーが出るので， http://stackoverflow.com/questions/11344547/how-do-i-compile-a-static-library 上記を参考に．"ar rcs"と修正．

../generatorがnot found と言われるので，mylib/MY_TIFF/generatorをbuild/mylibへこぴコピー．

これでmakeが成功しopticalflowの実行ファイルが出来る．

README.txtには terminal> opticalFlow syntheticData/rawData/TM00000 syntheticData/rawData/TM00001 syntheticData/rawData/test_backgroundPredictionIlastik_00001 syntheticData/opticalFlowOutput.bin 25 と書いてあるが，これだとSegmentaion faultが起きる． 正しくは， terminal> opticalFlow syntheticData/rawData/TM00000/test_00000 syntheticData/rawData/TM00001/test_00001 syntheticData/rawData/TM00001/test_backgroundPredictionIlastik_00001 syntheticData/opticalFlowOutput.bin 25 と実行する．

求めたフローはsyntheticData/opticalFlowOutput.binに入っているがバイナリのため，そのままでは開けない． 「The optical flow from TM00001 to TM00000 should be stored in binary file syntheticData/opticalFlowOutput.bin. The output binary file is stored in floats and has three channels (x,y,z optical float), with each channel beign the same size as the input images.」とのこと．

statgenetJimu commented 8 years ago

瀬尾先生：色々あったけれど、インストールはうまく行って、出力もあった、ということですね。 時間追跡して同一物体がどれかを決める処理がうまく行くのは、プロジェクトの肝なので、よい方法ならぜひ使いたいです。「追跡」がうまくいった後、それを「整数値」分離して4次元ボクセルに直せると、そのあと、別のアプリのアルゴリズムとかに乗せられるのですが、出力の変換も行けそうでしょうか？　山田

shigetosenoo commented 8 years ago

山田先生：はい，とりあえず出力までいきました．これの中身を見るとどうやら32bitのfloatになっているようです．xyzの各座標にxyz方向のフローが入っているはずなのですが，一本の数列になっているような状態です．これから並びを特定して可視化してみたいと思います．

shigetosenoo commented 8 years ago

test_00000.tifとtest_00001.tifをそれぞれ3D volumeで表示してgifアニメにしたもの．ガンマ補正してノイズを抑制．3Dイメージの大きさは96x128x26，グレースケール．

sample

上の動画像のオプティカルフロー．黒の点がフローの原点で赤の線がフロー（displacement）． of

フローを描くためのRのコード．恥を恐れずfor文を使う．

#バイナリファイルの読み込み（32bit float）
to.read<- file("opticalFlowOutput.bin","rb")
vec<-readBin(to.read, numeric(),128*96*26*3,size=4)
close(to.read)

#フローのプロット
array<-array(vec,dim=c(96,128,26,3))
carray<-apply(abs(array),c(1,2,3),sum)

xarray<-array[,,,1]
yarray<-array[,,,2]
zarray<-array[,,,3]

coords<-NULL
arrows<-NULL
for(x in 1:96){
  for(y in 1:128){
    for(z in 1:26){
      if(carray[x,y,z]!=0){
        coords<-rbind(coords,c(x,y,z))
        arrows<-rbind(arrows,rbind(c(x,y,z),
                                   c(x+xarray[x,y,z],y+yarray[x,y,z],z+zarray[x,y,z])))
      }
    }
  }
}
library(rgl)
plot3d(coords,xlab="x",ylab="y",zlab="z",xlim=c(1,96),ylim=c(1,128),zlim=c(1,26))
segments3d(arrows,xlab = "x", ylab = "y", zlab = "z",col="red")
statgenetJimu commented 8 years ago

いいですねぇ。

３D+時間=4Dですが、空間２次元(ｘ、ｙ）の解像度と空間１次元(ｚ)と時間次元との解像度に大きな格差があるので、おそらく （A) 4Dで攻める (B) 2D空間特徴量の時間変化を追う のように「３次元物体映画」と「２次元物体映画」との両方を解析の対象とすることになると思います。 そのときに、２次元と３次元とを別立てて解析処理することにすると、初期プロトタイプ段階は楽だとしても、２年後くらいには嫌になると思っています。というわけで、ある意味任意次元にしようと思っていて、そのうえで、以下のような処理が必要になるかなーと思っています(トポロジーについては今のところ不要)。

以下で言うところの「クラスタ」は「時間次元、コミ」の連結成分のつもりです。

１１月５日にお会いするときに、詳しく議論できればと思っています。

ラベリングしてボクセルのクラスタを作る

時系列ラベリングの場合には、空間的に全時刻で連結であることをもって

連結体(クラスタ)とみなす、などの制約を入れることも考慮するが

今は単純に全次元での連結とする

クラスタごとに特徴量抽出をする

1 IDを振る

2 全体積(４次元空間で連結体となっているものの全ボクセル数をカウント)

3 空間体積の時系列ベクトル

4 全表面積(４次元空間で連結体となっているものの周囲面積計算→メッシュ化してから…)

5 空間表面積の時系列ベクトル

6 全重心座標

7 空間重心の時系列ベクトル

8 空間三角化メッシュの時系列

9 空間曲率分布の時系列(頂点空間曲率の頂点重み付き分布)

10 全トポロジー

11 空間トポロジーの時系列

12 全PCA分解

13 空間PCA分解の時系列

14 全PCA超楕球からの残差

15 空間PCA分解超楕球からの残差の時系列

クラスターごとの抽出特徴量から「意味を与えた」解析項目を作る

クラスターごとの「意味」

1 軌道 (速度・向き・それらの変化)

2 拡大縮小　体積の時系列変化

3 回転　空間PCA分解の時系列変化

4 形の変化　空間PCA分解超楕球からの残差の時系列変化

5 曲率の分布、曲率の空間微分の分布、さらにその微分の分布

6 曲率の時間微分については、2D空間＋1D時間の3Dにすることで、いわゆる3D画像の離散外積部分幾何(DEC)を使うことができる

クラスター集合の分布に関する「意味」

1 情報量の多い視野という指標を作る

3 クラスターが作る分布の時系列変化

4 クラスターが作る分布の位置別比較

上記のうち、少しわかりにくいのは、立体形状全体をDEC扱いすることとその後の処理

３次元オブジェクトの形全体を扱う・その時間変化を扱うには、以下のようにする

Voxel x time = 4D data

3D x 1T とする

3Dごとに周縁をメッシュ化する

Fiaringアルゴリズムにて、正球にまで変形する

球面スカラー場とする

球面調和級数関数分解する

スペクトルの時間微分をとる

かたちの比較：スペクトル比較 動きの比較　：スペクトル時間微分の比較