1 / 93
Date:Wed, 21 Sep 2016 20:12:18 +0900

From:

Tsukasa NAKANO

To:

Kentaro UESUGICc:

MATSUNO Junya, Aiko Nakato, Sugimoto Miyama, 中村隆太, Akira Kitayama,⼟`⼭ 明, Hashiguchi Minako, ⽵内, 星野, 上椙Subject: oct_1

うえすぎさま、
GSJ/AIST のなかのです。昨⽇までの SPring-8 実験ではお世話になりました。実験の場で話したように、
offset CT のサンプル回転⾓ごとの投影画像の位置決めや画像再構成などの⼀連の処理を⾏うためのプログ
ラム群がとりあえず完成しました。そのソースコードなどの⼀式を以下の２個の書庫ファイル（同内容です）に⼊れておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zip
これらには GNU C-Compiler（gcc）でコンパイルした Windows ⽤の実⾏ファイル "*.exe" も⼊れてあります。ファイル名が "*_64.exe" のものは 64 bit Windows ⽤の実⾏ファイル、それ以外は 32 bit ⽤です。また、書庫ファイルを展開したディレクトリ oct/cuda/ には CUDA GPU（CBP engine for CUDA）を⽤いて offset
CT の画像再構成を⾼速に⾏うプログラム oct_sg2tg_g のソースコードの類も⼊れてあります。ただし、ぼくはその Windows の上でのコンパイルの仕⽅を知りません。Linux ⽤に書いた oct/cuda/Makefile の記述などを参考にして Windows ⽤の実⾏ファイルをコンパイルして下さい。
プログラム oct_sg2tg_g は CBP engine for CUDA の実⾏と再構成画像のファイル書き込みを別スレッ
ドで⾏います。oct_sg2tg_g のコードは WinAPI（から利⽤可能な Windows Thread）と pthread（POSIX
thread）のどちらを使ってその処理を⾏うかをコンパイル時に⾃分で決めるようになっていますが、
pthread を使う場合にはコンパイル時にそのライブラリ（のファイル）を指定する必要があります。
oct/cuda/Makefile の以下の部分がその設定で、WinAPI を使う場合には不要です。
LP= -lpthread

なお、こちらで使っている MinGW（Minimalist GNU for Windows）の gcc によるコンパイルでは WinAPIのライブラリ指定も不要です。
書庫ファイル oct.taz などに⼊れたプログラム群のインストール法をファイル oct/install.txt に記しておきま
した。プログラムそれぞれの使⽤法については説明書（oct/usage.txt）が完成した後に改めて紹介します。
ここではそれらのプログラムを⽤いた offset CT の位置決めおよび画像再構成のための必要最⼩限の処理の⼿順と、それに要した処理時間について書きます。

2 / 93(1) offset CT のペアの投影画像の位置決め
offset CT の測定データから得られるサンプルを 360°回転しながら撮影した投影画像 F(θ,x,y)から通常の CT の測定データから得られるものと等価なサンプルを 180°回転しながら撮影した投影画像 G(θ,x,y)を合成すれば、通常の CT の時と同じ⽅法で画像 G に対して画像再構成処理を⾏うことにより画像 F に対する CT 画像を得ることができます。画像 F からの画像 G の合成法は以下の通りです（注：ここでは Ox≧0 を仮定しています）。
G(θ,x,y)= F(θ,x,y)０≦ x ＜Ox= { F(θ,x,y) + H(θ,x - Ox,y) } / 2

Ox≦ x ＜Nx

= H(θ,x - Ox,y)Nx≦ x ＜Nx+Ox

θ：サンプル回転⾓[x,y]：投影画像の座標値
Nx：測定した投影画像 F の横画素数
Ox：180°異なる回転⾓の投影画像 F のペアの x ⽅向のオフセットH(θ,x,y) = F(θ + 180°,Nx - 1 - x,y)：回転⾓が 180°以上の投影画像 F を x ⽅向に反転した画像
この画像 G の合成の式では Ox≦x＜Nx の範囲において画像 F と「画像 H を x ⽅向に Ox だけ平⾏移動した画像」が完全に⼀致すると仮定しています。
F(θ,x,y) == H(θ,x - Ox,y)
しかし、これは実際の測定データには不適当なので、２画像が重なり合っている領域において同じ位置にある画素の投影値の「差の⾃乗の平均値（mean square difference、MSD）」が最⼩になっている Ox の値を探し出すことにしました。また、この MSD を⽤いた２画像の位置決めの処理では y ⽅向のオフセットの値 Oyも探し出します。これにより２画像の「縦ズレ量」を検出できますが、厳密には Oy == 0 でないと offset CTの画像再構成はできません。
offset CT の位置決め処理ではまずサンプル回転⾓ごとの MSD の「マップ」
MSD(θ,Ox,Oy) = <{ F(θ,x,y) - H(θ,x - Ox,y-Oy) }^2>_x,yを作成します。これが最⼩値になっている Ox(θ)と Oy(θ)の値から回転⾓θごとの「サンプルの動き」を解析できます（その詳細は次便以降の E-mail で紹介します）。そのためにプログラム oct_cmp_180 を書きました が 、 こ れ は 使 い 勝 ⼿ が や や 悪 い の で 、 offset CT の 画 像 再 構 成 を と り あ え ず ⾏ な い た い 場 合 に はMSD(θ,Ox,Oy) を回転⾓について平均したマップ
MSD(Ox,Oy) = <MSD(θ,Ox,Oy)>_θ
を作成し、それが最⼩値になっている Ox と Oy の値を探し出すプログラム oct_xy をお使い下さい。

3 / 93以前に何度か紹介した画像ペアの位置決め処理では RMSD（root MSD；MSD の平⽅根）を使っていましたが、それらは単純に平均できないので、今回のプログラムでは MSD を採⽤しました。
プログラム oct_xy は複数の投影画像のペアの MSD のマップをマルチスレッドで並列に作成できます。その起動前に稼動するスレッドの個数を環境変数 THREADS に設定して下さい（この指定がない場合には１個のスレッドだけを稼動します）。
setenv THREADS

8# tcshexport THREADS=8# bashset THREADS=8# Windows のコマンドプロンプト

その後、以下の⾏⼊⼒によって oct_xy を起動することができます。
oct_xy HiPic/
→

Ox

Ox1

Ox2
Oy

Oy1

Oy2 {MSD.tif}

MSD> xy.log# 標準エラー出⼒

HiPic/
offset CT の測定データファイル（ログファイル output.log と HiPic 画像 dark.img および "q*.img"）が⼊っているディレクトリの名前
Ox1 と Ox2
Ox の探索範囲。これらは -Nx＜Ox1≦Ox2＜Nx を満たす整数値でなければならない。
Oy1 と Oy2
Oy の探索範囲。Ny を offset CT で測定した画像 F の縦画素数として、これらは -Ny＜Oy1≦Oy2＜Ny を満たす整数値でなければならない。
MSD.tif
MSD(Ox,Oy) の マ ッ プ を 格 納 す る float TIFF の 画 像 フ ァ イ ル の 名 前 。 指 定 を 省 略 す る と
MSD(Ox,Oy) のマップの画像を出⼒しない。
xy.log標準出⼒に書き出される「サンプル回転⾓それぞれに対する下記の４個の値を並べた⾏」を書き込むテキストファイルの名前。
[1] サンプル回転⾓θに対応する投影番号（０〜）[2] MSD(θ,Ox,Oy) が最⼩値になっている Ox(θ) の値[3] MSD(θ,Ox,Oy) が最⼩値になっている Oy(θ) の値[4] MSD(θ,Ox,Oy) の最⼩値
Ox、Oy と MSD
標 準 エ ラ ー 出 ⼒ に 出 ⼒ さ れ る MSD(Ox,Oy) が 最 ⼩ 値 に な っ て い る Ox と Oy の 値 お よ び
MSD(Ox,Oy) の最⼩値。
７⽉に SPring-8 BL47XU で⾏った FZP-CT 実験で得た offset CT の４測定のデータを oct_xy で処理してみました。ただし、下記の Nθはサンプル回転の総ステップ数（投影数）、Dr は nm 単位の画素の辺⻑です。

4 / 93
測定

Nx

Ny

Nθ

Dr

note

160723[k,l]

1292

1032

3600

35.9

7 keV

160723[n,o] 964

1044

3600

31.1 41.1

8 keV

Ox と Oy の探索範囲として以下の値（適当に決めた値です）を指定しました。
160723[k,l]160723[n,o]Ox1 = Nx - 255

1067

739

Ox2 = Nx - 25

1267

939

Oy1 = -50
Oy2 = +50
oct_xy が標準エラー出⼒に書き出した Ox と Oy の値は以下の通りでした。
測定

Ox

Oy

160723k1169-4160723l1160-5160723n870-3160723o881-2

なお、oct_xy が出⼒した Ox と Oy が指定した範囲の端に⼀致する値の場合は MSD(Ox,Oy) の最⼩値の検索に失敗しています。Ox と Oy の範囲の指定を変えて再挑戦して下さい。
(2) offset CT の画像の再構成
前記のオフセットの値 Ox と Oy に従ってサンプル回転⾓が 180°異なる投影画像 F から合成した画像 G の画素数とその総数は以下の値になります。
画像 G の x ⽅向の画素数：nx = Nx + | Ox |y ⽅向の画素数：ny = Ny - | Oy |総数（投影数）：nθ = Nθ / 2
また、少し考えればわかるように、サンプル回転軸は回転⾓によらず画像 G の x ⽅向中央（x = (nx - 1)/2）に位置するはずです。ここで紹介する３個の画像再構成プログラムはこれを仮定しているので、「センター値」の指定は不要です。
offset CT ⽤に以下の３個の画像再構成プログラムを書きました。これらはいずれも測定画像を読み込みながら投影画像 F の作成と画像 G の合成を⾏うので、以下ではこの処理を「投影画像の読み込み」と呼ぶことにします。

5 / 93
プログラム oct_tg
サンプル回転⾓ごとの投影画像を読み込むたびにそれに対する畳み込みと逆投影演算を⾏う。その処理には再構成画像と同量の作業⽤の計算機メモリが必要。現在の設定では再構成画像の各画素の値を単精度浮動⼩数点数（float；４バイト⻑）として保持しているので、作業⽤のメモリとして最⼩限必要な量は「４× nx^2 × ny」バイト。
oct_sg2tg_[t,g]oct_sg2tg_t はマルチスレッド処理により、また、oct_sg2tg_g は CUDA GPU によってスライスごとに画像再構成を⾏う。投影画像すべてを読み込んだ後に CBP engine にスライスごとの投影値
sinogram のデータを渡すので、それに応じた作業⽤メモリが必要。投影値を float として保持している現在の設定ではその量は「４× nx × nθ × ny」バイト。
つまり、offset CT ⽤の画像再構成プログラムは多量の計算機メモリを要します。上記のバイト数のメモリを載せていない計算機を使う場合は再構成するスライスの範囲を限定して下さい（これは ny の値を⼩さくすることに相当します）。
プログラム oct_tg はサンプル回転⾓ごとの投影画像の読み込みと畳み込み・逆投影演算をマルチスレッドで並列に実⾏できます。その起動前に oct_xy と同様にして稼動するスレッド数を環境変数 THREADS に設定して下さい（その指定を省略すると畳み込み・逆投影演算を実⾏するスレッドは１個になります）。また、
oct_sg2tg_t では CBP engine（for CPU）がマルチスレッドで処理を⾏うので、そのためのスレッド数を事
前に環境変数 CBP_THREADS に設定して下さい（その設定法は THREADS の場合と同じです）。なお、
oct_sg2tg_[t,g] は画像再構成処理と再構成画像のファイル書き込み処理を別スレッドで⾏いますが、それに関するスレッド数の指定は不要です。
さて、offset CT ⽤の画像再構成プログラムの起動法ですが、それは３個のプログラムで概ね同じです（oct_tgと oct_sg2tg_[t,g] で起動パラメータが少しだけ異なりますが）。２通りの起動法があります。まず、offset
CT の測定画像のファイルが⼊ったディレクトリの名前（HiPic/）と投影画像のペアのオフセット Ox および
Oy の値だけを指定して起動すると、標準エラー出⼒に合成した投影画像 G の画素数などを書き出します。
oct_tg HiPic/ Ox Oy
oct_sg2tg_[t,g]→nxHiPic/ Ox
ny

Oy
nθ

# 標準エラー出⼒

これらの値から処理に必要な作業⽤メモリの量を計算し、使⽤する計算機の搭載メモリの量に応じて再構成するスライスの範囲を決めて下さい。その後、以下の⼊⼒によって offset CT の画像再構成を実⾏できます。
oct_tg HiPic/Ox Oy Dr RA0

TG_format > tg.log

oct_tg HiPic/ Ox Oy layer1 layer2 Dr RA0

TG_format > tg.log

oct_sg2tg_[t,g]HiPic/ Ox

Oy rangeList

Dr RA0

TG_format > tg.log

oct_sg2tg_[t,g]HiPic/ Ox

Oy rangeList

Dr RA0BPS TG_format > tg.log

oct_sg2tg_[t,g]HiPic/ Ox

Oy rangeList

Dr RA0base step

→

nx

ny

nθ

# 標準エラー出⼒

BPS TG_format > tg.log

6 / 93
HiPic/、 Ox と Oy
説明済み
layer1 と layer2（oct_tg のみ）
再構成を⾏いたいスライスの番号（０〜 ny-1）の最⼩値と最⼤値。これらの値の指定を省略するとすべてのスライスを再構成する。
rangeList（oct_sg2tg_[t,g] のみ）
再構成すべきスライスの番号の範囲を指定する「range list」。"-" を指定するとすべてのスライスを再構成する。
Dr
投影画像の画素（正⽅形と仮定）の辺⻑。再構成画像の画素（⽴⽅体と仮定）の辺⻑もこの値になる。なお、再構成した CT 値の単位はこの値の単位の逆数（cm 単位の辺⻑を指定すれば CT 値の単位は 1/cm）。
RA0
サンプル回転⾓のバイアス（単位は度）。通常は０を指定すれば良い。０でない値を指定するとスライス画像を画質の劣化なしに回転できる。
base と step および BPS（oct_sg2tg_[t,g] のみ）
再構成画像の各画素の float の CT 値 f を整数画素値 i に変換する際の下記の関係式の２個の係数値と i の上限の値を決めるビット数。
f = base + step × i
i = 0 〜 2^BPS - 1base と step の値の指定を省略するとスライスごとに CT 値 f の最⼩値と最⼤値で正規化した画素
値 i になるような base と step の値を使う。また、これら３個のパラメータ値の指定を省略するとCT 値 f をそのまま画素値として float TIFF の画像ファイルに書き込む。
TG_format
各スライスの再構成画像を書き込むファイルのパス名の決定に使⽤する C ⾔語の標準ライブラリ
関数 printf ⽤の format の⽂字列。スライス番号 z = 0 〜 ny - 1（整数値）に応じたファイル名にするには "%04d"（４桁以内のスライス番号の場合）を含む⽂字列を指定すれば良い。なお、ディレクトリ名を含むパス名にした場合には事前にそのディレクトリを作成しておくこと。
tg.log標準出⼒に書き出される「再構成したスライス画像それぞれに対する下記の３個の値を並べた⾏」
を書き込むテキストファイルの名前。
[1] スライス番号 z = 0 〜 ny - 1[2,3] そのスライスに含まれる CT 値の最⼩値と最⼤値
nx、ny と nθ
説明済み
なお、上の説明には陽に書きませんでしたが、oct_tg は再構成画像の各画素の CT 値をそのまま float TIFFの画像ファイルに書き込みます（oct_sg2tg_[t,g] とは異なり、oct_tg は普通の整数画素値の integer TIFF の画像ファイルを作成できません）。

7 / 93(3) 実測した処理時間
２台の Linux 計算機で７⽉の実験で得た４個の offset CT の測定の画像のそれぞれに対する位置決めと画像
再構成の処理を⾏いました。それに使⽤した４個のプログラムそれぞれの処理時間は以下の通りです。
使⽤した計算機の諸元
vrm
CPU：Intel Xeon E5-2609 v2（8 cores、64 GB memory）
GPU：Nvidia GeForce GTX 750 Ti （２ GB memory）
firmware：Scientific Linux 6.8 + gcc 4.4.7 + CUDA 7.5
test01CPU Intel Core i3 530（4 cores、16 GB memory）
GPU：Nvidia Tesla C2075（５ GB memory）
firmware：Scientific Linux 6.8 + gcc 4.4.7 + CUDA 6.0
処理時間（中 or 左列が vrm、右列が test01 による処理時間；単位は秒）
測定 160723k (1292×1032×3600 / 2461×1028×1800)
oct_xy324.198599

453.565869

oct_tg3214.311764

9049.451679

oct_sg2tg_t3523.193628

10204.277470

oct_sg2tg_g

937.934966

1093.643824

測定 160723l (1292×1032×3600 / 2452×1027×1800)
oct_xy378.926121

455.483427

oct_tg3427.715060

8958.444803

oct_sg2tg_t3617.878824

10100.010755

oct_sg2tg_g1004.150321

1065.524604

測定 160723n (964×1044×3600 / 1834×1041×1800)
oct_xy252.022873

373.306258

oct_tg2063.673255

4983.110979

oct_sg2tg_t2100.376438

5751.708006

oct_sg2tg_g

653.065814

646.277122

測定 160723o (964×1044×3600 / 1845×1042×1800)
oct_xy251.760916

371.965420

oct_tg1916.112358

5043.500624

oct_sg2tg_t2160.073730

5822.809195

oct_sg2tg_g693.666390

675.303937

注１２台の計算機のどちらでもすべての処理の前に環境変数 THREADS と CBP_THREADS の両⽅に使⽤するスレッド数として値８を設定しました。

8 / 93
注２
計算機 test01 ではメモリ不⾜のため各測定のスライスすべてを⼀度に再構成できません。oct_tg とoct_sg2tg_[t,g] のいずれによる画像再構成でも測定ごとに全体の半数づつのスライスを指定した
合計２回の処理を⾏いました。そのため、test01 による画像再構成の処理時間は offset CT の投影
画像の２回⽬の読み込み処理の分だけ増えています。なお、上に記した test01 の画像再構成の処理
時間は２回起動したプログラムの処理時間の和です。
注３
上記の測定番号の後の括弧の中に記した "/" の前後の⽂字列はそれぞれ offset CT の投影画像 F とそれから合成した画像 G の「横画素数×縦画素数×投影数（[N,n]x×[N,n]y×[N,n]θ）」です。
これらより、ある程度⾼性能な GPU 搭載の計算機であれば、７⽉の測定のものと同程度のサイズの offset
CT の画像の再構成処理を oct_xy と oct_sg2tg_g を⽤いて 30 分以内に実⾏可能なことがわかりました。
⻑い E-mail になりました。とりあえず以上です。

Date:

Fri, 23 Sep 2016 17:20:16 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya , Aiko Nakato , Sugimoto Miyama , 中村隆太 , Akira Kitayama ,⼟`⼭ 明 , Hashiguchi Minako , ⽵内, 星野, 上椙Subject: oct_2

うえすぎさま、
GSJ/AIST のなかのです。offset CT に関する⼀昨⽇の E-mail の補⾜です。
(1)
前の E-mail ではプログラム oct_xy の使⽤法について不親切なことを書きました：
On Wed, 21 Sep 2016 20:12:18 +0900 Tsukasa NAKANO wrote:> なお、oct_xy が出⼒した Ox と Oy が指定した範囲の端に⼀致する値の場合は> MSD(Ox,Oy) の最⼩値の検索に失敗しています。Ox と Oy の範囲の指定を変えて> 再挑戦して下さい。
oct_xy による検索が失敗した場合は oct_xy で MSD(Ox,Oy) の画像 MSD.tif を作成し、それを以下のようにして検討して下さい：例えば、この E-mail に添付した 1607_oct_xy.pdf を御覧下さい。７⽉の測定
160723[k,l,n,o] のそれぞれに対する oct_xy で作成した画像 MSD.tif の内容をグレースケールと等⾼線で⽰しました。MSD のマップはこのような卵型のプロファイルなので、その中央部をカバーする Ox[1,2]Oy[1,2] の値を oct_xy に再指定して下さい。

や9 / 93(2)
ぼくの主⼒計算機 gsjgix による offset CT の画像再構成の処理時間（単位は秒）を調べ、この E-mail に添付した 1607_time.txt に前回の E-mail で紹介した結果と⼀緒に⽰しました。
この処理を⾏っている時に嫌なことに気づきました。最新版（ver. 17.0.0）の Intel C Compiler（icc）に最適化のオプション指定（-fast や "-O[1-3]"）を付けてコンパイルすると画像再構成プログラム oct_sg2tg_t に組み込んである cbp_thread_int.c（整数演算による逆投影の処理をマルチスレッドで実⾏する CBP engine
for CPU）のコードが誤動作します。これは cbp_int.c（整数演算による逆投影をシングルスレッドで実⾏）のコードでも同様です。最新版の icc に最適化なしのオプション（-O0）を指定した場合や、以前の版の iccを使った場合にこれらのコードは正常に動作するので、この問題の原因は icc 17.0.0 の最適化処理にあると思われます。ただし、誤動作している箇所を特定できないので、Intel にクレームをつけることができません。
そのため、今回の gsjgix による処理では gcc でコンパイルしたプログラムを使いました。
"icc -fast" と⼊⼒してコンパイルしたプログラム oct_xy と oct_tg の処理時間は 1607_oct_time.txt に記した秒数の半分以下でした。
(3)前の E-mail で紹介した offset CT の画像再構成処理の「精度」について書くのを忘れていました。
まず、oct_xy で推定したオフセットの値 Ox と Oy の妥当性についてですが、この E-mail に添付した
1607_oct_sg.pdf を御覧下さい。測定 160723[k,l,n,o] のそれぞれに対して以下の４個の sinograms（SG）を並べてあります。
上段の SG（投影画像 G の SG）
ct_xy で推定した値を指定して（未説明の）プログラム oct_sg で合成したスライス番号 500 の SG。
下段左の SG（投影画像 F の SG）
以前に紹介したプログラム sg_hg を使って offset CT の測定データから得たスライス番号 500 の
SG のうちのサンプル回転⾓が 180°未満の領域。表⽰輝度と投影値の対応関係を G のものと同じにしてある。
下段右の SG（投影画像 H の SG）
offset CT の測定データから得たスライス番号 500 - Oy の SG を左右反転した SG。表⽰輝度と投影値の対応関係は G や F のものと同じ。
下段中央の SG（差分 SG、F - H）
投影画像 F と H が重複している部分の投影値の差分値を格納した SG。表⽰輝度と差分値の関係を
上記の３個の SGs の５倍に強調してある。
測定 160723[l,n] の差分 SG にはそれら以外のものよりも⼤きな違いがありますが、合成した画像 G の SGと「元の画像 F と H の SG」を⽬視で⽐較しただけではそれらの違いを識別できません。160723[k,l,n,o] では oct_xy が（⾃動的に）検出した Ox と Oy の値に問題はないように思われます。

10 / 93
offset CT ⽤の３個の画像再構成プログラム（oct_tg と oct_sg2tg_[t,g]）のそれぞれが出⼒した CT 画像には⽬視で識別できる違いはない感じです。これらの３次元画像全体の画素値ヒストグラムを調べてみました。
この E-mail に添付した 1607_oct_hg.pdf です。⻘、緑、⾚⾊の３本のヒストグラムの線がほぼ重なっている
ので、３個の画像再構成プログラムの計算精度は概ね同じだと⾔えそうです。
とりあえず、以上です。
添付ファイル 1607_oct_xy.pdf

添付ファイル 1607_time.txt
9/26 の E-mail に添付したものと重複するので、ここには掲載しません。

11 / 93
添付ファイル 1607_oct_sg.pdf

添付ファイル 1607_oct_hg.pdf
測定 160723[n,o] の画像再構成に間違った検出器間隔の値を使っていました。
それを正したものを 10/8の E-mail に添付したので、この E-mail の添付ファイルはここには掲載しません。

12 / 93
Date:

Mon, 26 Sep 2016 11:03:33 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya , Aiko Nakato , Sugimoto Miyama , 中村隆太 , Akira Kitayama ,⼟`⼭ 明 , Hashiguchi Minako , ⽵内, 星野, 上椙Subject: oct_3

うえすぎさま、
GSJ/AIST のなかのです。前回の E-mail に「最新版の Intel C Compiler（icc 17.0.0）に最適化オプションを指定してコンパイルした oct_sg2tg_t は誤動作する」と書きました。その後、先週⼟曜⽇に配布開始のアナウンスがあった icc 16.0.4（ひとつ前の major revision の最新版の icc）で試すと、oct_sg2tg_t の誤動作の原因だった cbp_thread_int.c が正常に動くことがわかりました。そこで、icc 16.0.4 に最適化オプション "-fast"を指定してコンパイルした offset CT ⽤のプログラム oct_xy、oct_tg と oct_sg2tg_t の計算機 gsjgix の上での処理時間を再計測してみました。その結果を以前のものと⼀緒にこの E-mail に添付した 1607_time.txt に記しておきましたので、どうぞ御覧下さい。
前回の E-mail に書いた「"icc -fast" としてコンパイルした oct_xy と oct_tg の処理時間は gcc でコンパイルしたものの半分になる」は間違いでした。oct_xy の処理時間は約半分になっていますが、oct_tg の
それは icc と gcc のどちらを使っても⼤差なしです。なお、icc を使った oct_sg2tg_t の処理時間は gccを使ったものの値の 2/3 程度になっています。
とり急ぎ、
添付ファイル 1607_time.txt
gsjgixCPU

Intel Xeon E5-2687W (16 cores, 128GB memory)GPU

Nvidia Tesla K20c (5GB memory)

F/W

CentOS 6.8 + icc 16.0.4 and/or gcc 4.4.7 + CUDA 5.0

CPUIntel Xeon E5-2609 v2 (8 cores, 64GB memory)GPU

Nvidia GeForce GTX 750 Ti (2GB memory)

F/W

Scientific Linux 6.8 + gcc 4.4.7 + CUDA 7.5

vrm

test01CPUIntel Core i3 530 (4 cores, 16GB memory)GPU

Nvidia Tesla C2075 (5GB memory)

F/W

Scientific Linux 6.8 + gcc 4.4.7 + CUDA 6.0

13 / 93160723k (1292×1032×3600 / 2461×1028×1800)
gsjgix (icc)

gsjgix (gcc)

vrm

test01

oct_xy136.169601

224.782030

324.198599

453.565869

oct_tg2491.025488

2488.979275

3214.311764

9049.451679

oct_sg2tg_t2315.625867

3011.798583

3523.193628

10204.277470

582.626244

937.934966

1093.643824

oct_sg2tg_g160723l (1292×1032×3600 / 2452×1027×1800)
gsjgix (icc)

gsjgix (gcc)

vrm

test01

oct_xy111.662926

226.722810

378.926121

455.483427

oct_tg2742.630453

2477.232121

3427.715060

8958.444803

oct_sg2tg_t2225.590112

3049.520787

3617.878824

10100.010755

658.687303

1004.150321

1065.524604

oct_sg2tg_g160723n (964×1044×3600 / 1834×1041×1800)
gsjgix (icc)

gsjgix (gcc)

vrm

test01

oct_xy52.073374

178.812962

252.022873

373.306258

oct_tg1347.201550

1388.172050

2063.673255

4983.110979

oct_sg2tg_t1292.401548

1772.655100

2100.376438

5751.708006

361.263332

653.065814

646.277122

oct_sg2tg_g160723o (964×1044×3600 / 1845×1042×1800)
gsjgix (icc)

gsjgix (gcc)

vrm

test01

oct_xy51.878586

195.631270

251.760916

371.965420

oct_tg1378.125639

1410.539692

1916.112358

5043.500624

oct_sg2tg_t1255.359127

1793.865064

2160.073730

5822.809195

361.086667

693.666390

675.303937

oct_sg2tg_g

Date:Sat, 08 Oct 2016 15:05:04 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya , Aiko Nakato , Sugimoto Miyama , 中村隆太 , Akira Kitayama ,⼟`⼭ 明 , Hashiguchi Minako , ⽵内, 星野, 上椙Subject: oct_4

うえすぎさま、
GSJ/AIST のなかのです。9/21 の E-mail に記した（７⽉の offset CT 実験で得た）測定 160723 [n,o] の検出器間隔の値が間違っていました。
On Wed, 21 Sep 2016 20:12:18 +0900 Tsukasa NAKANO wrote:
>

測定

Nx

Ny

Nθ

Dr

note

>

160723[k,l]

1292

1032

3600

35.9

7 keV

14 / 93>160723[n,o]

964

1044

3600

31.1

8 keV

160723[n,o] の Dr の値は 31.1 nm ではなく 41.1 nm でした。この誤りにより再構成画像の上の estimated
LAC が真の値の 41.1 / 31.1 倍になるだけですが、念のため oct_tg、oct_sg2tg_t と oct_sg2tg_g を使った
160723[n,o] の画像再構成をやり直しました（ただし、処理時間の再計測はしませんでした）。また、9/23 のE-mail に添付した画素値ヒストグラムの図（1607_oct.pdf）も作り直しました。この E-mail にそれを添付します。とり急ぎ、添付ファイル 1607_oct_hg.pdf

15 / 93
Date:

Mon, 17 Oct 2016 14:32:15 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya , Aiko Nakato , Sugimoto Miyama , 中村隆太 , Akira Kitayama ,⼟`⼭ 明 , Hashiguchi Minako , ⽵内, 星野, 上椙Subject: oct_5

うえすぎさま、
GSJ/AIST のなかのです。9/21 の E-mail で紹介した oct_xy より以前に書いた Offset CT のペアの投影画像の位置決め⽤のプログラムなどを紹介します。特に断らない限り、それらのソースコードや Windows の実⾏
ファイルは Offset CT ⽤の書庫ファイルに⼊れてあります。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zip(1) サンプル回転⾓が０度と 360 度の投影画像の位置決めプログラム oct_cmp_0 を使えば Offset CT の測定の最初と最後に撮影した０度と 360 度の投影画像のペアの位置決めを容易に⾏うことができます。
注
oct_cmp_0 のものと等価な処理は下記の既存のプログラムで⾏えます。
hp2xp と t2t_float
SPring-8 の X 線 CT 実験の測定データから指定した投影画像を float TIFF のファイルとして取り出すプログラムとそれを普通の integer TIFF のファイルに変換するプログラム。
rmsd_2d もしくは rmsd_2d.CC
ペアの２次元画像が重複している領域に位置する画素の値の RMSD のマップを作成する
ためのプログラム。これらのソースコードなどは以下の書庫ファイルに⼊れてある。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rmsd.tar
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rmsd.zip
ぼくは oct_cmp_0 の処理結果をこれらを使って確認しました。
oct_cmp_0 の起動法は以下の通りです。
oct_cmp_0 HiPic/ {y1y2}

Ox1 Ox2

Oy1 Oy2{MSD.tif}

HiPic/
offset CT の測定データファイル（測定ログファイル output.log と HiPic 画像 dark.img および"q*.img"）が⼊っているディレクトリの名前

16 / 93y1 と y2
ここでの処理に⽤いる投影画像上の領域の y ⽅向の範囲。これらの値の指定を省略すると投影画像の全領域を処理対象とする。
Ox1 と Ox2 および Oy1 と Oy2
０度と 360 度の投影画像のペアの x および y ⽅向の位置ズレ Ox と Oy の探索範囲。投影画像の[x,y] ⽅向の画素数を N[x,y] とすると、これらは -N[x,y]＜O[x,y]1 ≦O[x,y]2 ＜N[x,y] を満たす
整数値でなければならない。
MSD.tif
投影画像のペアの重複部分の画素の投影値の MSD のマップを格納する float TIFF の画像ファイルの名前。指定省略時は画像を出⼒しない。
oct_cmp_0 は「MSD が最⼩になっている位置ズレの値（Ox と Oy；０度の画像に付随した座標系における360 度の画像の原点の画素の座標値）」と「その MSD の最⼩値」の合計３個の数値を標準出⼒に１⾏にまとめて書き出します。
７⽉の実験の測定 160723[k,l,n,o] に対してパラメータ値 Ox1= -200、Ox2=200、Oy1= -100、Oy2=100 を指定して得た oct_cmp_0 の処理結果は以下の通りです。
測定番号

Ox

Oy

MSD(Ox,Oy)

160723k-8-66.161827e-04

160723l-2-68.971565e-04

160723n1-53.031774e-04

160723o1-52.660161e-04

つまり、７⽉のどの測定でも０度と 360 度の画像に位置ズレがあるようです。その確認のため、プログラム
hp2xp、sub_float と t2t_float を使ってこれらの画像の差分画像を作ってみました。
# 測定 160723k の０度と 360 度の投影画像の抽出# ただし、０度の画像の投影番号は０、360 度のそれは 3600
hp2xp

160723k/raw

0,3600 /tmp/xp_%d.tif

0-0.142479

0.611193← 投影値の最⼩値と最⼤値

3600-0.137838

0.584854← 〃# 位置ズレなし（Ox=Oy=0）とした差分画像の作成
sub_float /tmp/xp_0.tif /tmp/xp_3600.tif-0.260605

0.362302

0 0 /tmp/sub_0.tif← 差分値の最⼩値と最⼤値# oct_cmp_0 で推定した位置ズレを与えた場合
sub_float /tmp/xp_0.tif-0.148427/tmp/xp_3600.tif -8 -6

0.165524/tmp/sub_1.tif← 差分値の分布域が狭い

17 / 93
# ２個の差分画像を差分値の分布域の幅で正規化した８ビット画像に変換# ただし、( 0.362302 - -0.260605 ) / 255 = 0.00244277254901960784
t2t_float /tmp/sub_0.tif/tmp/sub_0.tif-0.260605

t2t_float /tmp/sub_1.tif/tmp/sub_1.tif-0.260605 0.00244277254901960784

8/tmp/sub_0.tif

8/tmp/sub_1.tif

2.442773e-03-0.260605 0.00244277254901960784-0.260605

2.442773e-03

この E-mail に添付した 1607_oct_0.pdf に上記の処理で得た８ビット画素値の２個の差分画像を並べました。
また、それらの上部に oct_cmp_0 で得た０度と 360 度の重複部分の投影値の MSD のマップの画像も並べました。これらより、７⽉のどの測定でも０度と 360 度の投影画像に位置ズレがあることは明らかです。
(2) サンプル回転⾓が 180 度異なる投影画像のペアすべての位置決め
9/21 の E-mail で説明した oct_xy とほぼ同じ処理（サンプル回転⾓が 180 度異なる投影画像のペアのすべての位置決め）を⾏うプログラム oct_cmp_180 を書きました。両者の起動法も概ね同じです（ただし、
oct_cmp_180 はシングルスレッドのプログラムなのでスレッド数の指定は不要です）。
On Wed, 21 Sep 2016 20:12:18 +0900 Tsukasa NAKANO wrote:>>

oct_xy HiPic/→Ox

oct_cmp_180
→

Ox1
Oy

HiPic/ {y1
Ox

Oy

Ox2

Oy1

Oy2 {MSD.tif} > xy.log

MSD
y2}

Ox1

# 標準エラー出⼒
Ox2 Oy1 Oy2 {MSD/} > cmp_180.log

MSD# 標準エラー出⼒

HiPic/：oct_xy や oct_cmp_0 の起動パラメータと同じ意味
y1 と y2：oct_cmp_0 のものと同じ
Ox1 と Ox2 および Oy1 と Oy2：oct_xy や oct_cmp_0 のものと同じ
MSD/
oct_xy と同じ⽅法で計算した 180 度異なる投影画像のペアのそれぞれの重複部分の MSD のマップの画像を格納する「ディレクトリ」の名前。ただし、そこに書き込む画像ファイルの名前は「ペアの投影画像のうちの 180 度未満の画像の投影番号 +".tif"」。
cmp_180.log
oct_cmp_180 が標準出⼒に書き出す下記の「投影画像のペアに対する５個の値（最初の４個は
oct_xy の出⼒と同じ値で、５個⽬の値だけが異なる）」の⾏を書き込むテキストファイルの名前。
[1] 回転⾓θが 180 度未満の画像の投影番号（０〜）[2,3] MSD(θ,Ox,Oy) が最⼩になっている O[x,y](θ) の値[4] MSD(θ,Ox,Oy) の最⼩値
[5] MSD(θ,Ox,Oy) が最⼩になっている状態の画像ペアの重複部分にある 2 画素の「投影値の積（product）の平均値」

18 / 93
Ox、Oy と MSD
標準エラー出⼒に書き出される値は oct_xy のものと同じ。
Ox
サンプル回転⾓が 180 度未満の投影画像と「それと 180 度異なる回転⾓の投影画像を左右反転した画像」の重複部分の投影値の MSD が最⼩になっている x ⽅向のオフセット。
Oy：上と同じ２画像の MSD が最⼩の配置の y ⽅向のオフセット
MSD：上記の MSD の最⼩値
起動パラメータ y1 と y2 による領域指定をせずに oct_cmp_180 で７⽉の測定のデータを処理すると、標準
エラー出⼒に書き出されたオフセット Ox と Oy（および MSD）の値は oct_xy が出⼒した値とまったく同じになりました。
>

測定

Ox

Oy

>

160723k

1169

-4

>

160723l

1160

-5

>

160723n

870

-3

>

160723o881-2oct_xy とは異なり oct_cmp_180 はサンプル回転⾓が 180 度異なる「投影画像のペアのそれぞれの MSD のマップ」の画像を作成します。プログラム t2g_float を使えば、これらの画像を投影番号の順にまとめた GIFの動画を作成できます。
# それぞれのマップを最⼩値と最⼤値で正規化した表⽰輝度の動画の作成法
t2g_float

160723k_msd/- 160723k_msd.gif# カラーマップを指定することにより擬似カラー画像にすることもできる。
t2g_float

160723k_msd/- CM.txt 160723k_msd.gif

カラーマップのデータファイル CM.txt
各⾏に GIF 画像の上の画素値と表⽰⾊の R、G、B 成分の値の合計４個の値（いずれも０〜255 の整数値）を書き並べたテキストファイル
また、逆に、プログラム add_float を使えば、oct_cmp_180 が出⼒した複数個のマップの画像から oct_xy のものと等価な「回転⾓について平均した MSD のマップ」の画像を作成することもできます。
ls 160723k_msd/* |add_float

0 160723k_msd.tif# ただし、add_float のパラメータの値「０」は今の場合は「ダミー」です。
oct_cmp_180 や oct_xy が標準出⼒に書き出した（ファイル cmp_180.log や xy.log に書き込んだ）データか
ら画像ペアのオフセットの値 Ox(θ) と Oy(θ) を近似した「サンプルの回転⾓θに対応した投影番号 v に関する１次式」の係数値を計算するプログラム oct_fit_180 を書きました。

19 / 93
oct_fit_180 org.log rangeList {new.log} > fit_180.log
org.log
oct_xy が出⼒した xy.log もしくは oct_cmp_180 の cmp_180.log
rangeList
最⼩⾃乗法の処理対象とするオフセット値の投影番号 v の範囲指定⽤の range list。すべてを処理対象とする場合は "-" を指定すれば良い。
new.log
rangeList で指定した回転⾓の範囲に該当する org.log のデータ⾏だけを新しいファイル new.log にコピーする。この指定を省略すると new.log を作成しない。
fit_180.log
投影番号 v のオフセット値 Ox(v) と Oy(v) を近似する１次式
Ox(v) = Ax * v + Bx ± σx
Oy(v) = Ay * v + By ± σyの合計６個の係数値を２⾏に分けて書き込むテキストファイルの名前。ただし、σ[x,y] は O[x,y]の近似の推定誤差。
通常の X 線 CT で使っている detector offset（DO）の推定プログラム hp2DO と同様に、oct_fit_180 は org.logに書き込まれている投影値の積の値（== 投影画像の「影の濃さ」）を「重み」として最⼩⾃乗法を⾏います。
ただし、この積の値が書き込まれていない（oct_xy のログファイル xy.log を処理する）場合には重みの値を１とするようにしてあります。
さらに、oct_cmp_180 が出⼒した cmp_180.log や oct_fit_180 の new.log のデータを GNU plot を使⽤してグラフ表⽰する C-shell script "plot.oct_180" を書きました（GNU plot に加えて UNIX の標準コマンド群を使っているので、この script は Linux or MacOSX ⽤です）。
plot.oct_180 cmp_180.log {Ox Oy} {fit_180.log} > oct_180.ps
cmp_180.log
oct_cmp_180 が出⼒した cmp_180.log や oct_fit_180 の new.log
Ox と Oy
Ox(θ)と Oy(θ)のグラフの上に⾚線で⽰すオフセット値の「基準値」。oct_cmp_180 が標準エラー
出⼒に書き出したオフセット値を指定すれば良い。これらの指定を省略すると⾚線を描かない。
fit_180.log
Ox(θ)と Oy(θ)のグラフの上に⻘線で⽰す１次式の係数の値を書き込んだファイルの名前。
oct_fit_180 の標準出⼒をリダイレクトしたファイルを指定すれば良い。指定を省略すると⻘線を描かない。
oct_180.ps
グラフを記述する PostScript のコードを書き込むファイルの名前。
plot.oct_180 で描いた７⽉の Offset CT の４個の測定それぞれに対するグラフをこの E-mail に添付した

20 / 93
1607_oct_180.pdf にまとめて⽰しました。ただし、各測定の４個のグラフの横軸はいずれもサンプル回転⾓
θに対応する投影番号 v（= 0〜1799）です。また、測定 160723l では v = 1004〜1033（サンプル回転⾓に換算するとθ= 100.4〜103.3 度）のオフセットの値が「異常」だったので、その部分のデータを取り除いたグラフ（最⼩⾃乗法の計算でもその部分のデータを使わなかった）を 1607_oct_180.pdf の２ページ⽬に載せました。これは具体的には、以下の⼊⼒によって oct_fit_180 で作成したファイル new.log と fit.log を
plot.oct_180 の起動パラメータとして指定して描きました。
oct_fit_180 140723l_cmp_180.txt -1003,1034- new.log> fit.log(3) 測定中のサンプルの動き
Offset CT ⽤のプログラムを書き始めた当初は 1607_oct_180.pdf のような投影番号 v に対するオフセット値
Ox(v) と Oy(v) のグラフを作れば測定中のサンプルの動きがわかると思っていたのですが、それは以下の理由で難しいです。
[1] oct_cmp_180 や oct_xy で得た Ox(v) や Oy(v) の推定精度が低い。
[2] 測定中のサンプルの動きを表現する適当なモデルがない。
これらについて考えがまとまっていないので、今⽇のところはここまでとします。
10/11 にお送りした７⽉の FZP-CT 実験の報告書で紹介した「browse 画像」のうちの Offset CT の分だけを抜き出し、この E-mail に添付した 1607_oct_xyz.pdf に並べました。実測した２×２＝４測定の分と、それらの画像をつなぎ合わせた（仮想的な）測定 160723kl と 160723on の合計６セットの browse
画像です。1607_oct_180.pdf に⽰した「mean value of product（画像が重複している部分の投影値の積）」のグラフの変化の様相はサンプルの形状によるものであることがわかります。
添付ファイル 1607_oct_0.pdf

21 / 93
添付ファイル 1607_oct_180.pdf（１ページ⽬）

22 / 93
添付ファイル 1607_oct_180.pdf（２ページ⽬）

23 / 93
添付ファイル 1607_oct_xyz.pdf

24 / 93
Date:

Mon, 17 Oct 2016 19:11:00 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya , Aiko Nakato , Sugimoto Miyama , 中村隆太 , Akira Kitayama ,⼟`⼭ 明 , Hashiguchi Minako , ⽵内, 星野, 上椙Subject: oct_6

うえすぎさま、
GSJ/AIST のなかのです。たびたびですみません。Offset CT ⽤の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zipに⼊れたプログラムのうち、これまでに説明していないものの概要を書きます。
(1) Offset CT ⽤のユーティリティプログラム
oct_xt HiPic/Ox Oy {rangeList XT_format}

oct_xp HiPic/Ox Oy {rangeList XP_format}

ディレクトリ HiPic/ の下の Offset CT の測定データから rangeList で指定した投影番号の X 線透過率（XT）や投影（XP）画像を抽出するプログラム。rangeList と出⼒画像のファイル名のフォー
マット指定を省略すると oct_tg と同様に出⼒画像の画素数や枚数だけを出⼒する。
oct_sg HiPic/ Ox Oy {rangeList

SG_format}

rangeList で指定した番号のスライスの sinogram のバイナリデータを取り出すためのプログラム。
書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.zipに⼊っているプログラム sg2tg などを使えば取り出したバイナリファイルから画像再構成を⾏うことができる。
また、それらの書庫ファイルに⼊っているプログラム a2tiff を使えば sinogram を TIFF
画像（通常の integer TIFF 画像）に変換できる。
(2) float-TIFF の画像⽤のユーティリティプログラム
xp2xt_float XP.tif
xp2xt_float XP/

XT.tifnameFile XT/

xt2xp_float XT.tif XP.tif
xt2xp_float XT/

nameFile XP/

XP の画像を XT の画像に、もしくは、その逆の変換を⾏うプログラム。ディレクトリ名と nameFileを指定すると複数ファイルの変換を⼀度に⾏うことができる。

25 / 93(3) その他のユーティリティプログラム
hp2xt、hp2xp と trim_float および
t2t_float と t2g_float（これらについては少しだけ説明しました）
SPring-8 の X 線 CT 測定のデータファイル群の読み込みのための汎⽤サブルーチン "rhp.[h,c]" の利⽤例として書いたプログラム。これらの使⽤法については説明書
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rhp.pdf
を御覧下さい。
(4) Offset CT の画像再構成の⼿順を記したファイル
これはプログラムではありませんが、上杉君が SPring-8 の X 線 CT のページ
http://www-bl20.spring8.or.jp/xct/にリンクしていた「⽊材？」を撮影した Offset CT の測定データに対する画像再構成の実際の⼿順をテキス
トファイル wood.txt に記しました。この E-mail に添付したファイルと同じものを前記の Offset CT ⽤の書
庫ファイルにも⼊れておきました。ただし、これを Windows で試す場合は以下のようにして下さい。
プログラムとサンプルデータのインストール
Offset CT ⽤の書庫ファイルやサンプルデータのダウンロードと解凍・展開は Windows ⽤の適当なインターネットもしくはファイルブラウザで⾏って下さい。なお、Offset CT ⽤の書庫ファイルには Windows ⽤の実⾏ファイルを⼊れてあるので、プログラムのコンパイルは不要です。
プログラムの実⾏時の⼊⼒[1] 使⽤するスレッド数の設定法については wood.txt に記した通りです。
[2] プ ロ グ ラ ム が 標 準 出 ⼒ に 書 き 出 し た デ ー タ を 捨 て る た め に 使 っ て い る "/dev/null" はWindows では "nul"（"l" は１個です）として下さい。
[3] UNIX の標準コマンド “cut” や “bc” などを⽤いて estimated LAC の値（CT 値）をいじっている⾏は実⾏できなくても問題ないので無視して下さい。
とりあえず以上です。
添付ファイル wood.txt# installation of programs
wget http://www-bl20.spring8.or.jp/\~sp8ct/tmp/oct.taz...
tar xzf oct.taz
cd oct
make
...
cd ..

26 / 93
set path=(${cwd}/oct $path)# execution path# installation of sample data
wget http://www-bl20.spring8.or.jp/xct/soft/lzh/offsetCT.zip...
unzip offsetCT.zip...
cat offsetCT/readme.txt...
sample: wood
pixel size: 5.4um# == 5.4e-4 cm

energy: 15keV
beamline: BL20B2
number of projection: 1800 for 360degree
image size: 1000 x 32 pixels
rotation axis: 909...
# setting of number of threads for oct_xy and oct_tg
setenv THREADS

8# tcshexport THREADS=8# bashset THREADS=8# Windows# estimation of offset-values
oct_xy offsetCT 729 909 -16 16819

0> /dev/null

3.105794e-03# Ox, Oy and MSD# reconstruction of tomograms with float- and byte-values
mkdir

tg

oct_tg offsetCT 819 0
1819
cut

32

5.4e-4 0 tg/%02d.tif > tg.log
1800# N[x,y], Nz and views-f3 tg.log | sort -n | tail -1
13.905724# maximum value of estimated LAC in 1/cm

echo 13.905724 / 255| bc -l.05453225098039215686
mkdir# upper limit of estimated LAC / 8 bit PV

byte

t2t_float tg -0 0.025 8

byte > /dev/null27 / 93
Date:Thu, 03 Nov 2016 14:44:21 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya , Aiko Nakato , Sugimoto Miyama , 中村隆太 , Akira Kitayama ,⼟`⼭ 明 , Hashiguchi Minako , ⽵内, 星野, 上椙Subject: oct_7

うえすぎさま、
GSJ/AIST のなかのです。Offset CT ⽤の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zipに「⽊材？」を撮影した Offset CT の測定データに対して Windows 上で画像再構成処理を⾏う以下の３個のバッチファイルを⼊れておきました。
wood_xy.bat
プログラム oct_xy を使ってオフセット値を推定する。
wood_tg.bat および wood_sg2tg_t.bat
oct_tg および oct_sg2tg_t を使って画像再構成を⾏う。
本来は oct_xy で推定したオフセットの値を oct_tg や oct_sg2tg_t の起動時に指定して画像再構成を⾏いますが、それを埋め込んであるので、これらのバッチファイルの起動の順番は任意です。上記の順番でバッチ
ファイルを起動した時のログファイル wood.log をこの E-mail に添付しました。
バッチファイル “wood_*.bat” は実⾏中のコマンド⾏を標準出⼒にエコーバックする設定にしてあります。ただし、ぼくがバッチファイルに記したコマンドの標準出⼒のファイルへのリダイレクションの書式が以下のように改変されていました（どちらでも⽂法的には問題なしです）。
ぼくが指定した書式：command > file
エコーバックの書式：command 1> file
また、逆に、再構成したスライス画像のファイル名のフォーマットはバッチファイル⽤の書式から端末
⼊⼒⽤のものに変更されたエコーバックになっています（バッチファイルの上では⽂字 “%” は “%%”に「エスケープ」しなければならない）。
バッチファイル⽤のファイル名のフォーマットの書式：directory/%%02d.tif
端末⼊⼒を⾏う際のファイル名のフォーマットの書式：directory/%02d.tif
なお、wood.log に記されているように、以前に紹介したプログラム stop_watch を使ってそれぞれのバッチ
ファイルの実⾏時間を測定しました。

28 / 93
使⽤した Windows 機（東芝 dynabook R73/W6M）
CPUIntel Core i7-4710MQ（2.50 GHz）
4 cores → hyper threading on なので⾒かけ上は 8 cores
memory and SSD：16 GB and 1 TB
OS：Windows 8.1（64 bit 版）
実⾏時間
wood_xy.bat：6.4337449 秒
wood_tg.bat：126.0409685 秒
wood_sg2tg_t.bat：83.2228138 秒
とりあえず以上です。
添付ファイル wood.log
cd Desktop\OffsetCT\woodpath=..\oct;%path%
stop_watch wood_xy.bat 2>&1
set THREADS=8
oct_xy offsetCT 729 909 -16 16819

06.4337449

xy.tif1> xy.log

3.105794e-003← stop_watch の出⼒

stop_watch wood_tg.bat 2>&1
set THREADS=8
mkdir

tg

oct_tg_32 offsetCT 819 0
1819

32

5.4e-4 0 tg/%02d.tif 1> tg.log

1800

t2t_float tg -0 0.025 8

126.0409685← stop_watch の出⼒

stop_watch wood_sg2tg_t.bat

tg 1>nul

2>&1

set CBP_THREADS=8
mkdir

sg2tg_toct_sg2tg_t_32

offsetCT

1819

32

83.2228138

819

0

-

5.4e-41800← stop_watch の出⼒

0

0

0.025

8

sg2tg_t/%02d.tif

1>sg2tg_t.log

29 / 93
Date:

Wed, 23 Nov 2016 14:37:04 +0900

From:

Tsukasa NAKANO

To:

Akira TSUCHIYAMA

Cc:

Kentaro UESUGI, 三宅 亮, matsuno j, Aiko Nakato, Sugimoto Miyama, kitayama.akira.46a,
中村隆太, tkadokawa, Takashi Sakurama, Minako HASHIGUCHI, take, hoshino, uesugi

Subject: oct_lm+sa+sbs

みなさま、
GSJ/AIST のなかのです。投影画像の「つなぎ合わせ⽅」に関する Offset CT ⽤の以下の６個のプログラムの改修作業が終わりました。
oct_x[t,p]サンプル回転⾓が 180 度異なる X 線透過率（XT）もしくは投影値（XP）画像のそれぞれをつなぎ合わせる。
oct_sg
投影値 sinogram（SG）をつなぎ合わせる。
oct_tg投影画像をつなぎ合わせた後に以前に紹介したプログラム cbp や fdk で使っている（投影画像ごとに convolution と back-projection 処理を⾏う）⼿法で画像再構成を⾏う。
oct_sg2tg_[t,g]SG をつなぎ合わせた後に oct_sg2tg_t は CPU ⽤の CBP engine で、また、oct_sg2tg_g は CBPengine for CUDA を使って画像再構成を⾏う。
Offset CT ではサンプルを 360 度回転（full rotation；FR）した投影画像からサンプルを 180 度回転（halfrotation；HR）した通常の CT で撮影したものと等価な投影画像 C を合成した後に画像再構成を⾏います。
先⽇の実験の場でも話したように、その際に⾏うサンプルの回転⾓が 180 度異なる FR の投影画像のペア Lと R の重複している部分（Ox ≦ x < Nx）の投影値の処理の仕⽅には以下の３通りが考えられます（記号の意味の説明は省略しますが、これでわかりますよね？）。
simple average（SA）
for Ox ≦ x < Nx

C[x] = (L[x] + R[x]) / 2

side by side（SBS）
for Ox ≦ x < xC

C[x] = L[x]for x = xC

C[x] = (L[x] + R[x]) / 2

for xC < x < Nx

C[x] = R[x]where xC = (Ox + Nx - 1) / 2← サンプル回転軸の位置

linear mixing（LM）
for Ox ≦ x < Nx

C[x] = { (Nx - x) * L[x] + (x - Ox + 1) * R[x] } / (Nx - Ox + 1)

30 / 93
注
先⽇の実験の場で計算機 MUSASHIMARU にインストールした linear mixinig のテスト⽤の画像再構成プログラム oct_sg2tg_g_LM（だったと思う）では上記の式とは若⼲異なる式を使っています：
for Ox ≦ x < Nx

C[x] = { (Nx - 1 - x) * L[x] + (x - Ox) * R[x] } / (Nx ‒ Ox - 1)

どちらの式を使っても⼤差はないと思いますが、...。
以前の Offset CT ⽤のプログラムでは simple average で Offset CT の画像をつなぎ合わせていましたが、
linear mixing を使うように改修しました。また、それ以外の⼿法も拡張⼦ "_sa" と "_sbs" のプログラムで実⾏できます。
⼿法 LM（linear mixing）で Offset CT の投影画像をつなぎ合わせる
oct_x[t,p]、oct_sg、oct_tg、oct_sg2tg_[t,g]SA（simple average）
oct_x[t,p]_sa、oct_sg_sa、oct_tg_sa、oct_sg2tg_sa_[t,g]SBS（side by side）
oct_x[t,p]_sbs、oct_sg_sbs、oct_tg_sbs、oct_sg2tg_sbs_[t,g]
改修したプログラムのソースコードと Windows ⽤実⾏ファイルを以前と同じ Offset CT ⽤の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zipに⼊れておきました。以前と同様にしてコンパイル・インストールして下さい。
先⽇の SPring-8 実験の測定 161119[q,r,s] のデータを⽤いて LM、SA と SBS による処理結果の違いを調べてみました。この E-mail に添付した 19qrs.pdf をまず御覧下さい。３ページそれぞれの上段の横列は同じ図で、つなぎ合わせる前の投影値 sinograms とそれらの重複部分の差分 sinogram です（ただし、差分値の表
⽰輝度は投影値の５倍に強調しています）。そして、中段の横列は合成した sinogram、また、下段の横列はそれらから再構成した画像です（再構成画像の CT 値と表⽰輝度の対応関係はすべて同じにしてあります）。
これらより SA を使った再構成画像のリング状の偽像が顕著ですが、LM と SBS の画像の良否は良くわかりません。そこで、再構成したスライス画像それぞれの中央部の同じ画素数の領域を切り出し、この E-mail に添付した 19qrs.jpg に並べました（サイズが⼤きいので 100% quality ではない JPEG にしましたが、画質の劣化はほとんどない感じです）。これらを眺めても SBS と LM の画像の差異はわかりません。と⾔うことで、
Offset CT ⽤のプログラムのデフォルトを LM にしましたが、SBS でも良かったようです。
とりあえず以上です。

31 / 93
添付ファイル 19qrs.pdf（１ページ⽬）

添付ファイル 19qrs.pdf（２ページ⽬）

32 / 93
添付ファイル 19qrs.pdf（３ページ⽬）

添付ファイル 19qrs.jpg

33 / 93
Date:

Wed, 30 Nov 2016 10:13:10 +0900

From:

Tsukasa NAKANO

To:

三宅 亮, Kentaro UESUGI, Akira TSUCHIYAMA

Cc:

matsuno j, Aiko Nakato, Sugimoto Miyama, kitayama.akira.46a, 中村隆太, tkadokawa,Takashi Sakurama, Minako HASHIGUCHI, take, hoshino, uesugi

Subject: oct_frct_hrct

みなさま、
GSJ/AIST のなかのです。
(1)
昨⽇の三宅さんからの E-mail で「先週の SPring-8 BL47XU での FZP-CT 実験の測定 161122n は⽋測」が確定したので、結局、Offset CT（OCT）のテスト⽤の測定は下記の３つのグループのもの（合計 3 + 4 + 5 =12 測定）だけですね。
グループ名

測定サンプル

測定番号

19qrs

MIL090657

161119[q,r,s]

21dfhi

MIL090657

161121[d,f,h,i]

22pqrstforsterite

161122[p,q,r,s,t]

11/23 の E-mail に測定グループ 19qrs の再構成画像の画質などについて書いたので、この E-mail では残りの測定グループに関する話をします。
(2)
その処理中に OCT ⽤の画像再構成プログラム "oct_sg2tg_*" のバグに気づきました。致命的なバグではありません（"oct_sg2tg_*" の起動パラメータの base と step の値を指定して integer TIFF の画像ファイルに格納した再構成画像を既存のプログラム tg2tg を使って処理する場合にのみ問題になるバグなので、ぼく以外のユーザに影響なし？）。修正したソースコード oct_sg2tg.c とそれをコンパイルした 32 および 64 bit
Windows ⽤の実⾏ファイルを以下の OCT ⽤の書庫ファイルに⼊れておきました（注：これらを展開した
oct/oct_sg2tg.c と oct/cuda/oct_sg2tg.c はハードリンクされている同⼀のファイルです）。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zip(3)測定グループ 21dfhi に対しても 19qrs の場合と同様に「再構成したスライスの中央部の概ね同じ広がりの領域を切り出した画像」を並べてみました。この E-mail に添付した 21dfhi.pdf を御覧下さい。これらは 21dfhiの４測定それぞれの８ビット画素値の３次元 CT 画像の COA（吸収重⼼）を通るスライスで、サンプル回転
軸が同じ位置になる領域を切り出しました。そして、その右下に記したように 21dfhi.pdf を拡⼤・縮⼩せず

34 / 93にそのままプリンタ出⼒すると、実⼨の 2000 倍になります（出⼒した紙⾯上の各画像の横幅は約 8 cm なので、その実⼨は 40μm 程度）。これらの４測定の画像の撮影条件は以下の通りです。
161121d

7 keV かつ 1800 投影の half rotation CT で撮影

161121f

8 keV、1800 投影、half rotation CT

161121h

8 keV、900（×２）投影、full rotation CT

161121i

8 keV、1800（×２）投影、full rotation CT

８keV で撮影した画像のリングアーチファクトは上記の順で減っているように⾒えます。つまり、FZP-CT の画像の画質向上に full rotation CT が有効だと⾔えそうです。しかしながら、投影数を２倍にしても画質は概ね同じ感じです。
(4)
測定グループ 22pqrst の画像群はサンプルの下部に「台座」があったため COA で「位置決め」できませんでした。それで、このグループの画像に対しては画素値ヒストグラムを調べるだけにしました。この E-mail に添付した 22pqrst.pdf を御覧下さい。７keV と８keV の画像のヒストグラムの図を別々に描きました。それ
ぞれの再構成画像に関するラベルや、forsterite（密度 3.226 g/cc、組成 Mg2SiO4 と仮定）の LAC の理論値も⽰しました。これらより forsterite だと思われるヒストグラムのピークはどれも概ね同じ位置にあることがわかります。ただし、測定 16112p（８keV の full rotation CT）では他のものよりやや低い値の位置にピー
クがあります。また、これらのピークの CT 値と forsterite の LAC の理論値の⽐は７keV と８keV で違っています。これは何故だ？
(5)
こちらで処理した Offset CT のテスト⽤の測定に関するリスト 1611_test.txt をこの E-mail に添付します。
MS-Excel などで開いてみて下さい。その１⾏⽬に記したように、テスト⽤の測定に関するファイルはすべて
SPring-8 の計算機 vrm のバードディスクに置いてあります。
vrm.spring8.or.jp:/media/disk/tsukasa/1611/test/
raw/測定番号
測定データのディレクトリへのシンボリックリンク
測定グループ名.log
各グループの画像再構成のログファイル（1611_test.txt と同内容）
測定グループ名/測定番号_ext/
８ビット画素値の３次元再構成画像が⼊っているディレクトリ
測定グループ名/測定番号_ext.log
再構成画像の各スライスに出現した CT 値の最⼩値と最⼤値
測定グループ名/測定番号_ext.csv
８ビット画素値の再構成画像上の CT 値のヒストグラム
なお、上記の ext を含む 1611_test.txt で⽤いた記号の意味は以下の通りです。

35 / 93ext（extension の略のつもり）
full rotation（FR）CT や half rotation（HR）CT、および、「linear mixing（LM）もしくは side byside（SBS）でつなぎ合わせた投影画像を使⽤した OCT」のそれぞれを区別するための記号。
set
測定番号
nm
測定した投影画像の正⽅形画素の nm 単位の辺⻑。再構成画像の⽴⽅体画素の辺⻑もこの値になる。
Ox と Oy もしくは DO_base と DO_step
FRCT と OCT の画像再構成に使う投影画像の x と y ⽅向のオフセットの値、もしくは、HRCT の画像再構成で使うサンプル回転軸の位置（DO）のスライス番号に対する切⽚（base）と傾き（step）の値。
Nr、Nz と views
再構成に使う投影画像の横・縦画素数とそれらの枚数。ただし、FRCT と OCT ではつなぎ合わせの後の投影画像の値。再構成画像の画素数は Nr^2（スライス⾯内）×Nz（スライス数）。
CTV/PV
８ビット画素値の再構成画像（byte 画像）の画素値１階調ごとの CT 値（estimated LAC）の増分。
実験の場で決めた測定ごとにまちまちな byte 画像の CTV/PV の値を平均した各組ごとに単⼀の値を使⽤した。
添付ファイル 21dfhi.pdf

36 / 93
添付ファイル 22pqrst.pdf

添付ファイル 1611_text.txt
vrm.spring8.or.jp:/media/disk/tsukasa/1611/test/
19qrs
ext.

set

lm

nm

Ox

Oy

Nr

Nz

views

CTV/PV

161119q 39.18

644

-1

2115

690

1800

2.90233333333333333333

sa161119q 39.18

644

-1

2115

690

1800

2.90233333333333333333

sbs

161119q 39.18

644

-1

2115

690

1800

2.90233333333333333333

lm

161119r 39.18

899

-1

2370

690

1800

2.90233333333333333333

sa161119r 39.18

899

-1

2370

690

1800

2.90233333333333333333

sbs

161119r 39.18

899

-1

2370

690

1800

2.90233333333333333333

lm

161119s 39.18

1153

-1

2624

690

1800

2.90233333333333333333

sa161119s 39.18

1153

-1

2624

690

1800

2.90233333333333333333

sbs

161119s 39.18

1153

-1

2624

690

1800

2.90233333333333333333

set

Ox or

Oy or

Nr

Nz

views

CTV/PV

DO_baseDO_step

21dfhi
ext.

nm

hr

161121d 70.24-684.937360 -0.000051

1348

731

1800

1.83425hr

161121f 80.70

-497

0

916

677

1800

1.83425fr

161121h 80.70

82

0

998

677

900

1.83425fr

161121i

83

-1

999

676

1800

1.8342580.70

37 / 93
22pqrst
ext.

set

nm

Ox or

Oy or

DO_baseDO_step

Nr

Nz

views

CTV/PVfr

161122p 80.70

46

-2

1480

809

900

1.1414fr

161122q 70.24

99

-1

1533

810

900

1.1414lm

161122r 70.24

-180

-2

639

584

1800

1.1414sbs

161122r 70.24

-180

-2

639

584

1800

1.1414lm

161122s 80.70

159

-2

618

584

1800

1.1414sbs

161122s 80.70

159

-2

618

584

1800

1.1414hr

161122t 80.70-417.648958 0.001928

743

565

1800

1.1414Date:

Wed, 30 Nov 2016 15:52:34 +0900

From:

Tsukasa NAKANOTo:

MiyakeCc:

Kentaro UESUGI, Akira TSUCHIYAMA, matsuno j, Aiko Nakato, Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, tkadokawa, "Takashi Sakurama", Minako HASHIGUCHI,
take, hoshino, uesugi

Subject: do-not-use-oct_sbs

みやけさま、
みなさま、
GSJ/AIST のなかのです。三宅さんの質問に答えるべく測定グループ 22pqrst の画像を調べていてすごいことに気づきました。結論から⾔うと
Offset CT（OCT）の画像再構成において SBS（side by side）で投影画像をつなぎ合わせるとその「つなぎ⽬」に相当するサンプル回転軸の位置（== 再構成画像のスライスの中央）に無視できない偽像が発⽣することがあるです。つまり、OCT の画像再構成では SBS ではなく LM（linear mixing）で投影画像をつなぎ合わせないとダメなようです。以下では今回⾏った処理の内容を紹介します。
今回の処理の結果をまとめたファイル 22pqrst.txt をこの E-mail に添付しました。それを⾒ながら以下をお読み下さい。
実は、ここでの話は再構成画像を⽬で⾒て観察すればすぐにわかることなのですが、ぼくが理想としている「画像を⽬で⾒ずに解析する」を実⾏してしまったので少し⻑い話になります（すみません）。
(1)
測定グループ 19qrs や 21dfhi の場合と同様に 22pqrst の画像についてもおおよその位置決めを⾏おうと思いました。22pqrst で撮影したサンプルの下部には測定ごとに縦幅が異なる台座があるので、とりあえずはサンプルの上端の位置を決めようと（19qrs や 21dfh のスライス画像のトリミングの場合と同様にして）画像を２値化し、それによって識別した「物体像の画素」のすべてを内部に含む最⼩の直⽅体（bounding box）
をまず調べました。

38 / 93

FRCT の画像 161122p_fr の画素値 50 以上の画素の bounding box の調査
si_area 22pqrst/161122p_fr -

50-

→

511

296

34

1032

1339

808

#

x1

y1

z1

x2

y2

z2

#

x1 〜 z2 は bounding box の対⾓点の座標値

ただし、今回の処理で採⽤した「画素値 50 以上」に深い意味はありません。19qrs や 21dfhi のスライスのトリミングではこの criterion で OK でした。
(2)
22pqrst.txt の「area] がその処理結果です。投影画像を⾒ると、サンプルの上部は空気なのに値 z1（boundingbox の上端の z 座標値）が０の画像があります。ただし、通常の half rotation CT の画像 161122t_hr でも値
z1 は０なので、これは今回の CT 画像が全体的に noizy だということかもしれません。
(3)
そこで、２値化に加えて erosion-dilation（ED）処理を⾏ってみました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/osp.pdf#page=2
これにより、指定した画素幅単位の層厚（layers）以下の幅の「ゴミ」や突起を２値画像上の物体像から除去できるはずです。プログラム si_ed を使えばその結果の物体像の bounding box を調べることができます。
layers = 1 の ED 処理（「ed_1」）
si_ed

22pqrst/161122p_fr

→

- 50-

1

511

296

34

1032

1339

808

642

463

241

973

1075

807

641

462

240

974

1076

808

#

x1

y1

z1

x2

y2

z2

#

１⾏⽬は２値化後の画像上の物体像の bounding box#

２⾏⽬は erosion の後の物体像の bounding box#

３⾏⽬がここで知りたい ED の後の物体像の bounding box(4)
22pqrst.txt の「ed_1」からわかるように、161122s_sbs を除く画像では layers = 1 の ED 処理によって値 z1
が０でなくなっています。しかしながら、layers = 2 の ED 処理（「ed_2」）を⾏っても 161122s_sbs の z1 は０のままです。これは何故かと（この時点で初めて）再構成画像を眺めると、その上端のスライス画像 000.tifの中央部（== サンプル回転軸の位置）に数画素幅以上の明るい部分があることに気づきました。これは明らかに偽像です。
(5)
そこで、161122s_sbs の画像再構成に⽤いた SBS でつなぎ合わせた投影画像を作ってみました。以下の⼊⼒によって（例えば）回転⾓０の投影画像を作成できます。

39 / 93

oct_xp_sbs raw/161122s 159 -2 0
→

618

584

0

-0.075998

xp.tif1800
2.980101

この画像を viwer で強調表⽰して眺めると、つなぎ合わせた左右の投影画像の空気の投影値にかなり⼤きな違いがあることがわかります。つまりは、OCT ⽤に撮影したサンプル回転⾓が 180°異なっているペアの画
像のノイズレベルは⼤きく異なっており、それらの画像のつなぎ⽬（の段差）が OCT_SBS の再構成画像の上の偽像の原因になっているようです（12/7 の E-mail に添付した 22s.jpg の右側最下段の OCT_SBS の再構成画像を御覧下さい；その中央に投影画像のつなぎ⽬の段差に起因する偽像が発⽣しています）。
(6)という訳で、OCT において SBS で投影画像をつなぎ合わせるのは NG です。
(7)
最後に、三宅さんの質問に答えます。
> 1122p (8kev, 360degree/1800)だけが⼤きくずれているようにみえるのは、> バックグラウンドが⾼くなっている影響はあるのでしょうか？
そうかもしれません。ただ、forsterite のピークの「幅」から考えると、同じ full rotation（FR）CT を⾏った
測定 161122[p,q] では７keV の 161122q の⽅が noizy なはずですが、...。
> そもそもなぜ⾼くなっているのかも気になりますけど。
FRCT の測定 161122[p,q] の視野のサイズ（== 投影画像の画素数）は他のものよりも⼤きいので CT 値が０付近の部分のヒストグラムが⾼い（== 画素数が多い）のは当然ですが、forsterite のピーク付近までそれが⾼い理由は不明です。
> あと、1122s で、内側の重なっている部分と、外側の部分では違うのでしょうか？
これは 161122s_lm と 161122s_sbs の再構成画像のスライスの中央部（== サンプル回転⾓が 180°異なっている投影画像の重複部分）とそれ以外の部分の CT 値の違いに関する話ですかね？ それについてはまだ調べていません。もしくは、forsterite のピークを挟んだ左右の部分のヒストグラムの⾼さの違いについての話ですか？ ぼく的には両者の違いはわずかなので気になりませんが、...。
とりあえず、以上です。
On Wed, 30 Nov 2016 10:39:57 +0900 Miyake wrote:> 中野様> ありがとうございます。

40 / 93>
> （３）> >８ keV で撮影した画像のリングアーチファクトは上記の順で減っているように> >⾒えます。つまり、FZP-CT の画像の画質向上に full rotation CT が有効だと> >⾔えそうです。しかしながら、投影数を２倍にしても画質は概ね同じ感じです。
> 確かに。
> 贔屓⽬に⾒ると、1121i がリングアーチファクトによると思われる、むらが少ないように思いますね。
>> （４）> 1122p (8kev, 360degree/1800)だけが⼤きくずれているようにみえるのは、> バックグラウンドが⾼くなっている影響はあるのでしょうか？
> そもそもなぜ⾼くなっているのかも気になりますけど。
>> あと、1122s で、内側の重なっている部分と、外側の部分では違うのでしょうか？
添付ファイル 22pqrst.txt
area
group

set_ext.

nm

x1, y1, z1

x2, y2, z2

x2-x1+1, y2-y1+1, z2-z1+1

22pqrst 161122p_fr

80.70

511, 296, 34

1032,1339, 808

522, 1044, 775

22pqrst 161122q_fr

70.24

0, 1, 01532, 1531, 809

1533, 1531, 810

22pqrst 161122r_lm

70.24

1, 4, 3637, 637, 583

637, 634, 581

22pqrst 161122r_sbs

70.24

1, 1, 0638, 638, 583

638, 638, 584

22pqrst 161122s_lm

80.70

186, 73, 107

599, 478, 583

414, 406, 477

22pqrst 161122s_sbs

80.70

159, 3, 0613, 478, 583

455, 476, 584

22pqrst 161122t_hr

80.70

251, 57, 0661, 593, 564

411, 537, 565

group

nm

x1, y1, z1

x2, y2, z2

x2-x1+1, y2-y1+1, z2-z1+1

22pqrst 161122p_fr

80.70

641, 462, 240

974, 1076, 808

334, 615, 569

22pqrst 161122q_fr

70.24

617, 400, 210

1112, 1066, 809

496, 667, 600

ed_1
set_ext.

22pqrst 161122r_lm

70.24

125, 36, 78

636, 591, 583

512, 556, 506

22pqrst 161122r_sbs

70.24

81, 16, 77

636, 591, 583

556, 576, 507

22pqrst 161122s_lm

80.70

248, 119, 109

568, 405, 583

321, 287, 475

22pqrst 161122s_sbs

80.70

201, 104, 0574, 405, 583

374, 302, 584

22pqrst 161122t_hr

80.70

310, 109, 111

649, 484, 564

340, 376, 454

group

nm

x1, y1, z1

x2, y2, z2

x2-x1+1, y2-y1+1, z2-z1+1

80.70

658, 526, 243

972, 1008, 808

315, 483, 566

ed_2
set_ext.

22pqrst 161122p_fr
22pqrst 161122q_fr

70.24

663, 518, 212

1072, 969, 809

410, 452, 598

22pqrst 161122r_lm

70.24

166, 76, 78

609, 588, 583

444, 513, 506

22pqrst 161122r_sbs

70.24

183, 26, 78

610, 588, 583

428, 563, 506

22pqrst 161122s_lm

80.70

256, 124, 110

549, 381, 583

294, 258, 474

22pqrst 161122s_sbs

80.70

210, 114, 0548, 382, 583

339, 269, 584

22pqrst 161122t_hr

80.70

323, 137, 111

613, 439, 564

291, 303, 454

41 / 93
Date:Thu, 01 Dec 2016 14:26:13 +0900

From:

Tsukasa NAKANOTo:

MiyakeCc:

Kentaro UESUGI, Akira TSUCHIYAMA, matsuno j, "Aiko Nakato", Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, tkadokawa, "Takashi Sakurama", Minako HASHIGUCHI,
take, hoshino, uesugi

Subject: oct_sg

みやけさま、
みなさま、
GSJ/AIST のなかのです。
On Wed, 30 Nov 2016 17:57:38 +0900 Miyake wrote:> >> あと、1122s で、内側の重なっている部分と、外側の部分では違うのでしょうか？
>>> >これは 161122s_lm と 161122s_sbs の再構成画像のスライスの中央部（== サン> >プル回転⾓が 180° 異なっている投影画像の重複部分）とそれ以外の部分の CT> >値の違いに関する話ですかね？ それについてはまだ調べていません。
>> こちらになります。
> あと、sa でも、違うのかどうか気になりますけど。
OCT や FRCT では投影画像のデータを「混ぜ合わせて」画像再構成を⾏うので、その結果の CT 画像の上の値の空間分布の違いを⾒るのは難しいです。そこで、11/23 の E-mail で紹介した測定グループ 19qrs の場合と同様に、つなぎ合わせた sinograms（SGs）を調べてみました。この E-mail に添付した以下の４個のファイルを御覧下さい。
21hi.pdf
ついでに調べた測定グループ 21dfhi の FRCT の測定 161121[h,i] の SA でつなぎ合わせた SGs を並べてある。
22pq.jpg
測定グループ 22pqrst の FRCT の測定 161122[p,q] の SGs。
22[r,s].jpg
22pqrst の OCT の測定 161122[r,s]のそれぞれの LM および SBS でつなぎ合わせた SGs を並べた。
注E-mail の添付ファイルとしてはサイズが⼤き過ぎる感じだったので、測定グループ 22pqrst の図の
PDF を JPEG に変換しました。計算機 vrm のディレクトリ/media/disk/tsukasa/1611/test/ の下にオリジナルの PDF の 22pqrs.pdf を置いてあるので、必要ならそれを御覧下さい。

42 / 9321hi.pdf に⽰した測定グループ 21dfhi の FRCT の SGs には特に問題はない感じです。あえて⾔うなら、差分 SG の⽐較から、投影数が少ない 161121h に⽐べて投影数が多い 161121i の⽅が投影画像の重複した部分の違いが⼤きいように思われます。測定時間の違いがその原因かもしれません。
問題は測定グループ 22pqrst の SGs です。まず、22pq.jpg に⽰した FRCT の測定 161122[p,q]の SGs ですが、これらのどちらにも「つなぎ⽬」は⾒えず、特に問題がないと思えるかもしれません。しかしながら、これらから再構成した画像はサンプル回転軸の位置が間違っています（それに伴う偽像が発⽣しています）。
実は 22pq.jpg の差分 SGs からもその兆候を読み取れます。つまり、差分 SGs の物体像の部分の両端に⽩もしくは⿊の「輪郭線」が発⽣しており、そのため、投影画像のつなぎ合わせが上⼿く⾏っていない。これはMSD による位置決めの⼿法が原因ではなく、測定画像（投影画像や SG）に含まれていた物体像と空気の境界に出現している「屈折？」のために⽣じたのではないかと思われます。
22[r,s].jpg に⽰した OCT の測定 161122[r,s] の SGs ではこの「屈折？」の影響がさらに顕著です。特に
161122r では投影画像のつなぎ合わせがダメで、つなぎ合わせた SGs につなぎ⽬が⾒えます。⼀⽅、161122s
では「屈折？」があるにも関わらずつなぎ合わせは上⼿く⾏っています。ただし、昨⽇の２通めの E-mail で⾔及した SBS でつなぎ合わせた SGs から明らかなように、サンプル回転⾓が 180°未満とそれ以上の投影
画像のノイズレベルの違いが⼤問題です。
なお、Offset CT のテスト⽤に回転軸から⼤きく離れた位置にサンプルを置いた測定 161122[r,s] においてより顕著なので、ぼくはこの「屈折？」をサンプルによる屈折だと考えましたが、これは FZP のフォーカスの影響かもしれません。いずれにせよ、OCT ではサンプルを回転軸の直上に置くことが重要ですね。

以下は今回⾏った処理に関する補⾜です。
(1)
再構成画像の上のサンプル像の位置決めが出来ていなかった測定グループ 22pqrst の SG のスライスの位置の決定には昨⽇の２通⽬の E-mail で紹介した ED 処理の結果を使いました。
ED 処理で決めたサンプルの上端に相当するスライスの位置 z0
測定番

z0

スライスの層厚 dz（単位は nm）

161122p

240

80.70

161122q

210

70.24

161122r

7770.24

161122s

109

80.70

161122t

110

80.70

通常の HRCT の測定 161122t の投影画像を調べるとサンプルの上下中央のスライスの位置が 275 だったので、上記の z0 と dz の値を⽤いてそれ以外の測定の再構成画像の上のサンプルの上下中央のスライスの位置
Z を計算しました：Z = z0 + ( 275 - 110 ) * 80.70 / dz

43 / 93(2)
この E-mail で紹介した OCT と HRCT の SG の図は OCT ⽤の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zipに⼊っている C-shell script "plot.oct_sg*" を使って作成しました。
plot.oct_sg RAW OX OY

Z {title}> SG.psplot.oct_sg_sa RAW OXOY Z {title} > SG.ps

plot.oct_sg_sbs RAW OXOY Z {title} > SG.ps

ただし、これらはそれぞれ LM、SA もしくは SBS で投影画像をつなぎ合わせます。また、起動パラメータの意味は以下の通りです。
RAW
HiPic 形式の測定画像などが⼊っているディレクトリの名前
OX と OY
投影画像のつなぎ合わせに使うオフセットの値Z作成する SG のスライス番号
SG.ps
SG の図の PostScript コードを書き込むファイルの名前
title
図の上部に付加するタイトル⽤の⽂字列。指定を省略すると RAW として指定した⽂字列を使う。
最終的には図の上部に以下の⽂字列が記される。
title /

Ox = OX / Oy = OY / z = Z / NV views

ただし、NV は OCT の画像再構成で使⽤する投影数。
とりあえず以上です。
添付ファイル 21hi.pdf、22pq.jpg、22r.jpg および 22s.jpg
後に⽰す 12/7 の E-mail の添付ファイル 21hi.jpg、22pq.jpg、22r.jpg および 22s.jpg を御覧下さい。

44 / 93
Date:Thu, 01 Dec 2016 19:06:49 +0900

From:

Tsukasa NAKANO

To:

⼟`⼭ 明

Cc:

Miyake, Kentaro UESUGI, Akira TSUCHIYAMA, matsuno j, "Aiko Nakato", Sugimoto Miyama,
北⼭晃, 中村隆太, kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Akihisa TAKEUCHI, Masato Hoshino, Masayuki Uesugi

Subject: Re: oct_frct_hrct

つちやまさま、みなさま、なかのです。E-mail ありがとうございます。
On Thu, 1 Dec 2016 17:19:13 +0900 ⼟`⼭ 明 wrote:> > （３）> > ８ keV で撮影した画像のリングアーチファクトは上記の順で減っているように> >> ⾒えます。つまり、FZP-CT の画像の画質向上に full rotation CT が有効だと> >> ⾔えそうです。しかしながら、投影数を２倍にしても画質は概ね同じ感じです。
>> 確かに。
> 贔屓⽬に⾒ると、1121i がリングアーチファクトによると思われる、むらが少ないように思いますね。
>> FRCT にするとリングアーチファクトが軽減（激減？）するということが、よくわかりますね。
>> 161121h と 161121i の違いですが、僕にはリングアーチファクトよりもコントラ> ストが違うように思えます。すなわち、投影数の少ない 161121h の⽅が、画像が> シャープな（コントラスト差がはっきりしている）ように思えます。ヒストグラム> どうなっていますかね？
昨⽇の１通⽬の E-mail で紹介したように FRCT の測定 161121[h,i] の再構成画像の画素値ヒストグラムのデータファイル "161121[h,i]_fr.csv" は計算機 vrm に置いてあるので「それを⾃分で⾒てくれ」と⾔いたいところですが、今回は⼤サービスでその図を描きました。この E-mail に添付した 21hi_hg.pdf です。上段の図は３次元画像全体の、下段の図は昨⽇の E-mail に添付した 21dfhi.pdf に⽰したスライスのヒストグラム（ただし、トリミングしていないスライス画像のヒストグラム）です。画像全体とスライスのヒストグラムのピークの位置は概ね合っていますが、全体のヒストグラムは 161121i に⽐べて 161121h のピークの CT 値がやや低く、スライスではそれは逆にやや⾼くなっています（これらの違いは誤差の範囲内かもしれません）。
これらの図からはつちやまさんが⾔われているようなコントラストの違いはわかりません。
それよりも気になるのは投影数が少ない 161121h の空気の部分のヒストグラムの⽅が 161121i のものよりも⾼いことです。これは今⽇の１通⽬の E-mail に記した差分 sinograms の⽐較に基づくぼくの予想（投影45 / 93
数が多い == ⻑い測定時間の 161121i の⽅が 161121h よりも noizy）とは逆の様相を⽰しています。
On Thu, 01 Dec 2016 14:26:13 +0900 Tsukasa NAKANO wrote:> 21hi.pdf に⽰した測定グループ 21dfhi の FRCT の SGs には特に問題はない感じです。
> あえて⾔うなら、差分 SG の⽐較から、投影数が少ない 161121h に⽐べて投影数が多い> 161121i の⽅が投影画像の重複した部分の違いが⼤きいように思われます。測定時間の> 違いがその原因かもしれません。
とり急ぎ、
添付ファイル 21hi_hg.pdf
後に⽰す 12/7 の E-mail の添付ファイル 21fhi_hg.pdf を御覧下さい。

Date:Thu, 01 Dec 2016 19:14:20 +0900

From:

Tsukasa NAKANO

To:

⼟`⼭ 明

Cc:

Miyake, Kentaro UESUGI, Akira TSUCHIYAMA, matsuno j, Aiko Nakato, Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, tkadokawa, "Takashi Sakurama", Minako HASHIGUCHI,
Akihisa TAKEUCHI, hoshino, Masayuki Uesugi

Subject: Re: do-not-use-oct_sbs

つちやまさま、みなさま、なかのです。先程の E-mail にも書いたように、ここでの議論に関連するデータファイルはすべて計算機 vrmの上にコピーしてあるので、FRCT の測定 161122p の CT 画像の観察やその画素値ヒストグラムのチェックなどはそちらで⾏って下さい。よろしくお願いします。とり急ぎ、
On Thu, 1 Dec 2016 18:32:12 +0900 ⼟`⼭ 明 wrote:> >> 1122p (8kev, 360degree/1800)だけが⼤きくずれているようにみえるのは、> >> バックグラウンドが⾼くなっている影響はあるのでしょうか？
>>> > そうかもしれません。ただ、forsterite のピークの「幅」から考えると、同じ> > full rotation（FR）CT を⾏った測定 161122[p,q] では７keV の 161122q の⽅が> > noizy なはずですが、...。
>> 161122p の視野のサイズを他のものと同じになるように狭めた時に、> ヒストグラムのピークトップの位置がどうなるかをみれば、「バック> グラウンドが⾼くなっている影響」をチェックすることができるはずです。
>

46 / 93> >> そもそもなぜ⾼くなっているのかも気になりますけど。
>>> > FRCT の測定 161122[p,q] の視野のサイズ（== 投影画像の画素数）は他のもの> > よりも⼤きいので CT 値が０付近の部分のヒストグラムが⾼い（画素数が多い）> > のは当然ですが、forsterite のピーク付近までそれが⾼い理由は不明です。
>> これは、不思議ですね。
> 僕が上に書いたように、161122p の視野のサイズを他のものと同じになるように> 狭めた時のヒストグラムではどうなっているのか？
>> CT 像の⽐較も⾒たいのですが、お願いします。

Date:Fri, 02 Dec 2016 10:22:35 +0900

From:

Tsukasa NAKANOTo:"TAKEUCHI, Akihisa"

Cc:

MIYAKE, Akira TSUCHIYAMA, Kentaro UESUGI, matsuno j, Aiko Nakato, Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Masato Hoshino, Masayuki Uesugi

Subject: Re: 180degCT と 360degCT

たけうちさま、なかのです。E-mail ありがとうございます。
> 360degCT が効果があったのは、offsetCT のときだけだったでしょうか、> それとも通常の CT でも屈折の効果はキャンセルされてたでしょうか。
先⽉の実験ではすべての FRCT の画像再構成を Offset CT ⽤のプログラムを⽤いて⾏いました。その場合に「投影画像の⽔平⽅向のオフセットの値（Ox）」が０なら 360 °回転で通常の CT をしたことになるので、これから FRCT の画像再構成で使⽤した Ox の値を調べてみます。とり急ぎ、P.S.
先⽉の実験の設定は「光軸の位置 == 検出器の中央」になっていたのでしょうか？ HiPic の画像ファイ
ルにコメントとして記録されている撮影時の ROI のデータを使えば画像再構成に使⽤した投影画像上の光軸の位置を知ることができる？

47 / 93
On Fri, 2 Dec 2016 09:29:15 +0900 "TAKEUCHI, Akihisa" wrote:> 関係者みなさま>> ⽵内です。
> 先⽇の⼟'⼭さんのビームタイムで、> 360degCT は屈折の効果がキャンセルされて 180degCT よりきれい？
> というお題を頂き、現在進⾏中の僕のビームタイムの中でそれについて検証しました。
> 結果を添付します。試料がナニかについては追求しないでください。
>> 結果についてですが、おそらく> ・普通のセットアップ（回転軸が光軸上）ではそんなことはない（添付⼀枚⽬）> ・でも光軸をずらすととても効果あり（添付２枚め）> ということのようです。
>> 添付２枚めはようするに OffsetCT と同じ条件になるのですが> この場合は試料回転軸のオフセットに併せて検出器も移動してるので> 検出⾯上の回転中⼼は検出⾯の中⼼のままです。
>> 当初僕がでっちあげた仮説は、>デフォーカスの効果をキャンセルするんじゃ？

> というものでしたが、そうでなく、>

FZP の光軸から外れる（視野の外側）につれて出て来るフリンジをキャンセル> しているようです。おそらくですが、でもこちらの考えの⽅がよりしっくりします。
>> そこで、先⽇⼟'⼭さん実験で参加されていたみなさんにお聞きしたいのですが、> 360degCT が効果があったのは、offsetCT のときだけだったでしょうか、> それとも通常の CT でも屈折の効果はキャンセルされてたでしょうか。
> 上記の結論だと、通常の CT では 360degCT の効果はあまりなかったはずですが、> もし効果はやっぱりあった、というのであれば、ひょっとしたら通常の CT でも> ひょっとして光軸がずれていたのかもしれません。まあ、結果的にそれがよかったのですが。
>> この⽅法が使い物になるのかどうかは、また次のビームタイムで余裕があれば試してみてください。
> 論⽂にできるかどうかは、うーんもうちょっと考えてみないといけないかも。
>> 以上です。

48 / 93
Date:

Fri, 02 Dec 2016 12:59:52 +0900

From:

Tsukasa NAKANO <tsukasa.nakano@aist.go.jp>

To:"TAKEUCHI, Akihisa"

Cc:

MIYAKE, Akira TSUCHIYAMA, Kentaro UESUGI, "matsuno j", Aiko Nakato, Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, "kadokawa tokiyuki", Takashi Sakurama, Minako HASHIGUCHI,
Masato Hoshino, Masayuki Uesugi

Subject: Re: 180degCT と 360degCT

たけうちさま、
GSJ/AIST のなかのです。先⽉の SPring-8 BL47XU での FZP-CT 実験においてほぼ連続して⾏った４測定の投影画像と再構成画像の「幾何学」を調べました。
180°回転の通常の CT（half rotation CT、HRCT）
測定 161121d

X 線エネルギー 7 keV、投影数 1800 views

161121f

8 keV、1800 views

360°回転の CT（full rotaion CT、FRCT）
161121h

8 keV、1800 views → 900 views にして再構成

161121i

8 keV、3600 views → 1800 views

まず、測定で得た HiPic 画像のファイルにコメントとして書き込まれていた、測定時に指定した ROI のsubarray の左端の座標（x0）と横幅（Nx）を調べました（値 Nx が測定した HiPic 画像の横画素数と⼀致することを確認しました）。ROI の設定前の検出器の横幅が 2048 だとすると、測定した HiPic 画像の上での光軸の位置は LC = 2048 / 2 - x0 のはずです。これらの値と画像再構成⽤に推定したサンプル回転軸の（HiPic
画像の上での）位置 RC は以下の通りです。ただし、Offset CT ⽤の画像再構成プログラムを使った FRCT の測定に対してはペアの投影画像の x ⽅向のオフセット値 Ox（と Nx の値）から RC を算出しました（これはプログラム "oct_sg2tg*" でも⾏っています）
。
測定

x0

Nx

LC

RC

161121d

397

1348

627

〜 685

161121f

566

916

458

497

161121h

566

916

458

498.5 = (916+82-1)/2← Ox = 82

161121i

566

916

458

499.5 = (916+83-1)/2← Ox = 83と⾔う訳で、先⽉の測定では光軸と回転軸の位置は⼀致していないようですね。
ついでながら、これらの測定の再構成画像のスライス画像を並べた（11/30 の E-mail に添付した PDF を少しだけ改変した）21dfhi_rc.pdf をこの E-mail に添付します。161121d のスライス画像の上にサンプル回転
軸の位置を⽰しました。これ以外の測定でもそれぞれのスライス画像の上の同じ位置に回転軸があるはずです（これらはいずれも、回転軸からの距離が同じになるようにトリミングしたスライス画像です）。これらの49 / 93
画像からわかるように回転軸を中⼼としたリングアーチファクトは HRCT よりも FRCT において、また、
FRCT 同⼠では投影数が多い⽅で少なくなっています。とり急ぎ、On Fri, 2 Dec 2016 11:00:01 +0900 "TAKEUCHI, Akihisa" wrote:> 中野様>> ⽵内です。PS のご質問について。
> 先⽉の実験も含めて毎回、光軸＝検出器中央、になっているはずです（ROI なしの場合で）。
> それは毎回、そのように慎重にアライメントしているので、間違ってはないと思います。
>> ただ、先⽇の実験データで Ox = 0 でも屈折の効果が消えていたとなれば、> アライメントがやっぱりいつのまにかずれていたか、> それとも、僕の２つ⽬の仮説もまちがいで光軸上でも効果がでる条件がちゃんとあるのか・・・。
添付ファイル 21dfhi_rc.pdf

50 / 93
Date:

Fri, 02 Dec 2016 13:40:09 +0900

From:

Tsukasa NAKANOTo:"TAKEUCHI, Akihisa"

Cc:

MIYAKE, Akira TSUCHIYAMA, Kentaro UESUGI, "matsuno j", Aiko Nakato, Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, "kadokawa tokiyuki", Takashi Sakurama, Minako HASHIGUCHI,
Masato Hoshino, Masayuki Uesugi

Subject: Fw: Re: 180degCT と 360degCT

たけうちさま、なかのです。先程の E-mail の訂正と追加です。
訂正
>

測定

x0

Nx

LC

RC

>

161121d

397

1348

627

〜 685

>

161121f

566

916

458

497

>

161121h

566

916

458

498.5 = (916+82-1)/2← Ox = 82>

161121i

566

916

458

499.5 = (916+83-1)/2← Ox = 83

上記の 161121i の RC の値 499.5 は間違いで 499.0 でした。
追加
先⽉の実験の FRCT で画像再構成に指定した Ox の値を調べたら、すべて０ではありませんでした。つまり、先⽉の実験の FRCT は Offset CT として処理せざるをえません。

Date:Fri, 2 Dec 2016 14:36:38 +0900

From:"TAKEUCHI, Akihisa"

To:

Tsukasa NAKANOCc:

MIYAKE, Akira TSUCHIYAMA, Kentaro UESUGI, "matsuno j", Aiko Nakato, Sugimoto Miyama,
kitayama.akira.46a, 中村隆太, "kadokawa tokiyuki", Takashi Sakurama, Minako HASHIGUCHI,
Masato Hoshino, Masayuki Uesugi

Subject: Re: 180degCT と 360degCT

中野様
検証ありがとうございました。
検出器上での LC と RC の違いが５０ピクセル程度
これは、試料⾯に直すと２ミクロン弱、この程度であれば回転軸はほぼ光軸上に乗っていたと考えてよい範囲です。
あとは FZP と検出器中⼼が本当に LC 光軸上にあったかどうかですが、これはもう確かめようがないですね。

51 / 93
ただ、アライメントの⼿順を思い返しても、光軸がずれていたようには思えません。
頂いた再構成像からは、HR と FR で、屈折の出⽅の違いまではよく分かりませんね。
リングアーティファクトが FR で少なくなる件については、どうしてなのかまだちょっとノーアイデアです。SN が増えたから、くらいしか。
引き続き、通常の CT で HR と FR の屈折の効果の違いが出ていたかどうかについて、
ご報告お待ちしてます。

Date:Sat, 03 Dec 2016 13:22:09 +0900

From:

Tsukasa NAKANO

To:

⼟`⼭ 明

Cc:

Miyake, Kentaro UESUGI, Akira TSUCHIYAMA, matsuno j, Aiko Nakato, Sugimoto Miyama,
北⼭晃, 中村隆太, kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Akihisa TAKEUCHI, Masato Hoshino, Masayuki Uesugi

Subject: Re: oct_frct_hrct

つちやまさま、なかのです。
> 161121i、161121h を⽐べると> (1) ピーク値は⼀緒（だと思う）だが、ピーク⾼さと幅が異なる> （投影数の多い 161121i の⽅がシャープ）> (2) ピークの CT 値、ピーク⾼さ、幅が異なる> （投影数の多い 161121i の⽅がシャープ、CT 値がやや⾼い）> (3) 両者ほぼ変わらないが、投影数の多い 161121i の⽅がややピークが⾼く幅が狭いように⾒える> （つまりシャープ）
画像を⾒ただけでこれらに気づいたのは「⼼眼 or 神眼」ですか？
> ピークのシャープさだけを考えると、投影数の多い⽅がシャープなので、> これは予想通りということなのでしょう。
予想通りと⾔うよりも、⾃分の印象に合う都合の良い証拠を探した気がする。
> 161121i の⽅が 161121h について CT 像の⾒かけは、そんなに差はない、> ということでまあいいかと思います。
これが普通の感じ⽅ですね。とり急ぎ、

52 / 93
On Fri, 2 Dec 2016 17:12:03 +0900 ⼟`⼭ 明 wrote:> 中野さん> ありがとうございます。
>> 撮影⽅法が違うと、ヒストグラムが微妙に変わっています。
>> このサンプルには、⼤まかに３つのピークがあります。
> (1) 空気 CT 値=0 cm-1
> (2) Olivine or enstatite 粒⼦

CT 値〜105 cm-1

> (3) Matrix CT 値〜230-240 cm-1>
> 161121i、161121h を⽐べると> (1) ピーク値は⼀緒（だと思う）だが、ピーク⾼さと幅が異なる （投影数の多い 161121i の⽅がシャープ）> (2) ピークの CT 値、ピーク⾼さ、幅が異なる（投影数の多い 161121i の⽅がシャープ、CT 値がやや⾼い）> (3) 両者ほぼ変わらないが、投影数の多い 161121i の⽅がややピークが⾼く幅が狭いように⾒える> （つまりシャープ）>> ピークのシャープさだけを考えると、投影数の多い⽅がシャープなので、> これは予想通りということなのでしょう。
> (2)のピークの CT 値が異なるのは、気になります。なぜなのだろう？中野さんもご指摘のように、> １スライスのヒストグラムでは、さらに⼤⼩関係が逆になっているのも、不思議ですね。
> そもそも、CT 値にはこれくらいの誤差があるということなのでしょうか？
>> 161121i の⽅が 161121h について CT 像の⾒かけは、そんなに差はない、> ということでまあいいかと思います。

Date:

Wed, 07 Dec 2016 11:25:48 +0900

From:

Tsukasa NAKANOTo:

Miyake, Akira TSUCHIYAMA

Cc:

Kentaro UESUGI, matsuno j, Aiko Nakato, Sugimoto Miyama, 北⼭晃, 中村隆太,
kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Akihisa TAKEUCHI, Masato Hoshino, Masayuki Uesugi

Subject: oct_sg+tg+hg_rev.figs

みなさま、
GSJ/AIST のなかのです。先⽉の SPring-8 BL47XU FZP-CT 実験で撮影した Offset CT のテスト⽤の画像に関するこれまでの E-mails に添付した図の改定版を送ります。E-mail に添付するためサイズが⼤きい PDF
を低分解能の JPEG 画像に変換しました。そのオリジナルのものを含む PDFs はすべて計算機 vrm のディレ
クトリ/media/disk/tsukasa/1611/test/に置いてあります。

53 / 9321hi.jpg
12/1 の E-mail（Subject: oct_sg）に添付した「測定 161121[h,i] の simple average（SA）でつなぎ合わせた full rotation（FR）CT ⽤のものなどの sinograms（SGs）を並べた図」の下部に HRCTの CT 画像を付加した図。vrm の上のオリジナルの PDF のファイル名は 21hi.pdf。
22pq.jpg、22r.jpg および 22s.jpg
12/1 および 12/2（Re: oct_sg）に添付した「測定 161122[p,q,r,s] の FRCT や Offset CT ⽤のものなどの SGs を並べた図」の下部に CT 画像を付加。これら３画像は vrm の 22pqrs.pdf の各ページに対応。
21fhi_hg.pdf
12/1 の E-mail（Re: oct_frct_hrct）に添付した「測定 161121[h,i] の FRCT の３次元 CT 画像全体
および１枚だけのスライス画像の画素値ヒストグラムの図」に測定 161121f の half rotation（HR）
CT の画像のヒストグラムを追加した図。
とりあえず以上です。
添付ファイル 21hi.jpg

54 / 93
添付ファイル 22pq.jpg

添付ファイル 22r.jpg

55 / 93
添付ファイル 22s.jpg

添付ファイル 21fhi_hg.pdf

56 / 93
Date:

Wed, 07 Dec 2016 17:21:53 +0900

From:

Tsukasa NAKANO

To:

Kentaro UESUGICc:matsuno j, Aiko Nakato, Sugimoto Miyama, 北⼭晃, 中村隆太, kadokawa tokiyuki,
Takashi Sakurama, Minako HASHIGUCHI, Akihisa TAKEUCHI, Masato Hoshino,
Masayuki Uesugi, Miyake, Akira TSUCHIYAMA

Subject: oct_f2h

みなさま、
GSJ/AIST のなかのです。SPring-8 でサンプルを 360°回転して撮影した full rotation（FR）CT の測定デー
タからサンプル回転⾓が180°未満（０ 〜 180°）の「normal view」の CT ⽤のデータ
180°以上（180 〜 360°）の「reverse view」の CT ⽤のデータを取り出し、それらの２セットの half rotation（HR）CT の画像再構成を⾏うための３種類（合計６個）のC-shell scripts と C ⾔語のコードを書きました。
rc_f2h_[n,r].csh
FRCT の測定データから normal もしくは reverse view のデータを取り出し、通常の HRCT ⽤の Cshell script（check.rc；HASPET ⽤の center.csh も check.rc を使⽤している）と同じ「０°と 180°の投影画像の位置決め」によってそれそれの画像再構成に必要なサンプル回転軸の位置（いわゆる「センター値」）などを推定する。
tg_f2h_[n,r].csh
FRCT の測定データから normal もしくは reverse view のそれぞれの HRCT のすべてのスライス画
像を再構成する。
sg_f2h_[n,r]
FRCT の測定データから hp_sg などの HRCT ⽤の既存のプログラムで取り出したサンプル回転⾓が０〜360°の sinogram のバイナリデータの
「前半」
と「後半」
を分離して、normal もしくは reverseview の sinogram のそれぞれのデータを抽出する（C ⾔語コードをコンパイルした）プログラム。
注
上記の C-shell scripts "*.csh" などの名前に含まれている⽂字列 "[n,r]" は normal と reverse の頭⽂字である「"n" もしくは "r"」を表しています。
これらのファイル（および sg_f2h_[n,r] の Windows ⽤の実⾏ファイル）を Offset CT ⽤の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zipに⼊れておきました。
それぞれの起動法などは以下の通りです。

57 / 93(0) 実⾏に必要なプログラム
C-shell scripts "rc_f2h_[n,r].csh"と"tg_f2h_[n,r].csh" はぼくが以前に書いた下記のプログラムを使っています。また、プログラム"sg_f2h_[n,r]"を⽤いた「画像再構成のテスト」でもプログラム sg2tgs を使います。
これらはすべて HASPET ⽤のプログラムのパッケージに含まれています。
HASPET ⽤のプログラムのインストール法を含む説明書
http://www-bl20.spring8.or.jp/~sp8ct/tmp/haspet.pdf
実⾏に必要なプログラム
dit_tbl と pig および CPU ⽤の sg2tgs と hp2tg_t
説明書
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.pdf
書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.zip--- 以下のプログラムの説明書はないので、書庫ファイルだけを記した。
CUDA GPU ⽤の sg2tgs と hp2tg_t
http://www-bl20.spring8.or.jp/~sp8ct/tmp/cuda_cbp.taz
rar_gray（２次元画像の直⾓回転⽤プログラム）
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rar.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rar.zip
rmsd_2d.CC（２次元画像ペアの RMSD 計算プログラム）
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rmsd.tar
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rmsd.zip(1) C-shell script "rc_f2h_[n,r].csh"
起動法
rc_f2h_[n,r].csh

HiPic/ {x0} {y1 y2} Rx Ry

起動パラメータの意味
HiPic/
HiPic 画像などの測定データが⼊っているディレクトリの名前
x0
assumed rotation center の位置。測定画像の横画素数を N として、指定省略時は x0 = (N -1 )/2だと⾒なされる。
y1 と y2
サンプル回転⾓が０°と 180°の投影画像の位置決めに使う測定画像の縦⽅向の範囲。測定画像の縦画素数を L として、これらの値の指定省略時は y1 = 0 かつ y2 = L - 1 とする。
Rx と Ry投影画像の位置決めの際の横と縦⽅向の探索範囲を指定する０〜１の値。２枚の投影画像の横と縦の画素数を N と L = y2 - y1 + 1 だとすると、０°の画像に付随する座標系において 180°の画像

58 / 93の原点 Ox と Oy を以下の範囲の整数座標値とした配置における２画像の重複部分の投影値の差の⾃乗の平均の平⽅根（RMSD）を計算して位置決めを⾏う。
Ox = Dx0 - Dx 〜 Dx0 + Dx
Oy = - Dy 〜 Dy
ただし、
Dx0 = ２・x0 - (N - 1)
Dx = Rx・N ← 探索範囲の横幅の半値は Rx・N
Dy = Ry・L← 探索範囲の縦幅の半値は Ry・L = Ry・(y2-y1+1)

なお、HASPET ⽤の center.csh は check.rc に Rx = 0.5 と Ry = 0.25 を指定して HRCT のサンプル回転軸を推定するように設定してある。
出⼒
rc_f2h_[n,r].csh は以下の６個の値が並んだ⾏を標準出⼒に書き出す。
[1] Cx1 = (N - 1 + Dx0 - Dx) / 2 = x0 - Dx / 2 = x0 - Rx・N / 2[2] Cx2 = (N - 1 + Dx0 + Dx) / 2 = x0 + Dx / 2 = x0 + Rx・N / 2[3] Cx0：サンプル回転軸の位置（いわゆる「センター値」）の推定値[4] Oy1 = - Dy = - Ry・(y2-y1+1)[5] Oy2 = Dy = Ry・(y2-y1+1)[6] Oy0：２枚の投影画像の「縦ズレ」の推定値
注
Cx0 == Cx[1,2] もしくは Oy0 == Oy[1,2] になっている場合、rc_f2h_[n,r] によるサンプル回転軸の位置の推定は失敗です。
(2) C-shell script "tg_f2h_[n,r].csh"
起動法
tg_f2h_[n,r].csh

HiPic/ DO_base DO_step dr {f0 df} BPS

TG/ > TG.log

起動パラメータの意味
HiPic/
HiPic 画像などの測定データが⼊っているディレクトリの名前
DO_base と DO_step
スライス番号０における detector offset（DO）とスライス番号ごとの DO の増分の値。スライス番号に依存しない単⼀のセンター値 C を指定する場合は DO_base = -C かつ DO_step = 0 とすれば良い。
dr
測定画像の正⽅形画素の辺⻑ == 再構成画像の⽴⽅体画素の辺⻑
f0 と df
再構成画像上の各画素の整数画素値 i と CT 値 f の関係式の係数値。
f = f0 + df × i
これらの指定を省略すると、再構成したスライス画像ごとにそこに出現した CT 値の最⼩値と最⼤値で正規化した係数値を使う。

59 / 93
BPS
再構成画像の画素値のビット数（Bit Per Sample）
TG/
再構成画像のスライス画像のファイルを格納するディレクトリの名前
TG.log標準出⼒に書き出される「再構成したスライスごとに以下の４個の値が並んだ⾏」のデータをリダ
イレクトするテキストファイルの名前。
[1] スライス番号[2] そのスライスの再構成に⽤いた DO の値[3,4] そのスライスに出現した CT 値の最⼩値と最⼤値(3) プログラム "sg_f2h_[n,r]"
起動法
sg_f2h_[n,r]SG_F.bin SG_H.bin

起動パラメータの意味
SG_F.bin
FRCT の測定データから抽出したサンプル回転⾓が０〜360°の sinogram のバイナリデータが⼊っているファイルの名前。"-"（負の符号）を指定すると標準⼊⼒からバイナリデータを読む。
SG_H.bin
SG_F.bin の sinogram データの前半（normal view の場合）もしくは後半（reverse view）を書き込む HRCT ⽤のバイナリファイルの名前。"-" を指定するとバイナリデータを標準出⼒に書き出す。
(4) sg_f2h_[n,r] を⽤いた画像再構成のテスト
FRCT ⽤の測定データから取り出した normal もしくは reverse view の HRCT において「サンプル回転軸の位置を様々に変えた画像再構成のテスト」を⾏う場合の処理⼿順は以下の通りです。まず、プログラム hp_sg
を⽤いてサンプルの像があるスライスの sinogram のバイナリデータを取得します。
hp_sg

161122p 405 sg_f.bin

ここでは、sg_f2h_r を使ってその reverse view の sinogram のデータを抽出してみます。
sg_f2h_rsg_f.bin sg_r.bin

プログラム sg2tgs を使えば指定した⼀連の DO の値のそれぞれを使った複数の画像の再構成（のテスト）を⼀度に⾏うことができます。
注
通常 CPU のシングルコア⽤やマルチスレッド⽤、および CUDA GPU ⽤などプログラム sg2tgs には複数のものがありますが、それらのどれも起動法は同じです。

60 / 93

mkdir

tg_r

sg2tgs sg_r.bin 80.70e-7 -738 -0.5 28 180

16 tg_r > tg_r.log

ただし、ここで指定したパラメータの意味は以下の通りです。
80.70e-7
cm 単位の検出器間隔 == 再構成画像の⽴⽅体画素の辺⻑-738、-0.5、28
再構成に使う DO の値の初期値、刻み幅と総個数。j 番⽬の DO の値はDO（j = 0 〜 27） = -738 - 0.5・jとなる。
180
度単位のサンプル回転⾓の初期値。180°を指定した reverse view の再構成画像は０°を指定した
normal view の画像と同じ向きになる。
16
再構成画像の画素値のビット数。再構成画像上の画素値 i と CT 値 f の関係式 f = f0 + df・i の係数
値 f0 と df の指定を省略したので、ここでは出現した CT 値の最⼩値と最⼤値で正規化した画素値になる。なお、f0 と df の値はコマンド⾏でビット数の直前に指定できる。
tg_r
DO の値の⼩数点 "." を下線 "_" に変えた⽂字列を名前とする再構成画像のファイルを書き込むディレクトリの名前。上に⽰したように、このディレクトリを事前に⾃分で作成しておく必要がある。
tg_r.log
再構成画像ごとに以下の３個の値を並べた⾏を書き込むファイルの名前。
[1] その画像の再構成に⽤いた DO の値[2,3] その再構成画像に出現した CT 値の最⼩値と最⼤値
なお、Linux や MacOSX などでは作業⽤のバイナリデータのファイルを作らずに上記の⼀連の処理をパイプ
ラインで連結して実⾏できます（Windows では不可；MSYS2 なら実⾏可能）
。
mkdir

tg_r

hp_sg 161122p 405 - |sg_f2h_r -

\

- |\

sg2tgs - 80.70e-7 -738 -0.5 28 180
とりあえず以上です。

16

tg_r > tg_r.log

61 / 93
Date:Thu, 08 Dec 2016 13:18:43 +0900

From:

Tsukasa NAKANO

To:

Kentaro UESUGI, Miyake, Akira TSUCHIYAMA

Cc:

Aiko Nakato, Sugimoto Miyama, 北⼭晃, 中村隆太, kadokawa tokiyuki, Takashi Sakurama,
Minako HASHIGUCHI, Akihisa TAKEUCHI, Masato Hoshino, Masayuki Uesugi, matsuno j

Subject: f2h_21hi+22pq

みなさま、
GSJ/AIST のなかのです。昨⽇の E-mail で紹介したプログラムやスクリプトを⽤いて先⽉の SPring-8
BL47XU FZP-CT 実験で得た FRCT の測定 161121[h,i] と 161122 [p,q] のそれぞれの測定データから normal
および reverse view の HRCT の画像（HRCT_[n,r] 画像）を再構成しました。
そのためにまず、プログラム oct_xy で推定したオフセット値 Ox から FRCT の測定画像上でのサンプル回転
軸の位置（RC；いわゆる「センター値」）RC_F を計算し、同じ測定画像を⽤いた HRCT の normal および
reverse view ⽤に推定した RC の値（RC_[n,r]）と⽐較してみました。
FRCT の測定画像の横画素数を N とすると、
FRCT のつなぎ合わせた投影画像の横画素数 NO= N + Ox

when Ox ≧ 0= N - Ox

when Ox ≦ 0

つまり、Ox の正負に関わらず NO = N + | Ox |FRCT の RC はつなぎ合わせた画像のちょうど中央にあるはずなので、その測定画像の上での位置 RC_F
= (NO - 1) / 2

when Ox ≧ 0= Ox + (NO - 1) / 2

when Ox ≦ 0

つまり、Ox の正負に関わらず RC_F = (N - 1 + Ox) / 2
値 RC_[n,r] は昨⽇の E-mail で紹介した rc_f2h_[n,r].csh で推定しました。この E-mail に添付した
rc_21hi+22pq.txt が測定 161121[h,i]と 161122[p,q]の３種類の RC の推定処理のログファイルです。各⾏の⾏頭は測定番号ですが、それ以外は⾏ごとに書式の差異があるので、測定 161122p の５⾏のログを例にして各⾏に並んでいる値の意味を説明します。
測定 161122p のログの１⾏⽬に並んでいる値
161122p

46

1480

1434

739.5

-2

１番⽬の値（161122p）：測定番号２番⽬（46）：プログラム oct_xy で推定したオフセット値、Ox
３番⽬（1480）：つなぎ合わせた投影画像の横画素数、NO４番⽬（1434）：測定画像の横画素数、N５番⽬（739.5）：前記の式で計算した RC の推定値、RC_F

62 / 93
６番⽬（-2）：プログラム oct_xy で推定したオフセット値、Oy161122p のログの２⾏⽬
161122p

n

358

1075

740.5

-202

202

-2

１番⽬（161122p）：測定番号２番⽬（n）：normal もしくは reverse view を指す頭⽂字
３番⽬（383）：RC_[n,r] の探索範囲の最⼩値、RC_[n,r]_1
４番⽬（1075）：RC_[n,r] の探索範囲の最⼤値、RC_[n,r]_2
５番⽬（740.5）：rc_f2h_[n,r].csh で推定した RC の値、RC_[n,r]６番⽬（-202）：投影画像ペアの縦ズレの探索範囲の最⼩値、Oy1
７番⽬（202）：投影画像ペアの縦ズレの探索範囲の最⼤値、Oy2
８番⽬（-2）：rc_f2h_[n,r].csh で推定した縦ズレ、Oy0161122p のログの３⾏⽬
161122p

n

358

1075

740.5

-142

142

-2

240

810

１〜８：値の違いはあるが、それぞれの意味は２⾏⽬のものと同じ
９（240）： RMSD 計算に使った投影画像の領域の上端の座標値、y1--- 注：y1 は投影画像の上のサンプル像の上端の位置
10 （810）：RMSD 計算に使った投影画像の領域の下端の座標値、y2--- 注：y2 は投影画像の最下端 == サンプル像の下端の位置
161122p のログの４〜５⾏⽬
161122p

r

358

1075

750

-202

202

-1

161122p

r

358

1075

750

-142

142

-1

240

810

それぞれ２〜３⾏⽬の値と同意。
という訳で、注⽬すべきは rc_21hi+22pq.txt の各⾏の５番⽬の RC の値です。
161121h161121i

161122p

161122q

RC_F

498.5

499

739.5

766

RC_n

498.5

499

740.5

768.5

RC_r

498.5

499

750

772.5

測定 161121[h,i] では RCs は測定ごとに同じ値です。⼀⽅、161122[p,q] の RC_[F,n,r] はすべて異なる値になっています（サンプルの像がある領域だけを使っても全体を使った場合と同じ値になるので、RC_[n,r]は信頼できる推定値？）。つまり、普通に考えると、161122[p,q] では測定中にサンプルが動いた可能性が⾼いです。しかしながら、昨⽇の E-mail で上杉くんが⾔及した sinograms の図にはこれに相当する異常は⾒当たりません。why ?
推定した RCs が同じ値だった測定 161121[h,i] の normal および reverse view の HRCT 画像（HRCT_[n,r]画像）の再構成を昨⽇の E-mail で紹介した tg_f2h_[n,r].csh を⽤いて⾏いました。この E-mail に添付した
21hi.tif にそれらの代表的なスライスの画像を FRCT のものと⼀緒に並べました。いずれもサンプル回転軸が同じ位置になるようトリミングした画像です。また、それらの FRCT および HRCT_[n,r] の３次元画像全体と代表的なスライス画像のそれぞれの画素値ヒストグラムの図をこの E-mail に添付した 21hi_hg.pdf に⽰し63 / 93
ました。画像 21hi.tif を⾒る限り FRCT の画像にはリングアーチファクトが少ないので、HRCT の画像よりも⾼画質だと⾔えます。また、21hi_hg.pdf の HRCT_[n,r] 画像のヒストグラムの CT 値 100（1/cm）付近のピークは FRCT 画像のものよりもわずかに低くかつ幅広になっています。つまり、FRCT 画像の⽅がHRCT_[n,r] 画像よりもコントラストが⾼いと⾔えそうです。これらの様相は予想もしくは期待通りです。
測定 161122[p,q] に対しては推定した相異なる RC_[F,n,r] をカバーする値域の様々な RC の値を指定した
HRCT_[n,r] 画像の再構成のテストを⾏ないました。昨⽇の E-mail で紹介した処理⼿順で再構成した複数のスライス画像をサンプル回転軸との位置関係が同じにようにトリミングし、それぞれをマトリックス状に並べました。この E-mail に添付した以下の４個の画像ファイルを御覧下さい。
161121p_[n,r]_405.tif
測定 161121p のスライス 405 の HRCT_[n,r] 画像を並べた画像
161121q_[n,r]_400.tif
測定 161121q のスライス 400 の HRCT_[n,r] 画像を並べた画像
注
ここで⽤いた 405 と 400 はそれぞれサンプル像の鉛直⽅向中央の概ね同じ位置にあるスライスです。
マトリックスの左上隅に並べた左上隅に「FRCT」と記したものは各測定の FRCT のスライス画像です（これもまた、サンプル回転軸が同じ位置になるようにトリミングした画像です）。それ以外は左上隅に記した数値（ただし、"_" を "." に読み替えて下さい）を RC として指定して再構成した HRCT_[n,r] の画像で、右下隅の「FRCT」、「normal」と「reverse」はそれぞれ値 RC_[F,n,r] を指定した画像であることを意味します。
これらのテストの結果の画像から推定値 RC_[n,r] は妥当である（それらの値を指定して再構成した画像は他の値のものよりもそれなりに⾼画質である）ことがわかります。念のため、上記のスライスから鉛直⽅向に±10μm 離れた位置にあるスライスに対する合計 2（位置）×２（[n,r]）×２（161122[p,q]）=８セットのテストも⾏ってみましたが、この結論は変わりません。RC_F の値も妥当な感じなので、これは 161122[p,q]で撮影した同じサンプルが測定中に⼀定の速度でクリープしていたためかもしれません。why ?
とりあえず推定値 RC_[n,r] を指定して、測定 161122[p,q] の HRCT_[n,r] の画像すべてを再構成しました。
この E-mail に添付した 22pqrst_nr.pdf がそれらの３次元 CT 画像全体の画素値ヒストグラムの図です。測定161121[h,i] の場合と同様に、HRCT_[n,r] 画像のヒストグラムのサンプル像（forsterite）のピークは FRCTのものよりもわずかに低くて幅広になっているように⾒えます。
以上のように FRCT の画像は HRCT_[n,r] のものよりも⾼画質（⾼コントラスト）になっているようです。
ただし、これは投影数の違いに伴う「スタッキング効果」かもしれません。
「FRCT 画像が⾼画質になることに理論的な根拠がないこと」や「FRCT は測定時のアクシデントに弱いこと」を考え合わせると、現時点では
FRCT ではなく「投影数を多くした HRCT」を⾏う⽅が良いかもしれません。
な お 、 今 回 再 構 成 し た HRCT_[n,r] 画 像 な ど も す べ て SPring-8 の 計 算 機 vrm の デ ィ レ ク ト リ/media/disk/tsukasa/1611/test/ にコピーしてあります。11/30 の E-mail で紹介した測定グループ 21dfhiと 22pqrst のディレクトリに「ext（extension）」として "[n,r]" を付けてそれぞれの測定の HRCT_[n,r] 画

64 / 93像などを格納しました。ただし、161121[h,i] と 161122[p,q] の HRCT_[n,r] 画像の再構成の全容を記した
ログは「測定グループ.log」ではなく「21hi.log」と「22pq.log」に書き込まれています（これらのログの書式はすべて同じです）。これまでの E-mails に添付したファイルも追加してあるので、必要なら vrm のディ
レクトリ/media/disk/tsukasa/1611/test/ の全体をダウンロードし直す⽅が良いかもしれません。その総容量は 18 GB です。とり急ぎ、On Wed, 30 Nov 2016 10:13:10 +0900 Tsukasa NAKANO wrote:> vrm.spring8.or.jp:/media/disk/tsukasa/1611/test/
>
>
>
>
>
>
>
>
>
>

raw/測定番号
測定データのディレクトリへのシンボリックリンク
測定グループ名.log
各グループの画像再構成のログファイル（1611_test.txt と同内容）
測定グループ名/測定番号_ext/
８ビット画素値の３次元再構成画像が⼊っているディレクトリ
測定グループ名/測定番号_ext.log
再構成画像の各スライスに出現した CT 値の最⼩値と最⼤値
測定グループ名/測定番号_ext.csv
８ビット画素値の再構成画像上の CT 値のヒストグラム>> なお、上記の ext を含む 1611_test.txt で⽤いた記号の意味は以下の通りです。
>>

ext（extension の略のつもり）>

full rotation（FR）CT や half rotation（HR）CT、および、「linear mixing（LM）>

もしくは side by side（SBS）でつなぎ合わせた投影画像を使⽤した OCT」のそれ>

ぞれを区別するための記号。

添付ファイル rc_21hi+22pq.txt
161121h

82

998

916

498.5

0

161121h

n

228.5

686.5

498.5

-169

169

0

161121h

r

228.5

686.5

498.5

-169

169

-1

161121i

83

999

916

499

-1

161121i

n

228.5

686.5

499

-169

169

-1

161121i

r

228.5

686.5

499

-169

169

-1

161122p

46

1480

1434

739.5

-2

161122p

n

358

1075

740.5

-202

202

-2

161122p

n

358

1075

740.5

-142

142

-2

161122p

r

358

1075

750

-202

202

-1

161122p

r

358

1075

750

-142

142

-1

161122q

99

1533

1434

766

-1

161122q

n

358

1075

768.5

-202

202

-1

161122q

n

358

1075

768.5

-150

150

-1

161122q

r

358

1075

772.5

-202

202

-1

161122q

r

358

1075

772.5

-150

150

-1

240

810

240

810

210

810

210

810

65 / 93
添付ファイル 21hi.tif

添付ファイル 21hi_hg.pdf

66 / 93
添付ファイル 161122p_n_405.tif

添付ファイル 161122p_r_405.tif

67 / 93
添付ファイル 161122q_n_400.tif

添付ファイル 161122q_r_400.tif

68 / 93
添付ファイル 22pqrst_nr.pdf

Date:Thu, 08 Dec 2016 16:18:24 +0900

From:

Tsukasa NAKANOTo:

MiyakeCc:

Kentaro UESUGI, Akira TSUCHIYAMA, Aiko Nakato, Sugimoto Miyama, 北⼭晃, 中村隆太,
kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Akihisa TAKEUCHI, Masato Hoshino, Masayuki Uesugi, matsuno j

Subject: Re: f2h_21hi+22pq

みやけさま、
なかのです。E-mail ありがとうございます。
> 161121[h,i]の n,r で同じで、161122[p,q] の n,r で異なるのがわからないのですが、> １８０度のときに、サンプルを逃して、white を取りに⾏くと思います。
> この影響はないのでしょうか？
測定のログファイル output.log を⾒ればわかりますが、3600 投影の 161121i ではサンプル回転⾓が 180°の時にも透過 X 線強度（I0）画像を撮影していますが、それ以外の 1800 投影の 161121h、161122[p,q]では測定前（０°）と後（360°）でだけ I0 画像を撮影しています。161122[p,q] では 180°でサンプルを脇に寄せていないようです。

69 / 93

> ⾮常に連続的なうごきなので、sinograms の図には、気づかないのですが、> 半分づつに分けると、RC の値がずれるとか。
多分、そうなのでしょうね。とり急ぎ、
On Thu, 8 Dec 2016 15:55:52 +0900 Miyake wrote:> 中野様 上杉様 ⽵内様> ありがとうございます。理解できていないのですが、>> >これは 161122[p,q] で撮影した同じサンプルが測定中に⼀定の速度でクリープ> >していたためかもしれません。why ?
>> 161121[h,i]の n,r で同じで、161122[p,q] の n,r で異なるのがわからないのですが、> １８０度のときに、サンプルを逃して、white を取りに⾏くと思います。
> この影響はないのでしょうか？
>> 今回の撮影で気になったのは、ステージを回転させ始めた際に、ほぼ必ず上下に> 動きました。このため、助⾛があったと思います。
> 1121h,i を取るときには、まだ⼤丈夫だったけど、こうした影響が蓄積してきて、> ほぼ最後の測定に顕著に出てきたということはないのでしょうか？
> ⾮常に連続的なうごきなので、sinograms の図には、気づかないのですが、> 半分づつに分けると、RC の値がずれるとか。

Date:Sat, 10 Dec 2016 15:37:58 +0900

From:

Tsukasa NAKANO

To:

⼟`⼭ 明

Cc:

Kentaro UESUGI, Miyake, Akira TSUCHIYAMA, Aiko Nakato, Sugimoto Miyama,
北⼭晃, 中村隆太, kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Akihisa TAKEUCHI, Masato Hoshino, Masayuki Uesugi, matsuno j

Subject: Re: f2h_21hi+22pq

つちやまさま、なかのです。E-mail ありがとうございます。ぼくが⾔っている「測定時のアクシデント」はつちやまさんが⾔われている「測定時に起こったアクシデント」ではありません。HRCT であろうが FRCT であろうが、投影数が多くなるにつれて測定時にアクシデントが発⽣するリスクが⾼まることは⾔うまでもありません。問題はそれを検知して軽減する⼿法（と⾔うか⼿順）を確⽴していないことです。通常の CT（HRCT）では以下の４通りの⼿法でそれを⾏うことができます。

70 / 93[1] プログラム hp2xta を使って測定データから total absorption（TA）の画像を作成し検討する。
[2] プログラム hp2DO（もしくは HASPET ⽤の C-shell script "DO.csh"）を使って０°と 180°の投影画像からスライスごとのサンプル回転軸の位置を推定し、その「傾き」やバラツキを検討する。
[3] C-shell script "check.rc"（もしくは HASPET ⽤の "center.csh"）を使って０°と 180°の投影画
像の位置ズレを調べる。
[4] サンプル回転軸の位置を様々に変えた「画像再構成のテスト」を⾏ない、それらの画像を⽐較する。
FRCT でもこれらの処理を⾏うことは可能です。ただし、[1] ぼくは時々TA を調べていますが、その重要性を理解していないためか、ぼくをまねてこの⼿法を実⾏する⼈は皆無です。
[2] これは [3] とも関連することですが、
FRCT の場合は０°と 180 °の２枚の投影画像だけでなく、
サンプル回転⾓が 180°づつ異なる画像のペアのそれぞれに対してスライスごとの回転軸の位置（の変化）を調べることになります。この処理は⼤変で、かつ、得られた情報をどのように特徴付
けて理解するか、...。
[3] 投影画像のペアのそれぞれの位置ズレの調査は Offset CT ⽤のプログラム oct_xy で調べることができます。また、それをもっと詳しく⾏うためのプログラム oct_cmp_180 も書きました。しかしながら、[2] と同様に、それらで得た情報をどのようにしてコンパイルするかが問題。
[4] Offset CT や FRCT の「画像再構成のテスト」は現時点でも可能ですが、通常の CT の場合の１個のパラメータ（センター値や DO）ではなく、２個のパラメータ（投影画像の x と y ⽅向のオフセ
ット値）を様々に変える必要があります。
このように、通常の CT の場合とは異なり、FRCT を含む Offset CT では測定時に発⽣したアクシデントに対応するのはかなり⼤変です。まずは「（以前に紹介した C-shell script "plot.oct_sg*" によって描いた）Offset
CT の sinograms を毎回チェックする」ですかね。とり急ぎ、On Sat, 10 Dec 2016 12:18:27 +0900 ⼟`⼭ 明 wrote:> 中野さん、皆さん>> ようやく、追いつきました。
> 確認の質問です（思い違いしていないといいのですが）。
>> 161121[h,i]（隕⽯サンプル）では、サンプル固定やステージが安定していて、> FRCT がうまくいっている（リングアーチファクトの減少など）。
> ⼀⽅、161122[q,r,p,s,t]（フォルステライト）では、おそらくサンプル固定や> ステージが安定していなくて("測定時のアクシデント")、FRCT では回転中⼼の> 推定がうまくいかないので、画質が悪いが、HRCT だとまあまの絵が出る。
>> という理解かと思いますが、これでいいでしょうか？
>> > 「FRCT 画像が⾼画質になることに理論的な根拠がないこと」

71 / 93> > や「FRCT は測定時のアクシデントに弱いこと」を考え合わせると、現時点では> > FRCT ではなく「投影数を多くした HRCT」を⾏う⽅が良いかもしれません。
>> 上と同じことかもしれま s んが、投影数が同じ場合に"測定時のアクシデント”は、> FRCT に⼤きな影響を与えるが、HRCT では影響はそれほどでもない、という> 理解でいいでしょうか？（"測定時のアクシデント”の起こる確率は、おおよそ> 投影数によると思うのですが）Date:Sat, 10 Dec 2016 16:12:45 +0900

From:

Tsukasa NAKANOTo:

MiyakeCc:"TAKEUCHI, Akihisa", Kentaro UESUGI, Akira TSUCHIYAMA, matsuno j, Aiko Nakato,
Sugimoto Miyama, 北⼭晃, 中村隆太, kadokawa tokiyuki, Takashi Sakurama,
Minako HASHIGUCHI, Masato Hoshino, Masayuki Uesugi

Subject: Re: SP8-BL47XU-FZP-CT-ring-again

みやけさま、
なかのです。Offset CT のテスト⽤の 12 測定の CT 画像の縦断⾯の画像などを並べた図を作りました。このE-mail に添付した 1611_xyz_test.pdf です。なお、これは昨⽇の E-mail に添付した 1611_xyz_Mt.pdf と1611_xyz_AV-3.pdf でも同様ですが、PostScript ファイルとして作った図を Adobe acrobat XI で変換する際に「サイズが縮⼩された PDF（⽇本語として変ですが、そう記されている）」にしました。また、計算機 vrmの/media/disk/tsukasa/1611/ には「縮⼩なし」の PDFs も置いてありますので、必要ならお使い下さい。
とり急ぎ、
On Fri, 09 Dec 2016 19:53:38 +0900 Tsukasa NAKANO wrote:> みやけさま、> なかのです。測定 161122[p,q,r,s,t] で撮影した forsterite の画像は昨⽇の E-mail に書いた
> サンプル回転軸の問題があるので低画質です。ぼくはそれらの Offset CT のテスト⽤の画像> を別扱いしているので、縦断⾯の作成は少々お待ち下さい。とり急ぎ、>> On Fri, 9 Dec 2016 19:26:02 +0900 Miyake wrote:> > 中野さま ⽵内様> > あとは、Fo（Forsterite）が均質試料ですね。
>>> > AV ⾙、magnetite ともに測定条件は違わないはずですけど・・・。
>>> > 三宅72 / 93
添付ファイル 1611_xyz_test.pdf（１ページ⽬）

73 / 93
添付ファイル 1611_xyz_test.pdf（２ページ⽬）

74 / 93
Date:

Tue, 20 Dec 2016 15:29:54 +0900

From:

Tsukasa NAKANO

To:

Kentaro UESUGI, MiyakeCc:matsuno j, Aiko Nakato, Sugimoto Miyama, 北⼭晃, 中村隆太, kadokawa tokiyuki,
Takashi Sakurama, Minako HASHIGUCHI, Masato Hoshino, Masayuki Uesugi,"TAKEUCHI, Akihisa", Akira TSUCHIYAMA

Subject: Ox_and_RC

みなさま、
GSJ/AIST のなかのです。
On Mon, 12 Dec 2016 07:55:16 +0900 Miyake wrote:> 中野様 ⼟`⼭様
> > > 161121[h,i]（隕⽯サンプル）では、サンプル固定やステージが安定していて、> > > FRCT がうまくいっている（リングアーチファクトの減少など）。
> > > ⼀⽅、161122[q,r,p,s,t]（フォルステライト）では、おそらくサンプル固定や> > > ステージが安定していなくて("測定時のアクシデント")、FRCT では回転中⼼の> > > 推定がうまくいかないので、画質が悪いが、HRCT だとまあまの絵が出る。
>> ふと思ったのですが、ステージの影響はどれも⼀緒だとすると、> 161121[h,i]（隕⽯サンプル）は、FIB から TEM-tomo まで考えて作ったホルダーを使っています。
> ⼀⽅ 161122[q,r,p,s,t]（フォルステライト）、AV ⾙は、W 針につけたものを、> SP8 で、真鍮台にシリコンチューブに固定しています。
> この固定⽅法のしかたが影響しているかもしれませんね。
先⽉の SPring-8 BL47XU FZP-CT 実験で撮影した測定 161122[p,q] の投影画像を⾊々と調べた結果、以下のことがわかりました。
[1] 161122[p,q] のサンプルもしくはそれを載せていたステージが測定中に動いていたことは確か。特に後者の場合にはサンプル回転軸が測定中に動いていたことになる。
[2] その回転軸の位置 RC(θ) が回転⾓θ =０〜 360°に関する２次関数
RC(θ) = A × θ^2 + B × θ + C
であれば、θが 180°異なる投影画像のペアのそれぞれに対して推定した⽔平（x）⽅向の画像のオフセット値の「時系列データ」Ox(θ) を説明できる（ただし、Ox(θ) のθはそれぞれのペアの画像のうちの 180°未満の回転⾓で撮影した投影画像の回転⾓の値；θ = ０〜180°）。
これらの結論に⾄った処理内容を以下に記します。

75 / 93(1) 測定 161122[p,q] の FRCT 画像の問題点
まず復習です。測定 161122[p,q] の CT 画像には以下の難点がありました[a] FRCT の再構成画像の画質が悪い。
[b] その測定の０〜 180°のサンプル回転⾓で撮影した「normal view」の投影画像と 180〜360°の「reverse view」の投影画像のそれぞれを使って再構成した HRCT_[n,r] 画像の画質も悪い。
[c] これら３画像の画素値ヒストグラムの上のサンプル像のピークの位置や形状が⼀致しない。
これらの画像の再構成に適切なサンプル回転軸の位置が異なるので、161122[p,q] では測定中に投影画像上のサンプルの像が動いていたようです。しかしながら、投影値 sinograms に特に問題はない感じです。その確認⽤に以前の E-mail で紹介した投影値 sinograms などの図 21hi.pdf（対⽐⽤）と 22pq.pdf をこの E-mailに再度添付しました。測定ごとに以下の５個の図が載っています。
上段左normal view の sinogram、SG_n

上段右reverse view の sinogram を左右反転したもの、SG_r

上段中央

SG_[n,r] の重複部分の差分 sinogram

中段

SG_[n,r] を simple average（SA）で⾜し合わせたもの、SG_SA

下段

SG_SA から再構成した tomogram、TG_SA

22pq.pdf をよくよく⾒ると、161122[p,q] の TG_SA の上のサンプル像の左右の端に画像再構成に⽤いたサ
ンプル回転軸の位置が不適切だったことを⽰す偽像が発⽣していますね。
(2) ０°と 360°の投影画像の⽐較
161122[p,q] の 測 定 中 に サ ン プ ル 像 が 動 い て い た こ と を 確 認 す る た め 、 以 前 に 紹 介 し た プ ロ グ ラ ム
oct_cmp_0 を⽤いてサンプル回転⾓が０と 360°の投影画像のペアの位置決めを⾏いました。この E-mail に添付した 0.pdf を御覧下さい。
最上段に記した測定ごとのラベル
oct_cmp_0 で推定した０と 360°の画像ペアの位置ズレ (Ox, Oy)
上段の画像
その位置ズレの配置で重ね合わせた０と 360°の投影画像の差分画像
中段の画像
位置ズレがない場合のペアの投影画像の差分画像
下段の画像（オマケ）
サンプル回転⾓が０の投影画像。その左上隅のラベルは画像上の投影値の値域、また、外側下のものはその画像の画素の辺⻑（Drz）および横と縦の画素数（Nr と Nz）。
最下段に記したラベル（オマケ）
FRCT 画像の再構成の際に指定したオフセット値 (Ox, Oy)

76 / 93
0.pdf の差分画像からわかるように、oct_cmp_0 によるサンプル回転⾓が０と 360°の投影画像のペアの位
置ズレの推定値は正しいようです。そして、正常な測定 161121[h,i] の場合と⽐較すると、161122[p,q]では測定中に投影画像の上のサンプルの像が⽔平（x）⽅向に⼤きく動いていたことは確かです。
(3) 180°異なる投影画像のペアのオフセット値の時系列
161122[p,q]の測定中のサンプル像の動きを知るため、サンプル回転⾓が 180°異なる投影画像のペアのオフ
セット値などの「時系列」を調べました。ここでは以前に紹介したプログラム oct_cmp_180 で測定データか
ら新たに抽出した時系列を（これも以前に紹介した）C-shell script "plot.oct_180" で描画しました。この Email に添付した 180.pdf を御覧下さい。各測定の以下の４個の時系列のグラフを並べました。
１番上
回転⾓θの投影画像と「θ+180°の投影画像を左右反転した画像」の x ⽅向のオフセット値の推定
値 Ox(θ) の時系列のグラフ
２番⽬
上記の画像ペアの鉛直（y）⽅向のオフセット値の推定値 Oy(θ)
３番⽬
上記のオフセット値の配置にした画像ペアの重複した領域における２画像の同じ位置の画素の投影値の差の⾃乗の平均値（mean square difference)、MSD(θ)
４番⽬
上記のオフセット値の配置にした画像ペアの重複した領域における２画像の同じ位置の画素の投影
値の積の平均値（mean value of product)、MVP(θ)
これらのグラフの横軸はいずれも回転⾓θ（正確には、投影画像のペアのうち 180°未満の⾓度で撮影した画像の回転⾓の値；θ =０〜180°）に対応した投影番号（０〜）です。そして、Ox と Oy の値の単位は投
影画像の画素幅、MSD と MVP は投影値が無次元なので無次元値です。また、Ox(θ) と Oy(θ) のグラフ上部のラベルに記した「回転⾓ごとのマップを積算した MSD(Ox,Oy) が最⼩になっている Ox と Oy の値（HRCTの画像再構成で指定したオフセット値）」を⾚横線で⽰しました。なお、すべての測定の Oy、161121[h,i] のOx および 161122[p,q] の MVP のグラフの縦軸の値域が⾮常に狭いことに注意して下さい。つまり、Oy や
161121[h,i] の Ox の時系列は⼤変動しているように⾒えますが、実質的には⼀定値です。また、161122 [p,q]の MVP もほぼ⼀定の値で、これは撮影した forsterite が均質な円断⾯のサンプルであるためです。
180.pdf から 161122[p,q]の Ox(θ) がθの２次関数的な変動パターンを⽰していることがわかります。測定
中にサンプルがズルズルと動いたクリープなら１次関数なので、これはぼくの想定外の変動パターンです。
そこで、このような変動が先⽉の実験の他の FRCT（および Offset CT）の測定で発⽣していないかどうかを調べてみました。この E-mail に添付した xy.pdf を御覧下さい。
実は、oct_cmp_180 のものだけでなく先⽉の実験の場で使っていたプログラム oct_xy が出⼒した時系
列も描画できるように plot.oct_180 を改造しました。以下のようにして実験の場で oct_xy を使って作成したファイル xy.log（および xy.txt に書き込んだ値）を指定してやれば plot.oct_180 で 180.pdf のものと概ね同じグラフ（oct_xy の出⼒にない MVP 以外の３つの時系列のグラフ）を描くことができます。

77 / 93
# 念のために plot.oct_180 の起動法を表⽰
plot.oct_180→

usage : plot.oct_180 cmp_180.txt {OxOy}{fit_180.txt}# oct_xy で推定した Offset CT の画像再構成⽤の Ox と Oy の値を表⽰
cat→161119q/xy.txt
644-1

8.171594e-03# ３つの時系列のグラフの PostScript コードのファイル xy.ps を作成
plot.oct_180 161119q/xy.log

644 -1 > 161119q/xy.ps

xy.pdf から時系列 Ox(θ) もしくは Oy(θ) が⼀定値でないのは以下の５個の測定だけのようでした（その判定の「しきい値」をどうするかにもよりますが）。
161119s（撮影したサンプルは杉本さんの MIL090657 House5）
Offset CT のテスト⽤の投影画像の重複部分の幅を最も狭くした測定。この測定では Ox(θ) とOy(θ) がどちらも⼀定値ではないが、これはペアの投影画像の重複部分が狭過ぎるために⽣じた
推定ミスに起因した⾒かけの変動だと思われる。そして、幸運なことに、グラフ上に横⾚線で⽰した画像再構成に使⽤した Ox と Oy はこの⾒かけの変動の影響を受けていない値になっていた。
161122[p,q,r,s]（forsterite）
Offset CT や FRCT のテスト⽤に forsterite を撮影した測定すべてで Ox(θ) は２次関数的な時系列になっていた。
つまり、先⽉の実験の Offset CT や FRCT の測定のうちで Ox(θ) が「特異な」θの２次関数的なパターンになっていたのは forsterite を撮影したテスト⽤の４個の測定 161122[p,q,r,s] だけだったと⾔えそうです。
(4) オフセット値とサンプル回転軸の位置の関係式
以上のような Ox(θ) とサンプル回転軸の位置 RC(θ) の間の関係式を考えてみます。Ox(θ) は oct_xy で計算した回転⾓が 180°異なる投影画像のペアの重複部分の MSD(θ) が最⼩になっている場合のオフセット値ですが、これは通常の画像の位置決めの場合と同様に、ペアの２画像のそれぞれの上のサンプル像の「同⼀点」の x 座標値の差に⼀致します。この同⼀点としてサンプル回転軸の直上の点を考えます。oct_xy は「回転⾓が 180°未満の投影画像」と「それと 180°だけ異なる回転⾓の画像を左右反転した画像」のペアの重複部分の MSD を計算しているので、結局、Ox(θ) と RC(θ) の間の関係式は以下の通りです。
θの単位を radian、Nx を測定した投影画像の横画素数とすると、Ox(θ) = RC(θ) - {Nx-1 - RC(θ + π)} = RC(θ) + RC(θ + π) - (Nx-1)
数学的には、測定データから推定した Ox(θ)の時系列に上式を適⽤して RC(θ)を⼀意的に決定することはできません（例えば、θ≧180°で R(θ)≡ 0 の場合もあり得るから）。そこで、Ox(θ)を説明できる RC(θ)
について考えることにします。RC(θ) がθの２次関数なら Ox(θ) も（測定 161122[p,q,r,s] の推定値の特異な変動パターンを説明できる）２次関数になります。

78 / 93
RC(θ) = A × θ^2 + B × θ + C なら、
Ox(θ) = RC(θ) + RC(θ+π) - (Nx-1) = a × θ^2 + b × θ + c
ただし、
a=2×A
b = 2 × (A × π + B)
c = A × π^2 + B × π + 2 × C - (Nx-1)
この E-mail に添付した図 ox+rc.pdf を御覧下さい。測定 161122[p,q,r,s] のそれぞれに対する以下の２個のグラフが並んでいます。
グラフ "horizontal offset"
⿊の折れ線（IPR_MSD；"image pair registration using MSD" の略のつもり）
回転⾓が 180°異なる投影画像のペアの MSD の計算による位置決めで推定した Ox(θ) の時
系列データ。
⾚の曲線（by LSM；"by least square method"）
最⼩⾃乗法（LSM）で Ox(θ) の時系列のデータに当てはめた２次関数
緑の横線（FRCT；"full rotaion CT"）
FRCT 画像の再構成で指定した（θに関して⼀定の）オフセット値
⻘の縦線
Ox(θ) に当てはめた２次関数の頂点の位置"rotation center"
⾚の曲線
Ox(θ) の変動パターンを説明できる２次関数の RC(θ)。ただし、その係数値 A、B、C は Ox(θ)の時系列のデータに対して LSM で決定した係数値 a、b、c から前記の式で計算した。
緑の横線
FRCT 画像の再構成で指定したオフセット値を換算した回転軸の位置。以前の E-mail に記したようにこの値は RC_F = ( Nx - 1 + Ox ) / 2
シアンもしくはマジェンタの横線（RC_normal もしくは RC_reverse）
FRCT や Offset CT の normal もしくは reverse view の投影画像から HRCT 画像を再構成する
場合に最適な回転軸の位置。ただし、Offset CT の測定 161122[r,s] ではこれらは不完全な
HRCT 画像になる。
⻘の縦線
RC(θ) の２次関数の頂点の位置
これらのグラフの横軸は回転⾓θ（グラフの下端に度単位の値を記した）もしくはそれに対応する投影番号（上端の整数値）です。そして、グラフの左端には投影画像の画素幅単位の縦軸の値を、また、右端には⽔平⽅向のオフセット値やサンプル回転軸の位置のμm 単位の実⼨の値を記しておきました。
ox+rc.pdf のグラフからわかるように、161122[p,q,r,s,t] では測定中にサンプル回転軸が数μm の範囲を⾏ったり来たりしたようです。この２次関数的な変動は回転⾓θを時間に読み替えると等加速度運動なので、測定中ずっと、サンプル回転軸（およびそれが載っていたサンプルステージ全体）に投影画像の x ⽅向正の

79 / 93
向きの⼀定の⼒が働いていたことになります。このような⼒が働いていた可能性はあるでしょうか？ なお、
FZP-CT で撮影したものが X 線の上流側から⾒たサンプルの正⽴像なら、この⼒の向きはリングからホールの⽅向です。
それが⽣じた原因はさておき、以上のようにしてサンプル回転軸の位置の時系列 RC(θ)を得ることができました。それを⽤いて normal および reverse view の投影画像のつなぎ合わせを⾏えば、測定 161122[p,q,r,s]の⾼画質な FRCT や Offset CT 画像を再構成できるはずです。ただ、２次関数的変動パターンが⼀般的かど
うかわからないので、それ⽤の実⽤的なプログラムを書くかどうかは、...。
それよりも重要なのことは、Ox(θ) などの時系列を⾒れば FRCT や Offset CT の測定中に発⽣した「アクシデント」を検知できることです。今後の Offset CT や FRCT の測定では plot.oct_180 を⽤いて前記の⼿順で描いた時系列のグラフを毎回確認することにしましょう。
⻑い E-mail になりました。とりあえず以上です。

オマケ：０と 360°の投影画像の位置ズレ Ox0 と Ox(θ) の関係
前述のプログラム oct_cmp_0 で推定したθ=０と２×π（360°）の投影画像のペアの位置ズレ Ox0 は通常の画像ペアに対するものと同じ単純な２画像の位置決めの結果なので、RC(θ)との関係は以下のようになる。
Ox0 = RC(0) - RC(２×π)
そして、前記の Ox(θ)と RC(θ)の関係式から、Ox0 と「Ox(0)および Ox(π)」の関係式は以下のようになる。
０と 180°の画像のオフセット値：Ox1 = Ox(0) = RC(0) + RC(π) - (Nx-1)
180°と 360°の画像：Ox2 = Ox(π) = RC(π) + RC(２× π) - (Nx-1)
これらの差をとると、Ox0 = Ox1 - Ox2
実際の値を代⼊してみると、測定 161121[h,i] と 161122[p,q] ではこの関係式（Ox0 = Ox1 ‒ Ox2）がほぼ完全に成り⽴っていることがわかった。ただし、180°と 360°の投影画像のペアのオフセットの推定値 Ox2がないので、その代わりに 180°と 360°のそれぞれの１ステップ前の投影番号の画像ペアの値を使った。
Ox0

Ox1

Ox2

Ox1 - Ox2

161121h

0

82

82

0

161121i

0

83

83

0

161122p

-19

48

67

-19

161122q

-8

104

111

-7

80 / 93
添付ファイル 21hi.pdf と 22pq.pdf
これらはそれぞれ 12/7 の E-mail の添付ファイル 21hi.jpg と 22pq.jpg を御覧下さい。
添付ファイル 0.pdf

81 / 93
添付ファイル 180.pdf

82 / 93
添付ファイル xy.pdf（１ページ⽬）

83 / 93
添付ファイル xy.pdf（２ページ⽬）

84 / 93
添付ファイル ox+rc.pdf

Date:Wed, 21 Dec 2016 16:38:47 +0900

From:

Tsukasa NAKANO

To:

⼟`⼭ 明

Cc:

Kentaro UESUGI, Miyake, matsuno j, Aiko Nakato, Sugimoto Miyama, 北⼭晃, 中村隆太,
kadokawa tokiyuki, Takashi Sakurama, Minako HASHIGUCHI,
Masato Hoshino, Masayuki Uesugi, Akihisa TAKEUCHI, Akira TSUCHIYAMA

Subject: Re: Ox_and_RC

みやけさま、
つちやまさま、
GSJ/AIST のなかのです。E-mail ありがとうございます。昨⽇の E-mail の補⾜です。ぼくが⾒つけた測定161122[p,q,r,s] のサンプル回転軸の⽔平⽅向の「等加速度運動」についてですが、それと同様なサンプルの⾒かけの動きは測定中にサンプルが徐々に傾いた場合にも⽣じる可能性があることに気づきました。傾斜運
動と 360°回転の結果、投影画像上をサンプル像が⾏ったり来たりしたのかもしれない（ただし、平⾏移動を仮定した MSD 計算による画像ペアの位置決め処理でその「回転+傾斜運動」がどのように検知されるかは不明ですが）。

85 / 93
回転⾓が０と 360°の（投影画像と等価な）X 線透過率画像のペアの観察から、測定 161122p ではそのようなサンプルの傾斜運動が起こっていないことを確認しました。この E-mail に添付した 22p_0+360.pdf を御覧下さい。上段の０°の X 線透過率画像の上の「サンプルの台（真鍮？）の左右の端を⽰す四辺形」を下段の 360°の X 線透過率画像の上の同じ位置にコピーしました。これより、360°の画像の上のサンプル像は０°の画像のものから傾いておらず、それが単純に⽔平⽅向に平⾏移動しただけであることがわかります。
> > それが⽣じた原因はさておき、以上のようにしてサンプル回転軸の位置の時系列> > RC(θ) を得ることができました。それを⽤いて normal および reverse view
> > の投影画像のつなぎ合わせを⾏えば、測定 161122[p,q,r,s] の⾼画質な FRCT
> > や Offset CT 画像を再構成できるはずです。ただ、２次関数的変動パターンが> > ⼀般的かどうかわからないので、それ⽤の実⽤的なプログラムを書くかどうかは、…。
>> 再構成画像、格段に良くなるのでしょうね。
> 修正したもの、是⾮⾒てみたい気はしますが・・・。
161122[p,q,r,s,t] の forsterite の円断⾯の像はそもそもが⾒栄えしないので、...。
とり急ぎ、
On Wed, 21 Dec 2016 09:08:24 +0900 ⼟`⼭ 明 wrote:> 中野さん>> 詳細な解析ありがとうございます。
>> > ox+rc.pdf のグラフからわかるように、測定 161122[p,q,r,s,t] では測定中に> > サンプル回転軸が数 μm の範囲を⾏ったり来たりしたようです。この２次関数> > 的な変動は回転⾓ θ を時間に読み替えると等加速度運動なので、測定中ずっと、> > サンプル回転軸（およびそれが載っていたサンプルステージ全体）に投影画像の> > x ⽅向正の向きの⼀定の⼒が働いていたことになります。このような⼒が働いて> > いた可能性はあるでしょうか？ なお、FZP-CT で撮影したものが X 線の上流側> > から⾒たサンプルの正⽴像なら、この⼒の向きはリングからホールの⽅向です。
>> 三宅さんからん回答> > 161122[p,q,r,s]

なんか通常と違うことしたかなぁ。

> > 最終⽇の朝で、フォルステライトを取るはずが AV ⾙を間違って撮っていて、> > 慌てて、フォルステライトにしたのです。
> > なので、その時の⼿順ミスや試料のグリップの⽢さとかあったのか・・・。
> > ちょっと、今となっては思い出せません。
>> サンプルステージ全体（？）の⽔平⽅向のドリフトというのは、単なる⼿順ミス> のようなものではないはずですよね。サンプルステージが押されていたような86 / 93
> ことも考えにくいし。
> サンプルホルダーの取り付け⽅の問題というのもないですよね（？）
。
>> > それが⽣じた原因はさておき、以上のようにしてサンプル回転軸の位置の時系列> > RC(θ) を得ることができました。それを⽤いて normal および reverse view
> > の投影画像のつなぎ合わせを⾏えば、測定 161122[p,q,r,s] の⾼画質な FRCT
> > や Offset CT 画像を再構成できるはずです。ただ、２次関数的変動パターンが> > ⼀般的かどうかわからないので、それ⽤の実⽤的なプログラムを書くかどうかは、…。
>> 再構成画像、格段に良くなるのでしょうね。
> 修正したもの、是⾮⾒てみたい気はしますが・・・。
>> > それよりも重要なのことは、Ox(θ) などの時系列を⾒れば FRCT や Offset CT> > の測定中に発⽣した「アクシデント」を検知できることです。今後の Offset CT> > や FRCT の測定では前記の⼿順で plot.oct_180 によって描いた時系列のグラフ> > を毎回確認することにしましょう。
>> そうしましょう。
添付ファイル 22p_0+360.pdf

87 / 93
Date:Thu, 22 Dec 2016 15:08:26 +0900

From:

Tsukasa NAKANO

To:

Kentaro UESUGI

Cc:

tsuchiyama.akira.4x, miya, jmatsuno, nakato, sugimoto.miyama.63x, kitayama.akira.46a,ryuta0609, tkadokawa, 櫻間, hashiguchi.minako.123, hoshino, uesugi, take, atsuchi

Subject: Re: Ox_and_RC

うえすぎさま、
GSJ/AIST のなかのです。E-mail ありがとうございます。
(1)
何度もすみません。⼀昨⽇の E-mail（Subject：Ox_and_RC）の補⾜です。
On Tue, 20 Dec 2016 15:29:54 +0900 Tsukasa NAKANO wrote:> xy.pdf から時系列 Ox(θ) もしくは Oy(θ) が⼀定値でないのは以下の５個の> 測定だけのようでした（その判定の「しきい値」をどうするかにもよりますが）。
先⽉の FZP-CT 実験の時にプログラム oct_xy を使って得た FRCT や Offset CT の Ox(θ) などの時系列の図（xy.pdf）を再度吟味したところ、先の E-mail に記した 161119s と 161122[p,q,r,s] だけでなく、イトカワ
粒⼦ RB-QD04-0040 を７および８keV で撮影した FRCT の測定 161121[m,n] の Ox(θ) も⼀定値とは⾒なせない変動パターンになっていました（ただし、こちらは⼀次関数的な変動パターンです）。そこで、それらの sinograms を調べてみました。この E-mail に添付した 21mn.pdf を⾒る限り、測定 161121[m,n] のsinograms には特段の問題は無いように思えます。
(2)
このように、ぼくが書いたプログラム oct_xy および oct_cmp_180 で推定したオフセット値の時系列は完全には信頼できません。そして、上杉くんのコメント「測定中に回転ステージが勝⼿に動くことはあり得ない」
や三宅さんのコメント「サンプル設定の際の不⼿際は今となっては思い出せない」を考え合わせると、
161122[p,q,r,s] の「等加速度運動」は再現できない現象のように思いました。それゆえ、現時点ではこの現象の原因や対応策を考えないことにします。つまり、[1] FRCT や Offset CT では測定中にオフセット値が変化しないと仮定する。
[2] オフセット値の変動が⼤きい測定は「失敗」だと⾒なして取り扱わない。
[3] FRCT や Offset CT ⽤の画像再構成プログラムもオフセット値が⼀定値であることを仮定し、その変動に対応したプログラムは作成しない。
という訳で、先⽉の FRCT や Offset CT のテスト⽤の測定 161122[p,q,r,s] は失敗なので、今後それらに関することは議論しないつもりです。
とり急ぎ、

88 / 93
On Thu, 22 Dec 2016 09:26:15 +0900 Kentaro UESUGI wrote:> 上杉です>> 忙しくてちゃんと読んでいませんが、> ステージは触らなければスキャン前後、あるいは複数回のスキャン中には動きません。
> これは１１⽉末の⽵内さんの実験の時に確認しています。
> （というか、⾃分のスタディの時はいつも）>> 試料交換による取付時の操作と試料そのものの問題のような気がします。
> とりあえず、試料を取り付け時にどのように動くか測⻑器でも付けてみますか。
添付ファイル 21mn.pdf

89 / 93
Date:

Tue, 15 Aug 2017 19:21:20 +0900

From:

Tsukasa NAKANO

To:

Kentaro UESUGICc:"TAKEUCHI, Akihisa", "TSUCHIYAMA, Akira", Miyake, MATSUNO Junya

Subject: bug_of_oct_*_sbs

うえすぎさま、なかのです。
「180 度異なる回転⾓の投影画像を side by side（SBS）でつなぎ合わせる Offset CT ⽤のプログラム」に間違いを⾒つけました。具体的には下記の４種類・６個のプログラムです。
ソースコード実⾏ファイル

oct_xt.coct_xt_sbs と oct_xp_sbs

oct_sg.coct_sg_sbs

oct_tg.coct_tg_sbs

oct_sg2tg.c

oct_sg2tg_t と oct_sg2tg_g

修正済みのコードとそれらをコンパイルした Windows ⽤実⾏ファイル（32 bit OS ⽤の "*_32.exe" と 64 bit
OS ⽤の "*_64.exe"）を Offset CT ⽤の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zipに⼊れておきました。
今回のバグは SBS で処理した Offset CT 画像が[1] 横⽅向のオフセット値 Ox が負の値
[2] ２枚の投影画像が重なっている部分の横画素数 W が奇数のいずれかの場合に顕在化します。そして、[2] の画像に⽣じる乱れはごく僅かですが、[1] の画像の乱れはかなりひどいです。これらは昨年の 11 ⽉にテスト⽤に撮影した下記の２測定の Offset CT with SBS 画像で発⽣していました：
161119q（Ox = 644：正、W = 2115 - 644 * 2 = 827：奇数）[2] の画像。どこが乱れているかわからない。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.pdf#page=32
161122r（Ox = -180：負、 W = 639 - 180 * 2 = 279：奇数）[1] と [2] の画像。⾔われてみると画像のつなぎ⽬などが変。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.pdf#page=54
とりあえず以上です。

90 / 93
Date:

Fri, 01 Sep 2017 11:31:14 +0900

From:

Tsukasa NAKANO

To:

Kentaro Uesugi

Cc:

MATSUNO Junya, Miyake, "TSUCHIYAMA, Akira", Masayuki Uesugi, ⽵内, 星野

Subject: bug_of_rhp.c+oct_*うえすぎさま、
GSJ/AIST のなかのです。
(1)
SPring-8 の通常の CT 実験で得た測定データから X 線透過率画像を取り出す C ⾔語コード rhp.c に間違いを⾒つけました。⼊射 X 線強度（I0）の測定画像が１個だけの場合に「コアダンプ」が発⽣します。以下の書庫ファイルに修正したコードやそれを組み込んだ Windows ⽤の実⾏ファイルを⼊れておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rhp.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rhp.zip
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.zip(2)
8/15 の E-mail で報告した SBS（side by side）で⾏うものに加えて「180 度異なる回転⾓の投影画像を linear
mixing（LM）でつなぎ合わせる Offset CT ⽤プログラム」にも間違いを⾒つけました。横⽅向のオフセット
値 Ox が負の場合の LM による画像のつなぎ合わせ⽅が間違っていました。修正したコードや Windows ⽤の実⾏ファイルを前記の書庫ファイル oct.taz と oct.zip に⼊れておきました。また、今回の LM と 8/15 の SBSのバグの修正に伴う Offset CT の処理結果の違いを⽰す以下の２枚の画像をこの E-mail に添付しました。
161207_22r.jpg
2016/12/07 の E-mail に添付した「バグあり」の画像
170830_22r.jpg
修正したプログラムで得た画像。SBS の sinogram がまだ変？
上記のものからレイアウトを変えた画像を後に貼り付けました。
(3)
SPring-8 の CT 実験などで得た「サンプルを 180 度回転（half rotation、HR）した測定データ」から Offset
CT ⽤の「サンプルを 360 度回転（full rotation、FR）した測定データ」を合成するプログラム hp_h2f を書きました。このデータを Offset CT ⽤のプログラム oct_tg などで処理すれば、HR-CT の測定データから得られるものと概ね同じ再構成画像が得られるはずです（これらの CT 画像の画質には画像再構成の際の数値計算の誤差に起因する若⼲の差異が⽣じます）。

91 / 93
プログラム hp_h2f のソースコードや Windows ⽤の実⾏ファイルを前記の書庫ファイル oct.taz と oct.zip に⼊れておきました。hp_h2f の起動法や起動時に指定すべきパラメータの意味は以下の通りです。
起動法
hp_h2f HRCT/ RC

Ox

Oy {XT_min XT_max} FRCT/ >H2F.txt

起動パラメータ
HRCT/
HR-CT の測定データファイル（output.log、dark.img と"q*.img"）が⼊っているディレクトリ
の名前
RC
HR-CT のサンプル回転軸の位置（いわゆる「センター値」）。HR-CT の測定画像の横画素数よりも⼩さい正の値でなければならない。
Ox と Oy
FR- or Offset CT の測定データの合成に⽤いる x（⽔平）と y（鉛直）⽅向のオフセットの値。
これらはどちらも整数値で、Ox はその絶対値 | Ox | が HR-CT の測定画像の横画素数の半分
より⼩さい値でなければならない。これらのオフセット値を Offset CT ⽤のプログラム "oct_*"に指定してやれば、合成した測定データを使った画像再構成などの処理を実⾏できる。
XT_min と XT_max
HR-CT の測定データから合成した Offset CT ⽤の「浮動⼩数点数画素値の X 線透過率画像」
を「整数画素値の X 線強度画像」に変換する際に⽤いる「X 線透過率の最⼩値と最⼤値」。合成した X 線強度の測定画像から計算が可能な X 線透過率はここで指定した範囲の値になる。
これらの値の指定を省略するとそれぞれを０と１と⾒なして処理を⾏う。後述するように
hp_h2f は処理に⽤いた浮動少数点数の X 線透過率の値域を出⼒するので、必要ならそれらの値を指定して hp_h2f を再実⾏すれば良い。
FRCT/合成した Offset CT ⽤の測定データファイル（output.log、dark.img と"q*.img"）を格納する
ディレクトリの名前。ただし、このディレクトリは事前に⾃分で作成しておく必要がある。
H2F.txt
hp_h2f が標準出⼒に書き出す以下の２⾏・５個の値のテキストデータをリダイレクトするためのファイルの名前。
１⾏⽬（タブコード区切りで３個の値が並んでいる）[1] 合成した測定画像の横画素数 = HR-CT の測定画像の横画素数 - | Ox |[2] 合成した測定画像の横画素数 = HR-CT の測定画像の縦画素数 + | Oy |[3] Offset CT の投影数 = HR-CT の投影数の２倍の値
２⾏⽬（２個の値）[1,2] 合成処理に⽤いた X 線透過率の最⼩値と最⼤値
上記のパラメータ Oy の値域には制限はありませんが、それを使った hp_h2f の y ⽅向の「オフセット処理」
はいい加減です。つまり、hp_h2f はサンプル回転⾓に応じて X 線強度画像（サンプル像）を y ⽅向に徐々に移動するのではなく、回転⾓が 180 度異なっている normal view と reverse view の測定画像の⼀⽅を他⽅に

92 / 93
対して y ⽅向に⼀定値 Oy だけ異なる配置にします。これは⾒た⽬は変ですが、現状の Offset CT ⽤のプログラム群で問題なく処理できます。
と⾔うか、現状の Offset CT ⽤のプログラム群は回転中に Oy の値が徐々に変わる場合を想定していま
せん：http://www-bl20.spring8.or.jp/~sp8ct/tmp/oct.pdf#page=87(4)
既存の HR-CT の測定データから hp_h2f に正・負および奇・偶数のオフセット値 Ox と Oy を指定して合成した FR-CT の合計 16 セットの測定データを⽤いて、下記の Offset CT ⽤の合計 18 個のプログラムの動作
テストを⾏いました。
oct_xt、oct_xt_sa、oct_xt_sbs、
oct_xp、oct_xp_sa と oct_xp_sbs
与えられた normal view と reverse view のものを LM、SA（simple average）と SBS のそれぞれでつなぎ合わせた X 線透過率（XT）もしくは X 線投影値（XP）の画像を作成する。
oct_sg、oct_sg_sa と oct_sg_sbs
LM、SA と SBS のそれぞれでつなぎ合わせた X 線投影値の sinogram の画像を作成する。
oct_tg、oct_tg_sa、oct_tg_sbs、
oct_sg2tg_t、oct_sg2tg_sa_t、oct_sg2tg_sbs_t、
oct_sg2tg_g、oct_sg2tg_sa_g と oct_sg2tg_sbs_g
LM、SA と SBS のそれぞれで画像のつなぎ合わせを⾏なった後に Offset CT 画像を再構成する。
これらで作成した画像に LM、SA と SBS の違いに伴う差異はありませんでした。つまり、
hp_h2f による HR-CT のものから FR-CT の測定データの合成の⽅法と上記の Offset CT ⽤プログラム "oct_*" のそれぞれに組み込んだ normal view と reverse view の画像のつなぎ合わせの仕⽅には間違いはないようです（すべてを系統的に間違っていなければ）。
とりあえず以上です。

93 / 93
LM と SBS のバグの修正に伴う測定 161122r の sinograms や tomograms の画質の変化：