サルにもわかる Houdini: dicom import

今更ながら、かなり昔 (2014年) の Dicom Import に関して

使用したビルド:

Houdini 18.5.499 Python 2.7 (Windows 10)
Pydicom (非常に古い) 2.7 対応のもの。ここから pydicom-0.9.8.zip をダウンロードしたのだと記憶している。
HDA: dicom_hdas.zip

1. Pydicom 設定および動作確認

ここでは、環境変数 PYTHONPATH を C:\Python27\Lib\site-packages に設定している。

dicom.zip を解凍してできる dicom ディレクトリを site-packages に配置する。
dicom_hdas.zip を解凍し、 awong_copfromdicom.otl と awong_volumefromcop.otl の二つのファイルを HDA のパスに配置する (例: ~/houdini18.5/otls)。
Houdini を起動、ネットワークエディタで COPs (/img/comp1 など) に行き、Tab-> Digital Assets -> COP from DICOM を実行、ノードを配置してエラーが出なければ、pydicom モジュールが認識ていると考えられる。以下のようなエラーが出た場合にはモジュールが認識されていない。
```
Error 
Python error: Traceback (most recent call last):
File "<stdin>", line 2, in <module>
ImportError: No module named dicom
```
以下のエラーの場合には入力画像が指定されていないということで、次 (以下) に進む。
```
Error 
Python error: Error generated by Python node.
Unable to read dataset. 
```

2. DICOM ファイルの読み込み

ここで扱う DICOM ファイルは2次元の画像ファイルフォーマットなので、まず COPs に読み込んでから3次元 (ボリューム) 化する。ここでは、 vhm.[1001-1245].dcm という連番ファイルを読む。

COPs で COP from DICOM ノードを作成 (copfromdicom1)、DICOM Files にファイルを指定。
- ブラウザで最初のファイル vhm.1001.dcm を選択
- その後、vhm.1001.dcm の部分を vhm.`1000+$F`.dcm と書き換える。
- Composite View でタイムラインをスクラブすれば、1フレーム目から245フレーム目の間で連番画像がすべて読み込まれているのがわかる。
  
  上の画像は 122 フレーム目。左上で解像度が 512x512、fr (フレーム) 122、A (アルファ) チャネル画像が表示されていることがわかる。
  Composite View 上で 'I' キーを押せば、カーソル直下のピクセル値を数値で表示できる。このピクセルでは Alpha の値が 0.0372 とわかる。
- 左上の Histogram ボタンを押せば入力画像がヒストムグラム表示される。
  
  タイムラインをスクラブすれば各フレームでのアルファ値の分がわかる。ここでは、ほぼすべてフレームにおいて 0 から 0.05 の間で、0.1 を超えることがないことがわかる。
Null ノードを作成し、copfromdicom1 の出力に接続、名前を OUT として、表示とレンダーフラグの両方を指定。

copfromdicom1 と OUT の間に画像処理ノードを挿入してチャネル補正やらその他編集を、後からであっても、行っても良い。

3. 2D から 3D

ネットワークエディタで /obj に戻り、Geometry ノードを作成、その中に入る。
TAB-> Digital Asset -> Volume from COP を実行し、ノード (volumefromcop1) を配置。
```
Error 
Error generated by Python node.
Invalid COP Path 
```
とエラーが出るが、 COP パスが指定されていないだけなので問題ない。よって、これを指定する。
COP Node に上で作成した OUT ノードを指定する (例: /img/comp1/OUT)。
Channels 脇の '+' を一度クリックしてチャネルを追加。
- Target に density
- Source に A を指定する。
- 残りはデフォルトのまま。
しばらくすると、3D ビューポートにスモークボリュームによる画像シーケンスをまとめた (スタックした) ものが表示される。

なお、上記の設定をしてもボリュームが表示/計算されない場合は、volumefromcop1 を一度ロック(赤フラグ)してアンロックすると実行される(はず)。

なお、背景は Dark (d キーを押して Display Options を表示、 Background タブの一番上にある Color Scheme を Dark にする) と良い。
Z Up でない場合は、volumefromcop1 の出力に Transform (transform1) を追加、Rotate X を -90 とすると正しい向きになる。
transform1 の出力に Volume Mix (volumemix1) を追加
- Post Mult の値を 0.025
- 残りはそのまま
すると、骸骨の部分だけが表示される。
volumemix1 上で RMB->Save->Geometry... として .bgeo で保存すれば、File ノードで別途読み込むことができる。そうすれば Python を何らかの拍子で走らせる可能性がなくなる。このデータの場合
- .bgeo で保存すると 140MB
- .bgeo.sc で保存すると 81.6MB
.bgeo は SOP で扱えるデータタイプをすべて保存できる。

必ずしも Volume Mix を使わなければならないわけではなく、他のノードでもできそうだが、これで2D画像からの3D化ができたことになる。

4. データビジュアリゼーション (Volume Slice)

3D ビューポートでデータを可視化するには、以下の方法をとることも可能。

volumemix1 の出力に Volume Slice (volumeslice1) ノードを接続。
- Method を Mesh に
- Attribute は density
- Visualization Ramp を Infra-Red (いろいろ試すと良い)
- Visualization Range を 0 から 0.001 ぐらいにする
Offset の値をスライダで -1 から 1 の間で動かせばスライス面での density の値が色表示される。

上の画像は Keep Origianl Geometry をオンにした状態。
Visualization Ramp と Visualization Range を調整して様々な可視化を試すのも良い。

5. メッシュ化

上記の可視化のプロセスで分かったように、density の値は、0 ~ 0.03 など非常に小さい範囲にある。ボリュームデータからポリゴンに変換する前に、スケールするのが良い。volumemix1 に戻って Post-Multi の値を上げても良いが、ここでは Volume Wrangle を使う。

volumemix1 の出力から分岐する形で Volume Wrangle (volumewrangle1) を追加、VEXpression 部分に
```
f@density*=2.0;
```
とする。
個別の値を入力するのではなく、スライダで調整したい場合は、
```
f@density*=chf("scale");
```
とすれば volumewrangle1 のパラメータペーンに scale スライダが追加される。
volumewrangle1 の出力に Convert Volume Node を追加
- Isovalue を 0.00098 などと小さい値にする。
するとポリゴンに変換される。
ノイズになっているところは値を調整して出力されないようにするか、選択して削除するかいずれかが考えらる。

6. HDA に関して

ここで使っている二つの HDA はどちらも Python で記述されている。よってそれぞれのノード上で RMB->Type Properties... で Type Properties ウィンドウを開き、Code タブに行けば、実際のコードが確認できる。

COP from DICOM は単一ファイルに連番が含まれているタイプの .dcm には未対応。

最終更新: 2021-03-02

サルにもわかる Houdini

2021-02-26

dicom import