haspet.doc : 1

画像再構成と基本的な画像処理の手順
目次(1) ソフトウェアのインストール(1-1) 書庫ファイルのダウンロードと解凍・展開(1-2) ディレクトリ local/ の下の実行ファイル群のインストール(1-3) image/、cuda/ および etc/ の下の実行ファイル群のインストール(1-4) 実行パスの設定(2) 作成するディレクトリとファイル(3) 画像再構成の手順(3-1) 測定データのコピー(3-2) X 線透過率画像の動画の作成(3-3) X 線投影画像の作成と観察(3-4) DO の推定(3-5) 画素の辺長の設定(3-6) 画像再構成のテスト(3-7) 本番の画像再構成(4) 基本的な画像処理の手順(4-1) 再構成した画像の画素値ヒストグラムの調査(4-2) 8 ビット画素値の画像の作成(4-3) 画像のトリミングと縦断面の画像の作成(4-4) 縮小画像の作成(4-5) 鳥瞰図画像の作成(5) 実例のログ（040711j/log.txt）

haspet.doc : 2(1) ソフトウェアのインストール(1-1) 書庫ファイルのダウンロードと解凍・展開
ユーザ haspet のホームディレクトリ ~/ で以下の入力を行う。
wget http://www-bl20.spring8.or.jp/\~sp8ct/tmp/haspet.taz
tar xzf haspet.taz
rm haspet.taz
この結果、以下のディレクトリとファイルが ~/ に展開されるはずである。
local/
image/
cuda/
etc/
run/
PATH
Cshrc

# 中野が書いたものではないソフトウェア用# 中野が書いた画像再構成や画像処理用のソフトウェア# CUDA-GPU を使って画像再構成などを行うプログラム群# 中野が日常的に使っているユーティリティ群# SPring-8 で撮った X 線 CT 画像の再構成と基本的な画像処理用# 実行パスに設定すべきディレクトリ名を記したテキストファイル# 実行パスの設定を記したテキストファイル(1-2) ディレクトリ local/ の下の実行ファイル群のインストール
まず、ディレクトリ local/ の下にある中野が書いたものではないソフトウェアをインストールする。以下の入力によってそれぞれをインストールできるはずだが、処理中にエラーが発生するかもしれない。
TIFF のライブラリ libtiff（slice* だけが使う）と TIFF 用のツール一式のインストール
pushd local/src/tiff
csh INSTALL-3.6.1
popd
HiPic 形式画像を取り扱えるよう改造した画像閲覧ツール xv（もしくは、xv-hipic）
pushd local/src/xv
csh INSTALL-hipic
popd
動画表示用のソフトウェア xanim
pushd local/src/xanim
csh INSTALL2801
popdhaspet.doc : 3

グラフの修正などに使える「お絵かき」用ソフトウェア tgif（不要かもしれない）
pushd ~/local/src/tgif
csh INSTALL-2.14-p8
csh INSTALL-3.0.18
csh INSTALL-ln-s
popd# ver. 2.14 patch level 8 のインストール# ver. 3.0.18# 普段は ver. 3.0.18 の tgif を使うようにする設定(1-3) image/、cuda/ と etc/ の下の実行ファイル群のインストール
ディレクトリ image/、cuda/ と etc/ の下にある中野が書いたプログラム群のインストール中にエラーは発生しないと思われるので、以下のような入力によってすべてを機械的に処理すればよい。
foreach src (image/src/*
pushd $src
csh INSTALL
popdendcuda/src/*

etc/src)

ただし、cuda/ の下のプログラム群のインストールは CUDA toolkit 2.3 の環境を想定しており、自動でダウンロードした書庫ファイルに入っているものではなく cuda/src/*/ の下に用意した Makefile を使用する。これらの Makefile の先頭付近を書き換えれば 2.3 以降の版の CUDA toolkiｔ の環境にも対応可能である。
(1-4) 実行パスの設定
ファイル Cshrc に記されている下記の内容を ~/.cshrc（なければ作成する）の最後に付加する。
setpath=(. `sed s@\^@~/@g ~/PATH`$path)

その後、新しい端末を開くか以下の入力を行うかすれば上記の実行パスの設定が有効になる。
source~/.cshrc

haspet.doc : 4(2) 作成するディレクトリとファイル
以下ではディレクトリ run/ の下にある C-shell scripts "*.csh" などを用いて、X線 CT測定ごとのディレクトリの中に以下のディレクトリやファイルを作成する。ただし、テキストファイル "*.log" の内容はそれを書き出した
C-shell script のそれぞれが端末に表示するものと同一である。
raw/
測定データのファイル（output.log、dark.img と q*.img）を入れるディレクトリ
xts.gif
透過 X 線強度画像の動画のファイル。作成しないかもしれない。
xp/
０度と 180 度方向から透視した X 線投影画像のファイルを入れるディレクトリ
xp.log
各行に作成した X 線投影画像に関する以下の３個の値が記されているファイル[1] サンプルを透視した方向（== 度単位のサンプルの回転角）[2-3] その画像に出現した X 線投影の値の最小値と最大値
DO.tbl と DO.log
DO とは Detector Offset の略である。これは測定画像上でサンプル回転軸の位置を横方向の原点とした座標系における画像の左端にある画素の位置を指す座標値であり、いわゆる「センター値」の符号を変えた（それに負の符号を付けた）値と一致する。
DO.tbl には０度と 180 度の方向から透視した２枚の X 線投影画像を使って推定したDO の値に関する以下の情報（それぞれに複数の値が並んだ行）が書き込まれている。
１行目（２個の値が並んでいる）[1] 測定画像の横画素数[2] 縦画素数
２行目以降（４個の値）[1] スライス番号（０～）[2] そのスライスに対する DO の推定値
[3] そのスライスの投影値に含まれる「ノイズの強度」[4] そのスライスの「影の濃さ（投影値の自乗和）」
最後の行（３もしくは５個の値）[1] スライスごとの DO の推定値を最小自乗近似して得た回転軸の傾き[2] 上記の最小自乗近似によって得たスライス番号０の位置の DO の値[3] 上記の最小自乗近似の推定誤差（１以下の値であることが望ましい）[4-5] 上記の最小自乗近似に用いたスライスの範囲
また、DO.log には DO.tbl の最後の行と同じ内容が転記されている。
center.log
０度と 180 度方向から透視した２枚の X 線投影画像の２次元パターンマッチング（０度の方向から透視した画像と左右反転した 180 度の画像を様々に配置した場合の、重複した部分の投影値の RMSD

haspet.doc : 5の最小値の計算）によって推定したセンター値とそれらの画像の縦方向の位置ズレに関する以下の６個の値が記されている。
[1-2] 吟味したセンター値の範囲[3] センター値の推定値[4-5] 吟味した２枚の画像の縦方向の位置ズレの範囲[6] 縦ズレの推定値（単位は画素幅）
test/
DO の値を様々に変えて再構成した画像用のディレクトリ。この下に設けたスライスごとのディレクトリの中に「センター値.tif」というファイル名で再構成画像が格納されている。
test.log
test/ の下の再構成画像に関する以下の４個の値が各行に書き込まれている。
[1] スライス番号[2] 再構成の際に指定した DO の値[3-4] その再構成画像に出現した CT 値（LAC の推定値）の最小値と最大値
tg/
本番の画像再構成で得た 16 ビット画素値のスライス画像を入れるディレクトリ。ただし、それらの画像上の CT 値と画素値はスライスごとに異なった対応関係になっている。
tg.tbl と tg.log
tg.tbl には tg/ の下に格納されている再構成したスライス画像それぞれに関する以下の４個の値が各行に書き込まれている。
[1] スライス番号[2] そのスライスの画像再構成に用いた DO の値[3-4] そのスライスに出現した CT 値の最小値と最大値
また、tg.log にはすべてのスライスに関する以下の６個の値が書き込まれている。
[1] 画素（立方体）の辺長[2] スライス画像（正方形）の横および縦の画素数[3-4] 再構成したスライスの番号の最小値と最大値[5-6] すべてのスライス画像に出現した CT 値の最小値と最大値

haspet.doc : 6

ｗord/
tg/ の下にある再構成画像から変換した、それらのスライスすべてに出現した CT 値の値域（最小値と最大値）で正規化した 16 ビット画素値の画像を入れるディレクトリ
word.log
tg/ の下のすべてのスライス画像に出現した CT 値の値域が書き込まれている。
word.csv
word/ の画像に出現した CT 値の出現頻度分布（ヒストグラム）に関する以下の４個の値がコンマ区切りで各行に書き込まれている。
[1] ヒストグラムのビン番号（０～65535）[2-3] そのビンに対応する CT 値の下限と上限の値[4] そのビンに相当する CT 値を保持している画素の総数
byte/
tg/ の下にある再構成画像から変換した、指定した CT 値と画素値の対応関係（ただし、CT 値０に対応する画素値は０）に従う ８ ビット画素値の画像を入れるディレクトリ
byte.csv
byte/ の画像上の CT 値のヒストグラム情報。書き込み形式は word.csv と同じ。
trim.log
指定した画素値のしきい値と、byte/ の下の３次元画像上でそのしきい値以上の画素値を持つ画素が占める直方体領域の対角点 (x1,y1,z1) と (x2,y2,z2) の６個の座標値がこの順に記されている。
ｔrim/word/ の下の画像から切り出した３次元画像を入れるディレクトリ
trim_xy/ および trim_xz/ と trim_zy/
byte/ の画像から切り出した３次元画像と、それを切り直した相互に直交する２セットの縦断面のスライスからなる３次元画像のそれぞれを入れるディレクトリ
画素数_xy/、画素数_xz/ と 画素数_zy/
trim_xy/、trim_xz/ と trim_zy/ の下の３次元画像を縮小した画像を入れるディレクトリ。それらの名前の一部となっている、各辺沿いの画素数が指定した値の立方体領域に含まれる画素群の値の平均が新しい画素値となる。
bev_word_画素数_画素値.log…
bev_word_画素数_画素値/…
bev_byte_画素数_画素値.log…
bev_byte_画素数_画素値/…

haspet.doc : 7(3) 画像再構成の手順(3-1) 測定データのコピー
測定データ用にディレクトリ raw/ を作成し、そこに以下のファイルをコピーする。
output.log
X 線 CT の測定画像 "q*.img" の測定シーケンスが記録されているファイル
dark.img
暗電流強度画像の HiPic 形式ファイル
q*.img
入射もしくは透過 X 線強度画像の HiPic ファイル
ただし、これら以外のファイルやディレクトリが raw/ の下にあってもかまわない。 raw/ の下に画像再構成に必要なファイルすべてがそろっているかどうかは以下の入力によってチェックできる。
dit_tblraw > /dev/null(3-2) X 線透過率画像の動画の作成
以下を入力すれば X 線透過率画像の動画のファイルを作成できる。
xts.csh# 起動パラメータなし

これにより、測定画像から得られるサンプル回転角が０～180 度のX線透過率画像を以下のように間引きしたフレームから構成される GIF 形式の動画のファイル xts.gif ができるはずである。
フレーム数
総フレーム数が 150 以下になるように等間隔の回転角の画像を抜き出す。
それぞれのフレームの画像の画素数横画素数が 500 以下になるように隣接した複数の画素をまとめる（平均する）。
以下のように入力すれば作成した動画を閲覧できる。
動画の各フレームをスチール画像としてじっくりと見たい場合
xvxts.gif# PGUP や PGDN キーなどで表示フレームを変更できる。

動画として表示したい場合は以下のいずれかを入力すればよい。
xanim xts.gif
firefox xts.gif# メニューボタンによって動画停止や逆順表示を指示できる。
# GIF 形式動画は普通の WEB ブラウザで表示可能なはず。

haspet.doc : 8(3-3) X 線投影画像の作成と観察
まず、以下の入力により、サンプルを回転角が０度の方向から透視した X 線投影の画像 000.tif と 180 度の方向から透視した画像 180.tif をディレクトリ xp/ の下に作成する。
xp.csh# 起動パラメータなし

ただし、xp/ の下の２個の画像の８ビット画素値は画像それぞれに出現した投影値の値域（最小値と最大値）で正規化した値である。これらの値域はテキストファイル xp.log に記されている。
その後、以下のように入力してこれらの画像を表示する。
xvxp/*

表示された投影画像をじっくりと観察して以下の情報を調べる。
Nr と Nz
画像の横と縦の画素数。Nr は再構成画像（正方形）の各辺に沿った画素数、Nz は再構成可能なスライスの総数となる。これらの値は xv で表示した画像上でマウス右ボタンをクリックすると出現する xvの処理メニューに示されているはずである。
z1_DO と z2_DO
DO の値の推定の際に指定する、画像上で投影値が特に高い部分（== サンプルの像の影が濃いスライス）の縦座標値の範囲。これらは画像の上でマウス左ボタンのクリックによって表示される座標値を読み取って決める。
z_test
画像再構成のテストを行うスライスの位置。影が濃いスライスを指す縦座標値を選ぶ。
z1_tg と z2_tg
本番の画像再構成を行うスライスの範囲。

haspet.doc : 9(3-4) DO の推定
以下のように入力して位置 z のスライスにおける DO の値 DO(z) を推定し、さらに、最小自乗法でそれらを近似する式 DO(z) = α×z + β±σ の係数値（σは推定誤差）を計算する。
すべてのスライスの DO(z) を用いて最小自乗近似を行う場合
DO.csh
X 線投影画像の観察で決めた範囲のスライスの DO(z) だけを使う場合
DO.csh

z1_DO

z2_DO

いずれを行った場合も、テキストファイル DO.tbl の２行目以降にそれぞれのスライスの DO の推定値 DO(z)が記録される。
端末表示された最後の行（== DO.tbl の最終行やテキストファイル DO.log に記されている内容）にα、βおよびσの値がこの順に並んでいる。サンプル回転軸の傾きを表すαが 1/Nz 以上、もしくは、σが１以上の値な
らここで推定した DO を画像再構成に使えない。その場合は、
rm

DO.tblDO.logと入力して DO.tbl と DO.log を消去した後、z1_DO と z2_DO を別の値に変えて DO.csh を再試行すべきである。ただし、その前に、以下の入力によって測定中にサンプルが動いていないか否かを調べた方が良い（DO.csh による推定の結果に問題がない場合でも、念のために以下の入力を行った方が良い）。
center.csh# 起動パラメータなし

これにより center.log に記録される６個の値が端末にも表示される。それらのうちの最後の値が指す０度と180 度の方向から透視した２個の画像の縦方向の位置ズレが大きいなら、X 線 CT の測定中にサンプルが動いた可能性がある。

haspet.doc : 10(3-5) 画素の辺長の設定
画像再構成の前に以下の入力を行って、測定画像の正方形画素の辺長（もしくは、再構成画像の立方体画素の辺長）を環境変数 PS_CM に設定しておく必要がある。
setenv PS_CM

画素の辺長の値

ただし、単位が 1/cm の LAC を推定（再構成）するためには cm の単位で画素の辺長の値を設定しておく必要がある。なお、PS_CM に設定した値は以下の入力によって端末に表示できる。
echo $PS_CM(3-6) 画像再構成のテスト
X 線投影画像の観察によって決めたサンプル像の影が濃いスライス z_test に対して、DO の値を系統的に変えた画像再構成のテストを行う。ここでは以下のようにして指定した値を中央にした -0.5 刻み、11 通りのDO の値を用いた 11 枚の画像の再構成をそれぞれ行う。
スライスごとの DO の推定値を最小自乗近似した式で中央の DO の値を計算する場合
test.cshz_test

陽に指定した適当な DO を中央の値とする場合
test.cshz_testDO

この結果、ディレクトリ test/z_test/ の下に「センター値.tif」という名の再構成画像のファイルができる。また、それらに関する情報が端末とテキストファイル test.log に出力される。
このようにして得た test/z_test/ の下の一連の再構成画像の画質を観察する。
xv test/*/*もしくは、不適切な DO の値を指定した再構成画像に生じる偽像（artifact）を際立たせるため、画像の表示輝度をヒストグラム平滑化（histogram equalization）で強調した方が良いかもしれない。
xv-histtest/*/*

いずれにせよ、最も画質が良い（偽像が生じていない）再構成画像の DO の値を確認し、それを本番の画像再構成で使う。

haspet.doc : 11(3-7) 本番の画像再構成
以下のように入力して本番の画像再構成を行う。
すべてのスライスを再構成する場合
tg.cshortg.cshDO

X 線投影画像の観察によって決めた範囲のスライスだけを再構成する場合
tg.csh

z1_tg

z2_tg

tg.csh

z1_tg

z2_tg

orDO

ただし、DO.csh や center.csh で DO の値を推定できなかった場合は、画像再構成のテストによって決めたDO の値（再構成するスライスすべてに対して同じ値）を上のようにして指定する。そして、その指定がない場
合は、tg.csh は DO.tbl の最終行に記されている係数値（αとβ）を調べて、それらを用いて計算したスライスごとに異なる DO の値を使って画像再構成を行う。
なお、本番の画像再構成には長い処理時間を要するはずである。以下の入力によってスライス画像１枚当たりの再構成の処理時間を得ることができる。
time.csh# 起動パラメータなし

ただし、これによって端末に表示される３個の値の意味は以下の通りである。
[1] その時点で画像再構成が完了したスライスの枚数[2] それらの１枚当たりの画像再構成の処理時間（単位は秒）[3] 処理時間の標準偏差（単位は秒）haspet.doc : 12(4) 基本的な画像処理の手順(4-1) 再構成した画像の画素値ヒストグラムの調査
word.csh # 起動パラメータなし(4-2) 8 ビット画素値の画像の作成
byte.csh

LAC_step(4-3) 画像のトリミングと縦断面の画像の作成
trim.csh

PV_threshold(4-4) 縮小画像の作成
reduce.csh

bin_size(4-5) 鳥瞰図画像の作成
bev_word.csh bin_size LAC_threshold
bev_byte.csh bin_size PV_threshold(5) 実例のログ（040711j/log.txt）

haspet/040711j/log.txt : 1% cd 040711j% ls -Ftr
raw.csh* memo.txt log.txt% raw.csh...
raw/d01.img : time data not found (warning).
raw/d02.img : time data not found (warning).
% dit_tbl raw >/dev/null% xts.csh% xanim xts.gif% xp.cshd.img = raw/d.img
i.img = raw/q001.img
t.img = raw/q002.img at angle = 000.5000
d.img = raw/d.img
i.img = raw/q398.img
t.img = raw/q399.img at angle = 180.5000
000-0.077514
1.882011
180-0.080043
1.892347% xv xp/*% DO.csh 240 480
D000 = raw/d.img
I000 = raw/q001.img
T000 = raw/q002.img at angle = 000.5000
D180 = raw/d.img
I180 = raw/q398.img
T180 = raw/q399.img at angle = 180.5000
0.000000-498.500000
0.000028% center.csh
using z1 and z2 at the last line of "DO.tbl".
d.img = raw/d.img
i.img = raw/q001.img
t.img = raw/q002.img at angle = 000.5000
d.img = raw/d.img
i.img = raw/q398.img
t.img = raw/q399.img at angle = 180.5000
249.5 749.5 498.5 -60
60
-1

240

480

haspet/040711j/log.txt : 2% test.csh 360
using DO at the last line of "DO.tbl".
using pixel size (cm) of 0.0001 as the dummy value (warning).
360-496-14.172132
47.400288
360-496.5 -13.166102
47.025486
360-497-11.543563
45.441139
360-497.5 -10.387976
44.337341
360-498-8.775000
43.321228
360-498.5 -8.668720
42.373047
360-499-7.891128
42.976723
360-499.5 -10.095286
44.542908360-500-12.719671
45.441109360-500.5 -14.557191
46.236580
360-501-16.099836
46.932644% xv -e 0.5 -hist test/*/*% setenv PS_CM 0.000583% tg.csh0.000583

1000% time.csh
720
0.158554% word.csh-4.375421

000719-4.375421

35.706387

0.365259

35.706387% byte.csh 0.02930016% trim.csh 125
125
102
132

0876

865

705

% reduce.csh 2% bev_word.csh 2 3.662520
23.662520
13142% bev_byte.csh 2 125
2125% ls -Ftr
raw.csh* DO.tbl
memo.txt DO.log
log.txt center.log
raw/
test/
xts.gif test.log
xp.logword.log
xp/tg/

tg.tbltg.logword/
word.csv
byte/
byte.csv
trim.log

trim/
trim_xy/
trim_xz/
trim_zy/
2_xy/
2_xz/
2_zy

bev_word_2_13142.log
bev_word_2_13142/
bev_byte_2_125.log
bev_byte_2_125/

haspet̲040711j̲input.rtf : 1(3-1) 測定データのコピー
raw.cshcd ~/hdd
mkdir 040711j
ln -s ~/sdd/040711j raw
dit_tbl raw >/dev/null(3-2) X 線透過率画像の動画の作成
xts.cshxanim xts.gif(3-3) X 線投影画像の作成と観察
xp.cshxv xp/*(3-4) DO の推定
DO.csh 240 480
center.csh(3-5) 画素の辺長の設定
setenv PS_CM 0.000583(3-6) 画像再構成のテスト
test.csh 360
test.csh 360 -498.5
xv -hist test/*/*(3-7) 本番の画像再構成
tg.cshtg.csh -498.5
tg.csh 0 719
tg.csh 0 719 -498.5
time.cshhaspet̲040711j̲input.rtf : 2(4-1) 再構成した画像の画素値ヒストグラムの調査
word.csh(4-2) ８ビット画素値の画像の作成
byte.csh 0.02930016(4-3) 画像のトリミングと縦断面の画像の作成
trim.csh 125(4-4) 縮小画像の作成
reduce.csh 2(4-5) 鳥瞰図画像の作成
bev_word.csh 2 3.662520
bev_byte.csh 2 125