6mct/recipe.txt : Tue Jan 23 11:13:28 JST 2018
SIXM の画像再構成などの手順[0] プログラムのインストール法[1] 作業用ディレクトリと測定データの準備[2] 測定データのチェックなど[3] MF のための RP 画像の観察[4] R[I,R] 画像の作成[5] 作成した R[I,R] 画像に対するオプションの処理[6] R[I,R] 画像上の物体像を囲む bounding box の選択[7] スライスごとの値の「時間変動」や鉛直方向のサンプルの「動き」のチェック[8] サンプル回転軸の位置の推定[9] サンプル回転軸の位置を変えた画像再構成のテスト[10] 吸収および位相 CT 画像の再構成とそれらの確認[11] 吸収 CT 画像に対する後処理[12] 位相 CT 画像に対する後処理[13] 長さと表示輝度のスケールバーを付けた吸収と位相 CT 画像の作成（オマケ）[0] プログラムのインストール法
書庫ファイルのダウンロードと展開% cd ~/Desktop% wget http://www-bl20.spring8.or.jp/\~sp8ct/tmp/6mct.taz% tar xzf 6mct.taz% rm 6mct.taz
CPU 用のプログラムのコンパイル% cd 6mct% makeGPU 用のプログラムのインストール% cd cuda% vi Makefile
NVCC =nvcc -O3 -arch=sm̲50
CUDA =/usr/local/cuda% make% mv *̲g ..
% cd ..
ImageJ のランチャーとマクロのインストール% cp -p ImageJ.csh ~/Desktop/ImageJ% cp -p tlrb.txt ~/Desktop/ImageJ/macros
C-shell の設定ファイルへの追加% vi ~/.cshrc
set path=(~/Desktop/6mct $path)alias IJ '~/Desktop/ImageJ/ImageJ.csh \!*'
[1] 作業用ディレクトリと測定データの準備
測定ごとのディレクトリを作成し、その中に移動してすべての作業を行う：% mkdir 170516c% cd 170516c/dev/shm/170516c

SIXM の測定で得た２個のファイルをリンクしておく（コピーでも良い）：% ln -s /work1/tsukasa/1705/170516c/a.log .
% ln -s /work1/tsukasa/1705/170516c/a.HIS .
[2] 測定データのチェックなど
SIXM の測定データファイルのバイト数などを調べる：% ls -trlL
total 41275208-rwx------ 1 tsukasa tsukasa
14074443 May 16 18:52 a.log-rwx------ 1 tsukasa tsukasa 42251730944 May 16 18:52 a.HIS
ファイル a.log に記されている SIXM の測定のパラメータ値を調べておく：% head -1 a.log
551,601,100,1,0
これらのコンマ区切りの数値のうちの最初の３個の意味は以下の通りで、それらを後述するプログラム his2raw の起動パラメータとしてそのまま指定する：
scans = 551
各 view の最初の scan はゴミなので、R[I,R] 画像の横画素数は 550
views = 601
CT に使用する 180 度回転のステップ数は 600 views
darks = 100
データファイル a.HIS の先頭に格納されている暗電流 RP 画像の個数
データファイル a.HIS に格納されている RP 画像のサイズなどを調べる：% check̲his a.HIS | uniq -c
42251730350
331251 128
498
2
ただし、check̲his などが出力した数値の意味は以下の通り：
42251730350 ≦ 42251730944
有効な RP 画像の総バイト数。
この値はファイル a.HIS のバイト数以下のはず。
331251 == 551（scans）* 601（views）+ 100（darks）
有効な RP 画像の個数。
ファイル a.log に記されている値から計算した値と一致するはず。
128
RP 画像の横画素数。
正確には、16 個づつの binning 後の画素数。
498 = Nz
RP 画像の縦画素数。
この値と等しい枚数の CT 画像のスライス画像を再構成できる。
2
RP 画像のタイプ（2 == 16 ビット整数画素値の画像）。
後述する his2raw はこのタイプの RP 画像しか処理できない。
これら以外に以下の SIXM の以下のパラメータ値を知っておく必要がある：
Dz = 108.1 (nm)RP 画像の画素の鉛直方向の辺長。
この値は以下の処理では陽に指定しないが、RP 画像から作成する R[I,R] 画像の画素の鉛直方向の辺長、および、それらから再構成する３次元 CT 画像のスライスの厚さとなる。
Dr = 100 (nm) = 100e-7 (cm)
SIXM のスキャンのステップ幅。

この値は位相 CT 画像の再構成処理では陽に指定しないが、それを含めた３次元 CT 画像のスライス画像の正方形画素の辺長となる。
SDD = 6220 (mm) = 6220e-3 (m)サンプル回転軸と検出器面の間の距離
Dp = 6.5 (um) * 16 (binning) = 104 (um)RP 画像の画素の水平方向の辺長。
RP 画像が正立像なので、Dp として負の符号を付けた -104 を指定する。
RA0 = 0 (degree)
サンプル回転角の初期値。
非０の値を指定すれば軸回転した画像を画質の低下なしで再構成できる。
他のパラメータ値と区別するため、以下では RA0 として +0 を指定する。
なお、後述する画像再構成の処理では Dr、SDD および Dp として再構成する値に応じた上記の単位の値を指定する必要がある：
吸収 CT 画像
単位が 1/cm の LAC の値を推定するので cm 単位の Dr の値を指定する。
位相 CT 画像
単位が 1e-6（100 万分の１）の RID の値を推定するので、SDD の単位に対して 100 万分の１の単位の Dp の値を指定する。
[3] MF（mountain flattening；山の平坦化）のための RP 画像の観察
R[I,R] 画像の作成時に RP 画像に対して MF の処理を行うかどうかを判断する。
そのためには複数の RP 画像を吟味することが望ましいが、ここではとりあえずR[I,R] 画像の「中央」の画素の値の計算に使う１個の RP 画像だけを調べる：
R[I,R] 画像の「中央」の画素の位置
view = views / 2 = 300
scan = scans / 2 = 275
layer = Nz（RP 画像の縦画素数）/ 2 = 249
これらとファイル a.log に記されていた値を用いてファイル a.HIS から１個のRP 画像を抽出し、その画像（もとのRP 画像）とそれに対してパラメータの値
Wm (mountain width) = Ws (samplineg width) = 5 および 10を指定して MF を行った結果の RP 画像などを ImageJ で表示してみる：% rp+mf̲ij.csh a 300 275 249 5,5 10,10-380.329987
63949.281250-380.329987
63949.281250[5] 作成した R[I,R] 画像に対するオプションの処理必要なら以下のような R[I,R] 画像の画素値のビット数の削減処理を行っても良い。
float（32 ビット画素値）の R[I,R] 画像を 16 ビットのものに変換% t2t̲float ri - 16 ri > /dev/null% t2t̲float rr - 16 rr > /dev/null
float（32 ビット画素値）の R[I,R] 画像を８ビットのものに変換% t2t̲float ri - 8 ri > /dev/null% t2t̲float rr - 8 rr > /dev/null
これらの処理を行うと R[I,R] 画像の容量が激減するが、それを用いた画像再構成の結果が float の画像を使ったものと完全に同じになるかどうかは定かではない。
[6] R[I,R] 画像上の物体像を囲む bounding box の選択
ImageJ のマクロ tlrb を使って R[I,R] 画像のいずれかを表示し、それらの上の物体像を完全にカバーする bounding box のひろがりを調べる：% IJ -macro tlrb rr 1 10 tlrb.log
ただし、ここでは物体像の境界の識別が容易な RR 画像（ディレクトリ rr）を指定した（RI 画像を使いたければディレクトリを ri に置き換えれば良い）。また、処理時間短縮のため RR 画像のうち１番目から 10 枚おきの view だけを表示するよう指定した（すべての view を表示したければ「1 10」を削除すれば良い）。
画像下部のスライダ使えば指定した複数の view の画像のそれぞれを表示できる。
マウス左ボタンをドラッグしてそれらの画像の上に bounding box の長方形を描く。
その後、マウスの右ボタンを押すと bounding box の情報を示すウィンドウが開くので、その下部に並んでいる以下の３つのボタンのいずれかをクリックする：
YesNo

これらの画像や同時に表示されるラインプロファイルを観察してもとの RP 画像の「山」の高さが 2^16 = 65536 に近い「飽和した値」か？
Wm と Ws の値ごとの MF の結果のプロファイルに「不自然さ」がないか？
を確認し、MF の処理の実行の是非やその処理で用いる Wm と Ws の値を決める。
なお、rp+mf̲ij.csh は MF の処理を行うたびにもとの RP 画像が保持している画素値の値域を上記のように（標準エラー出力に）書き出す。これらのうちの最大値から「山」の高さが「飽和した値」になっているか否かを判断できる。
[4] R[I,R] 画像の作成
新しい２個のディレクトリをまず作成し、それらの中に a.log の値を指定してプログラム his2raw で抽出した a.HIS の上の RP 画像、もしくはそれに MFを行った結果の RP 画像のいずれかから計算した R[I,R] 画像を書き込む：% mkdir ri rr% his2raw a.HIS 551 601 100
ri rr > raw.log% his2raw a.HIS 551 601 100 5 5 ri rr > raw.log% his2raw a.HIS 551 601 100 10 10 ri rr > raw.log

した測定 170516c ではもとの RP 画像上の「山」が「飽和した値」だったので、
Wm = Ws = 10 を指定した上記の３番目の行の入力によって MF を行った R[I,R]
画像を作成した。

# MF を行わない場合# Wm = Ws = 5# Wm = Ws = 10

言うまでもないことだが、上記の３通りのいずれか、もしくは Wm と Ws に任意の値を指定した起動法で his2raw を１度だけ実行すれば良い。なお、ここで例と表示されている bounding box の情報を記録した後に処理を続行する。

表示されている bounding box の情報を記録せずに処理を続行する。
Cancelマクロ tlrb による bounding box の選択の処理を終了する。
上記の Yes により bounding box の情報がファイル tlrb.log にも書き込まれる。
Cancel（および ImageJ の終了）を行った後にそれを確認すれば良い：% tail -1 tlrb.log
100
50

120

470

これら４個の数値は R[I,R] 画像の上の物体像に関する以下の情報を示しており、後述するプログラム "r[i,r]*" のそれぞれに起動パラメータとして指定する：
top = 100
物体像の上部の空気の領域の幅、もしくは、物体像の上端の座標値
left = 50
物体像の左側の空気の領域の幅
right = 120
物体像の右側の空気の領域の幅
bottom = 470
物体像の下端の座標値
なお、以下のようにすれば選択した bounding box を確認することができる：

% IJ -macro tlrb ri 1 10 100 50 120 470% IJ -macro tlrb rr 1 10 100 50 120 470[7] スライスごとの値の「時間変動」や鉛直方向のサンプルの「動き」のチェックプログラム r[i,r]2m[p,r] に top、left と right の値を指定して R[I,R] 画像
から平均投影値（MP）と平均屈折量（MR）の画像を作成する：% ri2mp ri - 100 50 120 mp.tif-0.019094
0.119755% rr2mr rr - 100 50 120 mr.tif-1.014936
0.893499
M[P,R] 画像の横軸は view（左右端がそれぞれサンプル回転角０度と 180 度）に、また、縦軸は R[I,R] 画像の縦軸に対応している。そして、M[P,R] 画像の画素値はtop、left と right の値が指す空気の領域の画素の値を用いて
R[I,R] 画像のそれぞれの画素の値から計算した
X 線投影値（XP）と X 線屈折量（XR）の画像のそれぞれの上の横一列の画素の値の平均値である。これらの値が view に関して大きく変動している（厳密には一定値になっていない）と高画質の CT 画像を再構成することができないので、M[P,R] を眺めてその様相を調べる：% xv mp.tif mr.tif
M[P,R] 画像の上の「横方向に伸びている一定値の縞模様」が傾斜している場合は
SIXM の測定中にサンプル（もしくはサンプルが載っているステージ）が鉛直方向に動いた可能性が高い。view に関して定常的（steady vertical creep）だと仮定
できるなら、その動きを画像再構成の際に補正できる。ただし、そのためには回転角が０度と 180 度の X[P,R] 画像の上のサンプル像の鉛直方向の変位量（縦ズレ）の値 Dz を指定する必要がある。
後述するプログラム r[i,r]2rc を使えば値 Dz を推定できる。また、オプションの起動パラメータ（/dev/null）の指定により r[i,r]2m[p,r] でも自身が作成した
M[P,R] 画像から値 Dz を推定できる。下記の r[i,r]2m[p,r] のそれぞれの出力の２行目の最初の値（今の場合はどちらも -2）が Dz の推定値である：% ri2mp ri - 100 50 120 mp.tif /dev/null-0.019094
0.119755-2
4.369533e-03% rr2mr rr - 100 50 120 mr.tif /dev/null-1.014936
0.893499-2
3.253222e-02[8] サンプル回転軸の位置の推定
サンプル回転角が０度と 180 度の X[P,R] 画像の位置関係を変えながらそれらに含まれているサンプル像を比較することにより、SIXM の測定時のサンプル回転軸の位置（いわゆる「センター値」）を推定する：% ri2rc ri - 100 50 120 100 470 0.25 0.25 ri2rc.tif
206
343
241
-124
124
-3
% rr2rc rr - 100 50 120 100 470 0.25 0.25 rr2rc.tif
206
343
241
-124
124
-6
物体像の領域だけに限定すれば０度と 180 度の画像の比較の精度が向上するので、ここでは空気の領域の拡がりを指す top、left と right の値に続けて物体像が占めている縦方向の範囲を指す top と bottom の値をプログラム r[i,r]2rc の起動時に指定した。その後の２個の値「0.25 0.25」は２画像の位置関係を変えるために行う横と縦方向の並進移動のそれぞれの総量の範囲を指す、画像の横と縦の幅を１とした相対値である。
上記の r[i,r]2rc のそれぞれの２行の出力の最初の２個の値が指定した相対値に対応した回転軸の位置と縦ズレの調査範囲の実際の値域である。そして、それらの後の各行の最後の値が回転軸の位置（１行目；センター値）と SIXM の測定中に生じたサンプルの縦ズレの推定値（２行目；前述の Dz）である。
測定 170516c では X[P,R] 画像から推定した回転軸の位置はどちらも同じ値 241であり、それを吸収と位相 CT の両方の画像再構成で指定すれば良いはずである。
なお、後述する画像再構成のテストの結果にもよるが、原理的には吸収と位相 CTの画像再構成の際に指定するサンプル回転軸の位置は同じ値にすべきである。
測定 170516c では r[i,r]2rc のそれぞれで推定した２個の縦ズレ Dz は同じ値になっていない。さらに、それらは前述の r[i,r]2m[p,r] で推定した Dz の値とも異なっている。つまり、ここで例として処理した 170516c のサンプルは測定中に鉛直方向に少しだけ動いたようだが、その動きを steady vertical creep だと仮定して補正するのは難しい。
なお、r[i,r]2rc で推定した回転軸の位置や縦ズレの値が指定した範囲の両端の値に一致している場合、それらの推定値は信頼できない。その場合には物体像が占める縦方向の範囲を変えて r[i,r]2rc を再実行してみるべきである。また、内部構造が特異（例えば、内部が一様）な物体像に対して回転軸の位置や縦ズレの値の推定処理を正常に行えないこともある。その場合は r[i,r]2rc の起動時に指定したファイル r[i,r]2rc.tif に格納されている「RMSD のマップ」の画像を観察して対応策を考える必要があるが、ここではその詳細の説明は省略する：% xv
ri2rc.tif rr2rc.tif% xv -preset 4 ri2rc.tif rr2rc.tif# RMSD のマップを白黒で表示する# 擬似カラー表示（低い値が青色）[9] サンプル回転軸の位置を変えた画像再構成のテスト
念のため、先に推定した値（241）の前後のサンプル回転軸の位置を指定した画像
再構成のテストを行う。ここでは以下のパラメータを指定する：
再構成するスライスの位置 = ( top + bottom ) / 2 = 285
テストする回転軸の位置の最小値 = 241 -0.5 * 7 = 237.5
テストする回転軸の位置の間隔 = 0.5
テストする回転軸の位置の総数 = 1（推定値）+ 7*2（その前後）= 15
画像再構成のテストは CPU で走る r[i,r]̲stg̲t と GPU を使う r[i,r]̲stg̲g のどちらのプログラムで行っても良い（以下では両者の実行法を併記する）。ただし、
CPU 用を使う場合は実行前に稼働するスレッド数の設定を行っておく必要がある：% setenv CBP̲THREADS 8# ８個のスレッドで並列に画像再構成を行う

吸収と位相 CT の画像再構成のテスト用のディレクトリ test̲[a,p] のそれぞれを作成しておく。ここまでで説明したパラメータの値を指定して r[i,r]̲stg̲[t,g]を実行すれば、それぞれのディレクトリの中に指定したスライス位置の１個の sinogram と回転軸の位置を変えた複数個の tomogramsの２種類の画像を作成できる：% mkdir test̲a% ri̲stg̲t ri - 100 50 120 285 100e-7% ri̲stg̲g ri - 100 50 120 285 100e-7

237.5 0.5 15 +0 test̲a > test̲a.log
237.5 0.5 15 +0 test̲a > test̲a.log% mkdir test̲p% rr̲stg̲t rr - 100 50 120 285 6220e-3 -104 237.5 0.5 15 +0 test̲p > test̲p.log% rr̲stg̲g rr - 100 50 120 285 6220e-3 -104 237.5 0.5 15 +0 test̲p > test̲p.log

なお、steady vertical creep の補正を行う場合は先に推定した縦ズレの値 Dz
をサンプル回転角の初期値（+0）の直前に指定すれば良い（ただし、ここで例とした測定 170516c に対しては steady vertical creep を仮定できないので不要）：% mkdir test̲a% ri̲stg̲t ri - 100 50 120 285 100e-7% ri̲stg̲g ri - 100 50 120 285 100e-7

237.5 0.5 15 Dz +0 test̲a > test̲a.log
237.5 0.5 15 Dz +0 test̲a > test̲a.log% mkdir test̲p% rr̲stg̲t rr - 100 50 120 285 6220e-3 -104 237.5 0.5 15 Dz +0 test̲p > test̲p.log% rr̲stg̲g rr - 100 50 120 285 6220e-3 -104 237.5 0.5 15 Dz +0 test̲p > test̲p.log
その後、まずはファイル「test̲[a,p]/スライス番号.tif」の sinogram の画像を観察して測定中のサンプルの不穏な「動き」の有無を確認する：
xv test̲a/285.tif test̲p/285.tif# sinograms of layer 285

ファイル「test̲[a,p]/スライス番号̲回転軸の位置.tif」が回転軸の位置を変えて再構成した画像である。これらを個別に眺めて回転軸の位置の違いに伴う CT 画像の画質を評価しても良いが、それよりも回転軸の位置が異なる画像を並べたものを眺める方が効率的である。以下のようにすれば吸収と位相 CT のそれぞれのテスト用のディレクトリの下にある再構成画像すべてを並べた画像を作成できる：
横方向優先の順番で５×３のマトリックス状に並べる場合% pile̲h.csh test̲a 1 1 15 5 3 test̲a.tif% pile̲h.csh test̲p 1 1 15 5 3 test̲p.tif
縦方向優先の順番で３×５のマトリックス状に並べる場合% pile̲v.csh test̲a 1 1 15 3 5 test̲a.tif% pile̲v.csh test̲p 1 1 15 3 5 test̲p.tif
これらの画像を詳細に観察して、偽像のない高画質な CT 画像の再構成に最適なサンプル回転軸の位置を決定する：
xv test̲a.tif test̲p.tif

なお、steady vertical creep の補正を行う場合は先に推定した縦ズレの値 Dz
をサンプル回転角の初期値（+0）の直前に指定すれば良い（ただし、この補正はプログラム r[i,r]2tg では実行できない；また、ここで例とした測定 170516c
に対しては steady vertical creep を仮定できないので、この補正は不要）：% mkdir tg̲a% ri̲tg̲t ri - 100 50 120 100e-7% ri̲tg̲g ri - 100 50 120 100e-7

241.0 Dz +0 tg̲a > tg̲a.log
241.0 Dz +0 tg̲a > tg̲a.log% mkdir tg̲p
% rr̲tg̲t rr - 100 50 120 6220e-3 -104 241.0 Dz +0 tg̲p > tg̲p.log% rr̲tg̲g rr - 100 50 120 6220e-3 -104 241.0 Dz +0 tg̲p > tg̲p.log
再構成した画像をざっと見るため、テストのものと同様に吸収と位相 CT のそれ
ぞれの等間隔になるように選び出した 15 枚のスライスを５×３もしくは３×５のマトリックス状に並べた画像を作成してみる。
スライスの間隔 = ( bottom - top + 1 ) / 15 ≒ 24
最初のスライス = top + ( bottom - top + 1 - 24 * ( 15 - 1 )) / 2 ≒ 117
横方向優先の順番で５×３のマトリックス状に並べる場合% pile̲h.csh tg̲a 117 24 15 5 3 tg̲a.tif% pile̲h.csh tg̲p 117 24 15 5 3 tg̲p.tif
縦方向優先の順番で３×５のマトリックス状に並べる場合% pile̲v.csh tg̲a 117 24 15 3 5 tg̲a.tif% pile̲v.csh tg̲p 117 24 15 3 5 tg̲p.tif
これらの画像を眺めて画像再構成の結果に不備がないかを確認する：% xv tg̲a.tif tg̲p.tif[11] 吸収 CT 画像に対する後処理[10] 吸収および位相 CT 画像の再構成とそれらの確認
サンプル回転軸の位置を決定した後（ここでは吸収と位相 CT のどちらに対しても同じ値 241.0 とした）、本番の画像再構成を行う。テスト用の場合と同様に CPU
用と GPU 用の合計２×３個のプログラムのいずれかを利用できる（以下ではそれらの実行法を併記する）。ただし、CPU 用のプログラムを実行する場合には稼働するスレッド数（８とする）を事前に設定しておく必要がある：
プログラム r[i,r]2tg を実行する場合% setenv THREADS 8
プログラム r[i,r]̲tg̲t を実行する場合% setenv CBP̲THREADS 8
テストの場合と同様に吸収と位相 CT の再構成画像用のそれぞれのディレクトリを作成した後、これまでに紹介したパラメータの値を指定して画像再構成プログラム
r[i,r]2tg もしくは r[i,r]̲tg̲[t,g] を実行する：% mkdir tg̲a% ri2tgri - 100 50 120 100e-7% ri̲tg̲t ri - 100 50 120 100e-7% ri̲tg̲g ri - 100 50 120 100e-7% rr̲tg̲t rr - 100 50 120 6220e-3 -104 241.0 +0 tg̲p > tg̲p.log% rr̲tg̲g rr - 100 50 120 6220e-3 -104 241.0 +0 tg̲p > tg̲p.log

241.0 +0 tg̲a > tg̲a.log
241.0 +0 tg̲a > tg̲a.log
241.0 +0 tg̲a > tg̲a.log% mkdir tg̲p
% rr2tg
rr - 100 50 120 6220e-3 -104 241.0 +0 tg̲p > tg̲p.log

再構成した浮動小数点数画素値の CT 画像を画像処理に適した整数画素値の画像に変換する。そのためにまず、再構成した LAC の値（CT 値）の値域を調べる：% mm.csh tg̲a-375.003174

1699.017944

SIXM に限らずこれまでの X 線 CT の測定ではすべてのスライスで上記のようなCT 値の値域全体を 16 ビット画素値にマッピングした word 画像を作成していた。
その作成時に指定が必要な 16 ビット画素値１階調ごとの CT 値の増分の値(1699.017944-(-375.003174))/(2^16-1) = 0.0316475336537728
も先と同じ C-shell script "mm.csh" で計算できる：% mm.csh tg̲a 16-375.003174

1699.017944

0.0316475336537728

このようにして得た値を指定して吸収 CT 画像の word 画像とその画素値ヒスト
グラムのデータファイルを作成する：% mkdir word̲a% t2t̲f.csh tg̲a-375.003174 0.0316475336537728 16 word̲a > word̲a.txt

また、word 画像を作成せずにその画素値ヒストグラムだけを得ることもできる：% t2t̲float tg̲a - -375.003174 0.0316475336537728 16 -> word̲a.txt

その後、再構成した CT 値の全域のヒストグラムを表示してサンプルの大部分を占める「主要物質」の CT 値を調べる：% plot̲x11.csh word̲a.txt# CT 値 = 0 〜 500 の部分を表示表示した画素値ヒストグラムからサンプルの主要物質に相当するピークの CT 値
を読み取り、その 1/100（この値に理論的根拠なし）の値を後述する byte 画像の８ビット画素値１階調ごとの CT 値の増分とする。
測定 170516c の吸収 CT 画像のヒストグラムのピークの CT 値 ≒ 124 (1/cm)→ byte 画像の８ビット画素値１階調ごとの CT 値の増分 = 1.24 (1/cm)
このようにして決めた CT 値の増分の値を指定して CT 画像の byte 画像とその画素値ヒストグラムのデータファイルを作成する。
% mkdir byte̲a% t2t̲f.csh tg̲a 0 1.24 8 byte̲a > byte̲a.txt
先と同様に byte 画像の等間隔 15 枚のスライス画像を５×３もしくは３×５のマトリックス状に並べた画像を作成してみる：
横方向優先の順番で５×３のマトリックス状に並べる場合% pile̲h.csh byte̲a 117 24 15 5 3 byte̲a.tif
縦方向優先の順番で３×５のマトリックス状に並べる場合% pile̲v.csh byte̲a 117 24 15 3 5 byte̲a.tif
また、byte 画像の画素値ヒストグラムの PostScript の図も作成する：% plot̲a.csh ./byte̲a.txt > byte̲a.ps
なお、byte 画像の画素値０と 255 には「ゴミ」が混じっているので、画素値０に対応する CT 値の上限の値 = 1.24/2 = 0.62
画素値 255 に対応する CT 値の下限の値 = 1.24*254.5 = 315.58を除外した範囲だけを描けばヒストグラムの見栄えが良くなる：
plot̲a.csh ./byte̲a.txt 0.62 315.58 > byte̲a.ps
このようにして作成した画像と図の出来栄えを確認する：% xv byte̲a.tif% evince byte̲a.ps% plot̲x11.csh word̲p.txt 0 15

# CT 値 = 0 〜 150 の部分を表示

測定 170516c の位相 CT 画像のヒストグラムのピークの CT 値 ≒ 7.52 (1e-6)→ byte 画像の８ビット画素値１階調ごとの CT 値の増分 = 0.0752 (1e-6)% mkdir byte̲p% t2t̲f.csh tg̲p 0 0.0752 8 byte̲p > byte̲p.txt
横方向優先の順番で５×３のマトリックス状に並べる場合% pile̲h.csh byte̲p 117 24 15 5 3 byte̲p.tif
縦方向優先の順番で３×５のマトリックス状に並べる場合% pile̲v.csh byte̲p 117 24 15 3 5 byte̲p.tif% plot̲p.csh ./byte̲p.txt > byte̲p.ps
byte 画像の画素値０に対応する CT 値の上限 = 0.0752/2 = 0.0376
byte 画像の画素値 255 に対応する CT 値の下限 = 0.0752*254.5 = 19.1384% plot̲p.csh ./byte̲p.txt 0.0376 19.1384 > byte̲p.ps% xv byte̲p.tif% evince byte̲p.ps[13] 長さと表示輝度のスケールバーを付けた吸収と位相 CT 画像の作成（オマケ）
吸収と位相 CT 画像の両方に使える長さが 25 μm のスケールバーの画像の作成：% bar̲ls.csh 250 5 50 '' '25 um' ls.tif
吸収と位相 CT のそれぞれの byte 画像の値に応じた表示輝度のスケールバー：% bar̲gs.csh 750 25 50 0 "̀echo 1.24\*255 | bc -l̀ (1/cm)"bar̲a.tif% bar̲gs.csh 750 25 50 0 "̀echo 0.0752\*255 | bc -l̀ (1e-6)" bar̲p.tif
スケールバーの画像と吸収 CT の byte 画像を並べた画像を張り合わせる：% pile̲gray 2 1 -250 0
255 bar̲a.tifbar̲a.tif ls.tif% pile̲gray 1 2 0-50 255 browse̲a.tif byte̲a.tif bar̲a.tif
スケールバーの画像と位相 CT の byte 画像を並べた画像を張り合わせる：[12] 位相 CT 画像に対する後処理
位相 CT 画像に対する後処理の手順は吸収 CT 用のものと概ね同じである：% mm.csh tg̲p-23.171337> word̲p.txt% plot̲x11.csh word̲p.txt

その際に CT 値の範囲を限定したければ以下のように指定すれば良い：% plot̲x11.csh word̲a.txt 0 500% t2t̲float tg̲p - -23.171337 0.0007342873426413 16 -24.950184% pile̲gray 2 1 -250 0
255 bar̲p.tifbar̲p.tif ls.tif% pile̲gray 1 2 0-50 255 browse̲p.tif byte̲p.tif bar̲p.tif
作成した吸収と位相 CT の browse 画像の確認：% xv browse̲a.tif browse̲p.tif

16 ビット画素値１階調ごとの CT 値の増分の値= (24.950184-(-23.171337))/(2^16-1) = 0.0007342873426413% mm.csh tg̲p 16-23.171337

24.950184% mkdir word̲p% t2t̲f.csh tg̲p-23.171337 0.0007342873426413 16 word̲p > word̲p.txt

0.0007342873426413