1 / 124

SIXM_mails
date
2015/
4/28

subject

DPC-CT用再構成フィルタ

6/17
6/18
6/18
6/19
6/22

SIXM用の画像再構成プログラム
Re: SIXM用の画像再構成プログラムの訂正
Re: SIXM用の画像再構成プログラム
SIXM用の画像再構成プログラムの改造
Re: SIXM用の画像再構成プログラムの改造← Intel C++ Compiler for Windowsによるプログラムのコンパイルに関する話

6/23
6/24

Re: paper
位相 != 屈折

6/25
7/2
7/7
7/7

SIXM用画像再構成プログラムに関する話
SIXM用画像再構成プログラムに関する話の続き
SIXM用画像再構成プログラムに関する話の続きの続き
Re: SIXM用画像再構成プログラムに関する話の続きの続き

6/26
6/27
7/9
7/13

180度の方向から見た微分位相の投影画像
SIXM用センター値推定スクリプト
SIXM用センター値推定プログラム
センター値が決まらない場合

7/16
7/16
7/17
7/16
7/16
7/16
7/17

SP8_BL47XU_SIXM_150603c_a.HIS
SP8_BL47XU_SIXM_150603c_a.HIS_訂正
Re: SP8_BL47XU_SIXM_150603c_a.HIS_訂正
SP8_BL47XU_SIXM_幾何学量
Re: SP8_BL47XU_SIXM_幾何学量
Re: SP8_BL47XU_SIXM_幾何学量
SP8_BL47XU_SIXM_150602

8/1

SIXM_150602_オマケ

8/10
8/18
8/20

SIXM_empty_HIS
Re: SIXM_empty_HIS
Re: SIXM_empty_HIS

2 / 124

8/27

Re: SIXM_150602_オマケ

9/25
9/25
9/30
9/30
9/30

RIDの計算プログラム
RID ~= density
Re: RID ~= density
Re: RID ~= density
Re: RID ~= density

10/1
10/7
10/8
10/9
10/14
10/14
10/14

sum_of_projection(SP)の計算プログラム
SPの計算プログラムの改造
Re: SPの計算プログラムの改造
filtered_projection（FP）の計算プログラム
Re: filtered_projection（FP）の計算プログラム
Re: filtered_projection（FP）の計算プログラム
Re: filtered_projection（FP）の計算プログラム

10/14

xp2tiff

10/16
10/16
10/16

位相CTシミュレータ
Re: 位相CTシミュレータ
Re: 位相CTシミュレータ

10/19
10/20
10/20
10/20
10/20

SIXMの屈折プロファイルの重心
Re: SIXMの屈折プロファイルの重心
Re: SIXMの屈折プロファイルの重心
Re: SIXMの屈折プロファイルの重心
Re: SIXMの屈折プロファイルの重心

10/21
11/23
12/5
12/28

Martiono論文
SIXMの位相回復プログラム
SIXMの位相回復プログラムの別バージョン
Martino論文の問題点

2016/
1/4

Re: Martino論文の問題点

3 / 124

Date:
Tue, 28 Apr 2015 18:11:31 +0900
From: Tsukasa NAKANO
To:
Kentaro Uesugi
Subject: DPC-CT 用再構成フィルタ
うえすぎさま、なかのです。
(1) 以前に教えてくれた Pfeiffer et al. の論文（2007、PRL 98）に記されていた DPC（Differential Phase
Contrast）CT 用の再構成フィルタの計算機コードを書いてみました。Nyquist 周波数の制約を考慮した「窓関数」を変えた３種類のものを用意しました（再構成フィルタの欄のものは計算機コードの関数名）：
窓関数
box 関数sinc 関数
Hanning 窓

DPC-CT 用
DPC_R()
DPC_S()
DPC_C()

既存の吸収 CT 用再構成フィルタ
Ramachandran()
Shepp()
Chesler()(2) 関数 DPC_[R,S,C]() のコード dpc_[r,s,c].c とそれらのテスト用プログラムのコード "*_test.c"（および、gcc によるそれらのコンパイル用の Makefile）が以下の書庫ファイルに入っています。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/dpc.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/dpc.zip
また、これらには３種類の窓関数を用いた吸収 CT 用の再構成フィルタ gj と DPC_CT 用のフィルタ hjの計算式（単純な積分）に関するぼくの手書きノートをスキャンした PDF ファイル dpc.pdf も入っています
ので、どうぞ御覧下さい。
(3) 関数 DPC_[R,S,C]()は CBP engine（CUDA 版を含む）から呼び出して使います。そのコンパイル法は以下の通りです（CBP engine に DPC_R() を組み込む場合）。
gcc

...

-DFilter=DPC_R...

dpc_r.c

cbp.c

...

ただし、dpc_s.c だけは ANSI-C 非標準（C90 標準）の関数 sincos()を使っているので、gcc でそれをコンパイルするにはオプション指定 "-D_GNU_SOURCE" を追加する必要があります（icc ではこの指定は
不要）：
gcc

...

-DFilter=DPC_S-D_GNU_SOURCE...

dpc_s.c

cbp.c

...

(4) このようにしてコンパイルした（従来から使っている）プログラムに検出器の間隔δ（CBP engine のパラメ
ータ dr）として値１を指定すれば DCT-CT の画像再構成を正しく行うことができます。ちょっと長くなりますが、その理由を以下に書いておきます。

4 / 124(5) CBP engine と再構成フィルタの計算用の関数は検出器間隔δの値を１とした値をやり取りします（そうす
るよう CBP engine を設計しました）。そのため、吸収 CT 用の Ramachandran() などの関数は「理論値」に含まれている 1/δ^2 のファクタ（前記の dpc.pdf を御覧下さい）を１とした再構成フィルタの値を返すよ
うになっています。そして、CBP engine では投影値 p(r) と再構成フィルタ g(r) の畳み込み積分
q(r) = ∫ g(r - ρ) ・ p(ρ) dρ
を和分近似して filtered projection の値 q(r) を計算するので、結局、その等間隔の点δ・j （ただし、j = 0、±１、...）における値 Qj は下式で計算されることになります（ここで、Pk は等間隔に並んだ検出器 k における投影値）。
Qj= δ ・Σ_k { (1/δ^2)・Ramachandran(j - k)・Pk }= (1/δ) ・Σ_k { Ramachandran(j - k)・Pk }[1](6) 一方、DCT-CT 用の関数 DPC_R()なども再構成フィルタの理論値に含まれているファクタ 1/δ（dpc.pdf
を御覧下さい）を１とした値を返すので、本当はQj= δ・Σ_k { (1/δ)・DPC_R(j - k)・Pk }= Σ_k { DPC_R(j - k)・Pk }[2]とすべきですが、吸収 CT 用の畳み込み演算を想定した CBP engine を流用すると式 [1] にあるように
Qj に余分な 1/δ のファクタが掛かってしまいます。これを防ぐ（つまり、式 [1] と [2] を同一の内容にする）には陽にδ = 1 を指定してやれば良いわけです。
ぼくが書いた画像再構成プログラムではδの値を CBP engine のパラメータとして使うだけでなく、それを再構成画像のファイルに「コメント（TIFF の image description tag）」として埋め込んでいます。
それゆえ、既存のプログラム（のメインルーチン）を改造せずに DPC-CT に使うのはちょっと具合が悪い。これに対処するには CBP engine に指定するδ（dr） の値を１.0 にすれば OK です。
(7) 前記の書庫ファイルに入っている dpc_[r,s,c].c の C 言語関数のテストプログラムについても簡単に説明
しておきます。
プログラム ci_test（起動パラメータなし）
dpc_s.c の DPC_S() で使っている余弦積分関数（Cosine integral）の値を計算するための関数
Ci(double t) のテストプログラム。標準入力から読み込んだ０以外の値に対する余弦積分関数の近似値を標準出力に書き出す。なお、dpc_s.c にも記したように、Wikipedia に載っていた余弦積分
関数の近似式を用いた。
dpc_test {B S N}
３（窓関数の種類）×２（吸収もしくは DPC-CT 用）の再構成フィルタが返す値そのものを表示する
プログラム。整数値 j = B、B+S、...、B+S*(N-1) のそれぞれに対する以下の７個の値をタブコード
区切りの行として表示する。

5 / 124[1] 整数値 j[2,3] DPC_R(j) と - DPC_R(j) / Ramachandran(j)[4,5] DPC_S(j) と - DPC_S(j) / Shepp(j)[6,7] DPC_C(j) と -DPC_C(j) / Chesler(j)ただし、分母が０で割り算が不可能だった場合には記号"+"を表示する。なお、起動パラメータとして B、S、N の値を指定しなかった場合は B=0、S=1、N=16 とする（つまり、j = 0～15 に対する値を計算する）。
前記の dpc.pdf にも書いてあるように、そこで導出している吸収 CT 用再構成フィルタの値と DPC 用の
それとの間には以下の関係式が成り立っています。
窓関数が box 関数の再構成フィルタ
DPC_R(j) = - j ・Ramachandran(j)窓関数が Hanning 窓のフィルタ
j が奇数なら、DPC_C(j) = - j・Chesler(j)このような関係式が DPC_S(j) と Shepp(j) や偶数の j に対する Hanning 窓のフィルタの間で成り立っているか否かを調べるために dpc_test を書きました。j が大きな偶数値なら DPC_C(j) = - j・Chesler(j) が成り立っているようです。しかし、DPC_S(j) と Shepp(j) の間の関係式はこれらとは異なります。理由は不明ですが、j が大きな値の場合には以下の関係式が成り立っているようです（計算式導出の間違いではないことは何度もチェックしました）。
sinc 関数を窓関数とした再構成フィルタ
j が大きければ、DPC_S(j) = (-1)^j ・(2 / π) ・j・Shepp(j)とりあえず、以上です。

6 / 124

Date:
From:
To:
Cc:

Wed, 17 Jun 2015 17:21:16 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM 用の画像再構成プログラム
みなさま、
GSJ/AIST のなかのです。ぼくが書いた SIXM の吸収および屈折トモグラフィー用の画像再構成プログラム
について簡単に説明します。ただし、出来立てなのでバグがあるかもしれません。
画像再構成に必要最低限のプログラムしか用意していません。
Windows 上でも走りますが、大容量メモリ（～64 bit OS）が必要です。
(1) 書庫ファイル
以下の２個の書庫ファイル（同内容です）にぼくが書いたプログラムの C 言語のソースファイルなどと 64 bit
Windows 用の実行ファイルが入っています。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(2) 個々のプログラムの概要
プログラムはいずれも UNIX や Windows の端末環境で起動します。
check_his HIS
機能
HIS 形式ファイルをスキャンしてそれに含まれる複数の画像それぞれの横および縦画素数と「画素値のタイプ」を表す自然数を表示する。
使用例
check_his 150602d/a.HIS
→
128
400
2
← 横および縦画素数と画素値のタイプ（2 = 16 ビット画素値）…
← 各画像（今の場合は合計 180951 画像）の３つの値の出力
18540966032← 150602d/a.HIS の総バイト数（標準エラー出力に書き出される）。
his2tiff HIS rangeList format
機能
HIS 形式ファイルに含まれている画像のうちで文字列 rangeList によって指定した０から始まる番号のもの（複数も可）だけを抜き出し、文字列 format で名前の仕様を指定した TIFF 形式画像ファ
イルに書き込む。

7 / 124

使用例
mkdir dark
his2tiff 150602d/a.HIS 0-99 dark/%02d.tif← 100 枚の暗電流画像をディレクトリ dark/ の下に格納する。
his2xp HIS darks I0s scans views aXP rXP {aXP_180 rXP_180}
his2xp_top HIS darks I0s scans views aXP rXP {aXP_180 rXP_180}
機能
HIS 形式のファイルから吸収および屈折トモグラフィーの画像再構成に使うデータを抽出し、バイ
ナリデータファイル aXP および rXP に書き出す。また、指定があれば、サンプルを 180 度回転した
場合の投影値などのデータも抜き出してファイル aXP_180 や rXP_180 に書き出す。ファイル名以外のこれらの起動パラメータの意味は以下の通り。
darks
ファイル HIS に含まれている暗電流画像の枚数。
I0s
入射光強度 I0 として用いるサンプル回転角ごとの「投影画像」の左および右端からの列の数（スキャンの数；his2xp の場合）、もしくは I0 とするサンプル上部のスライス数（his2xp_top）。
scans
サンプル移動に伴う「空撮影」したスキャンの回数（１と仮定）を含むサンプル回転のステ
ップごとの SIXM のスキャン数。
ｖiews
180 度回転を含むサンプル回転のステップ数。
なお、これらのパラメータのうちの scans、views と darks に対しては SIXM の測定のログファイル
a.log の１行目にコンマ区切りで並んでいる値をそのまま指定すれば良い。
使用例
後回しにします。
xp2xt XP.bin {base step} XT.gif
機能
バイナリデータファイル XP.bin に入っている his2xp などで得た投影値から X 線透過率を計算し、GIF 形式の動画のファイル XT.gif に書き込む。base と step は動画の８ビット画素値に変換する際に用いる X 線透過率の最小値と画素値１諧調あたりの増分の値で、これらの指定を省略するとbase = 0、step = 1/255 と見なされる。
使用例
省略します。
xp2sg XP.bin z SG.bin
機能
his2xp などで得た投影値のバイナリデータファイル XP.bin から通常の１スライスごとの画像再構
成用プログラムで処理可能なスライス番号 z の sinogram のバイナリデータファイル SG.bin を抽出する。
使用例

8 / 124

省略します。
a2tiff array {base step} BPS TIFF
機能
his2xp などで得たバイナリデータに含まれる複数の「投影画像」の最初の１枚や、xp2sg で得た
sinogram のバイナリデータを TIFF 形式の画像ファイルに変換する。起動パラメータ array はこれらのバイナリデータのファイルの名前、BPS は TIFF 画像の画素値のビット数である。また、base とstep は TIFF 画像の上の画素値に変換する際の投影値などの最少値と画素値１諧調あたりの増分で、これらの指定を省略するとそこに出現した値を最少値と最大値で正規化した画素値に変換する。
使用例
省略します。
xp2tg aXP Dr Or Θ0 {base step} BPS TG/ > LOG.txt
xp2tg_dpc rXP Dr Or Θ0 {base step} BPS TG/ > LOG.txt
機能
his2xp などで得たバイナリデータから３次元画像のスライス画像すべてを再構成する。吸収トモグ
ラフィー用の xp2tg はファイル aXP のデータを読み込み、Chesler タイプの再構成フィルタを用いて画像再構成を行う。また、xp2tg_dpc はファイル rXP の微分位相のデータを Pfeiffer の論文（PRL.,
2007）に載っていた再構成フィルタを用いて処理する（微分位相の線積分を画像再構成と同時に行う）。ただし、これらの処理には以下の幾何学パラメータを使用する。
Dr
再構成画像のスライス面内の正方形画素の辺長Or
いわゆる「センター値」（単位は画素の辺長）
Θ0
サンプル回転角の初期値（単位は度）
その後、再構成した CT 値を BPS で指定したビット数の画素値に変換する。係数値 base と step を指定した場合には前者を画素値０に対応する CT 値、後者を画素値１諧調あたりの CT 値の増分とする。また、それらの値を指定しなければ、スライスごとに出現した CT 値すべてがちょうどおさまる base と step の値が自動的に選ばれる。さらに、BPS として負の値を指定するとその絶対値が画素値のビット数となり、かつ、３次元再構成画像上の CT 値すべてがちょうどおさまる base と stepの値がすべてのスライス画像に対して用いられる（これは新機能）。
最終的には、CT 値から変換した画素値を書き込んだスライスごとの TIFF 形式画像ファイルをディレクトリ TG の下に書き込む。また、それと同時に、それぞれのスライスに出現した CT 値の最小値と最大値を標準出力に書き出す（ただし、上の起動例ではそのテキストデータをファイル
LOG.txt にリダイレクトしている）。
使用例
後回しにします。
(3) SIXM の測定 150602d で撮影した画像の再構成
画像再構成処理に必要最小限の端末入力は以下の通りです。
a.log の１行目のチェック
head -1 a.log

9 / 124→
401,451,100,1,0
a.HIS から投影値などのバイナリデータを抽出
his2xp_top a.HIS 100 30 401 451 a.bin

r.bin

画像再構成プログラムが使うスレッド数（４とする）の指定
set THREADS=4← Windows のコマンドプロンプト
setenv THREADS 4← csh もしくは tcsh
export THREADS=4← bash吸収トモグラフィーの画像再構成
mkdir axp2tg a.bin 150e-7 200.5

0

-8

a

> a.log

屈折トモグラフィーの画像再構成
mkdir rxp2tg_dpc r.bin 150e-7
200.5

0

-8

r > r.bin

ただし、画像再構成プログラムには辺長として 150e-7 cm（= 150 nm）、センター値 200.5 を指定しました。
とりあえず以上です。
Date:
From:
To
Cc:

Thu, 18 Jun 2015 11:23:55 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM 用の画像再構成プログラムの訂正
みなさま、なかのです。すみません。昨日の E-mail の間違いを思い出しました。
On Wed, 17 Jun 2015 17:21:16 +0900 Tsukasa NAKANO wrote:> (3) 測定 150602d の SIXM 測定で得た画像の再構成>> 画像再構成処理に必要最小限の端末入力は以下の通りです。
> ...
>屈折トモグラフィーの画像再構成>
mkdir r>xp2tg_dpc r.bin 150e-7 200.5 0 -8 r > r.bin>> ただし、画像再構成プログラムには画素の辺長として 150e-7 cm（= 150 nm）、

10 / 124

> センター値 200.5 を指定しました。
上記の xp2tg_dpc の起動法の例に間違いが２つありました。ひとつ目はログのファイル名に関するケアレス
ミスで、"r.bin" ではなく "r.log" として下さい。
２番目の間違いは起動時に指定する再構成画像の正方形画素の辺長 Dr の値です。4/28 の E-mail で上杉くんに説明しましたが、微分位相の線積分を同時に行う再構成フィルタを組み込んだ吸収 CT 用プログラム
コードを流用している関係上、Dr には実際の辺長（150e-7）ではなく常に１を指定して下さい。
とり急ぎ、
Date:
From:
To:
Cc:

Thu, 18 Jun 2015 10:39:55 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: Re: SIXM 用の画像再構成プログラム
みなさま、なかのです。昨日の E-mail の補足です。
On Wed, 17 Jun 2015 17:21:16 +0900 Tsukasa NAKANO wrote:> GSJ/AIST のなかのです。ぼくが書いた SIXM の吸収および屈折トモグラフィー> 用の画像再構成プログラムについて簡単に説明します。
(1) 書庫ファイル his.taz と his.zip に入っているコードのコンパイル法はそれを展開したディレクトリで"make" と端末入力するだけですが、その際に[1] GNU C-compiler（gcc）でコンパイル[2] POSIX thread（pthread）ライブラリを使用する設定にしてあります。[1] の変更には"Makefile"の修正が必要です。また、Windows では [2] のために pthread-win32 などをインストールしておく必要があるかもしれません（ぼくが使っている MSYS2の環境では Windows 用 pthread ライブラリが標準でインストールされていました）。
(2) プログラム xp2tg と xp2tg_dpc は上記の pthread ライブラリを用いたマルチスレッド処理で画像再構成を行います。昨日の E-mail にも書いたように、その際に使うスレッドの個数を環境変数 THREADS で指定できます。使用する CPU のコア数に応じたスレッド数を適宜設定して下さい。
(3) プログラム xp2tg と xp2tg_dpc は３次元再構成画像上の CT 値のデータすべてをメモリに保持しつつ処理を行います。メモリの節約を目してこの値を４バイト長の単精度浮動小数点数として保持するので、結局、３次元画像の総画素数に４を乗じたバイト数のメモリを搭載した計算機を使う必要があります。
(4) プログラム xp2tg と xp2tg_dpc で再構成した画像上の CT 値と画素値の対応関係の変更やそれらのヒス
トグラムデータの取得は書庫ファイル

11 / 124

http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.zipに入っている通常 CT 用のプログラム tg2tg で行うことができます。
(5) プログラム his2xp や his2xp_top が出力する画像再構成用のバイナリファイル aXP や rXP のバイト数は
それらのもとになった HIS 形式ファイルのバイト数のおよそ 8/(２×128) = 1/32 です（ただし、これはSIXM のスキャンで撮影する屈折プロファイルの「幅」が 128 画素の場合の比率；この幅の画素数は固定値？）。
(6) バイナリデータファイル aXP や rXP などの仕様は以下の通りです。
１行目（ファイルの先頭から改行コードまで）
以下の３個の数値がタブコード区切りの ASCII テキストで書き込まれている。
[1] SIXM の「有効な」スキャン数、Ns[2] SIXM のスキャンの縦画素数（= 再構成可能なスライス数）、Nz[3] サンプル回転のステップ数、Nv
それ以降
Ns×Nz×Nv 個の IEEE 754 形式の倍精度浮動小数点数（それぞれ８バイト長）のバイナリデータが３次元のスキャンライン順で書き込まれている。
なお、このファイルの形式は以前から sinogram などに使っていたものと概ね互換です。以前のものには
１行目の３個目の値（Nv）の記載がありませんでした。
とりあえずここまで。
Date:
Fri, 19 Jun 2015 14:07:26 +0900
From: Tsukasa NAKANO
To:
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Cc:
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: SIXM 用の画像再構成プログラムの改造
みなさま、なかのです。SIXM 用の画像再構成プログラムを少しだけ改造しました（バグ・フィックスではありません）。
ここではその内容説明と、改造後のプログラムを 64 bit Windows で実行した時の端末入力の実例について書きます。
(1) 以前に紹介した書庫ファイル his.taz や his.zip に入れてあった実行ファイル his2xp_top.exe を 64 bit
Windows で走らせると HIS 形式ファイルのデータを２GB だけ読み込んだ時点でエラーが発生しました。
これは実行ファイルが利用している Microsoft Visual C の run time library の関数 fseek() が馬鹿なために生じたエラーだったので、それを回避するようプログラム his2xp と his2xp のコードを書き換えました。
(2) 上杉くんからのリクエストに応じて画像再構成プログラム xp2tg と xp2tg_dpc のマルチスレッド処理用のコードを pthread と Windows thread の両方に対応できるものに書き換えました。以前の書庫ファイルのものと差し替えた新しい xp2tg.exe と xp2tg_dpc.exe は Windows thread を使うようにしてあります。

12 / 124

(3) 画像再構成プログラム xp2tg や xp2tg_dpc の起動パラメータとして画素値と CT 値の対応関係を決める
係数値 base と step を指定せず、画素値のビット数 BPS に負の値を指定した場合、３次元再構成画像の各画素の CT 値はそれらの最小値と最大値に基づいて BPS の絶対値のビット数の画素値に正規化された後に TIFF 形式画像ファイルに書き込まれます（つまり、画像再構成の時点で「word 画像」を作成できます）。
この処理で使った CT 値の最小値と最大値は標準出力に書き出されるログ
それぞれの TIFF 形式画像ファイルのコメント（Image Description Tag）
から得ることができますが、それは手間なので、それらの値を標準エラー出力に書き出すよう xp2tg とxp2tg_dpc のコードを修正しました。
(4) 以上のような改造後のコードや実行ファイルを以前のものと同じ書庫ファイルに入れておきましたので、どうぞご利用ください。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(5) 改造後のプログラムを Windows のコマンドプロンプトの上で実行してみました。その記録を後に添付します。ただし、Windows の Desktop に先の書庫ファイルを展開したディレクトリ "his" と SIXM 実験の測定 150602d のデータが入ったディレクトリ "150602d" をコピーした後に処理を開始しました。
-------- ここから
Microsoft Windows [Version 6.3.9600](c) 2013 Microsoft Corporation. All rights reserved.
C:¥Users¥tsukasa>prompt #$S# cd Desktop# path=C:¥Users¥tsukasa¥Desktop¥his;%PATH%# cd 150602d# more a.log
401,451,100,1,0...
# his2xp_top
usage : his2xp HIS darks I0s scans views aXP rXP {aXP_180 rXP_180}# his2xp_top a.HIS 100 30 401 451 axp.bin rxp.bin# set THREADS=4# xp2tg
usage : xp2tg XP Dr Or t0 {base step} BPS TG/

13 / 124

# mkdir atg# xp2tg axp.bin 150e‐7 200.5 0 ‐8 atg >atg.txt‐1051.934937
922.537537# xp2tg_dpc
usage : xp2tg XP Dr Or t0 {base step} BPS TG/# mkdir rtg# xp2tg_dpc rxp.bin 1 200.5 0 ‐8 rtg >rtg.txt‐1.003630 0.917540-------- ここまで
とりあえず以上です。

14 / 124

Date:
Tue, 23 Jun 2015 19:41:59 +0900
From: Tsukasa NAKANO
To:
Kentaro UESUGISubject: Re: paper
うえすぎさま、なかのです。すみません。E-mail に今頃気づきました。論文のコピー、どうもありがとうございます。著者のMartino は Pfeiffer のお仲間のようですね。これから読んでみます。とり急ぎ、お礼まで。
On Mon, 22 Jun 2015 14:53:13 +0900 Kentaro UESUGI wrote:> 来ました。
>>On Mon, 22 Jun 2015 13:42:50 +1000 Marcus Kitchen wrote:> > Hi Uesugi-san,> > Please find the requested paper attached. This looks quite useful.
> > Best wishes,> > Marcus.
>>> > On 20 June 2015 at 00:45, Kentaro Uesugi wrote:> > > Hi Marcus,> > > Could you give me this article?
> > > https://www.osapublishing.org/ol/viewmedia.cfm?uri=ol-38-22-4813&seq=0&html=true> > > Thank you,> > > K. Uesugi
Date:
From:
To:
Cc:

Wed, 24 Jun 2015 14:01:43 +0900
Tsukasa NAKANO
Kentaro Uesugi"TSUCHIYAMA, Akira"Subject: 位相 != 屈折
うえすぎさま、なかのです。「位相トモグラフィーが可能 == 位相に関しても Radon の定理が成り立つ」について一応理解しました。これは「完全に成り立っている」のではなく「概ね成り立っていると仮定する」ですね。X 線の位相（シフト）は本当は屈折率と対応しているが、X 線の屈折の程度があまりに小さいので、その光路は Radonの定理で仮定している直線と見なせる。そして、物質による位相シフトの総量は光路に沿った積分（「物質の位相シフト」×「それを通過した光路長」の和）なので、その観測値は Radon の定理が対象としている「投影」と等価なものと見なせる。この論理のポイントは「屈折」ではなく「位相（シフト）」という語を用いることですね（と言うのも、地震波トモグラフィーでは地震波が屈折することが本質的なので「画像再構成」にはRadon の定理から即座に導出できる CBP や FBP 法ではなく Huygens-Fresnel の原理にまで戻った手法を使っており、このことを知っている人には「屈折」ではなく「位相」と言わなければ混乱してしまうので）。
とり急ぎ、

15 / 124

Date:
From:
To:
Cc:

Thu, 25 Jun 2015 17:19:47 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM 用画像再構成プログラムに関する話
みなさま、なかのです。これまでに複数の E-mails で説明したように、ぼくが書いたプログラム his2xp（もしくは
his2xp_top）と xp2tg（および xp2tg_dpc）を使えば SIXM で得た HIS 形式のデータから吸収および位相トモグ
ラフィーの CT 画像を再構成できます。ここではその処理に関連した以下の六つの話をします。
(1) プログラム xp2tg（および xp2tg_dpc）のコードの修正(2) バイナリデータファイル XP.bin の（画像再構成以外への）利用法(3) his2xp と his2xp_top の処理内容の違いについて(4) xp2tg_dpc の起動パラメータ Dr に関する話(5) xp2tg のベンチマークテスト(6) SIXM の画像再構成処理に対する TODO
なお、これらのうちで (3) と (4) は割と重要な話だと思われます。
(1) プログラム xp2tg（および xp2tg_dpc）のコードの修正
まず、さほど重要でない話から。6/22 に上杉くんから xp2tg（と xp2tg_dpc）の C 言語コードを Intel C++
Compiler for Windows でコンパイルできないと言われました。かなり手間取りましたが（何せ手元にない処理系への対応なので）、この問題は一応解決しました。修正したコードなどを以前と同じ書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zipに入れておきましたので、どうぞ御利用下さい。
上杉くん、これらの書庫ファイルに入れた xp2tg.c の最新版は 6/22 の最後の E-mail に添付したのものをさらに修正したコードです（「字面」だけを少し変えました）。できることならこの最新版のコードをお使い下さい。
(2) バイナリデータファイル XP.bin の（画像再構成以外への）利用法
6/17 の E-mail でも説明したように、プログラム his2xp と his2xp_top が出力する「投影画像（群）」のバイナリ
データファイル XP.bin に対して「xp2tg と xp2tg_dpc による画像再構成」以外の処理を行うための以下の３つのヘルパープログラムを前記の書庫ファイル his.taz と his.zip に入れておきました。
X 線透過率の動画を作成するプログラム
xp2xt XP.bin {base step} XT.gif

16 / 124

１枚のスライスの画像再構成用の sinogram のデータを抽出
xp2sg XP.bin スライス番号 SG.bin
XP.bin の最初の１枚の投影画像を TIFF 形式画像ファイルに変換
a2tiff XP.bin {base step} BPS XP.tif
xp2xt は通常の X 線 CT 実験の時の画像再構成処理において C-shell script "xts.csh" で作る "xts.gif" と同等な GIF 形式の動画の作成用プログラムです。パラメータ base と step の値を指定しなくてもそれなりの表示コントラストの動画を作成できるので、このプログラムの詳しい説明は不要だと思われます。
また、xp2sg も "test.csh" によって「画像再構成のテスト」を行う際に使っているものと同等な sinogram のバイナリデータファイルを作るためのプログラムです。通常 X 線 CT の「中野的画像再構成法」の説明書
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.pdf
で紹介しているプログラム sg2tg など（これに類する複数のプログラムがあります）を使えば、SG.bin から
スライス画像を再構成できます。
それから a2tiff も「中野的画像再構成」用の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.zipに入れてあるプログラムですが、前記の説明書 sp8ct.pdf にはそれについての記載がありません。このプログラムの具体的な使用例を後で紹介します。
(3) his2xp と his2xp_top の処理内容の違いについて
6/17 の E-mail にも少し書いたように、his2xp と his2xp_top は HIS 形式ファイルに書き込まれた測定データから「投影値」を計算する際に用いる I0 の値の合成（or 推定）法が異なります。話を簡単にするため SIXM
で幅が１画素の「屈折プロファイル」を測定した場合の吸収トモグラフィーだけを考えます。この時、同じサン
プル回転角のスキャンのそれぞれで得たプロファイルの縦列を横に並べたものがその回転角の投影画像に対応する「（X 線強度の）測定画像」になります（ただし、実際の SIXM の測定では回転角ごとに１列づつ「ゴミ」を撮影していますが、それはあらかじめ除外しておきます）。そして、測定画像の横（もしくは左右）方向の座標値を x = ０～X、縦（上下）方向のそれを y = 0～Y、そして位置 (x,y) の測定画像上の画素の暗電流（dark）補正済みの値を I(x,y)、投影画像上の投影値を P(x,y) とします。
his2xp による投影値の計算法の概要は以下の通りです。まず、測定画像の左右両端のそれぞれから指定した合計 I0s 個の縦列で I0 を測定したと考えて（つまり、これらの縦列の領域には物体像がないと仮定して）、
それらに含まれる横列ごとの値の平均値 I1(y) と I2(y) を計算しておきます。その後、これらを x 方向に線形補間した値を I0 として、横列の各画素の投影値を計算します：
his2xp による投影値の計算の擬似コード
i0 = I0s - 1
x0 = i0 / 2
dx = X - i0
for y = 0 ～ Y
I1(y) = { Σ_x=0～i0 I(x,y) } / I0s

17 / 124

I2(y) = { Σ_x=dx～X I(x,y) } / I0s
for x = 0 ～ X
r = ( x - x0 ) / dx
P(x,y) = Ln[ { ( 1 - r ) * I1(y) + r * I2(y) } / I(x,y) ]一方、his2xp_top では測定画像の上端から I0s 個の横列群で I0 を測定したと考えて、それらの縦列ごとの平均値 I3(x) をまず計算します。その後、その値を I0 として縦列の各画素の投影値を計算します：
his2xp_top の擬似コード
for x = 0 ～ X
I3(x) = { Σ_y=0～I0s-1 I(x,y) } / I0s
for y = 0 ～ Y
P(x,y) = Ln[ I3(x) / I(x,y) ]his2xp と his2xp_top で作った測定 150602d の吸収と位相トモグラフィー用の投影画像のそれぞれの最初のもの（サンプル回転角が０度のもの）だけを前記の a2tiff で TIFF 形式画像に変換してみました。Linux 上で実行した時の手順をこの E-mail に添付したテキストファイル 150602d.txt に記しました。
150602d.txt（および以前の E-mail で紹介した例）では I0 とする列の数 I0s として 30 を指定しましたが、この値に深い意味はありません。I0s を変えると投影画像の画質がどう変化するかを試してみるのも一興です。
a2tiff の使用法について少し書いておきます。このプログラムは投影値を TIFF 画像に変換すると同時に以下の２行・４個の値を標準出力に書き出します（ただし、150602d.txt では１行にまとめました）。
１行目
投影画像の横および縦方向の画素数
２行目そこに出現した投影値の最小値と最大値
his2xp と his2xp_top のそれぞれが出力した吸収もしくは位相の画像上の投影値と画素値の対応関係を揃えるため、150602d.txt ではまず a2tiff で出力ファイルのデータを "/dev/null"（Windows なら"nul"）に捨てて投影値の最小値と最大値を調べています。その後、比較すべき画像ペアの投影値す
べてが８ビット画素値にちょうどおさまるような base と step の値を指定して a2tiff を再実行しています。
ただし、そのための計算は UNIX の標準コマンド bc で行ない、a2tiff はその結果の値を格納した C-shell
変数を（"$bs" として）引用しています。
こうして作った２×２個の画像を GIF 形式のものに変換してこの E-mail に添付しておきましたので、吸収と位相のそれぞれの画像ペアを比較してみて下さい。
axp.gif と rxp.gif
his2xp で得た吸収と微分位相の投影画像
axp_top.gif と rxp_top.gif
his2xp_top で得た吸収と微分位相の画像

18 / 124

さて、これらのすべて画像に「縦線」が入っています。つまり、スキャンごとの入射 X 線強度が大きく変動している（入射 X 線強度に短周期の時間的な変動がある）ようです。上杉くんによるとこの原因はモノクロメータの振動らしいです。
his2xp_top はスキャンごとに推定した I0 を使っているので X 線強度の短周期の時間変動をある程度軽減できているはずですが、それでも不十分です。一方、竹内さんのプログラム「HIDEKI」に準じた方法で I0 を推定している his2xp は X 線強度の空間的変動（「縦方向」の変動）を軽減できるはずですが、短周期の時間変動に対しては無力です。そして、通常の X 線 CT 実験の時のように測定の前後に I0 の画像を独立に「まとめ撮り」してもこれらの変動を軽減できそうにありません（SIXM の測定シーケンスを I0 の「まとめ撮り」を追加したものに変更するよう竹内さんにリクエストしていましたが、それは取り下げました）。
残念ながら X 線強度の時間と空間の両方の変動を軽減できる測定法や補正法をちょっと思い付きません。
このあたりが SIXM 測定の最大の問題点でしょうね。
すみません。疲れてきたので、以下の話は次便で行います。とり急ぎ、
添付ファイル 150602d.txt# head ‐1 150602d/a.log
401,451,100,1,0# his2xp 150602d/a.HIS 100 30 401 451 axp.bin rxp.bin# his2xp_top 150602d/a.HIS 100 30 401 451 axp_top.bin rxp_top.bin# a2tiff axp.bin 8 /dev/null
400
400‐0.369277
0.762644# a2tiff axp_top.bin 8 /dev/null
400
400‐1.415264
0.739011# set min=‐1.415264# set max=0.762644# set bs=($min `echo "( $max ‐ $min )/255" | bc ‐l`)
# a2tiff axp.bin $bs 8 axp.tif
400
400‐0.369277
0.762644# a2tiff axp_top.bin $bs 8 axp_top.tif
400
400‐1.415264
0.739011# a2tiff rxp.bin 8 /dev/null
400
400‐3.574706
4.888005# a2tiff rxp_top.bin 8 /dev/null
400
400‐3.422092
4.795531# set min=‐3.574706# set max=4.888005# set bs=($min `echo "( $max ‐ $min )/255" | bc ‐l`)

19 / 124

# a2tiff rxp.bin $bs 8 rxp.tif
400
400‐3.574706
4.888005# a2tiff rxp_top.bin $bs 8 rxp_top.tif
400
400‐3.422092
4.795531# rm axp.bin rxp.bin axp_top.bin rxp_top.bin
添付ファイル axp.gif、axp_top.gif、rxp.gif および rxp_top.gif
これらは動画であり、各自で容易に作成できると思われるので、ここには掲載しません。
Date:
From:
To:
Cc:

Thu, 02 Jul 2015 13:41:26 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM 用画像再構成プログラムに関する話の続き
みなさま、なかのです。すみません。先程の E-mail に間違いを見つけたので、それを修正したものを再送します。先のものを破棄してこちらを御覧下さい。とり急ぎ、みなさま、なかのです。先日の E-mail の続きです。
On Thu, 25 Jun 2015 17:19:47 +0900 Tsukasa NAKANO wrote:> ... プログラム his2xp（もしくは his2xp_top）と xp2tg（および xp2tg_dpc）> を使えば SIXM で得た HIS 形式のデータから吸収および位相トモグラフィーの> 画像を再構成できます。ここではその処理に関連した以下の六つの話をします。
>(1) プログラム xp2tg（および xp2tg_dpc）のコードの修正>(2) バイナリデータファイル XP.bin の（画像再構成以外への）利用法>(3) his2xp と his2xp_top の処理内容の違いについて---- 先の E-mail ではここまで書いたところで力尽きた。
>(4) xp2tg_dpc の起動パラメータ Dr に関する話>(5) xp2tg のベンチマークテスト>(6) SIXM の画像再構成処理に対する TODO
(2) と (3) の補足
すみません。投影のバイナリデータを画像化するプログラム a2tiff の使い方をきちんと説明していませんでした。その起動法は先の E-mail に書いた通りです。
a2tiffXP.bin

{base

step}

BPS XP.tif

20 / 124

起動パラメータ XP.bin と XP.tif はそれぞれ投影のバイナリデータファイルとそれを変換した TIFF 形式画像
ファイルの名前です。また、base は画素値０に対応させる投影値、step は画素値１階調あたりの投影値の増分、そして BPS はその画像の画素値のビット数です。ただし、base と step の指定は省略可能で、その場合には画像上の投影値すべてが BPS で指定した画素値の範囲にちょうどおさまるような base と step の値を使って変換処理を行います。
なお、先の E-mai l にも書いたように、この処理の終了後に a2tiff は画像の横・縦画素数と画像上の投影値の最小・最大値の合計４個の値を「標準エラー出力」に書き出します（先の E-maiｌ に書いた「標準出力」は間違いです；すみません）。
(4) xp2tg_dpc の起動パラメータ Dr に関する話
6/18 の２通目の E-mail にプログラム xp2tg_dpc で位相（屈折）トモグラフィーの画像再構成を行う場合、その起動パラメータ Dr の値として１を指定すれば良いと書きましたが、それは間違いでした。結論から言うと、以下の値を Dr として指定してやればサンプル内部の位相シフトの「真の値」を得ることができます。
Dr =L / ( Δr・k ) = λ・L / ( ２・π・ Δr )

[1]

ただし、
L
サンプルと投影面の間の距離
Δr
SIXM の屈折角測定用の画像の画素の辺長k = ２・π / λ
X 線の波数λX 線の波長（E = h・c/λより、λ[meter] ～ 1240 / E [eV]）
この式の導出には SIXM と微分位相の画像再構成用プログラムの両方の原理？に関する知識が必要です。
ちょっと長くなりますが、それらを以下で説明します。
まず SIXM に関してですが、それについてぼくが理解していることをこの E-mail に添付したファイル sixm.gifと sixm_ra.txt に記しておきました。
上杉さま、
竹内さま、
SPring-8 の SIXM システムの「売り」である FZP などに関することを省略してしまいましたが、上記のファイルに記した SIXM 測定の「原理」に関するぼくの理解に大きな間違いはありませんよね？

21 / 124

これらに記したように、SIXM を使えば屈折角Θr を介してサンプル内部の微分位相 ∂φ/∂r の線積分の値（微分位相の投影値）を計測することができます：∫ (∂φ/∂r) ds ≒ k・sin Θr ≒ k・tan Θr ≒ k・( rc - r0 ) / L[2]

sixm_ra.txt では上記の線積分の値に負の符号がついています。これは屈折率と位相シフトの関係
式からそうなるはずですが、WEB で見つけた 高野秀和さんや篭島靖さんの文書ではそうなっていませんでした。なぜそうするのかわかりませんが、とりあえず負符号なしの式を使いました。
上記の r0 と rc はそれぞれサンプルへの入射 X 線とサンプル内部で屈折した透過 X 線を受けた
検出器の位置で、実際には空間的な「ばらつき」があります。プログラム his2xp や his2xp_top では
r 方向に並んだ検出器列で受けた X 線強度（光子数）の空間分布の「重心」を r0 や rc と見なしています。
次に屈折トモグラフィー用の画像再構成プログラムに関する話です。これはぼく的には込み入った話なので順を追って説明します。まず、微分位相の投影値から位相シフトの画像再構成を行うプログラム xp2tg_dpcの主要部（CBP 法）の計算機コードは吸収トモグラフィー用のものとまったく同じです。このコードに add-on
する「再構成フィルタ」だけを変えており、それらは以下のように X 線強度検出器の間隔 δr の２もしくは１乗に反比例した値を返す関数です。
吸収トモグラフィー用の再構成フィルタ、g(r) ∝ 1 / δr^2
位相トモグラフィー用の再構成フィルタ、h(r) ∝ 1 / δrh(r) は Pfeiffer 達（2007）が考案した「CBP 法による画像再構成の畳み込み積分の際に微分位相の投影値も積分して位相シフトの投影値に 変換する」再構成フィルタです。なお、Pfeiffer 達の論文（PRL、98、10815, 2007）に載っていた h(r) の計算式をこの E-mail に添付したファイル
pfeiffer.txt で導出しているので、よろしければ御覧下さい。
15 年以上前に書いた以下の文書（CBP 法による画像再構成の「仕様書」）の式 [13] が g(r) の具体的な計算式です。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/CBP.A4.pdf
http://www-bl20.spring8.or.jp/~sp8ct/tmp/CBP.A5.pdf
また、書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/dpc.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/dpc.zipに入れた、手書きノートを変換した PDF ファイル dpc.pdf の１ページ目で g(r) と h(r) の具体的な計算式を導出しています。よろしければこれらの PDF も御覧下さい。
CBP 法ではこれらの再構成フィルタと与えられた投影値 p(r) の畳み込み積分を行うことにより "filteredprojection" q(r) を計算します。
q(r) = ∫ p(ρ) ・[g,h] (r - ρ) dρ

22 / 124

ただし、CBP 法の計算機コードはこのρに関する積分を δr を「刻み幅」とした和分近似（∫F(ρ) dρ ≒ Σj
F(ρj)・δr）で計算するので、結局、位相トモグラフィーの filtered projection は δr に依存しない値になります。
吸収トモグラフィーの filtered projection ∝ (1/δr^2)・δr = 1/δr
位相トモグラフィーの filtered projection ∝ (1/δr)・δr = 1
吸収トモグラフィー用の画像再構成プログラム xp2tg では起動パラメータ Dr として指定した値で filteredprojection を割るようになっています。つまり、吸収トモグラフィーの場合には Dr = δr を指定すれば良いわけです。
ただし、SIXM の吸収トモグラフィーの画像再構成で指定する値 δr は式 [1] のΔr ではありません。
SIXM のスキャンのステップ幅（SIXM 実験の場で「松、竹、梅」と呼んでいた値）を指定して下さい。
一方、xp2tg の主要部分を流用したプログラム xp2tg_dpc で位相トモグラフィーの画像を再構成する場合、与えられた投影値が適切なものであれば Dr = 1 を指定すれば良いことになります（これが 6/18 の２通目のE-mail の記述の根拠です）。しかしながら、プログラム his2xp や his2xp_top が出力する位相トモグラフィー用の投影値 p は微分位相の投影値ではなく、SIXM の屈折角測定用の画像上の「画素の幅Δr を単位とする屈折量」です：p = ( rc - r0 ) / Δr[3]

それゆえ、xp2tg_dpc で位相シフトの値を再構成するには式 [2] と [3] を考慮した適切な起動パラメータの値 Dr を指定してやる必要があります。そして、式 [2] と [3] より
p = L / (Δr・k)・∫ (∂φ/∂r) ds
なので、結局、式 [1] で表される Dr の値を指定してやれば xp2tg_dpc が画像再構成に使う p/Dr の値は微分位相の投影値 ∫ (∂φ/∂r) ds に一致します。あぁ、ややこしかった。
すみません。またまた力尽きてしまったので、残りは次便で。とり急ぎ、

23 / 124

添付ファイル sixm.gif

添付ファイル sixm_ra.txt
SIXM で測定できる屈折角、Θr
SIXM の測定時の座標系 (s, r, z)srz

光源から投影面に向かう投影面に直交する座標軸
SIXM のスキャン方向と平行な投影面上の座標軸
s 軸と r 軸の両方に直交する投影面上の座標軸

サンプル透過時の X 線光路の幾何学は下記の光線方程式で決まる。
d ( n τ) / dt = grad(n)n と grad(n)サンプルの内部の屈折率の値とその空間勾配（ベクトル） 。ただし、λを X 線の波長、φ を
その位相シフトとすると、
k = 2 π / λ（X 線の波数）n = 1 - φ/ k
また、前記の座標系における grad(n) の成分値は、grad(n) = (∂n/∂s, ∂n/∂r, ∂n/∂z )24 / 124

τと t
光線の単位接線ベクトルとその光路に沿った座標値（長さ）。ただし、X 線の屈折は微小なので、サンプル内部でのその光路は s 軸と平行（t≒ s）と見なしても実用上は問題ない。しかしながら、内部での屈折が積算されるため、サンプルを透過直後のτの方向は s 軸の方向から有意に偏移する。
サンプルを透過直後の X 線の屈折角
光線方程式をサンプル内部の光路（≒ s 軸）に沿って積分すると、∫ { d (n τ
) / dt } dt = [ n τ]_<entrance and exit>= ∫ grad(n) dt ≒ ∫ grad(n) ds
サンプル表面の物質の屈折率が１で、サンプルへの X 線の入射方向が s 軸と平行と見なせる場
合、この式を使ってサンプルを透過した後の X 線の方向からサンプル内部の屈折率の空間勾配の線積分の値を推定できる。つまり、サンプルを透過直後の光路の接線ベクトルの [r,z] 方向への方位角（s 軸方向を０とする）をΘ[r,z] とするとτ
の [r,z] 成分値は sin(Θ[r,z]) なので、
sin(Θ[r,z]) ≒ ∫ (∂n/∂[r,z]) ds = - (1/k) ∫ (∂φ/∂[r,z]) ds
SIXM では r と z の２方向の屈折角のうちのΘr だけを測定できる。
添付ファイル pfeiffer.txt
微分投影が与えられた場合の再構成フィルタの式（Pfeiffer、2007）
CBP 法や FBP 法で使用している（普通の）投影 p(r) に対する再構成フィルタ
G(k) = |k| （波数空間の再構成フィルタ）g(r) = ∫ G(k)・exp(2・π・i・r・k) dk （実空間の再構成フィルタ）
filtered projection の計算式
FBP 法（filtering）
投影の Fourier 変換
P(k) =∫ p(r)・exp(-2・π・i・k・r)
filtered projection
q(r) = ∫ P(k)・G(k)・exp(2・π・i・r・k) dk
CBP 法（convolution）
q(r) = ∬ p(ρ)・G(k)・exp{2・π・i・(r - ρ)・k} dk dρ= ∫ p(ρ)・g(r - ρ) dρ

[1]
[2]

[3]

25 / 124p(r) の代わりに微分投影∂p/∂r が与えられた場合
P(k) の計算式 [1] を部分積分して p(r) を∂p/∂r に置き換える。
∫ x・(dy/dz) dz = [ x・y ] - ∫ (dx/dz)・y dz∫ exp(-2・π・i・k・r) dr = exp(-2・π・i・k・r) / (-2・π・i・k)
ただし、その際に以下を仮定する。
p(r) = 0

at r = ±∞上記の公式と仮定から、
P(k) = [ p(r)・exp(-2・π・i・k・r) / (-2・π・i・k)]_r=±∞-∫(∂p/∂r)・exp(-2・π・i・k・r) / (-2・π・i・k) dr= {-i / (2・π・k)}・∫(∂p/∂r)・exp(-2・π・i・k・r) dr
この P(k) の式を filtered projection の計算式 [2] に代入すると、
q(r) = ∬ (∂p/∂ρ)・{-i / (2・π・k)}・G(k)・exp{2・π・i・(r - ρ)・k} dk dρ= ∬ (∂p/∂ρ)・H(k)・exp{2・π・i・(r - ρ)・k} dk dρ= ∫ (∂p/∂ρ)・h(r - ρ) dρ
ただし、H(k) = {-i / (2・π・k)}・G(k)= {-i /(2・π)}・(|k| / k) (波数空間の再構成フィルタ）h(r) = ∫ H(k)・exp(2・π・i・r・k) dk (実空間の再構成フィルタ）
それゆえ、h(r) を g(r) の代わりに使えば、投影から計算する場合と同じ式 [3] で微分投影から
filtered projection を計算することができる。
Date:
Tue, 07 Jul 2015 15:26:59 +0900
From: Tsukasa NAKANO
To:
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Cc:
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: SIXM 用画像再構成プログラムに関する話の続きの続きみなさま、なかのです。6/25 と 7/2 の E-mail の続きです。
On Thu, 25 Jun 2015 17:19:47 +0900

Tsukasa NAKANO wrote:26 / 124

> ... プログラム his2xp（もしくは his2xp_top）と xp2tg（および xp2tg_dpc）> を使えば SIXM で得た HIS 形式のデータから吸収および位相トモグラフィーの> 画像を再構成できます。ここではその処理に関連した以下の六つの話をします。
>(1) プログラム xp2tg（および xp2tg_dpc）のコードの修正>(2) バイナリデータファイル XP.bin の（画像再構成以外への）利用法>(3) his2xp と his2xp_top の処理内容の違いについて---- 6/25 の E-mail ではここまで書いたところで力尽きた。
>(4) xp2tg_dpc の起動パラメータ Dr に関する話---- 7/2 の E-mail ではここまで書いたところで力尽きた。
>(5) xp2tg のベンチマークテスト>(6) SIXM の画像再構成処理に対する TODO
(5) xp2tg のベンチマークテスト
xp2tg は SIXM の吸収トモグラフィーの画像再構成用に書いたプログラムですが、これは SPring-8 で行っている「普通の」吸収トモグラフィーの画像再構成にも使えます。ただし、それを行うためにはプログラム
hp2xp を用いて普通の X 線 CT 測定のデータファイル（output.log と「HiPic 形式」の測定画像のファイル
dark.img や "q*.img"）から xp2tg が処理できる投影データファイル XP.bin を作成しておく必要があります。
hp2xp はそもそもはつちやまさんが今年度導入するμ-focus X 線 CT（MFX）の画像再構成のために書いた、３月末の E-mail で紹介した下記の MFX 用の書庫ファイルに入れてあるものとまったく同じプログラムです。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/fdk.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/fdk.zip
SIXM 用の書庫ファイル his.taz と his.zip にも hp2xp のソースファイルと Windows 用実行ファイルをコピーしておきました。さらに、3/3 の E-mail で紹介した hp2xp の起動法もこの E-mail に添付したファイル hp2xp.txtに転記しました（hp2xp の具体的な使用例は後で紹介します）。
さて、このような hp2xp と xp2tg による普通の X 線 CT 測定の画像再構成の処理時間やその画像の画質を
既存のプログラム"hp2tg*"の実行結果と比較してみました。ただし、"hp2tg*"は「hp2xp による測定データからの投影値の計算」と「xp2tg による（CBP 法を用いた）画像再構成」の処理を一括して行うプログラムで、
この E-mail に添付したファイル hp2tg.txt に記してあるように、合計すると 30 種類のもの（12 個の CPU 版と18 個の GPU 版）があります。それらのうちで「測定データを最初にまとめてメモリに読み込み」、「CPU のマ
ルチスレッド処理」を利して「サンプル回転角方向の補間なしの画像再構成」を行う以下の３個（GPU を用いて画像再構成を行う１個と CPU だけですべての処理を実行する２個）を今回の比較の対象としました。
hp2tg_tm
hp2tg_tint_OM
hp2tg_tnai_OM

GPU が単精度浮動小数点数演算で画像再構成を行う。
逆投影を整数演算で高速に実行する。
逆投影を倍精度浮動小数点数演算で実行する。

27 / 124

上記の並びはプログラムの処理速度の降順になっていて（hp2tg_tm が最も速い）、また、計算精度は
hp2tg_tnai_OM、hp2tg_tm、hp2tg_tint_OM の順です。そして、hp2xp + xp2tg の処理速度と計算精度はそれ
ぞれ hp2tg_tnai_OM と hp2tg_tm のそれと同程度のはずです。
xp2tg と hp2tg_tnai が同程度の計算速度になるのは同じスレッド数のマルチスレッド処理で実行した場合の話です。なお、これらのスレッド数の指定に使う「環境変数」が異なることに注意して下さい。
xp2tg
hp2tg_t*

THREADS
CBP_THREADS

過去の測定 040711j（サンプルは多数のビーズ球）と 020607a（金沢大学・森下くんの symplectite）の画像を
hp2xp + xp2tg と hp2tg_t* のそれぞれで再構成してみました。これらの測定データの概要はこの E-mail に添付したファイル xp2tg.txt に記した通りです。なお、後述する理由で「横画素数（X 線検出器の総数）N が奇数の測定画像」の画像再構成を試したかったので、020607a の測定画像として再構成画像を「X 線 CT
シミュレータ」で「再撮影」したものを使いました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct_demo.pdf
hp2tg_t* と hp2xp + xp2tg による画像再構成の具体的な実行例も xp2tg.txt に記しておきました。ただし、すべてのスライスを再構成する場合は hp2xp にそれらの範囲を指定する必要はありません。また、Windows
では hp2xp から xp2tg への投影値のバイナリデータの受け渡しにパイプラインを利用できません。
これら２個の実行例の最後の「dit_tbl」で始まる複数行の実行に要した総処理時間と「画像再構成の時定数」
は xp2tg.txt の後半に記した通りです。
画像再構成の時定数 = 総処理時間 / (再構成画像の総画素数 × 投影数 views）
これらの処理時間や時定数の値から処理速度は概ね当初の予想通りで、hp2xp +xp2tg は hp2tg_tnai_OM
よりも少しだけ高速なようです。なお、xp2tg.txt に示したように今回のベンチマークテストは３台のホスト計算機で実行しましたが、それらのハードウェア、OS や C-compiler などに関する話は割愛します。
ハードウェアに関することを一つだけ言っておくと、今回は SSD ではなく HDD 上の測定データファイルを処理したので、hp2tg_t* の処理は以前に行ったベンチマークテストの時よりもかなり低速でした（特に、 GPU による画像再構成の時定数は最速の場合の数倍大きい値でした）。
次に hp2xp + xp2tg と hp2tg_t* で再構成した CT 値（再構成画像の画質）についてですが、実は xp2tg が出力する再構成画像は（hp2tg_t* を含む）従来のプログラムのものとスライス面内の画素の数や「CT 値を計算する点の位置」が違う（ことがある）ので、それらを単純・直接的な方法で比較できません。
もう少し詳しく言うと、まず、従来のプログラムは測定画像の横画素数 N を各辺沿いの画素数とする正方形のスライス画像を再構成しますが、３次元画像の各画素の CT 値すべてをメモリ上に保持しつつ処理を行う xp2tg ではメモリの節約のため「センター値」に応じてその半径が変化する「画像再構成可
能な円形領域」がちょうどおさまるスライス画像を出力するようになっています（つまり、xp2tg の再構成

28 / 124

画像の画素数は常に従来の再構成画像の画素数以下の個数になっている）。さらに、xp2tg も従来のものと同様にサンプル回転軸を中央に置いた正方形のスライス画像を再構成しますが、その位置が「スライス中央の画素の中央」に来るようにスライス画像の各辺沿いの画素数を常に奇数にします。そして、ほとんどの測定画像の横画素数 N は偶数なので、従来の画像再構成プログラムではこの回転軸の配置を実現できません。それゆえ、xp2tg と従来の画像再構成プログラムでは x と y の両方向に半画素幅だけ異なった画像上の点における CT 値を計算していることになります。
このような微妙な違いがあるので、差分画像の作成などによる hp2xp + xp2tg と hp2tg_t* の再構成画像の直接の比較ができないことがあります。この違いは再構成画像上の CT 値の空間補間によって解消できま
すが、それは面倒なので、今回はそれぞれの再構成画像の CT 値の出現頻度分布（画素値ヒストグラム）を比較することにしました。この E-mail に添付したファイル hp2tg.pdf にその結果のグラフを示しました。
これらのグラフから以下のことがわかります。
測定 040711j
hp2tg_t* で得た３個の画像のヒストグラムの線はほぼ完全に重なっているが、hp2xp + xp2tg の線はそれらと若干ずれている。
020607a４個のヒストグラムの線はほぼ完全に重なっている。
測定 040711j のヒストグラムの hp2tg_t* と hp2xp + xp2tg の線のズレは前記の「CT 値を計算する点の位置の微妙な違い」によるものだと思われます（ただ、それにしては大きな違いなので、もう少し調べてみます）。
020607a では４個のヒストグラムが完全に一致しているので、とりあえずは hp2xp + xp2tg で再構成した画像の画質は従来のものとほぼ同じと言えそうです。
(6) SIXM の画像再構成処理に対する TODO
今後 SIXM に関わる以下のプログラムを書くつもりでいます。
センター値推定プログラム
C-shell script 版は書いたので、Windows でも走る C 言語版を書く。
位相画像のノイズ除去プログラム
6/22 に上杉くんにコピーしてもらった Martias 達の論文（Opt. Lett.、38、4813、2013）を一読しましたが、その中で仮定していることに納得できません（物理的な根拠なしの、解くためだけの仮定のように思える）。SIXM の測定原理を考慮した新しい手法を考案したいと思っています。
また、すでに「HIDEKI」で実行済みかもしれませんが、これまで紹介したプログラムを使って６月の SIXM 測定のデータすべての画像再構成を行うつもりです。
長い E-mail になりました。とりあえず以上です。

29 / 124

添付ファイル hp2xp.txt
On Tue, 31 Mar 2015 18:08:29 +0900 Tsukasa NAKANO wrote:> hp2xp {d.tbl i.tbl t.tbl} {y1 y2}|{x1 y1 x2 y2} XP.bin>> 機能>測定した HiPic 画像から投影値 Pj,l,m のデータを抽出する。
>> 起動パラメータ>
d.tbl、i.tbl および t.tbl> 暗電流（d.tbl）、入射（i.tbl）および透過 X 線強度（t.tbl）を測定> した HiPic のシーケンス画像のファイル名のそれぞれとそれらの測定> 時刻（単位は任意）の２個の値が各行に書き込まれたテキストファイル> の名前。これらの指定を省略すると標準入力から測定ファイルの名前と> 測定時刻の情報を読み込む（hp2sg の場合と同様で、dit_tbl の出力を
> hp2xp の標準入力に流し込んでやれば良い）。
>
y1 と y2 もしくは x1、y1、x2 と y2
> 出力する投影画像の y 座標（～ 再構成画像のスライス番号）もしくは> そのひろがりを表す「bounding box」の２個の対角点の座標値。これら> の指定を省略すると測定した HiPic 画像のものと同じひろがりの投影> 画像のバイナリデータを出力する。
>
XP.bin> 投影値 Pj,l,m を格納するバイナリデータファイルの名前。その１行目> にはタブコード区切りで Nu、Nw と Nθ の３個の自然数の値が、また、> その後（正確には、改行コードの直後）には倍精度浮動小数点数の内部> 形式で表現された投影値 Pj,l,m のポータブルではないデータが添え字> j、l および m に関するスキャンライン順で書き込まれる。XP として> "-"（負符号もしくはハイフン）を指定するとこのデータは標準出力に> 書き出される（ただし、Windows ではバイナリデータを標準出力に流す> ことは許されていないかもしれない）。
添付ファイル hp2tg.txt
hp2tg* の起動法
hp2tg* {d.tbl i.tbl t.tbl} slice1 slice2 dr DO_base DO_step angle {CTV_base CTV_step} BPS TG/
CPU ですべての処理を行うプログラム
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.zip
実行ファイルの命名則
マルチスレッド処理か否か（２種類）
hp2tg_t*
マルチスレッドで処理30 / 124

その他
シングルスレッドで処理
逆投影の計算法の違い（３種類）
hp2tg_*int*
角度方向の補間なしの整数演算で逆投影
hp2tg_*nai*
角度方向の補間なしの実数演算で逆投影
その他
角度方向の補間を行いつつ実数演算で逆投影
投影データファイルのアクセス法の違い（３種類）
hp2tg_*OM
処理の最初にファイル上のデータすべてをメモリに読み込む
hp2tg_*OF
データファイルすべてを開いたままにして処理する
その他
ファイルをその都度開いてデータを読む
実行ファイルの名前
マルチスレッドで処理
hp2tg_tint_OM, hp2tg_tint_OF, hp2tg_tint
hp2tg_tnai_OM, hp2tg_tnai_OF, hp2tg_tnai
hp2tg_t_OM, hp2tg_t_OF, hp2tg_t
シングルスレッドで処理
hp2tg_int_OM, hp2tg_int_OF, hp2tg_int
hp2tg_nai_OM, hp2tg_nai_OF, hp2tg_nai
hp2tg_OM, hp2tg_OF, hp2tg
CUDA GPU で画像再構成を行うプログラム
http://www-bl20.spring8.or.jp/~sp8ct/tmp/cuda_cbp.taz
実行ファイルの命名則
GPU と CPU の並列処理および CPU のマルチスレッド処理の如何（３種類）
hp2tg_t*
GPU の処理と CPU のマルチスレッド処理を並列に実行
hp2tg_o*
GPU の処理と CPU のシングルスレッド処理を逐次実行
その他
GPU の処理と CPU のシングルスレッド処理を並列に実行
投影データファイルのアクセス法の違い（２種類）
hp2tg_*m*
処理の最初にファイル上のデータすべてをメモリに読み込む
その他
ファイルをその都度開いてデータを読む
逆投影の計算法の違い（２種類）
hp2tg_*wai
逆投影の際に角度方向の補間を行う
その他
角度方向の補間を行わない。
実行ファイルの名前
マルチスレッドかつ並列処理
hp2tg_tmwai, hp2tg_twai, hp2tg_tm, hp2tg_t
シングルスレッドかつ逐次処理
hp2tg_omwai, hp2tg_owai, hp2tg_om, hp2tg_o
シングルスレッドかつ並列処理
hp2tg_mwai, hp2tg_wai, hp2tg_m, hp2tg
添付ファイル xp2tg.txt
040711jsize of projection images = 1000 * 720 pixels * 360 views
interval of detectors, dr = Dr = 0.000583 cm

31 / 124

position of rotation center, -DO = Or = 498.5
size of tomographic images = [1000,997]^2 pixels * 720 slices
http://www-bl20.spring8.or.jp/~sp8ct/tmp/040711j.gif
020607asize of projection images (simulated) : 481 * 509 pixels * 450 views
interval of detectors, dr = Dr = 0.00007184104862007412 (cm)
position of rotation center, -DO = Or = 240
size of tomographic images = 481^2 pixels * 509 slices
http://www-bl20.spring8.or.jp/~sp8ct/hd2/morishita_symplectite/020607a/byte_222.gif
image reconstruction by "hp2tg_t*" under C-shell environment# setenv CBP_THREADS 8# mkdir 040711j/tg# dit_tbl 040711j/raw | \hp2tg_t* 0 719 \
0.000583 -498.5 0 0 \0 0.02930016 8 040711j/tg > 040711j/tg.log
image reconstruction by hp2xp and xp2tg under C-shell environment# setenv THREADS 8# mkdir 020607a/tg# dit_tbl 020607a/raw | \hp2xp 0 508 - | \xp2tg - 0.00007184104862007412 240 0 0 1 8 020607a/tg > 020607a/tg.log
calculation time in sec.
host / GPUprogramgsjgix / K20c
hp2tg_tm/hp2tg_tint_OM/hp2xp + xp2tg/hp2tg_tnai_OM
gsjvix / C2070
hp2tg_tm/hp2tg_tint_OM/hp2xp + xp2tg/hp2tg_tnai_OM
vrm / GTX750Ti hp2tg_tm/hp2tg_tint_OM/hp2xp + xp2tg/hp2tg_tnai_OM

040711j37.000848
70.974262
130.093564
117.461060
36.611471
95.058632
145.413615
147.363683
38.605144
123.978912
181.122563
241.213953

020607a16.456110
11.492650
21.766453
16.968540
9.581409
18.520894
29.096570
23.545237
12.785267
27.870776
35.419727
51.583618

calculation time constant in nano sec.
host / GPUprogramgsjgix / K20c
hp2tg_tm/hp2tg_tint_OM

040711j0.142750
0.273820

020607a0.310532
0.216870

32 / 124//gsjvix / C2070///vrm / GTX750Ti///-

hp2xp + xp2tg
hp2tg_tnai_OM
hp2tg_tm
hp2tg_tint_OM
hp2xp + xp2tg
hp2tg_tnai_OM
hp2tg_tm
hp2tg_tint_OM
hp2xp + xp2tg
hp2tg_tnai_OM

0.504929
0.4531680.141248
0.366739
0.564391
0.568533
0.148940
0.478314
0.702987
0.930609

0.410740
0.320202
0.180804
0.349495
0.549062
0.444306
0.241262
0.525931
0.668382
0.973400

添付ファイル xp2tg.pdf
0

50

100

150

200

250

0

50

100

150

200

250

1e+07

020607a1e+06

log(pixels)

log(pixels)

hp2xp + xp2tg
hp2tg_tnai_OM

1e+06

10000010000

040711j

hp2tg_tm (K20c)
hp2tg_tint_OM
hp2xp + xp2tg

1000

hp2tg_tm (K20c)
hp2tg_tint_OM

100000

hp2tg_tnai_OM
100

0

1

2

3

4

5

CT-value (estimated LAC, 1/cm)

6

7

0

50

100

150

200

250

CT-value (estimated LAC, 1/cm)

Date:
From:
To:
Cc:

Tue, 07 Jul 2015 18:04:00 +0900
Tsukasa NAKANO
Kentaro UESUGI"TAKEUCHI, Akihisa",
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: Re: SIXM 用画像再構成プログラムに関する話の続きの続き
うえすぎさま、なかのです。E-mail ありがとうございます。あなたの理解の通りです。hp2xp のバイナリデータ（>8*N*L*M
バイト）は以下のようなコードで処理できます。
バイナリデータファイル "XP.bin" の書き込み#include <stdio.h>#include <stdlib.h> /* 関数 exit() を使うため */int main() {FILE *file;
int m,l,n,N=100,L=200,M=300;
double p;if ((file=fopen("XP.bin","wb"))==NULL) exit(1);(void)fprintf(file,"%d¥t%d¥t%d¥n",N,L,M);33 / 124

for (m=0; m<M; m++)for (l=0; l<L; l++)for (n=0; n<N; n++) {p=出力する値;(void)fwrite(&p,sizeof(double),1,file);}(void)fclose(file); return 0;}
バイナリデータファイル "XP.bin" の読み込み#include <stdio.h>#include <stdlib.h> /* 関数 exit() を使うため */#define LEN 1024/* １行目の文字列用 */int main() {FILE *file;char str[LEN];
int N,L,M,m,l,n;double p;if ((file=fopen("XP.bin","rb"))==NULL ||fgets(str,LEN,file)==NULL ||sscanf(str,"%d %d %d",&N,&L,&M)!=3 ||N<=0 || L<=0 || M<=0) exit(1);for (m=0; m<M; m++)for (l=0; l<L; l++)for (n=0; n<N; n++) {if (fread(&p,sizeof(double),1,file)!=1) exit(1);p の値を使った処理;}(void)fclose(file); return 0;}
とり急ぎ、
On Tue, 7 Jul 2015 17:18:21 +0900 Kentaro UESUGI wrote:> 中野さん>> 上杉です>> ちょっと質問です。hp2xp で生成するバイナリデータは> １行目テキスト：N(画像の横画素数) L(画像の縦画素数) m(投影数)> ２行目以降：倍精度浮動小数点数(64bit/pixel)> で良いでしょうか？２行目は 32bit/pixel ではなかったですかね？
> （で、XP.bin は数十 GB のデカイ単一ファイル？）34 / 124

Date
From:
To:
Cc:

Fri, 26 Jun 2015 17:14:04 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: 180 度の方向から見た微分位相の投影画像
みなさま、なかのです。昨日の E-mail の最後に書いた SIXM で撮影した吸収と微分位相の投影画像に関することですが、これらを使った画像再構成で指定する「センター値」の推定プログラムを書いていて以下のことに気づきました。
吸収の投影画像
これは通常の X 線 CT で経験済みのことですが、サンプルの回転角が 180 度の画像は０度のものの鏡像（左右反転した画像）になっている。
微分位相の投影画像
180 度の画像は左右反転しているだけでなく、値も正負反転している。
この微分位相の投影値の「二重反転」は言われてみれば当たり前のことですが、ちょっと面食らいました。
物体像に付随した座標値を x とすると、０度方向から見た投影画像上の座標値 r は、r0 = x
180 度方向から見た画像上のそれは r180 = -xなので、(∂φ/∂r)(r180) = - (∂φ/∂x)(r180) = - (∂φ/∂x)(- r0)
測定 150602d のデータでこれを確認しました。この E-mail に添付したファイル 150602d.txt（昨日の E-mailの添付ファイルに追記したものです）とファイル 150602d.pdf に並べた画像を御覧下さい。とり急ぎ、
添付ファイル 150602d.txt# head ‐1 150602d/a.log
401,451,100,1,0# set his=(150602d/a.HIS 100 30 401 451)
# his2xp$his axp.bin
rxp.binaxp_180.bin
rxp_180.bin# his2xp_top $his axp_top.bin rxp_top.bin axp_top_180.bin rxp_top_180.bin# a2tiff axp.bin
8 /dev/null
400
400‐0.369277
0.762644# a2tiff axp_180.bin
8 /dev/null
400
400‐0.165889
0.777727# a2tiff axp_top.bin
8 /dev/null
400
400‐1.415264
0.739011

35 / 124

# a2tiff axp_top_180.bin 8 /dev/null
400
400‐0.988178
0.751577# set min=0# set max=0.777727# set bs=($min `echo "( $max ‐ $min )/255" | bc ‐l`)
# a2tiff axp.bin$bs 8 axp.tif
400
400‐0.369277
0.762644# a2tiff axp_180.bin$bs 8 axp_180.tif
400
400‐0.165889
0.777727# a2tiff axp_top.bin$bs 8 axp_top.tif
400
400‐1.415264
0.739011# a2tiff axp_top_180.bin $bs 8 axp_top_180.tif
400 400‐0.988178
0.751577# a2tiff rxp.bin
8 /dev/null
400
400‐3.574706
4.888005# a2tiff rxp_180.bin
8 /dev/null
400
400‐3.594933
5.370930# a2tiff rxp_top.bin
8 /dev/null
400
400‐3.422092
4.795531# a2tiff rxp_top_180.bin 8 /dev/null
400
400‐3.685021
5.129389# set min=‐3.685021# set max=5.370930# set bs=($min `echo "( $max ‐ $min )/255" | bc ‐l`)
# a2tiff rxp.bin$bs 8 rxp.tif
400
400‐3.574706
4.888005# a2tiff rxp_180.bin$bs 8 rxp_180.tif
400
400‐3.594933
5.370930# a2tiff rxp_top.bin$bs 8 rxp_top.tif
400
400‐3.422092
4.795531# a2tiff rxp_top_180.bin $bs 8 rxp_top_180.tif
400
400‐3.685021
5.129389# rm axp*.bin rxp*.bin

36 / 124

添付ファイル 150602d.pdf

Date:
From:
To:
Cc:

Sat, 27 Jun 2015 20:27:00 +0900
Tsukasa NAKANO
Kentaro Uesugi"TAKEUCHI, Akihisa",
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM 用センター値推定スクリプト
うえすぎさま、なかのです。０度と 180 度方向の投影のバイナリデータを使ってサンプル回転軸の位置（rotation center；
いわゆる「センター値」）を推定するための C-shell scripts を書いてみました（C 言語版はこれから書きます）。
吸収の投影画像用の xp2rc.csh と微分位相（屈折）用の xp2rc_dpc.csh の２個があります。これらは SIXM の画像再構成プログラム用の書庫ファイルに入れてあります。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip
なお、今回書いた C-shell scripts は既存の画像処理プログラムを使って処理を行うので、事前にそれらをインストールしておく必要があります（そのために必要な書庫ファイルの名前も以下に併記しておきました）。
a2tiff
前記の his.taz もしくは hiz.zip
trim_gray（２次元画像のトリミングのためのプログラム）

37 / 124

http://www-bl20.spring8.or.jp/~sp8ct/tmp/trim.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/trim.zip
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（FFT による２次元画像の RMSD 計算のためのプログラム）
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rmsd.tar
http://www-bl20.spring8.or.jp/~sp8ct/tmp/rmsd.zip
xp2rc.csh はまず 180 度の画像を左右反転し、それと０度の画像の位置関係を平行移動によって変えながら両者が重なり合っている部分の「投影値の RMSD（root mean square differences）」が最小になっている配置を探し出します。また、xp2rc_dpc.csh では 180 度の画像に対して左右反転と同時に値の正負の反転を行った後に RMSD を計算します。こうして得た RMSD が最小になっている２画像の位置ズレ（０度の画像に付随した座標系における 180 度の画像を反転した画像の原点の座標値）を (dx, dy) とすると、２画像の横画素数が同じなら、センター値、x0 = ( 画像の横画素数 - １ + dx ) / 2となっているはずです。そして、「縦ズレ」の値 dy は SIXM の測定の間にサンプルが「動いた」かどうかを評価するために使えます。
xp2rc.csh と xp2rc_dpc.csh の起動法はまったく同じです。
xp2rc.csh 0.bin 180.bin {y1 y2} Rx Ry {RMSD.txt}
xp2rc_dpc.csh 0.bin 180.bin {y1 y2} Rx Ry {RMSD.txt}
ただし、これらの起動パラメータの意味は以下の通りです。
0.bin と 180.bin
０度と 180 度の投影のデータファイルの名前。ただし、これらの投影画像の横および縦画素数は同じでなければならない。
y1 と y2
RMSD の計算に用いる投影画像の縦（再構成画像のスライス）の範囲。指定を省略するとスライスの範囲すべてが対象となる。撮影した物体像の「影が濃い」部分だけを指定した方が推定の精度が向上するはず。
Rx と RyRMSD 計算の対象とする２画像の x および y 方向の位置ズレの相対的な範囲を指定するためのパラメータ。例えば Rx = 0.5 を指定すると、０度の画像に対して 180 度の画像を画像の横幅の半分だけ左に移動した状態から右に半分移動した状態までのすべての配置に対する RMSD を計算する。なお、xp2rc.csh と等価な処理を行う、HASPET 用の画像再構成プログラムのひとつである
C-shell script "center.csh" は Rx = 0.5、Ry = 0.25 を採用している。
RMSD.txt

38 / 124

位置ズレごとの RMSD の値すべてを書き込むテキストファイルの名前。デバッグ用なので普通は指定しなくても良い。
最終的には、xp2rc.csh と xp2rc_dpc.csh のどちらも以下の２行・６個の値を標準出力に書き出します。
１行目
２行目

センター値の探索範囲 x1～x2 とセンター値の推定値 x0
縦ズレの探索範囲 dy1～dy2 と縦ズレの推定値 dy

なお、dy の値の推定精度（？）は１画素幅ですが、前記の dx との関係式からわかるようにセンター値は半画素幅まで意味があります。そして、x0 = x1 or x0 = x2 や dy = dy1 or dy = dy2 の場合には RMSD が最小の配置が探索範囲の端なので、センター値や縦ズレの推定に失敗している可能性が高いです。
さて、これらの scripts で測定 150603j（および 150602d）のセンター値などを推定してみました。その具体的な手順はこの E-mail に添付した 150603j.txt に記した通りです。このように、測定 150603j の吸収と微分位相の投影画像のどちらに対しても同じセンター値（238）を推定できていますが、縦ズレの推定値は両画像でわずかに異なっています。微分位相の画像の方が noisy なので、それに対する縦ズレの推定値が間違っているのだと思われます。念のため、測定 150602d に関しても同じ処理を行ってみました。150603j.txtの後半がその記録です。こちらも吸収の画像に対しては問題なくセンター値と縦ズレを推定できていますが、微分位相の画像では x0 = x1 となっていて推定に失敗しています。
その原因を知るため、上記の処理で使った０度と 180 度のすべての投影画像を並べて表示してみました。
この E-mail に添付したファイル 150603j.pdf を御覧下さい。
この PDF 上の測定ごとの２枚の吸収の投影画像（axp と axp_180）と微分位相の画像（rxp と rxp_180）
はそれぞれで投影値と画素値の対応関係を揃えてあります。吸収の画像では投影値０が画素値０で、２画像に出現した最大の投影値を画素値 255 にマッピングしています。また、微分位相の画像では２画像に出現した値の絶対値の最大値を P として、-P～P の範囲の微分位相の値を画素値０～255 に割り当てました。
150603j.pdf を見てまず気になるのは 150602d のサンプルの上部にある白金の像です。０度と 180 度の微分位相の画像上でその部分の微分位相の値が正負反転していません。これはなぜだ？ また、150602j.txtの最後に記したように、xp2rc_dpc.csh にその部分を含まないスライスの範囲を指定しても測定 150602d の微分位相の投影画像のセンター値や縦ズレを推定できません。これもなぜだ？
とりあえず以上です。
On Fri, 26 Jun 2015 23:35:41 +0900 Kentaro Uesugi wrote:> 中野さん>> 上杉です>> 0deg と 180deg の相関を取るのに、微分位相像を積分するのはノイズの関係でアレなので、> 相関を取る用のデータを作るのが良いかもしれません。

39 / 124

> 具体的には 180deg のデータの場合には、重心の軸方向を反転させて> 同じ微分位相像になるようにする感じでしょうか。
添付ファイル 150603j.txt# head ‐1 150603j/a.log
501,451,100,1,0# his2xp_top 150603j/a.HIS 100 30 501 451 a r a180 r180
# xp2rc.csh
usage : xp2rc.csh 0.bin 180.bin {y1 y2} Rx Ry {RMSD.txt}# xp2rc.csh a a180 0.5 0.25
124.5 374.5 238
‐100
100
0
# xp2rc_dpc.csh
usage : xp2rc_dpc.csh 0.bin 180.bin {y1 y2} Rx Ry {RMSD.txt}# xp2rc_dpc.csh r r180 0.5 0.25
124.5 374.5 238
‐100
100
3
## 150603j seems to be almost OK as above, but 150602d is NG as follows.
# head ‐1 150602d/a.log
401,451,100,1,0# his2xp_top 150602d/a.HIS 100 30 401 451 a r a180 r180
# xp2rc.csh a a180 0.5 0.25
99.5
299.5 200
‐100
100
0
# xp2rc_dpc.csh r r180 0.5 0.25
99.5
299.5 99.5
‐100
100
‐49
# xp2rc_dpc.csh r r180 135 399 0.5 0.25
99.5
299.5 99.5
‐66
66
‐42

40 / 124

添付ファイル 150603j.pdf

Date:
From:
To:
Cc:

Thu, 09 Jul 2015 17:06:00 +0900
Tsukasa NAKANO
Kentaro Uesugi, "TAKEUCHI, Akihisa""TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM 用センター値推定プログラム
みなさま、なかのです。6/27 の E-mail で紹介した SIXM の再構成画像のセンター値を推定するための C-shell scripts"xp2rc*.csh" とほぼ完全に互換な C 言語プログラムを書きました。
xp2rc
0.bin 180.bin {y1 y2} Rx Ry {RMSD.txt}
xp2rc_dpc
0.bin 180.bin {y1 y2} Rx Ry {RMSD.txt}
新しいプログラム xp2rc と xp2rc_dpc のソースファイルと 64 bit Windows 用実行ファイルなどを以前と同じ
SIXM 用の書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zipに入れておきましたので、どうぞ御利用下さい。
なお、これらのプログラムのテスト中に「SIXM の測定データファイル "*.HIS"」から「吸収と微分位相の投影のバイナリデータファイル」を作成するプログラム his2xp と his2xp_top の不具合に気づきました。これらはX 線強度の測定値が０の画素に対して「NaN（not a number）」を出力します。以前の C-shell scripts はこのデータを TIFF 画像に変換して処理していたので問題なかったのですが、新しいプログラムはバイナリデータを直接処理するので NaN があると正しいセンター値を推定できません。そこで、NaN を出力しないように41 / 124

his2xp と his2xp_top のコードを書き換えました。修正したソースファイルや 64 bit Windows 用実行ファイルを
上記の書庫ファイルのものと差し替えておきました。
とりあえず以上です。
Date:
From:
To:
Cc:

Mon, 13 Jul 2015 20:02:32 +0900
Tsukasa NAKANO
Kentaro Uesugi"TAKEUCHI, Akihisa","TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: センター値が決まらない場合
みなさま、なかのです。6/27 の E-mail でも報告したように、xp2rc_dpc_csh に「いつも通り」のパラメータの値（具体的には、Rx = 0.5 と Ry = 0.25）を指定すると測定 150602d の位相トモグラフィー用のセンター値の推定に失敗します。これは xp2rc_dpc でも同様です。その原因を知るため xp2rc* にファイル名を指定して「０度と 180 度の投影画像の位置ズレ (dx, dy) ごとの RMSD のデータ」を出力し、それを画像化してみました。
プログラム xyv2t を使えば xp2rc* が起動時に指定されたファイルに書き込んだ RMSD(dx,dy) のテキストデータを TIFF 形式画像に変換できます。xyv2t の使用法や RMSD(dx,dy) のテキストデータの構成などについては添付したファイル rmsd.pdf を御覧下さい。
この E-mail に添付した xp2rc.pdf を御覧下さい。そこに並べた横４×縦３の画像のうちの上２段のものは6/27 の E-mail の添付ファイル 150603j.pdf のものと同様な測定 150602d と 150603j の吸収（AXP）もしくは微分位相（RXP）の０度と 180 度の投影画像です。ただし、２段目の 180 度の画像は左右反転し、さらに微分位相の画像では投影値の符号も反転してあります。そして、３段目が RMSD(dx, dy) の「マップ」で、RMSDの値の高低をマジェンタ（高）→赤→黄→緑→シアン→青（低）の順に変換する色のグラデーションで表現しています。
０度と 180 度の投影画像の横および縦画素数を Nx と Ny とすると、xp2rc.pdf の RMSD のマップ領域の座標値の範囲は以下の通りです。
-0.5 ≦ dx / Nx ≦ 0.5-0.25 ≦ dy / Ny ≦ 0.25
ただし、上式の定数の絶対値 0.5 と 0.25 は xp2rc* の起動パラメータ Rx と Ry の値に対応しています。
さて、6/27 の E-mail に記したように、センター値 x0 の推定値は RMSD が最小となっている dx の値 Dx を
以下の式に代入すれば計算できます。
センター値の推定値、x0 = ( Nx - 1 + Dx ) / 2

42 / 124

Rx = 0.5 と Ry = 0.25 を指定した場合に xp2rc* が出力した x0 の 値（および、その時の「縦ズレ」の値 Dy）
を RMSD のマップの下に赤で記しておきました。測定 150602d の RXP では x0 に対応する Dx はマップの左端（Dx/Nx = -Rx）の値であり、これは RMSD が最小の点とは言えないので x0 の値の後に "#" を付けてあります。そして、この場合は投影画像に含まれている「異常なノイズ」のために RMSD（のマップ）に「ローカルミニマ」が発生し、それが邪魔で「正しいセンター値」を推定できなかったようです。
この異常なノイズは０度と 180 度の画像でサンプル上部にある白金の微分位相の値が食い違っていることと関連しているようですが、真相は不明です。それはさておき、このローカルミニマを回避したセンター値の推定は容易です。つまり、サンプル回転軸は投影画像のほぼ中央にあると仮定できるので、xp2rc* のパラ
メータ Rx にもっと小さな値を指定してやれば良い。xp2rc* に Rx = 0.25 を指定して得た x0 と Dy の値を
xp2rc.pdf の最下段に青で記しておきました。
とりあえず以上です。
添付ファイル rmsd.pdf
このファイルの内容は SIXM の処理とは直接関係しないので、ここに掲載しませんでした。
添付ファイル xp2rc.pdf

43 / 124

Date:
From:
To:
Cc:

Thu, 16 Jul 2015 11:05:16 +0900
Tsukasa NAKANOMIYAKE, "中村隆太"
Akira TSUCHIYAMA, "MATSUNO Junya", Aki Takigawa, Takashi Sakurama,
Kentaro Uesugi, "TAKEUCHI, Akihisa"
Subject: SP8_BL47XU_SIXM_150603c_a.HIS
中村さま、
三宅さま、
GSJ/AIST のなかのです。今頃になって何ですが、上杉くんが SPring-8 にある計算機 vrm にアップロードしてくれていた６月の SIXM 実験で得た測定データを処理していたら、測定 150603c （石英サンプルを撮影；
ログノートなし）のファイル a.HIS が壊れていることに気がつきました。その概要はこの E-mail に添付したファ
イル 150603c.txt に記した通りです。端的に言うと、測定データファイル「150603c/a.HIS」の一番最後（271151 番目）の測定画像を読むことができません。
ファイル 150603c.txt では測定 150603c の前後の、ほぼ同じ撮影条件の測定 150603[b,d] の a.HIS もチェックしています。150506c/a.HIS に格納されている測定画像はこれらのものよりも１個少ない。
vrm の上の 150603c/a.HIS はデータ転送時に壊れたのかもしれません。そこで、お手数ですが、そちらにある 150603c/a.HIS を調べてもらえませんか？ 具体的には Linux 機で以下を入力してそのファイルのチェックサムを調べて下さい：
cksum

150603c/a.HIS

これを vrm の上で行うと以下の出力を得ました（150603c.txt の２行目）：
1103610525 27772972435

150603c/a.HIS

そちらで得た出力の最初と２番目の数値（CRC の値とファイルのバイト数）が上と同じなら 150603c/a.HIS は測定時に壊れたと思われるので、どうしようもありません。しかし、そうでない場合はそちらのファイルは大丈夫かもしれないので、それをぼくにコピーさせて下さい。よろしくお願いします。とり急ぎ、
添付ファイル 150603c.txt
このファイルの内容に間違いがありました。訂正版が次の E-mail に添付されています。
Date:
From:
To:
Cc:

Thu, 16 Jul 2015 15:51:06 +0900
Tsukasa NAKANO <tsukasa.nakano@aist.go.jp>
MIYAKE, "中村隆太"
Akira TSUCHIYAMA, "MATSUNO Junya", Aki Takigawa, Takashi Sakurama,
Kentaro Uesugi, "TAKEUCHI, Akihisa"
Subject: SP8_BL47XU_SIXM_150603c_a.HIS_訂正44 / 124

みなさま、なかのです。先程の E-mail に間違ったことを書きました。計算機 vrm の上にあるファイル 150603c/a.HIS の測定画像で壊れているのは「一番最後」のものだけではなく、最後から 101 枚（271051 ～ 171151 番目）の測定画像でした。そして、添付ファイル 150603c.txt の中程に記した引き算が間違っていた。
→ それを修正したファイルをこの E-mail に添付しておきました。
E-mail 本文の２箇所が間違っていた：[1] 一番最後（271151 番目）の測定画像を読むことができません。
→ 最後の 101 枚の測定画像を読むことができません。
[2] 150506c/a.HIS に格納されている測定画像はこれらのものよりも１個少ない。
→ 101 個少ない。
いずれにせよ vrm の 150603c/a.HIS は壊れているので、京都大学に持ち帰ったファイルのチェックサムを調べて下さい。よろしくお願いします。とり急ぎ、
添付ファイル 150603c.txt% cksum 150603b/a.HIS% cksum 150603c/a.HIS% cksum 150603d/a.HIS
#

‐>
‐>
‐>

1755732483 27783218835 150603b/a.HIS
1103610525 27772972435 150603c/a.HIS
3027004535 27783218832 150603d/a.HIS[b],[c],[d]% head ‐1 150603b/a.log ‐>% head ‐1 150603c/a.log ‐>% head ‐1 150603d/a.log ‐>601,451,100,1,0
601,451,100,1,0
601,451,100,1,0% expr 601 ¥* 451 + 100 ‐>271151# [e]% check_his 150603b/a.HIS | cat ‐n | tail ‐1‐> 27783218835# $1 == [b]271151 128400
2# $1 == [e]% check_his 150603c/a.HIS | cat ‐n | tail ‐1‐> 27772869971# $1 != [c] ($1 == [c] ‐ 102464)
271050 128
400
2# $1 != [e] ($1 == [e] ‐ 101)% check_his 150603d/a.HIS | cat ‐n | tail ‐1‐> 27783218832# $1 == [d]271151 128400
2# $1 == [e]
Date:
From:

Fri, 17 Jul 2015 16:34:28 +0900
Tsukasa NAKANO

# [b]
# [c]
# [d]

45 / 124

To:
Cc:

Takashi Sakurama
Akira TSUCHIYAMA, MATSUNO Junya, Aki Takigawa, MIYAKE, "中村隆太",
Kentaro Uesugi, "TAKEUCHI, Akihisa
Subject: Re: SP8_BL47XU_SIXM_150603c_a.HIS_訂正
さくらまさま、なかのです。連絡ありがとうございます。お手数をおかけしました。150603c/a.HIS は測定時に壊れていた
ようで、どうしようもないですね。とりあえずの対処策（次の E-mail に書きます）があるので何とかなります
が、気分が悪いので、今後は測定直後にきちんとチェックすべきだと思います。とり急ぎ、On Fri, 17 Jul 2015 16:10:05 +0900 (JST) Takashi Sakurama wrote:> 中野さま>> 櫻間です。
> 150603c/a.HIS のチェックサムの結果を添付します。
> 中野さんが vrm で得たのと同じ数値が得られました。
Date:Thu, 16 Jul 2015 11:20:22 +0900
From: Tsukasa NAKANO
To:
Kentaro Uesugi, "TAKEUCHI, Akihisa"
Cc:
Akira TSUCHIYAMA, MIYAKE, "MATSUNO Junya", Aki Takigawa, Takashi Sakurama, "中村隆太"
Subject SP8_BL47XU_SIXM_幾何学量
上杉さま、
竹内さま、
GSJ/AIST のなかのです。今頃になって何ですが、６月の実験で使った SIXM の装置の以下の２個の「幾何学量」の値を教えていただけますか？
サンプルと検出器の間の距離、L測定画像の横方向（屈折方向；画素数はすべて 128）の画素の辺長、Δr
7/2 の E-mail に記した位相画像の上の値の「スケールファクタ」
Dr = λ× L / ( 2 × π× Δr )の計算に使いたいと思っています。とり急ぎ、
Date:
From:
To:
Cc:

Thu, 16 Jul 2015 14:09:20 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa”
Kentaro Uesugi,
Akira TSUCHIYAMA, MIYAKE, MATSUNO Junya, Aki Takigawa, Takashi Sakurama, 中村隆太
Subject: Re: SP8_BL47XU_SIXM_幾何学量46 / 124

たけうちさま、なかのです。E-mail ありがとうございます。すみません。ぼくは「Hideki」を持ち帰っていないので「走査位相顕微鏡の再構成」では確認できませんでしたが、測定ごとのディレクトリの下にあったテキストファイル a.parの内容から確認できました。このファイルの９～12 行目の値の意味は以下の通りですか？
９
10
11
12

スキャンのステップ幅、Dr（単位は nm）
レイヤーの厚さ、Dz（nm）。どの測定でも 108.9 だった。
L（mm）。どの測定でも 6160 だった。
Δr / 16（um）。どの測定でも 6.4935 だった。

ただし、９行目の Dr の値には実験時のログノートに記録されている値やキャプチャー画像 a.bmp から読み取った値（から計算した値）とは違っているものがありました。具体的には、以下の測定の Dr が食い違っています。
測定
150605d150605f
150605i
150605ja.par の値
150
1000
50
100

ログノートや a.bmp の値
100 = 25 × 4
150 = 25 × 6
100 = 25 × 4
150 = 25 × 6

これらのどちらを信じたら良いのでしょうか？ また、実験のログノートには Dz = 108 と記録されていますが、その真の値は 108.9 ですか？ とり急ぎ、
On Thu, 16 Jul 2015 12:54:30 +0900 "TAKEUCHI, Akihisa" wrote:> 中野様>> 竹内です。
>> L = 6160 mm> Δr = 103.896 um (6.4935 um X 16 binning)> です。
> 或いは、Hideki => 「走査位相顕微鏡の再構成」で、> 目的のデータが入っているディレクトリから.par ファイルを選んでもらうと> これらのパラメータが「測定・計算パラメータ」タグのところに表示されると思います。
Date:
From:
To:
Cc:

Thu, 16 Jul 2015 15:30:42 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa"
Kentaro Uesugi ,Akira TSUCHIYAMA, MIYAKE, MATSUNO Junya, Aki Takigawa, Takashi Sakurama,中村隆太

Subject: Re: SP8_BL47XU_SIXM_幾何学量47 / 124

たけうちさま、なかのです。E-mail ありがとうございます。それでは、Dr としてログノートに記録されている値を、また、Dzには a.par の値 108.9 を使うことにします。お手数をおかけしました。とり急ぎ、
On Thu, 16 Jul 2015 14:46:36 +0900 "TAKEUCHI, Akihisa" wrote:
> 中野様
>
> はい。
> 9-12 行目の値は全て中野さんの推測で合っています。
> １２に関しては、37 行目に binning の値（16)が出てきます。
>> Dr に関しては、a.bmp から読み取れる値が正解です。
> a.bmp の値を基に測定するのに対して、> .par の方は、Hideki で再構成する際テキストボックスに入力されている> 値が記録されるので。
> 再構成する人がわざわざテキストボックスの値を> 書き換えない限りその前の再構成の際の値が記録されることになります。
>> Dz に関しては、ログノートの値よりも.par の 108.9 を信用してください。
Date:
From:
To:
Cc

Fri, 17 Jul 2015 18:51:58 +0900
Tsukasa NAKANO"TSUCHIYAMA, Akira"Kentaro Uesugi, "TAKEUCHI, Akihisa",
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SP8_BL47XU_SIXM_150602
みなさま、
GSJ/AIST のなかのです。これまでに紹介した書庫ファイル
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zipに入れた計算機プログラムを使って、６月に SPring-8 BL47XU の SIXM 実験で得た測定データから吸収お
よび位相トモグラフィーの画像再構成を行いました。測定データの量が膨大だった（約 1.7 TB）ので、今回はすべての処理を計算機 vrm のディレクトリ/media/disk/tsukasa/150602/で実行しました。
上杉くんが vrm の/media/disk/bl47xu/2015.06.02_tsuchiyama/の下に６月の SIXM 実験の測定データ
すべてをアップロードしてくれていたので、ぼくはそれら（のうちの必要なものだけ）を自分のディレクトリ
にハードリンクして処理を行いました。
この E-mail に添付したテキストファイル vrm.txt にぼくが使用・作成したファイルとディレクトリのリストが記されています。以下ではこのリストに沿って、ぼくが今回行った処理の内容を説明します。

48 / 124

ディレクトリ /media/disk/tsukasa/150602/for_download/ の下に「測定データファイルと vrm でのみ有効なプログラムの実行ファイル」 以外のもの（総量は約 19 GB）をハードリンクしておきました。ぼくの処
理結果が必要ならこれだけを丸ごとダウンロードすれば良いです。
(0) ディレクトリ 150602d.org/と 150603j.org/の下には実験時にコピーした測定 150602d と 150603j のデータ
ファイル一式が入っていますが、今回の処理ではこれらは使いません。これらのうちの必要なものは別
途ハードリンクしました。
(1) ファイル sample.xls と log.pdf はそれぞれ三宅さんと中村くんが実験直後に送ってくれた「測定ログ」で、もとのものから名前を付け替えました。そして、これらをもとにして画像再構成が可能だと思えた 29 の測定それぞれのデータと撮影サンプルの概要をファイル sixm.csv と sixm.ods に記載しました。これらのファイルの内容は概ね同じです。後者は OpenDocument Spreadsheet format のファイルで、MS-Excel
でも閲覧・編集などは可能なはずです。
(2) 処理に先立ち、上杉くんがアップロードしてくれた前記の測定データファイルのうちの a.HIS、a.log とa.par の３個だけを測定ごとのディレクトリにハードリンクしました。ただし、a.par は竹内さんのプログラム
HIDEKI が使っているパラメータファイルで、彼に質問した 7/16 の E-mails にも書いてあるように、位相トモグラフィーの画像再構成に必要な装置の幾何学量の値が記されています。
(3) ここからが実際の画像再構成処理です。まず、プログラム his2xp と a2tiff を利用しつつ処理を行う
C-shell script "run.xp" を起動して、それぞれの測定の０度方向の吸収および屈折の投影画像のTIFF 形式ファイル axp.tif および rxp.tif を測定データファイル a.HIS から作りました。これらは被写体の粗い観察のためだけに使う画像なので、画素値のビット数を８とし、また、his2xp の起動パラメータ I0sの値（ダイレクトビーム I0 を測定している投影画像の左端と右端の画素列の横画素数）もとりあえず10 としました。
(4) これらの画像を目視観察して測定ごとに以下の２個の値（画像上の座標値の差）を読み取り、それらをテキストファイル i0s.tbl に書き込みました。
I0s_bs（bs == both side）
I0 を測定していると見なせる投影画像の左右両端に位置する空気の部分の縦画素列のうちで横幅が狭い方の横画素数。ただし、サンプル回転軸が投影画像の横中央にない場合に備えて、画像から読み取った値から 20（この値に深い意味なし）を減じた値を I0s_bs とした。
I0s_top
I0 を測定していると見なせる投影画像の上端に位置する空気の部分の横画素列の縦画素数。
４個の測定で I0s_top は０だった（つまり、４測定の投影画像でサンプル上部に空気の部分がなかった）。
なお、ファイル i0s.tbl の I0s_bs の値には正負の符号を付けてあります。"+" は投影画像左端の空気の
縦列の横幅が右端のものよりも狭いことを、また "-" はその逆を意味しています（実際にはこの区別を活用しませんでしたが）。
(5) プログラム his2xp、his2xp_top、xp2rc および xp2rc_dpc を利用しつつ処理を行う C-shell script "run.rc"を起動して、各測定の吸収および位相トモグラフィー用のセンター値のそれぞれを推定しました。ただし、

49 / 124

その際に先のファイル i0s.tbl に記されている値に応じて０度と 180 度の投影画像（のバイナリデータファイル）を作成する手法（プログラム）を変えました：
I0s_top が０の場合
起動パラメータ I0s に I0s_bs の値（の絶対値）を指定してプログラム his2xp で０度と 180 度の吸収および屈折の投影画像を作成する。
そうでない場合
I0s に I0s_top の値を指定して his2xp_top で投影画像を作成する。
最終的には、run.rc はこのような投影画像を入力した xp2rc と xp2rc_dpc の２行づつの出力を測定ごとに１行にまとめて標準出力に書き出します。それをテキストファイル rc.log に書き込みました。
(6) run.rc の実行中に（7/16 の E-mail で紹介した 150603c.txt に記した）測定 150603c のファイル a.HIS の不備に気づきました。そのため、run.rc ではセンター値を推定できなかったので、150603c/a.HIS で壊れている 180 度の投影画像をその１枚前の 179.6 度の画像で代用する C-shell script "run.rc_150603c"
を実行し、その結果を rc_150603c.log に書き込みました。
(7) さて、run.rc では投影画像のサンプル上部に空気の部分がない場合に i0s.tbl の I0s_bs の値をそのまま
his2xp のパラメータ I0s として使います。しかし、his2xp は投影画像の左右両端の I0 の値を線形補間するので、あまりに大きな I0s の値では線形補間の意味がなくなります。そこで、上部に空気がない場合に his2xp に I0s_bs の値ではなく一定値 10（この値に根拠なし）を指定して実行する run.rc_bs_10 を走らせて、その結果を rc_bs_10.log に書き込みました。
(8) C-shell script "list.rc" は上記の rc.log、rc_150603c.log と rc_bs.log の内容を整理して標準出力に書き出します。それを書き込んだファイル rc.txt のセンター値や「縦ズレ」の値を検討して、吸収および位相
トモグラフィーの両方の画像再構成で指定する測定ごとの「単一のセンター値」を決定しました。
(9) 以上で SIXM の画像再構成の前処理は終わりです。後は画像再構成プログラム xp2tg と xp2tg_dpc を実行するだけですが、後者に指定する起動パラメータ Dr の値をあらかじめ計算しておく必要があります。
7/2 の E-mail で論じたように、以下の Dr の値を xp2tg_dpc に指定してやれば位相トモグラフィーの再構成画像上の CT 値は位相シフトφの推定値となります。
Dr = Dr_P = ( L / Δr ) × λ/ ( 2 π)
ただし、
L サンプルと検出器の間の距離
Δr 測定画像の横方向（屈折方向；画素数 128）の画素の辺長λ 測定に使用した 8 keV の X 線の波長
しかし、竹内さんは屈折率 n から「n = 1 - δ」で計算できるδの値（refractive index decrement）を再構成しているようでした。そこでぼくも 10^-6 を単位とするδの値を再構成することにします。そして、この場合には、
δ= φ× λ/ ( 2 π)なので、結局、xp2tg_dpc の Dr として以下の値を指定すれば良いはずです。
Dr = Dr_D = L / Δr
このようなφ用の Dr（Dr_P）とδ用の Dr （Dr_D）の値の具体的な計算に関することをテキストファイル
dr_pd.txt に書き込んでおきました。

50 / 124(10)吸収と位相トモグラフィーの画像再構成に必要なパラメータをテキストファイル tg.tbl に書き込みました。
これを読み込みつつプログラム xp2tg と xp2tg_dpc のそれぞれで各測定の吸収と位相トモグラフィーの
３次元 CT 画像を再構成する C-shell script が run.atg（吸収）と run.dtg（位相）です。これらはまずディ
レクトリ「測定番号/[a,d]tg」を作成します。その後、xp2tg や xp2tg_dpc を起動して指定されたスレッド数（指定をしなかった場合は８スレッド）のマルチスレッド処理で画像再構成を行い、それらが標準出力に書き出したスライスごとの CT 値の最小値と最大値をテキストファイル「測定番号/[a,d]tg.log」にリダイ
レクトします。さらに、run.[a,d]tg は測定ごとの再構成画像全体の CT 値の最小値と最大値を標準出力に書き出します（ぼくはそれらをテキストファイル [a,d]tg.log にリダイレクトしました）。
(11)C-shell script "list.tg"は測定それぞれの再構成画像の以下の情報を標準出力に書き出します（ぼくはそれをファイル tg.txt にリダイレクトしました）
set
測定番号
Nr
スライス画像（正方形）の各辺沿いの画素数
Dr_nm スライス画像の画素の各辺の長さ（単位は nm）
Nz
スライスの枚数（今回の SIXM 実験の測定ではすべて 400）
Dz_nm スライス厚（単位は nm；今回の測定はすべて 108.9 nm）
A_min._1/cm と A_max._1/cm
吸収トモグラフィーの CT 値の最小値と最大値（単位は 1/cm）
D_min._1E-6 と D_max._1E-6
位相トモグラフィーの CT 値の最小値と最大値（単位は 10^-6）(12)こうして得た再構成画像に関して以下のことにご注意下さい。
[1] 吸収トモグラフィーの CT 値の単位は 1/cm です。
[2] 位相トモグラフィーの CT 値の単位は 10^-6 です。
[3] これらの CT 値は３次元画像に出現した最小値と最大値で正規化した 16 ビット画素値に変換されています（つまり、従来は別途変換していた word 画像を再構成の時点で作成しました）。
[4] これらの CT 値と画素値の対応関係の変更や画像上の CT 値の出現頻度分布（ヒストグラム）の調査はプログラム tg2tg で行えます。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/sp8ct.zip(13)今回の画像再構成に要した正確な処理時間は不明です。まず、7/16 の 19 時頃に run.atg を起動しました。これは８スレッドのマルチスレッド処理を行っているはずでしたが、CPU の稼動率は 200 % 程度でした。そこで、run.atg が３個目の測定の画像再構成を始めた頃に run.dtg を起動しました。このような画像再構成処理の同時実行は CPU 的には問題ないと思いますが、すべてのデータファイルを HDD に置いていたので、そのアクセスがかなり遅くなっていたかもしれません。いずれにせよ、7/16 の 23 時頃にすべての処理が終了したようです（つまり、今回の SIXM の画像再構成の総処理時間は約４時間）。
長い E-mail になりました。とりあえず以上です。
添付ファイル vrm.txt
次の E-mail に内容を追加した同名のファイルを添付したので、こちらのファイルは掲載しません。

51 / 124

Date:
From:
To:
Cc:

Sat, 01 Aug 2015 17:53:19 +0900
Tsukasa NAKANO"TSUCHIYAMA, Akira"Kentaro Uesugi, "TAKEUCHI, Akihisa",
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: SIXM_150602_オマケ
つちやまさま、
GSJ/AIST のなかのです。7/17 の E-mail で紹介した６月の SPring-8 BL47XU での SIXM 実験で撮影した
LAC と RID（refractive index decrement；１から屈折率の値を引いた値）の再構成画像に対して以下の処理を行いました。
[1] LAC と RID の推定値を 16 ビット画素値として格納した３次元再構成画像（tg 画像）のそれぞれの画素値ヒストグラムデータの取得[2] LAC と RID のそれぞれの８ビット画素値の再構成画像（byte 画像）の作成と、それらの画素ヒスト
グラムデータの取得[3] LAC と RID の byte 画像それぞれの代表的な xy、xz および zy 断面の画像を並べた browse 画像の作成[4] LAC の値を横軸、RID の値を縦軸とした tg と byte 画像のペアのそれぞれの「２次元ヒストグラム」の作成
今回の処理は手元にある計算機で行いましたが、それで使用・作成した C-shell scripts を含むほぼすべてファイルを以前に紹介した SPring-8 の計算機 vrm のディレクトリ/media/disk/tsukasa/150602/にコピーしておきました。この E-mail に添付したファイル vrm.txt にそれらの内訳を追記してあります。これに沿ってぼくが今回行った処理の内容を説明します。
vrm のディレクトリ/media/disk/tsukasa/150602/for_download/の下に vrm 用の実行ファイル si_pair_hg
以外の今回の処理で使用・作成したファイルをハードリンクしておきました。なお、このディレクトリの総容量は以前の値から 5 GB 増えた 24 GB です。
(0) まず、プログラム tg2tg を使って処理を行う C-shell script "run.tg_csv" を走らせて、以前に再構成した
29 測定のそれぞれのディレクトリの下に LAC と RID の tg 画像（ディレクトリ atg/ と dtg/）の画素値ヒス
トグラムのデータファイル atg.csv と dtg.csv を作成しました。なお、これらは HASPET 用の C-shell script"hg_*.csh" でグラフとして表示することができます。
(1) プログラム tg2tg によって処理を行う C-shell script "run.byte" を実行して、tg 画像から LAC と RID の
byte 画像（abyte/ と dbyte/）を測定ごとのディレクトリの下に作成しました。ただし、その際に、正の値を持つ tg 画像の上の画素の 99.9 % のものの値だけを byte 画像が表現するように LAC や RID の値の下限と上限（正確には下限値０と画素値１階調あたりの値の増分）を指定しました。なお、これらのbyte 画像それぞれに用いた増分の値はテキストファイル abyte.log と dbyte.log に記録してありますが、うっかりしていてそれらの値を列記したファイルを作り忘れました。さらに、byte 画像それぞれの画素値
ヒストグラムのデータを abyte.csv と dbyte.csv に書き込みました。

52 / 124(2) 測定それぞれの画像の代表的な断面図を切り出すため、byte 画像の上の LAC と RID の値の空間分布の「重心」の位置を楕円体近似用のプログラム si_ofwd を使って計算しました。C-shell script "test.coa"
で得た測定ごとの LAC と RID の重心の座標値をテキストファイル coa.log に列記しました。さらに、このスクリプトは xy 平面内の重心の位置を赤線で示した LAC と RID の byte 画像の重心の位置を通るスラ
イスを並べた画像のファイルをディレクトリ coa/の下に書き込むようになっていますが、これは vrm にはコピーしていません。
うっかりしていました。coa は「center of absorption」の意味で、RID の値の重心には不適切な名称でした。
(3) ディレクトリ coa/ の画像の観察などから測定それぞれの byte 画像の代表的な断面が通る点の座標値を以下のように決め、それらの情報をテキストファイル xyz.tbl に列記しました。
測定 150603c と 150605a
再構成に失敗した画像もしくは「極小」の画像なので、断面の browse 画像を作成しない。
150605[b,c,d]
物体像がないなど LAC と RID の重心は断面図用には不適当だったので、byte 画像の目視観察で適切と思われる点の座標値を読み取った。
150603d、150603f、150605j
RID の重心を代表的な断面が通る点とする。
上記以外の 21 測定
LAC の重心を代表的な断面が通る点とする。
(4) xyz.tbl に記した点を通る測定それぞれの LAC と RID の byte 画像の xy、xz および zy 断面を切り出し、それらを並べた browse 画像の PostScript（PS）データを出力する C-shell scripts を実行しました。
print.xyz
起動パラメータとして指定した測定（複数も可）の LAC と RID の３枚の断面図を並べた browse
画像の PS データを出力する C-shell script。
print.xyz_all
print.xyz を順次起動して xyz.tbl に記されている 27 測定の browse 画像の PS データを取得し、それらを A4 の各ページに３個づつ並べた合計９ページのデータとして出力する。
print.xyz_all で出力した PS データを変換して PDF ファイル byte_xyz.pdf を作成しました。
(5) DEA 用のプログラム si_pair_hg などを流用して LAC の値を横軸、RID の値を縦軸とする tg および byte
画像のペアの２次元ヒストグラムを作成しました。C-shell script "run.tg_hg" と "run.byte_hg" は測定それぞれの tg や byte 画像のペアから２次元ヒストグラムのデータファイル "*_hg.txt" を作成します。
そして、ヒストグラムの bin に落ちたペアの個数の対数の値を８ビット画素値とした画像 "*_hg.tif" に格納した後にプログラム t_pcs によってそれらのすべての「ピーク・クラスタ」を検出し、その情報をテキスト
ファイル"*_hg.log" に書き込みます（ただし、このファイルは作っただけで活用していませんが）。

53 / 124

(6) 以下の C-shell scripts を用いて DEA の場合と同様な２次元ヒストグラムなどの図を描きました。
draw.tg_hg と draw.byte_hg
これらはそれぞれ起動パラメータとして指定した測定（複数も可）の２次元ヒストグラムの図のPS データを出力する C-shell script
print.hg
draw.tg_hg と draw.byte_hg を順次起動して xyz.tbl に記されている 27 測定の２次元ヒストグラ
ムの図の PS データを取得し、それらを A4 の各ページに３測定分づつ並べた合計９ページのデータとして出力する。
ただし、"draw.*_hg" はファイル cm_rgb_b.txt に記されている成分値の色で２次元ヒストグラムの高さを表現します。また、プログラム sr_map でヒストグラムの起伏に応じた陰影付けを行っていますが、これをそちらの Linux 計算機にインストールしたかどうか定かではありません。さらに、以前にインストールに失敗した図形編集用（お絵かき）プログラム「tgif」（これはぼくが書いたものではないプログラムで、テキストファイル hg.obj はそれ用の「テンプレート」です）を使うので、これらのスクリプトはそちらの計算機では実行できない可能性が大です。
(7) いずれにせよ、こちらで print.hg を実行して得た PS データを変換した PDF ファイル tg+byte_hg.pdf をコピーしておきましたので、前記の byte_xyz.pdf と一緒に御覧下さい。
とりあえず以上です。なお、ぼくは来週は夏休み（GSJ 全員の一斉休業）です。とり急ぎ、
添付ファイル vrm.txt
tsukasa@vrm.spring8.or.jp:/media/disk/tsukasa/150602/
150602d.org/, 150603j.org/
sample.xlsx
log.pdf
sixm.csv, sixm.ods++

150603c/, 150605a/a.HIS, a.log, a.paraxp.tif, rxp.tif
atg/*.tif, atg.log
dtg/*.tif, dtg.log
atg.csv, abyte.log, abyte/*.tif, abyte.csv
dtg.csv, dbyte.log, dbyte/*.tif, dbyte.csv
150602d/, 150602e/, 150602f/, 150602g/, 150603a/,
150603b/, 150603d/, 150603e/, 150603f/, 150603g/,
150603i/,
150603j/, 150603k/,

54 / 124++++

150604b/,150604c/,150604d/,150604e/, 150604f/, 150604g/,
150605b/, 150605c/,
150605d/, 150605e/, 150605f/, 150605h/, 150605i/,
150605j/a.HIS, a.log, a.paraxp.tif, rxp.tif
atg/*.tif, atg.log
dtg/*.tif, dtg.log
atg.csv, abyte.log, abyte/*.tif, abyte.csv
dtg.csv, dbyte.log, dbyte/*.tif, dbyte.csv
tg_hg.txt, tg_hg.tif, tg_hg.log
byte_hg.txt, byte_hg.tif, byte_hg.log
his2xp, a2tiff
run.xpi0s.tblhis2xp_top, xp2rc, xp2rc_dpc
run.rc, rc.log
150603c.txt
run.rc_150603c, rc_150603c.log
run.rc_bs_10, rc_bs_10.log
list.rc, rc.txt
dr_pd.txt
xp2tg, xp2tg_dpc
tg.tbl, run.atg, run.dtg, atg.log, dtg.log
list.tg, tg.txt++

run.tg_csv
run.byte++

test.coa, coa.log
xyz.tbl+print.xyz, print.xyz_all, byte_xyz.pdf

+
+
+
+

si_pair_hg
run.tg_hg, run.byte_hg
cm_rgb_b.txt, hg.obj
draw.tg_hg, draw.byte_hg, print.hg, tg+byte_hg.pdf
for_download/55 / 124

Date:
From:
To:
Cc:

Mon, 10 Aug 2015 14:54:53 +0900
Tsukasa NAKANO
Kentaro Uesugi, "TAKEUCHI, Akihisa""TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: SIXM_empty_HIS
うえすぎさま、
たけうちさま、
GSJ/AIST のなかのです。SIXM で「サンプルなし」の状態を測定したデータファイル（a.HIS）はありません
か？ CT 用のものが望ましいですが、そうでなくてもかまいません。とり急ぎ、
Date:
From:
To:
Cc:

Tue, 18 Aug 2015 18:55:19 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa"
Kentaro Uesugi,"TSUCHIYAMA, Akira", MIYAKE, "Aki Takigawa", Junya Matsuno, 中村隆太, Takashi Sakurama

Subject: Re: SIXM_empty_HIS
たけうちさま、なかのです。すみません。先程の E-mail に大きな誤りを見つけました（Pm と書いたつもりが Pg と書いてしまった；添付ファイル msd+gp.pdf では Pg は global peak を表しています）。先程の E-mail を破棄してこちらを御覧下さい。
たけうちさま、
GSJ/AIST のなかのです。応答が遅くなってすみません。そちらにサンプルなしの HIS のデータがないのは残念です。ぼくは以下で説明する６月の実験で得たサンプルありの HIS のデータから抽出した X 線プロファ
イルの解析で先に進めなくなったので、そのものズバリの SIXM の入射 X 線のプロファイルの様相を知りた
いと思って empty な HIS データがあるか否かをお尋ねしました。
On Tue, 11 Aug 2015 07:39:45 +0900 "TAKEUCHI, Akihisa" wrote:> 中野様>> 早い夏休みに入っており、返事が遅くなりました。
>> empty な his データは残念ながら持ち合わせてないです。
> 急ぎの場合は、視野上部にサンプルが無いデータ（三宅さんのサンプルは殆どそうか> と思いますが）を切り出してどうにかする、というのは可能でしょうか(0) ぼくがしたいのは SIXM の吸収と屈折（位相シフト）の両方の３次元再構成画像からのノイズ除去です。
ただし、以下の理由で、それを竹内さんが使われている Martino の論文（Opt. Lett., 38, 4813, 2013）の方法で行いたくない：

56 / 124[a] そもそも、この方法では吸収の画像のノイズ除去はできない。
[b] 位相回復（線積分）は測定した微分位相のデータのノイズ除去の後に Pfeiffer の論文（Phy. Rev.
Lett., 98, 108105, 2007）の手法で画像再構成と同時に行いたい。
[c] Martino の論文では理解に苦しむことを仮定している（これについては別便で話します）。
(1) SIXM の２種類の画像のノイズ除去のためにはそれに用いた測定データに含まれているノイズの様相を知ることが第一だと思われます。ただし、SIXM では吸収トモグラフィーに使う投影（もしくは透過率）のデータと屈折トモグラフィー用の微分位相の両方のデータを同じ「屈折プロファイル」から算出しているので、結局、実際の「屈折プロファイル」がどうなっているかを調べました。、(2) ここで言っている「屈折プロファイル」とは HIS 形式ファイルに入っている多数の画像が表している X 線
強度の空間分布 I[p,x,y,v] のことです。ただし、v はサンプル回転のステップ、y は鉛直方向の座標値で、スキャンのステップが x の時に水平方向に拡大された屈折光の検出位置が p = 0～ n です（６月の実験では n = 127）。そして、以前に紹介したプログラム his2xp* では以下の式で SIXM の２種類の画像の再構成に使う投影値のそれぞれを計算しています：
吸収用の投影値
A[x,y,v] = log(I0) - log( Σ_p=0～n I[p,x,y,v] )ただし、I0 はサンプルなしの I[p,x,y,v] から推定した入射光強度
屈折用の投影値
R[x,y,v] = p0 - Σ_p=0～n ( p * I[p,x,y,v] ) / Σ_p=0～n I[p,x,y,v]
ただし、p0 はサンプルなしの I[p,x,y,v] から推定した入射光の位置で、位相回復した値の符号が正常になるように R の符号を反転した。
なお、I[p,x,y,v] は HIS 形式ファイルの先頭に入っている複数枚の暗電流画像 D[p,y] の平均値を差し引いた X 線強度の値です。
この E-mail に添付した TIFF 形式画像 "*_dark.tif" が６月の SIXM 実験の測定 150602d と 150603jのそれぞれの 100 枚の暗電流画像の各画素の平均値の画像です。ただし、もとの 16 ビット画素値の値域は150602d1514.00 ～ 1750.07
150603j1514.45 ～ 1749.72
でしたが、ここではこれらの範囲を正規化した８ビット画像にしました。
(3) さて、今回の解析ではさほど活躍しませんでしたが、HIS 形式のデータファイルから吸収もしくは屈折の「生の投影画像」をそれぞれ抽出するプログラム his2ri（ri==raw intensity）と his2rr（rr==raw refraction）
を書きました：his2ri HIS darks scans views RI.bin
his2ri HIS darks scans views 0.bin 180.bin
X 線強度（屈折プロファイルの「面積」）
a[x,y,v] = Σ_p=0～n I[p,x,y,v] / (n+1)のバイナリ画像を作成するプログラム。
his2rr HIS darks scans views RR.bin
his2rr HIS darks scans views 0.bin 180.bin

57 / 124

屈折プロファイルの「重心」の位置
r[x,y,v] = - Σ_p=0～n ( p * I[p,x,y,v] ) / Σ_p=0～n I[p,x,y,v]のバイナリ画像を作成するプログラム
ただし、これらの起動パラメータ HIS、darks、scans と views の意味は以前に紹介した his2xp* のものと同じです。最初の起動法では CT 用の 180 度未満の回転ステップ v ごとの「投影画像」をまとめたバイ
ナリファイルを出力します。また、２番目は０度と 180 度の投影画像を別々のファイルに書き出します。
そして、これらの画像はプログラム a2tiff で TIFF 形式のものに変換可能です。
プログラム his2ri と his2rr のソースと Windows 用実行ファイルを以前に紹介した以下の書庫ファ
イルのそれぞれに追加しておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(4) ここからが今回の話の肝心な部分です。SIXM 実験の測定 150602d と 150603j の回転角０（v = 0）の投
影画像の左端のサンプルなしの部分のスキャン（x = 10）と、サンプルがある左右中央部のスキャン（x =199 or 254）の屈折プロファイル I[p,x,y,0] を調べてみました。この E-mail に添付した４個の（A4 横置きの）PDF ファイル ri+rp.pdf、rr+sd.pdf、rr+lp.pdf と rr+gp.pdf を御覧下さい。これらの上段には 150602d の、下段には 150603j の「生の投影画像」と２個の屈折プロファイルを示す図もしくは画像が並んでいます。
繰り返しになりますが、これらの屈折プロファイル I[p,x,y,0] の値はプログラム his2tiff で HIS 形式
ファイルから抽出した画像の画素値から暗電流画像 D[p,y] の各画素の平均値を引いた値です。
(5) ri+rp.pdf に示した投影画像は前記の his2ri で得た X 線強度画像で、極端に高い画素値を除外するような表示輝度にしてあります。そして、２×２個の屈折プロファイルはいずれも５おきの y 座標値のものだ
けを立体的に表示しており、それらの「高さ」と I[p,x,v] の対応関係はプロファイルごとにまちまちです。
通常の X 線 CT 用の測定画像と比較すると、これらの X 線強度画像では「照明」の空間的な不均質性は小さいように思われます。しかし、150602d では短波長の「縦縞」が、150603j ではそれより長波長の「横縞」が目に付きます。前者は屈折プロファイルの時間変動、後者は鉛直方向の空間変動に起因するはずですが、それは立体表示した屈折プロファイルの図でははっきりしません（プロファイルの高さが「絶対値」ではないためかもしれません）。
屈折プロファイルの全体的なかたちに関して以下のことが気になります。まず、これらのプロファイルは
いずれも「台地」の上に「山脈」が乗ったかたちをしており、この山脈の頂上の位置は概ねプロファイルの中央（p～n/2 = 63.5）です。また、台地のかたちは左右非対称で、その様相は 150602d（台地は右肩上がり）と 150603j（左肩上がり）で反転しています。そのため、次で議論する屈折プロファイルの重心の位置が山脈の頂上から左右にずれた位置になっているようです。
(6) 次に rr+sd.pdf、rr+lp.pdf と rr+gp.pdf の３つです。これらには his2rr で得た屈折プロファイルの重心の位
置（p = Pm）の符号反転した値（- Pm）を示す画像と、先の ri+rp.pdf のものと本質的に同じ２×２個の屈折プロファイルの画像が示されています。そして、これらのプロファイルの画像の上にはサンプルなしのスキャン位置（x = 10）とサンプルありのスキャン位置（x = 199 or 254）の重心の位置（p = Pm）をそれぞ58 / 124

れ青と赤色で示してあります。また、３個の PDF のそれぞれの緑色で示した位置の意味は以下の通りです。
rr+sd.pdf| p - Pm | = σ となる位置 p。ただし、σは屈折プロファイルから以下の式で計算した重心を中央とする高まりの裾野の幅：
σ^2 = Σ_p=0～n ( p^2 * I[p] ) / Σ_p=0～n I[p] - Pm^2
rr+lp.pdf
I[p] が最大値をとる屈折プロファイルの局所的頂点（local peak）。
rr+gp.pdf
この E-mail に添付した文書ファイル msd+gp.pdf に記した「２次関数の最小自乗法によるあて
はめ」によって得た屈折プロファイル全体のかたちを反映したピーク（global peak）の位置。
先にも書いたように、青もしくは赤色で示した屈折プロファイルの重心の位置 Pm は山脈の頂上の位置と一致しておらず、150602d ではその右側に、150603j では左側に分布しています。そして、山脈の頂上の位置に比べると Pm の位置はサンプルなしの場合でも多少揺らいでいます。しかしながら、（Pm を平均値と読み替えた場合に標準偏差に相当する）σの値は驚くほど一定で、具体的には、どの屈折プロフ
ァイルでもσ = 22 ～ 25 でした。
Pm と山脈の頂上の位置の食い違いが気になるので、local peak を探しました。これによって山脈の頂上の位置を知ることができましたが、データ解析の手法としてはノイズに弱いのでその結果を信頼できません。
そこで、屈折プロファイル全体のかたちを考慮した global peak を調べました。これは山脈の頂上と重心の間の位置になるようです（このことは理論的？にも証明できて、その詳細は msd+gp.pdf に記した通りです）。
(7) 以上のようにして屈折プロファイルのかたちを調べましたが、結局、どこに注目した解析をすれば良いかがわかりません。つまり、[a] なぜ屈折プロファイルは台地と山脈からなるかたちをしているのか？
[b] 台地のかたちは左右非対称で良いのか？
[c] 屈折プロファイルの台地の部分はゴミではないのか？
屈折プロファイルは理論的にはどういう「かたち」になっているべきですか？入射光に MTF を施したものなのでしょうけれども、それに関する理論や情報はありませんか？
最後はいい加減になりましたが、とりあえず以上です。
添付ファイル 150602d_dark.tif と 150603j_dark.tif
これらは暗電流画像でほとんど何も観察できないので、ここには掲載しません。

59 / 124

添付ファイル ri+rp.pdf

添付ファイル rr+sd.pdf

60 / 124

添付ファイル rr+lp.pdf

添付ファイル rr+gp.pdf

61 / 124

添付ファイル msd+gp.pdf

62 / 124

Date:
From:
To:
Cc:

Thu, 20 Aug 2015 12:44:34 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa"
Kentaro Uesugi,"TSUCHIYAMA, Akira", MIYAKE, "Aki Takigawa", Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: Re: SIXM_empty_HIS
たけうちさま、
GSJ/AIST のなかのです。丁寧な説明、ありがとうございます。SIXM の FZP の働きが良くわかりました（やはり、装置のことを理解せずにデータ解析するのは良くないですね）。「屈折プロファイル」の「山脈」は非常に目立つ上に、サンプルがある部分ではいくつかに分離していたり頂上の位置が揺らいでいたりするので信号だと思い込んでいましたが、これがゴミだとは、...。山脈の部分を除外した屈折プロファイルに２次関数をあてはめれば、より高精度に台地のかたちを反映した global peak を算出できそうですが、どのような基準でそれを除外すべきか？ また、分光器のΔθの揺らぎに起因するノイズは短波長（短周期）すぎてソフトウェ
アで濾過できそうにありませんが、これをどうするか？（露光時間を長くするなどのハードウェア的手法でこれを緩和できませんか？）。今後、このようなことを色々と考えて、さらなる質問をさせてもらうことになると思いますので、よろしくお願いします。とり急ぎ、お礼まで。
On Wed, 19 Aug 2015 18:55:07 +0900 "TAKEUCHI, Akihisa" wrote:> 中野様>> 竹内です。
> 返事が遅くなりました。
> 中野さんの説明のどこまで理解できているかはわかりませんが、> とりあえず最後の質問に対してお答えします。
> これもどの程度中野さんの満足いく回答になっているかはわかりませんが・・・。
>> [a] なぜ屈折プロファイルは台地と山脈からなるかたちをしているのか？
> [b] 台地のかたちは左右非対称で良いのか？
> [c] 屈折プロファイルの台地の部分はゴミではないのか？
>> ちょっと分りにくいかもしれませんが、簡単な絵を描いたので添付ファイルを参照ください。
> SIXM では二つの 1 次元 FZP を走査用と結像用に使っていますが、> ここでは簡単のため走査用だけで説明します。
> （結像用はただ単に縦に引き伸ばしているだけと考えて差し支えないので）>> まず、「台地」と「山脈」について。
> 答えから言うと、> 「台地」が信号：集光プローブが広がってカメラに入ったもの、> 「山脈」がごみ：FZP をそのまま抜けてきた光、> になります。
> 「山脈」は、FZP 特有にみられるものです。
>

63 / 124

> FZP の一次回折光を集光プローブとして試料をスキャンして、> プローブが広がった像を後ろのカメラで取り込んでいるので> この広がったプローブが「台地」のような像として検出されます。
> ただ、FZP はレンズとして理想的な素子ではなく、集光する光以外にもいろんな光が出てきます。
> この場合一番大きな影響があるのが非回折光、つまり FZP をそのまま透過してくる光です。
> これをできるだけ止めるために、試料面直前に OSA と呼んでいるスリットをおいてます。
> 「山脈」に見えるのは、このスリットを抜けてきた FZP の非回折光です。
> 非回折光はほとんど広がらないのに対して、回折光は約 60 倍に広がるので、> その分ピーク強度としては弱くなります。
> ですが、全体強度としては「台地」の方が圧倒的に大きく（だいたい 25:1）なので> 「山脈」=ノイズ は、無視してそのまま一緒に計算しています。
>> もちろん、本来は、「山脈」は分離して計算すべきです。
> ただ、光学的にセンタービームストップを置くなどして「山脈」を分離した実験も試みましたが、> 実際のところほとんど変わらない> （むしろ位相像では若干感度が落ちる。重心演算上では重心位置をキャンセルする向きに働くので）> ので、「山脈」は気にしないで計算しているのが実情です。
> 実際、.his のデータでは「山脈」の強度が飽和しているものも多くありますが、それも気にしません。
> それよりも「台地」の信号量を増やす方が結局 SN の向上につながるので。
>> 「山脈」が出ない素子としては、ミラーや屈折レンズなどがⅩ線領域でもありますが、> ミラーは装置として安定性が問題＋スペックルノイズで感度が落ちる> レンズは 8keV 前後でいい素子がない> ということで FZP を使っています。勿論、いいミラーやレンズが出てくれば> そちらに乗り換えるのがベストかと思います。
>> [a][c]についてはこれで回答になったかと思います。
>> [b]については、勿論左右対称であるに越したことはありません。
> データによって対称性が崩れてるのは、OSA、Diffuser、Δθなどのずれが考えられます。
> Δθのズレ自体は、本来左右対称性には関係ないはずですが、> 47XU での実験の場合、FE スリット開口を極端に狭めていてしかもその開口の形がいびつなため、> Δθがずれるとその影響が出てきてしまいます。
> 左右非対称であっても、測定の間に大きく変化なければ大した問題ではないと考えていますが> 中野さんも指摘されている短波長の縞状のノイズはΔθが原因かと考えています。
>> Martino の方法は僕自身もあれが最適な方法とは思っていません。現状可能な妥協案と考えています。
> Pfeiffer の方法ベースでノイズ除去も納得いく方法があればそれがベストと思います。
> その点は中野さんと完全に同意です。
> 単に僕がいい方法を見つけられてないだけという。
>> 以上、長くなりました。

64 / 124

竹内さんからの E-mail の添付ファイル

65 / 124

Date:
From:
To:
Cc:

Thu, 27 Aug 2015 17:40:54 +0900
Tsukasa NAKANO
Akira TSUCHIYAMA
Kentaro UESUGI, Akihisa TAKEUCHI,MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
Subject: Re: SIXM_150602_オマケ
つちやまさま、
GSJ/AIST のなかのです。
(1) ぼくが書いた SIXM に関連するほとんどすべての E-mails を１個のファイルにまとめて SPring-8 の WEB
サイトにアップロードしておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/SIXM_mails.pdf
http://www-bl20.spring8.or.jp/~sp8ct/tmp/SIXM_mails.docx(2) SIXM の画像再構成はぼくが直接説明しなくても、6/17～19 の E-mails に書いた画像再構成用計算機
プログラムの使用法と、6/27～7/13 の E-mails に書いたセンター値決定用プログラムの説明を読めば誰でも実行可能だと思います。ついでに言っておくと、7/16～17 の E-mails で報告した６月の SIXM 実験の画像の再構成は（画像のサイズが小さいので）ノートパソコンでも遂行可能です。
(3) ぼくが理解している SIXM の画像再構成の基礎は 6/25～7/7 の E-mails に書いた程度のものですから、
こちらも容易に理解できると思います。
(4) SIXM に関してぼくが一番の問題と思っていることは 8/10～20 の E-mail に書いた通りです。測定データから「真の信号」だけをきちんと取り出す方法を考えない限り SIXM は本当の意味で実用にならない気がします。
(5) また、8/1 の E-mail で報告した LAC と RID の２次元ヒストグラムの解釈には RID についても LAC の場
合と同様な「元素依存性などの物理」を勉強する必要がありますが、それも全然できていません。
(6) 以上のように、現時点では SIXM についてみなさんに講習するようなレベルの知識を持ち合わせていません。次回（11/27～12/2）の SIXM 実験には参加するつもりなので、直接の質疑応答はその時にお願
いします。
とり急ぎ、
On Sun, 23 Aug 2015 14:59:03 +0900 Akira TSUCHIYAMA wrote:> 中野さん> いろいろなものが返事遅れておりすみません。ようやく SIXM に対応できるようになりました。
> 再構成と画像処理ありがとうございました。8/24 からアルバイトの人が来るようになるので、> いろいろ遅れていることがある程度進むと思います（それでも結局は僕がボトルネックなのでしょうが）。
> 但し、再構成から含めて一連のやり方について、メイルを掘り起こして行こうとしたのですが、どうも簡単> にはいきそうにありません。前回のメイルでは（もっと広い意味で）講習会といいましたが、そんな大層な> ものでなくていいので、とりあえず SIXM についてだけでも教えにきていただけないでしょうか。
> 僕が筑波にお伺いしてもいいのですが、三宅さんや学生さんとも情報を共有し、京大でやれるようになる> ことを考えると、申し訳ありませんが京都までお越しいただくのが一番ありがたいです。
>…

66 / 124

Date:
From:
To:
Cc:

Fri, 25 Sep 2015 13:06:35 +0900
Tsukasa NAKANO"TSUCHIYAMA, Akira" Kentaro Uesugi
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama, Akihisa TAKEUCHI

Subject: RID の計算プログラム
みなさま、
GSJ/AIST のなかのです。
(1) 密度と化学組成が既知の物質の 10 eV～30 keV の範囲のエネルギーにおける RID（refractive index
decrement；１から屈折率の実数部の値を差し引いた、いわゆる「δ」の値）の理論値を計算するプログラ
ム rid_cm と rid_dcc などを書きました。
(2) これらの使い方は以前に紹介した X 線 LAC（linear absorption coefficients）の理論値の計算プログラム
lac_cm と lac_dcc の使い方とまったく同じです。それについては以下の文書の "lac" を "rid" に読み替えて下さい。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/h_s.pdf#page=3(3) 新しいプログラムのソースコードや 32 bit Windows 用実行ファイルなどを以下の書庫ファイルに入れておきましたので、どうぞ御利用下さい。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/henke.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/henke.zip(4) これらの書庫ファイルの内容は以下の通りです。
henke/
index.txt
原子番号 1～92 の元素それぞれに対する以下の４個の値が各行にタブコード区切りで書き込まれているテキストファイル。
[1] 原子番号（01～92）[2] 原子番号 / 質量数[3] 元素記号[4] 元素の名前
なお、プログラム lac_cm や lac_dcc を作成した時もこれとまったく同じファイルを使用しました。
url.txt
参考文献および RID などの計算に用いた atomic scattering factors のデータファイルの URLが書き込まれているテキストファイル。
sf.pdf
RID などの計算式を記した PDF ファイル。
run.sf
LBL（Lawrence Berkeley National Laboratory）のサイトから atomic scattering factors のデー
タすべてをダウンロードし、それらから eV 単位のエネルギーごとの元素の MRID（mass RID）
や MPAC（mass photo-absorption coefficients）を計算する C-shell script。

67 / 124

sf/*.nff と sf/read.me
LBL のサイトからダウンロードした書庫ファイルに入っていた atomic scattering factors のデー
タファイルとそれらの説明のファイル。
ev_mrid/*.txt と ev_mpac/*.txt
run.sf で得た 92 元素の MRID と MPAC のデータファイル。
run.lin_mrid と run.log_mpac
ev_mrid/*.txt や ev_mpac/*.txt のデータを C 言語コードに変換する C-shell scripts。
lin_mrid.h と log_mpac.hrun.lin_mrid や run.log_mpac が出力した長大な C 言語コード。
mrid.c、mrid.h、rid_cm.c と rid_dcc.c
rid_cm や rid_dcc の C 言語コード。コンパイル時にマクロ定数 LPAC を与えてやれば後述する
プログラム lpac_cm や lpac_dcc を生成することもできる。
Makefile"make" 入力で新しいプログラムすべてを gcc（GNU C++ compiler）でコンパイルする設定にしてある「Makefile」。
rid_cm.exe と rid_dcc.exe
rid_cm と rid_dcc の 32 bit Windows 用実行ファイル
lpac_cm.exe と lpac_dcc.exe
X 線 LPAC（linear photo-absorption coefficients）の理論値を計算するプログラムの 32 bit
Windows 用実行ファイル。
log_rid.pdf と lin_rid.pdf、および、log_lpac.pdf と lin_lpac.pdf
rid_cm や lpac_cm を使って計算した 5～25 keV における８種類の物質それぞれの対数もしく
はリニアスケールの RID や LPAC のプロファイルのグラフを並べた PDF ファイル。なお、LPACのグラフには lac_cm で計算した通常の LAC のプロファイルもプロットしてある。
log_mrid.pdf と lin_mrid.pdf、および、log_mpac.pdf と lin_mpac.pdf
rid_cm や lpac_cm を使って計算した 5～25 keV における８元素の対数もしくはリニアスケールの MRID や MPAC のプロファイルのグラフを並べた PDF ファイル。なお、MPAC のグラフには
lac_cm で計算した通常の atomic MAC のプロファイルもプロットしてある。
(5) 上に書いたように、今回は元素ごとの atomic scattering factors の虚数部の値から算出できる、吸収の素因として光電効果だけを考慮した LPAC の理論値の計算プログラム lpac_cm と lpac_dcc も書きました。
rid_cm や rid_dcc の場合と同様にこれらの使い方も lac_cm や lac_dcc のものとまったく同じです。なお、
LPAC（および MPAC）と通常の LAC（および MAC）の違いについては以下の文書の「物質による X 線の吸収」の章を御覧下さい。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/nakano_JGSJ.pdf#page=2
http://www-bl20.spring8.or.jp/~sp8ct/tmp/nakano_JGeolSocJpn.pdf#page=8
とりあえず以上です。

68 / 124

Date:
From:
To:
Cc:

Fri, 25 Sep 2015 15:54:05 +0900
Tsukasa NAKANO"TSUCHIYAMA, Akira", Kentaro Uesugi
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama, Akihisa TAKEUCHI

Subject: RID ~= density
みなさま、
GSJ/AIST のなかのです。
(1) 先程の E-mail で紹介した「RID の理論値の計算プログラム」の作成・テストの過程ですごいことに気づきました。結論から言うと以下の通りです。
SIXM で使用している 8 keV の X 線エネルギーでは元素ごとの MRID（mass refractive index
decrement）の値は桁で違わない。それゆえ、 物質の RID は概ねその密度（重力密度）で決まる。
つまり、SIXM で RID を可視化した位相画像はサンプル内部の密度分布を表している。
(2) この E-mail に添付した PDF ファイル"henke_sf.pdf"を御覧下さい。これは先程の E-mail で紹介した書庫
ファイル "henke.*" に入っている sf.pdf と同じファイルで、元素ごとの atomic scattering factors（の実数部）から RID の理論値を算出する際に用いた式が記されています。その先頭付近に記したように、これらは X-ray data booklet（その旧版は朱色の小冊子で、ぼくは SPring-8 の利用業務部で配布していたものを入手した）の 1.7（および 1.6）章に載っていたものを書き換えた式です。ただし、複合物質（composite material、CM）の RID（δ）の理論値の計算には 1.7 章の式 1 と 2 を組み合わせた式δ = ρ・na・re / ( 2・π)・λ^2・Σ_jw(j)・f1(j) / A(j)のままでは不便なので、LAC（μ）と MAC（τ= μ/ρ）の場合と同様な関係式に従う mass RID（MRID、Δ =δ/ρ）を導入しました。さらに、MAC の場合と同様に元素ごとの MRID（atomic MRID）を導入することにより、LAC の場合と同じ式で密度と化学組成が既知の物質の RID を計算することができます：
CM の LAC → μ = ρ・Σ_j w(j)・τ(j)
CM の RID → δ = ρ・Σ_j w(j)・Δ(j)← atomic MAC← atomic MRID

(3) さて、この E-mail に添付したテキストファイル mrid_8keV.txt を御覧下さい。そこには rid_cm に密度の値１を指定して得た 8 keV の X 線エネルギーにおける 92 元素それぞれの MRID の理論値が記されています。また、念のため、LBL のデータベース http://henke.lbl.gov/optical_constants/getdb2.html で得た
8 keV の MRID の値も記しました（92 元素すべてに対してではありませんが）。これらから以下のことがわかります。
[1] rid_cm で得た値と LBL のデータベースで得た値は数 % 程度違うだけ。つまり、rid_cm で計算した MRID の理論値は概ね正しい。なお、この MRID の理論値の違いは計算に使用した元素の質量数の値の違いによるのではないかと思われます。
[2] 8 keV における元素ごとの MRID は桁では違わない。

69 / 124

(4) このように 8 keV の atomic MRID がどの元素でもほぼ同じ値なら、それらに重量率を乗じて足し合わせ
た MRID はどんな物質でもほぼ同じ値になるはずです。この E-mail に添付したテキストファイル
rid_8keV.txt にそれを確かめた結果を記しました。ここでも rid_cm と LBL のデータベースの両方で RIDを計算し、それらを密度で割った MRID の値を示しました。これより８個の物質のどれでも MRID はほぼ同じ値（桁で違わない値）になっていることがわかります。
(5) 以上のことからこの E-mail の最初に書いた結論を得ました。つまり、RID の値に基づいた物質の化学組
成の推定は不可能です。そして、RID と LAC を縦と横軸にしてプロットした２次元ヒストグラムの解析も（RID が桁で変わらない密度を表しているだけなので）有意な情報をもたらしてくれそうにありません。
(6) 吸収トモグラフィーで得られる LAC の画像に比べて位相トモグラフィーの RID の画像が高「値」分解能に見えるのは、LAC に比べて RID の値の分布域の幅が狭いためではないか？
とりあえず以上です。
添付ファイル henke_sf.pdf
X-ray data booklet (http://xdb.lbl.gov/xdb-new.pdf)
Eric M. Gullikson
1.7 atomic scattering factors
1.6 mass absorption coefficientsindex of refractionn = 1 - δ - i・β= 1 - re / ( 2・π )・λ^2・Σ_j N(j)・f(j)δ refractive index decrement (RID)
β absorption index
re classical electron radius = 2.817940289458e-15 [m]λ wave lengthN(j) number of element "j" per unit volume
f(j) atomic scattering factor of element "j"---λ [m] = ( h・c / e ) / E [eV] = 1.239841576558181422983e-6 / E [eV]hceE [eV]

Planck constant = 6.6260689633e-34 [J・s]speed of light = 2.99792458e+8 [m / s]electron charge = 1.6021768740e-19 [C]photon energy in eV

N(j) = na・ρ・w(j) / A(j)na Avogradro's number = 6.0221417930e+23 [1 / mol]ρ density of material (including element "j")w(j) weight ratio of element "j"A(j) atomic weight of element "j"(1.7.1)

70 / 124

f(j) = f1(j) + i・f2(j)(1.7.2)

http://henke.lbl.gov/optical_constants/sf/sf.tar.gz---RIDδ = ρ・na・re / ( 2・π )・λ^2・Σ_j w(j)・f1(j) / A(j)atomic mass RID (MRID)
Δ(j) = na・re / ( 2・π )・λ^2・f1(j) / A(j)→
δ = ρ・Σ_j w(j)・Δ(j)absorption index
β = ρ・na・re / ( 2・π )・λ^2・Σ_j w(j)・f2(j) / A(j)atomic photo-absorption cross section
σ(j) = 2・re・λ・f2(j)(1.7.3)

atomic mass photo-absorption coefficients (MPAC)
τ(j) = na・σ(j) / A(j)= na・2・re・λ・f2(j) / A(j)(1.6.2)linear photo-absorption coefficients (LPAC)
μ = ρ・Σ_j w(j)・τ(j)= ( 4・π / λ )・β(1.6.3)添付ファイル mrid_8keV.txt

#
01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21

Z
H
He
Li
Be
B
C
N
O
F
Ne
Na
Mg
Al
Si
P
S
Cl
Ar
K
Ca
Sc

element namehttp://henke.lbl.gov/optical_constants/getdb2.html“echo 8 | rid_cm 1 Z:1”↓

HydrogenHeliumLithium
Beryllium
Boron
CarbonNitrogen
OxygenFluorine
Neon
SodiumMagnesium
Aluminum
SiliconPhosphorus
SulfurChlorine
Argon
Potassium
Calcium
Scandium

6.437122e‐06
3.242334e‐06
2.805434e‐06
2.882782e‐06
3.006378e‐06
3.251669e‐06
3.257848e‐06
3.265745e‐06
3.099623e‐06
3.249260e‐06
3.143735e‐06
3.251946e‐06
3.177650e‐06
3.293890e‐06
3.206130e‐06
3.305668e‐06
3.178893e‐06
2.989763e‐06
3.222023e‐06
3.299399e‐06
3.082658e‐06

6.43596513E‐06
3.24170946E‐06

3.25116821E‐06
3.25724636E‐06
3.26513623E‐06
3.09907796E‐06
3.14317708E‐06
3.25114024E‐06
3.17706713E‐06
3.29325189E‐06
3.20554068E‐06
3.30506236E‐06
3.17832087E‐06
3.22146593E‐06
3.29878389E‐06

71 / 124

22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

Ti
V
Cr
Mn
Fe
Co
Ni
Cu
Zn
Ga
Ge
As
Se
Br
Kr
Rb
Sr
Y
Zr
Nb
Mo
Tc
Ru
Rh
Pd
Ag
Cd
In
Sn
Sb
Te
I
Xe
Cs
Ba
La
Ce
Pr
Nd
Pm
Sm
Eu
Gd
Tb
Dy
Ho
Er
Tm
Yb
Lu
Hf
Ta
W
Re
Os

Titanium
Vanadium
Chromium
Manganese
Iron
CobaltNickelCopperZinc
Gallium
Germanium
Arsenic
Selenium
BromineKryptonRubidium
Strontium
Yttrium
Zirconium
Niobium
Molybdenum
Technetium
Ruthenium
Rhodium
Palladium
SilverCadmium
IndiumTin
Antimony
Tellurium
IodineXenon
CesiumBariumLanthanum
CeriumPraseodymium
Neodymium
Promethium
Samarium
Europium
Gadolinium
Terbium
Dysprosium
Holmium
ErbiumThulium
Ytterbium
Lutetium
Hafnium
Tantalum
Tungsten
Rhenium
Osmium

3.012809e‐06
2.937696e‐06
2.972503e‐06
2.884763e‐06
2.880300e‐06
2.687532e‐06
2.779224e‐06
2.763830e‐06
2.824569e‐06
2.766771e‐06
2.762639e‐06
2.777073e‐06
2.729634e‐06
2.788133e‐06
2.746693e‐06
2.776454e‐06
2.788059e‐06
2.826335e‐06
2.832497e‐06
2.857309e‐06
2.839109e‐06
2.851387e‐06
2.829369e‐06
2.845122e‐06
2.815352e‐06
2.837713e‐06
2.779523e‐06
2.776521e‐06
2.737230e‐06
2.716551e‐06
2.640242e‐06
2.699032e‐06
2.649563e‐06
2.650472e‐06
2.597848e‐06
2.597744e‐06
2.597941e‐06
2.598057e‐06
2.549963e‐06
2.548220e‐06
2.424541e‐06
2.329319e‐06
2.193281e‐06
2.278474e‐06
2.224571e‐06
2.079915e‐06
2.281706e‐06
2.343544e‐06
2.354834e‐06
2.385236e‐06
2.396681e‐06
2.422718e‐06
2.435712e‐06
2.449593e‐06
2.435614e‐06

3.01342402E‐06
2.97204656E‐06
2.88432807E‐06
2.88005117E‐06
2.68758595E‐06
2.77820436E‐06
2.76322953E‐06
2.82396422E‐06
2.76620176E‐06
2.76210562E‐06

2.77593949E‐06
2.78753942E‐06
2.82578526E‐06
2.8322836E‐06
2.85680335E‐06
2.83872146E‐06

2.83717463E‐06
2.77904724E‐06
2.77603544E‐06
2.73670639E‐06
2.71585373E‐06
2.69854831E‐06
2.6500245E‐06
2.59740773E‐06

2.42223177E‐06
2.43565523E‐06

72 / 124

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

Ir
Pt
Au
Hg
Tl
Pb
Bi
Po
At
Rn
Fr
Ra
Ac
Th
Pa
U

Iridium
PlatinumGold
Mercury
Thallium
Lead
BismuthPolonium
Astatine
Radon
Francium
RadiumActinium
Thorium
Protactinium
Uranium

2.451691e‐06
2.451074e‐06
2.470965e‐06
2.462452e‐06
2.451275e‐06
2.451554e‐06
2.463028e‐06
2.491795e‐06
2.514110e‐06
2.415146e‐06
2.437553e‐06
2.430562e‐06
2.447732e‐06
2.422310e‐06
2.451454e‐06
2.364970e‐06

2.45127012E‐06
2.45060232E‐06
2.47051344E‐06
2.46200534E‐06
2.45111232E‐06
2.46254604E‐06

2.41487305E‐06
2.43010049E‐06
2.42186047E‐06
2.36454139E‐06

添付ファイル rid_8keV.txt“echo 8 | rid_cm density chemical_formula”
density (g/cc)
chemical formulamaterialwater
1
H:2 O:1
quartz2.7
Si:1 O:2
calcite2.715
Ca:1 C:1 O:3
Ca:1 Mg:1 C:2 O:6
dolomite
2.86
Mg:2 Si:1 O:4
forsterite
3.226
Fe:2 Si:1 O:4
fayalite
4.4048
Mg:1 Si:1 O:3
enstatite
3.21
magnetite
5.2
Fe:3 O:4

RID at 8 keV

3.620619e‐06
8.853034e‐06
8.898499e‐06
9.350504e‐06
1.053804e‐05
1.347144e‐05
1.049759e‐05
1.553156e‐05

RID / density
3.620619e‐06
3.278901e‐06
3.277532e‐06
3.269407e‐06
3.266596e‐06
3.058355e‐06
3.270278e‐06
2.986838e‐06

http://henke.lbl.gov/optical_constants/getdb2.html
material
density (g/cc)
chemical formula
RID at 8000 eV
RID / density
3.61994717E‐06 3.61994717E‐06
water
1
H2O
8.85135159E‐06 3.27827843E‐06
quartz2.7
SiO2
8.89687453E‐06 3.27693374E‐06
calcite2.715
CaCO3
9.34872423E‐06 3.26878489E‐06
dolomite
2.86
CaMgC2O6
1.05358367E‐05 3.26591316E‐06
forsterite
3.226
Mg2SiO4
fayalite
4.4048Fe2SiO4
1.34696238E‐05 3.05794242E‐06
1.04954597E‐05 3.26961344E‐06
enstatite
3.21
MgSiO3magnetite
5.2
Fe3O4
1.55297585E‐05 2.98649184E‐06

73 / 124

Date:
From:
To:
Cc:

Wed, 30 Sep 2015 09:41:54 +0900
Tsukasa NAKANO
Akihisa TAKEUCHIMIYAKE, Akira TSUCHIYAMA, Kentaro UESUGI, Aki Takigawa,
Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI
Subject: Re: RID ~= density
みやけさま、
C.C. たけうちさま、
なかのです。竹内さんが書いておられるのは CT の前の元データ（投影画像）の「精度」ですね。それを画像再構成すると、少なくとも 数 % 程度の誤差（模擬データを用いたテストの結果）が CT 値に含まれます。ただし、これはサンプルの「端」で生じるエラーなので、均質なサンプルの中央の CT 値はもっと高精度に決まっているかもしれません。画像再構成で生じる誤差については以下の文献の図６を御覧下さい（ただし、SPring-8 の画像再構成では Chesler フィルタを使っています）。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/nakano_JGSJ.pdf#page=7
なお、長くなるのでここでは説明しませんが、SIXM の吸収と微分位相の投影画像は、測定した「refraction
profile」から投影画像の値を計算する処理がいい加減なので、10 % 程度の誤差が含まれていてもおかしくありません（つまり、現時点では SIXM の投影画像に含まれる誤差を評価できていません）。とり急ぎ、On Wed, 30 Sep 2015 09:09:07 +0900 "TAKEUCHI, Akihisa" wrote:> 三宅様>> 以前の CT 研究会でもちらっと紹介しましたが、感度は delta も beta も３シグマで 10^-7(8keV)程度。
> これを delta を密度、beta を LAC（8keV)に直すとそれぞれだいたい 40mg/cm^3、130cm^-1 くらいです。
> もちろん測定条件や試料で変わります。
> 現状使ってもらってる測定条件（「竹」レベル）だと上記よりも少し悪くなると思いますが、> それでも 0.2～0.3 程度の密度差であれば識別できると思います。
>> ---------------------------------------------------------------> On 2015/09/30 8:28, MIYAKE wrote:> > 中野様 竹内様> > この質問はどちらかと言うと、竹内さんなのかもしれませんけど、> > どれくらいの分解能で、密度の差がわかるのでしょうか？ 5%?, 1%? それ以下？？
>>> > 例として、良いかは置いておいてなのですが、昨年度の卒業生が流体包有物中の鉱物を> > 探索したさいのデータからですけど、>>>>
density
LAC(20keV)>>
FeCl2・(H20)
2.13
29.39
>>
FeCl2・2(H20)
2.39
29.54
>>
NaCl
2.16
15.44
>>
CaCl2・0.5(H20) 2.73
15.07

74 / 124

>>
CaSO4
2.96
17.26
>>
K2SO4
2.68
17.30
>>
> > とかだと、5%分解能があれば、区別がつきそうです。（もちろん、それ以外の手法で> > 区別が付くのかもしれませんけど。）> > 特に水和物とかには、有効そうな気がします。
>>> > ---------------------------------------------------------------> >> 中野さん、> >> 竹内さん、三宅さん、他皆さん> >> 土`山＠フランスです> >>> >> コメント遅くなってすみません。
> >> 中野さんの解析、ありがとうございました。何故、RID がほぼ密度に比例するのかよくわかりました。
> >>
> >>>
(5)
> >>>
以上のことからこの E-mail の最初に書いた結論を得ました。つまり、RID の値> >>>に基づいた物質の化学組成の推定は不可能です。
> >> これは、おっしゃるとおりです。
> >>
> >>>
そして、RID と LAC を縦と横> >>>軸にしてプロットした２次元ヒストグラムの解析も（RID が桁で変わらない密度> >>>を表しているだけなので）有意な情報をもたらしてくれそうにありません。
> >> これは、三宅さんや竹内さんのコメントと同様、以下のように同意できません。もっとも中野さんが> >> いう「有意な」がなになのかによって、いい方が違うだけかもしれませんが。
> >>> >> 勿論、２次元ヒストグラムの解析も化学組成の情報が増えるわけではありませんが、RID（~密度）が> >> 桁で変わらなくともその分解能さえあれば、LAC と組み合わせることにより、DET ではできない物質> >> の同定がある RID-LAC 空間領域では可能となります（正確には、同定ではなく、物質候補が> >> わかっていてこれをカンニングしているだけなのですが・・・、DET でも同じです）。
> >>> >> 従って、空気と水や有機物との区別ができます。また、水の塩分濃度もある程度は議論できる。
> >> 一方、有機物は候補の密度がわかっていないと、同定はできませんね。また、密度が１に近ければ> >> 水とも区別できない。
> >>> >> ---------------------------------------------------------------> >> 2015/09/29 16:34、TAKEUCHI, Akihisa のメール：> >>> 中野様
> >>>> >>> すみません、RID に関しては中野さんのご指摘のように、ほぼ密度に比例した関係があり、> >>> 原子番号 Z と質量数 A の間にだいたい Z/A~0.5 の近似が成り立つとして> >>> RID ~ 1.35 x 10^-6 ρ[g/cm^3] λ^2[A]> >>> という近似式がよく利用されてます。

75 / 124

> >>> なので、MRID はほぼ定数、という考え方もできるわけです。
> >>>> >>> また、LAC と RID を利用した例としては、たとえば> >>> 位相像と透過像の Log の比から実効原子番号 Z_eff を求める方法> >>> http://pfwww.kek.jp/publications/pfnews/32_1/saikin3.pdf> >>> あるいは
> >>> 位相 CT 像（RID ~ 密度分布）と吸収 CT 像（LAC)の比から> >>> 実効 MAC の分布を求める方法なんてのも報告があります。
> >>> これらは対象が単元素であれば威力はありますが、> >>> 複数の元素で構成される場合はその平均がでてきます。
> >>>> >>> 僕の 2D ヒストグラムの場合はすでにある報告例を真似てもしょうがないし、比にしてしまうと> >>> 情報が減って勿体無いので、土’山さん達が 7keV 8keV の 2D ヒストグラムを示していたのに> >>> 「インスピレーションを受けて」位相ー吸収の 2D ヒストグラムをやってみた次第です。
> >>> なので逆に、何らかの情報を引き出すには上記の方法を利用されてみても良いかもしれません。
> >>>> >>> 中野さんのご指摘のように、2D ヒストグラムから化学組成がわかるか、と言われたら> >>> ちょっとむり、と言うのが正直なところです。
> >>> 添付の文献にもあるように、RID と LAC の比と元素番号の間には一つの関係式が成り立つので> >>> 2D ヒストグラムの分布はその関係式が示す曲線の上にほぼ乗ってきます。
> >>> ただ、単純に、2 つの別個の情報を一瞥で見える形にすることで、内容物の分類にはある> >>> 程度役に立つのではと考えています。
> >>> 吸収、位相それぞれ単独のヒストグラムでは分離できないものが 2D ヒストグラムで> >>> 浮き出て見えてくるケースも実際にありますし。
> >>>> >>> 中野さんの（6）の考察の「分布の幅」がちょっと理解できなかったのですが、> >>> 組成がわかるかどうかはさておき、RID と LAC が同じ程度の精度で測定できた場合、> >>> （特に軽元素で）その断面積が RID が格段に大きいので本質的に位相の方が格段に感度が> >>> 上がります。
> >>> 実際、現状の SIXM は複素屈折率の delta（位相）と beta（吸収）はほぼ同じ感度で> >>> 測定できていますが、もともと delta が beta に比べて格段に大きいので、像としては位相像が> >>> 格段に感度が良く見え、2D ヒストグラムにおいては分布が横長（吸収の分解能が悪い）に> >>> なります。
> >>>> >>> 以上です。
> >>> 僕の理解が不十分だったかもしれないので> >>> 上記の回答もひょっとしたら頓珍漢かもしれませんが。
> >>>> >>> ---------------------------------------------------------------> >>>> 2015/09/29 19:00、MIYAKE のメール：> >>>> 中野様> >>>> 僕でもなんとなくは理解出来、> >>>> >RID の値に基づいた物質の化学組成の推定は不可能です。

76 / 124

> >>>> は、分かりました。
> >>>>> >>>>> そして、RID と LAC を縦と横軸にしてプロットした２次元ヒストグラムの解析も（RID が> >>>>> 桁で変わらない密度を表しているだけなので）有意な情報をもたらしてくれそうにありません。
> >>>>> >>>> これは、LAC 値に大きな差がみえなくて、密度に違うある場合、> >>>> または、その逆の時には使えるということですよね。
> >>>> この例として何がふさわしいのか分かりませんけど。
> >>>> また、位相トモグラフィーの RID の密度分解能がどれくらいかにもよるのかもしれませんけど。
> >>>>> >>>> 三宅
Date:
From:
To:
Cc:

Wed, 30 Sep 2015 12:44:18 +0900
Tsukasa NAKANO
Akihisa TAKEUCHI
Kentaro UESUGI, MIYAKE, Akira TSUCHIYAMA, Aki Takigawa,
Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI
Subject: Re: RID ~= density
たけうちさま、
C.C. みなさま、なかのです。
> 中野さんが指摘された事以外に、位相の場合は> 計算エラーで生じるバックグラウンドの緩い傾斜なんかもあり、「ぼくが指摘したこと」が「画像再構成演算に伴う計算誤差」を意味するなら、それは実質的には問題ないと思います（均質なサンプル像の端の CT 値以外は）。「計算エラーで生じるバックグラウンドの緩い傾斜」とは、ぼくが 8/18 の E-mail で指摘した「refraction profile」の左右非対称性のことでしょうか？
これは「X 線が投影面に垂直に入射していないこと」によって生じる幾何学的な効果ではないかと思っています。そのため、ある条件を満たしていない場合、refraction profile の「重心演算」では「屈折角」をきちんと推定できない。これについての説明は少々お待ち下さい。とり急ぎ、On Wed, 30 Sep 2015 12:25:44 +0900 "TAKEUCHI, Akihisa" wrote:> 0.2-0.3 程度の密度差が識別可能、というのは、その都度絶対値にふらつきはあるものの、> 同じデータ内 2 個間の相対値差がその程度であれば見分けがつきますよ、という意味ですね。
> すみません。
> 先ほど僕が示した値は、CT 像から算出した値です。
> ちなみに試料はガラスロッド、シーホスター、ポリスチレン球です。
> 再構成の３D データ全体から算出したので、中野さんが指摘される> CT 値の誤差は含まれていると考えるのが妥当です。
> ただ、均一な試料なので影響は比較的小さいと思いますが。
>

77 / 124

> 中野さんが指摘された事以外に、位相の場合は> 計算エラーで生じるバックグラウンドの緩い傾斜なんかもあり、> それが土'山さんの云う「ふらつき」につながるのかな、と> 考えてますが、ちょっとまだわかりません。
>> ---------------------------------------------------------------> On 2015/09/30 9:45, Kentaro UESUGI wrote:> >>それでも 0.2～0.3 程度の密度差であれば識別できると思います。
> > 土`山さんが研究会の時に示していた「絶対値のふらつき」が> > 最大の敵と言うことになるでしょうか。そして、もしかするとこのふらつきは> >>なお、長くなるのでここでは説明しませんが、SIXM の吸収と微分位相の投影画像は、> >>測定した「refraction profile」から投影画像の値を計算する処理がいい加減なので、> >>10 % 程度の誤差が含まれていてもおかしくありません（つまり、…
> > が原因？ (いい加減、というよりもきちんと投影画像を求められない状況のこと)
Date:
From:
To:
Cc:

Wed, 30 Sep 2015 19:23:33 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa"
MIYAKE, Akira TSUCHIYAMA, 上杉 健太朗, Aki Takigawa, "Junya Matsuno",
中村隆太, "Takashi Sakurama", Yohei IGAMI

Subject: Re: RID ~= density
たけうちさま、なかのです。昨日いただいた E-mail の以下の部分に関して応答するのを忘れていました。
On Tue, 29 Sep 2015 23:34:20 +0900 "TAKEUCHI, Akihisa” wrote:> 中野さんの（6）の考察の「分布の幅」がちょっと理解できなかったのですが、> 組成がわかるかどうかはさておき、RID と LAC が同じ程度の精度で測定できた場合、> （特に軽元素で）その断面積が RID が格段に大きいので本質的に位相の方が格段に感度が上がります。
> 実際、現状の SIXM は複素屈折率の delta（位相）と beta（吸収）はほぼ同じ感度で> 測定できていますが、もともと delta が beta に比べて格段に大きいので、像としては位相像が> 格段に感度が良く見え、2D ヒストグラムにおいては分布が横長（吸収の分解能が悪い）に> なります。
ぼくが言った「分布の幅」は上記の「2D ヒストグラムにおいては分布が横長」で言っている「分布の幅」と同様な意味です。２次元ヒストグラム（および、普通の１次元ヒストグラム）を作る時、ぼくは画像に含まれる
CT 値の最小値から最大値の範囲をヒストグラムの調査範囲にしています（多分、つちやまさんも）。ただし、このこと（CT 値の分布域をヒストグラムの調査範囲と一致させること）には実用的であること以上の深い意
味はありません。「２次元ヒストグラムで分布が横長」になっているのは見かけ上の話で、それが「吸収の分解能の悪さ」を意味するわけではないと思います。とり急ぎ、

78 / 124

Date:
From:
To:
Cc:

Thu, 01 Oct 2015 16:26:24 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", Kentaro Uesugi"TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI
Subject: sum_of_projection(SP)の計算プログラム

みなさま、
GSJ/AIST のなかのです。
(1) SIXM などで得た投影画像のスライスごとの投影値の総和（sum of projection；SP）の計算プログラム
xp2sp（吸収 CT 用）と xp2sp_dpc（屈折 CT 用）を書いてみました。
注
SP は従来の吸収 CT において total absorption（TA）と呼んでいた量ですが、屈折 CT に対しても違和感のない呼び名に変えました。
(2) SP は画像再構成の過程で保存される、以下の式で表される量です。
スライス z の sum of projection、SP(z)= 回転角θの投影画像の投影値 P(r,z,θ) の横一列の和、Dr・Σ_r P(r,z,θ)= スライス z の再構成画像の CT 値 f(x,y,z) の総和、Dr^2・Σ_x&y f(x,y,z)
ただし、Dr は投影画像の画素の横方向の辺長で、ぼくが書いた画像再構成プログラムでは再構成した
スライス画像の正方形画素の辺長もそれと同じ値です。
(3) 以前に紹介したプログラム his2xp や his2xp_top で得られる SIXM の屈折 CT 用の投影画像には「サンプル内部の位相シフト（RID）の r 方向の微分値（微分位相）の X 線光路に沿った積算値 dP/dr」が入っているので、それに対応したプログラム xp2sp_dpc では微分を差分近似した以下の式によって dP/drから投影値 P を得た後に SP を計算するようになっています。
投影画像の左端の外部では P(r,z,θ) = 0 と仮定
P(r,z,θ) = P(r - Dr,z,θ) + Dr・dP/dr(r,z,θ)(4) xp2sp と xp2sp_dpc のコード一式や 64 bit Windows 用の実行ファイルを以前に紹介した SIXM 用の２個の書庫ファイルの両方に入れておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(5) xp2sp と xp2sp_dpc の起動法は以下の通りです。

79 / 124

xp2sp XP.bin {Dr} SP.bin
xp2sp_dpc XP.bin {Dr} SP.bin
XP.binhis2xp、his2xp_top や hp2xp で作成した投影画像のバイナリファイルの名前。"-" を指定するとデータを標準入力から読み込む。
Dr
投影画像の画素の横方向の辺長。指定を省略すると１と見なす。
SP.bin
横軸をサンプル回転角θ（０～180 度）、縦軸をスライスの位置 z とする SP(θ,z) のローカルなバイナリ形式の画像を書き込むファイルの名前。"-" を指定するとデータを標準出力に書き出す。
なお、SP.bin に書き込んだバイナリ形式の画像は以前に紹介したプログラム a2tiff を使えば TIFF 形式の画像に変換することができます：# SP の最小値と最大値の範囲を正規化した８ビット画素値の画像に変換
a2tiff SP.bin 8 SP.tif(6) 実際の測定データから SP の TIFF 画像を作る具体的な手順は以下の通りです。
[1] 測定データから投影画像のバイナリファイルを作成する。
# 通常の X 線 CT の測定 150616e の投影画像 xp.bin を作る場合
dit_tbl 150616e/raw | hp2xp xp.bin# SIXM の測定 150602d の HIS 形式データ（暗電流画像の枚数 100、
# スキャン数 401、サンプル回転ステップ数 451）から投影画像の左右# 両端の 30 画素幅の領域の値が入射 X 線強度そのものだと仮定して# 吸収 CT 用の投影画像 a.bin と屈折 CT 用の画像 r.bin を作る場合
his2xp 150602d/a.HIS 100 30 401 451 a.bin r.bin[2] 投影画像から SP の画像のバイナリファイルを作成する。
xp2sp xp.bin sp.bin
xp2sp a.bin sp_a.bin
xp2sp_dpc r.bin sp_r.bin[3] SP の画像を TIFF 形式のものに変換する。
a2tiff sp.bin 8 sp.tif
a2tiff sp_a.bin 8 a.tif
a2tiff sp_r.bin 8 r.tif(7) SIXM の投影画像のシーケンスに後述する「おかしなこと」が生じていないかを見るため、それらの画像や通常の X 線 CT の投影画像の SP を画像化してみました。この E-mail に添付した PDF ファイル sp.pdf
を御覧下さい。そこには７測定のデータから作成した 13 組の「回転角０の投影画像」と「SP の画像」のペアを並べてあります。これらの測定・投影画像の詳細は以下の通りです。

80 / 124--- 以下の２個は SP の「基準？」となる測定のデータ
020607a
BL47XU で撮影した金沢大学・森下くんの simplectite の CT 画像
http://www-bl20.spring8.or.jp/~sp8ct/hd2/morishita_symplectite/020607a/byte_222.gifから X 線 CT シミュレータで合成した、SPring-8 の通常の測定のものと同形式になっている測定データ。これは合成したデータなので SP はサンプル回転角に関して一定の値をとる理想的なものになっているはず。
040711j
BL20B2 で撮影した直径が 1 mm 以下の均質な球形ビーズの CT 画像
http://www-bl20.spring8.or.jp/~sp8ct/tmp/040711j.gifの測定データ。上記のような良好な画質の画像を再構成できている。
--- 以下の３個は SIXM で同じサンプルを撮影した FZP-CT の測定データ
150616e
BL47XU の FZP-CT で撮影した吉田くんの石英サンプルの測定データ。これは SIXM の測定
150602d で撮影したものと同じサンプル。
150619a および 150619b
BL47XU の FZP-CT でそれぞれ７および８ keV の X 線エネルギーで撮影したクロタイラギ貝
サンプルのデータ。これは SIXM の測定 150603j で撮影したもの同じサンプル。
--- 以下の８個は２個の測定から手法を変えて抽出した投影画像
150602d_a、150602d_a_top、150602d_r、150602d_r_top
150603j_a、150603j_a_top、150603j_r、150603j_r_top
これらはそれぞれ SIXM の測定 150602d と 150603j で得た投影画像。記号 "_a" と "_r" はそれぞれ吸収と屈折 CT 用の画像。また、"_top"付きのものは his2xp_top で上端の 10 画素幅の領域の値が入射 X 線強度だと仮定して作成した画像で、それなしのものは his2xp で左右両端の 10 画素幅の領域の値を入射 X 線強度だとした画像。
(8) これらの 13 個の SP の画像の比較から以下のことがわかります。
[1] 020607a の SP の横方向の縞に変動はなく、理想的なものになっている。
[2] それ以外の通常の X 線 CT の 040711j や 15061?? の SP では横方向の縞模様に若干の変動がある。しかし、再構成画像に致命的な問題はないようなので、この程度の変動はあっても問題ないと思われる。
[3] SIXM の SP でも上記のものと同程度の縞模様の変動があるが、位相 CT 用の投影画像（"*_r*"）
で目立つ縦線以外は問題ない？全体的に見て、SIXM の SP は通常の X 線 CT のものとさほど違わないように思える。
[4] 問題は SIXM の吸収 CT 用投影画像の SP のコントラストが入射 X 線の受光を仮定した領域の選択（"_top" とそれ以外）で違っていること。これは何故だ？
以上のような SP の比較結果をまとめると、SP で見る限り SIXM の投影画像は通常の X 線 CT のものと大差ないといえます。実は、ぼくは SIXM のマイクロビーム形成用のスリットがサンプル内で顕著に屈折した光束の両端を遮ることがあるのではないかと疑っていました。それが生じているならサンプル回転角方向のSP の変動（横縞の途切れ、もしくは縦縞）があるはずですが（これが前述の「おかしなこと」です）、それを示
唆するものは見つかりませんでした。最後はちょっといい加減になりましたが、とりあえず以上です。

81 / 124

添付ファイル sp.pdf

82 / 124

Date:
From:
To:
Cc:

Wed, 07 Oct 2015 18:56:03 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", Kentaro Uesugi"TSUCHIYAMA, Akira, MIYAKE, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI
Subject: SP の計算プログラムの改造
みなさま、
GSJ/AIST のなかのです。
(1) 10/1 の E-mai で紹介した SIXM などで得た投影画像のスライスごとの投影値の総和（sum of
projection；SP）を計算する２個のプログラムのうちの屈折 CT 用の xp2sp_dpc を少しだけ改造しました。
(2) 以前の xp2sp_dpc では「サンプル内部の RID の r 方向の微分値の X 線光路に沿った積算値 dP/dr 」
から微分を「後方差分」で近似した式>>

投影画像の左端の外部で P(r,z,θ) = 0 と仮定
P(r,z,θ) = P(r - Dr ,z,θ) + Dr・dP/dr(r,z,θ)を使って r 方向の投影値それぞれを計算（線積分）していましたが、微分値のデータに含まれる誤差のため、この場合には投影画像の右端で P(r,z,θ)≠0 となる可能性があります。また、逆に、微分を「前方差分」で近似した式
投影画像の右端の外部で P(r,z,θ) = 0 と仮定
P(r,z,θ) = P(r + Dr ,z,θ) - Dr・dP/dr(r,z,θ)を使った場合も同様で、画像の左端で P(r,z,θ)≠0 となる可能性があります。

(3) 新しい xp2sp_dpc ではこれら２個の近似式で得た投影値の平均値を使いました。その結果、微分値のデータの誤差は r 方向の P(r,z,θ)の値それぞれに均等に割り振られるはずです。以下の式を用いて投影値や SP を計算しています。
D_i ≡ Dr・dP/dr(Dr・i, z,θ) （ただし、i=0 ～ N - 1）とする。
後方差分の近似式を使って推定した
投影値、A_i = Σ_j=0～i D_j
SP の値、A = Dr・Σ_i=0～N-1 A_i = Dr・Σ_i=0～N-1 (N - i)・D_i
前方差分の近似式を使って推定した
投影値、B_i = - Σ_j=i～N-1 D_j
SP の値、B = Dr・Σ_i=0～N-1 B_i = - Dr・Σ_i=0～N-1 (i + 1)・D_i
以上からSP(z,θ) ≒ (A + B)/2 = Dr・Σ_i=0～N-1 { (N - 1) / 2 - i }・D_i(4) 改造した xp2sp_dpc の使用法は以前のものとまったく同じです。また、SIXM 用の２個の書庫ファイルの中のコードや 64 bit Windows 用実行ファイルを改造後のものに差し替えておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip83 / 124

(5) ７測定・13 セットの X 線 CT の投影画像の SP を再度計算してみました。この E-mail に添付したファイル
sp.pdf を御覧下さい。赤枠で囲んだ４個の SP の画像だけが 10/1 の E-mail に添付した sp.pdf に示したものと違っています。改造した xp2sp_dpc は微分位相のデータに含まれる誤差をそれぞれの投影値に割り振っているので、それらから計算した SP の値域の幅が狭くなっています。
dataset150602d_r
150602d_r_top
150603j_r
150603j_r_top

old or new
old
new
old
new
old
new
old
new

range of SP-35022.811775 ~ 142931.442415-13654.272798 ~ 43565.716593-64800.654479 ~ 132946.347413-12800.121870 ~ 40799.714474-70118.723454 ~ 77576.755749-6986.372399 ~ 32745.017990-162153.547507 ~ 39503.043495-4068.297573 ~ 30833.674252

これにより画像の輝度のコントラストが向上したため、10/1 の sp.pdf では識別できなかった４個の SPの画像の「黒い縦線」を観察できるようになりました。
(6) 問題はこれらの黒い縦線が生じた原因です。SP の画像上での横座標値（== 投影画像のシーケンスの順番、０～）は測定ごとに概ね同じ値なので、これらの黒い縦線はデータ処理の過程で生じたものではありません。
測定番号150602d150603j

黒い縦線の横座標値
96 および 320 前後の値
59 と 75 前後の値、および 295 ± 10

そこで、8/18 の E-mail で紹介したプログラム xp2rr を使ってこれらの測定の HIS 形式ファイルから屈折の「生の投影画像（のシーケンス）」を抽出し、動画を作りました。それらのファイルはサイズが大きいので、計算機 vrm の以下のディレクトリにコピーしておきました：
vrm:/media/disk/bl47xu/150602_sp/
150602d_r_sp.tif、150602d_r_top_sp.tif、150603j_r_sp.tif と 150603j_r_top_sp.tif
この E-mail に添付した sp.pdf の上のものと同じ屈折の SP の画像
150602d_rr.gif と 150603j_rr.gif
屈折の生の投影画像の動画（GIF 形式）
注
このディレクトリには測定 150602d と 150603j の以下の画像や動画もコピーしてあります。
sp.pdf に載せたものと同じ吸収 CT 用の SP の画像
150602d_a_sp.tif、150602d_a_top_sp.tif、150603j_a_sp.tif と 150603j_a_top_sp.tif
プログラム his2ri で抽出した吸収の生の投影画像の動画
150602d_ri.gif と 150603j_ri.gif
これらの屈折 CT 用の生の投影画像の動画を観察すればわかるように、多面体状のサンプルの大きな
facets が X 線の照射方向と平行になっている撮影方向において SP の画像の黒い縦線が発生しています。これは言われてみれば当たり前のことですが、以下が屈折 CT の SP の画像の黒い縦線の原因だと思われます。

84 / 124

[1] RID の空間勾配∇δは facet を横切る時に非常に大きな値になる。
[2] その facet が X 線の照射方向と平行になっている投影画像でベクトル∇δの絶対値 |∇δ| を測定することになる。
[3] この |∇δ| は facet の位置で「デルタ関数」的に空間分布しているので、そのデータを数値的に線積分しても facet の位置で「階段関数」的に変動する投影値の空間分布を正確に推定することは不可能。
[4] この「デルタ関数」の「高さ」として真の値よりも低い値を測定することになるので、それを線積分した「階段関数」の「段差」も真の値より小さくなる。すなわち、この撮影方向の画像の投影値の推定値は（他の方向に比べて）全体的に低い値になり、SP の値も小さくなる（== SP の
画像上に黒い縦線が発生する）。
とりあえず以上です。

85 / 124

添付ファイル sp.pdf

86 / 124

Date:
From:
To:
Cc:

Thu, 08 Oct 2015 09:19:16 +0900
Tsukasa NAKANOMIYAKE"TAKEUCHI, Akihisa", 上杉 健太朗, "TSUCHIYAMA, Akira", Aki Takigawa,
Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI
Subject: Re: SP の計算プログラムの改造

みやけさま、
なかのです。E-mail ありがとうございます。
> そのとき、吸収の像は、明るい千線として出るのかなぁ、と思うのですが、屈折の場合は「facet を横切る RID の空間勾配（これがデルタ関数的になる）」を線積分して投影値を計算しているため、その方向の SP が低い値になります。しかし、吸収では投影値そのものを測定できるので、この現象は起こりません。また、facet による屈折・回折は空間的な輝線・暗線になりますが、それらの投影値の総和である SP では、正・負の値が打ち消しあうため、輝線・暗線にとして出てこないのだと思われます。
とり急ぎ、
On Thu, 8 Oct 2015 08:38:05 +0900 MIYAKE wrote:> 中野様 竹内様>>>
黒い縦線の原因に関する考察は概ねその通りだと考えています。
>>
このように像として表現するとよりわかりやすいですね。
>> 黒い縦線につきまして、なるほど、とおもいました。
> そのとき、吸収の像は、明るい千線として出るのかなぁ、と思うのですが、> 頂いた像を見る限り、そのようには見えません。これは、何故なのでしょうか？
> 単に顕著でないだけなのかもしれませんけど。
>> また、こうした facet があるといけないとすると、やはり円柱が良いということですかね？
>> ---------------------------------------------------------------> > 中野様>>> > 解析ありがとうございます。
> > 前方差分と後方差分の平均をとる方法は僕も試してみたことがあるんですが、> > 位相像のストリークノイズを消す目的ではあまり効果がなかったのでそれ以降> > しつこくやってなかったんですが、こういう解析ではしっかり違いが出てくるんですね。
>>> > 黒い縦線の原因に関する考察は概ねその通りだと考えています。
> > このように像として表現するとよりわかりやすいですね。
>>> > 中野さんがおっしゃる「デルタ関数」の高さの誤差がどうしても測定値は低く出てしまうので> > 位相像の測定値は実際より小さくなり、> > また、右と左でゼロに収束しない原因の一つにもなっているのではと考えています。
>>> > 実は、位相の測定値の定量性を評価するために最初の頃は、テストチャートを利用していたんですが、> > （素性も既知で空間周波数との関係も考慮しやすそうという理由から）> > これもエッジが立ち過ぎていて同様のエラーが生じやすいのでやめた経緯があります。

87 / 124

Date:
From:
To:
Cc:

Fri, 09 Oct 2015 11:54:28 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", Kentaro Uesugi"TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa,
Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI
Subject: filtered_projection（FP）の計算プログラム
みなさま、
GSJ/AIST のなかのです。
注
すみません。「外部からの標的型メール等によるネットワーク障害」のため、昨日書いたこの
E-mail を今日になって送信できました。以下の文章の「今朝」は「昨日の朝」、「昨日」は「一昨日」、...、です。それを考慮して以下の文章をお読み下さい。よろしくお願いします。
On Thu, 8 Oct 2015 06:35:20 +0900 "TAKEUCHI, Akihisa" wrote:> 前方差分と後方差分の平均をとる方法は僕も試してみたことがあるんですが、> 位相像のストリークノイズを消す目的ではあまり効果がなかったのでそれ以降> しつこくやってなかったんですが、こういう解析ではしっかり違いが出てくるんですね。
(1) 今朝早くの竹内さんの E-mail に書いてあった上記のことに関してですが、実は、ぼくが書いた屈折 CT
用画像再構成プログラム xp2tg_dpc は微分値から投影値を計算する際に「前方差分と後方差分の平均」
を使っていません。7/2 の E-mail で説明したように Pfeiffer 論文（PRL、98、10815、2007）に載っていた手法を用いて、CBP 法による画像再構成の畳み込み演算（convolution；CBP 法の「C」の処理です）と同時に部分積分で投影の微分値を投影値に変換（線積分）しています。それゆえ、昨日の E-mail に添付した sp.pdf に載せた屈折 CT 用の SP が xp2tg_dpc で再構成した画像上の CT 値の総和と一致する「保存量」になっているかどうかは定かではありません。
(2) ところで、10/1 の E-mail でぼくは SP について以下のように説明しましたが、それは不十分でした。
On Thu, 01 Oct 2015 16:26:24 +0900 Tsukasa NAKANO wrote:> SP は画像再構成の過程で保存される、以下の式で表される量です。
>>
スライス z の sum of projection、SP(z)>= 回転角θの投影画像の投影値 P(r,z,θ) の横一列の和、Dr・Σ_r P(r,z,θ)>= スライス z の再構成画像の CT 値 f(x,y,z) の総和、Dr^2・Σ_x&y f(x,y,z)
画像再構成の畳み込み演算で得られる filtered projection（FP）の r 方向の和も上記の２種類の値の和と同じ値（SP の値）になるはずです：
SP(z) = Dr・Σ_r Q(r,z,θ)← 投影値の和の場合と同じ式

畳み込み演算で得られる FP の値 Q(r,z,θ)= ∫ P(ρ,z,θ) ・G(ρ - r) dρ
← 吸収 CT の場合= ∫ ∂P / ∂ρ・H(ρ - r) dρ← 屈折 CT の場合
ただし、G(r) と H(r) は事前に解析的に計算した「再構成フィルタ」

88 / 124

(3) そこで、xp2tg_dpc（および、吸収 CT 用の xp2tg）が畳み込み演算で計算する FP の値を出力するプログ
ラム xp2fp_dpc（および、xp2fp）を書いてみました。
xp2fp XP.bin {Dr} FP.bin
プログラム his2xp や hp2xp などで得た吸収 CT 用の投影のバイナリデータのファイル XP.bin
を読み込み、それぞれの投影画像の横一列の投影値と再構成フィルタ G（正確には Chesler
filter）との畳み込み演算を行う。この結果の filtered projection の画像（のシーケンス）を
XP.bin と同じ構成のバイナリファイル FP.bin に書き込む。なお、検出器の間隔 Dr の指定を省略するとそれは値１であると見なされる。また、xp2tg や xp2tg_dpc の場合と同じ環境変数
THREADS で使用するスレッドの個数を指定してやれば、投影画像の各横列の畳み込み演算をマルチスレッド処理で並列に実行することもできる。
csh や tcsh 環境の場合
setenv THREADS スレッドの個数（１～）
sh や bash
THREADS=スレッドの個数 ; export THREADS
orexport THREADS=スレッドの個数 ← bash のみで有効？
Windows のコマンドプロンプト
set THREADS=スレッドの個数
xp2fp_dpc XP.bin {Dr} FP.bin
XP.bin から読み込んだ投影値の微分値の画像に対して、[1] 微分値からの投影値の線積分[2] その結果の投影値と再構成フィルタ G の畳み込みを同時に行える DPC（diiferential phase contrast；微分位相）用の再構成フィルタ H を用いて
xp2fp と同様な処理を行い、その結果を FP.bin に書き出す（この FP.bin の値は微分値ではな
いことに注意）。
これらのプログラムのソースコードや 64 bit Windows 用の実行ファイルなどを SIXM 用の２個の書庫ファ
イルの両方に入れておきました：
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(4) FP（filtered projection）の値から計算した SP（sum of projection）の値は投影値や CT 値から計算した値と同じになるはずですが、ここでは便宜上、それを SFP（sum of filtered projection）と呼ぶことにします。
プログラム xp2sp を使えば xp2fp と xp2fp_dpc で得た FP（FP.bin）から SFP を計算することができます。
その具体的な手法は、...、説明するまでもありませんね。
注
先にも書いたように xp2fp_dpc が出力した FP は微分値ではないので、xp2fp と xp2fp_dpc で計算した FP のどちらに対しても常に xp2sp を使って SFP を計算して下さい。
(5) 以前と同じ７測定・13 セットの X 線 CT 画像の SFP を計算し、昨日の E-mail の sp.pdf と同様にして FPの画像（正確には、シーケンスの１枚目の画像）と SFP の画像を並べてみました。この E-mail に添付した sfp.pdf を御覧下さい。
(6) まず、赤枠で囲った SIXM で撮影した屈折 CT の SFP の画像についてですが、これは sp.pdf のものと概ね同じになっています。「前方差分と後方差分の平均」の投影値の和（SP）は Pfeiffer の方法で計算した

89 / 124

FP の和（SFP）と概ね同じ空間分布になっていると言えますが、現時点ではその理論的な理由は不明です。とりあえずは、これらの２つの方法は同じ結果をもたらしてくれるようです。
(7) 問題は屈折 CT のもの以外の SP と SFP の画像の違いです。sp.pdf と sfp.pdf を比較すると理想的な投影データ（020607a）、良質な画像を再構成できた投影データ（040711j）と SIXM の吸収 CT 用の投影データ（"15060?_a*"）から計算した SP と SFP の画像は概ね同じになっています。しかし、青枠で囲った
FZP-CT の測定 150616e と 150619[a,b]の３個の SFP の画像では SP の画像にない以下のような「縦黒線」と「横黒縞」が発生しています。
縦黒線
３個の SFP の画像のそれぞれで発生しているが、同じサンプルを X 線エネルギーを変えて撮影した 150619[a,b] では同じ位置にある。その横位置（== 投影画像のシーケンスの順番、０～）は以下の通り：
150616e（SIXM の 150602d と同じ石英サンプルを撮影）
360～560 と 1260 ～1460 のそれぞれの位置に２本づつある。
150619[a,b]（SIXM の 150603j と同じ AV 貝を [7,8] keV で撮影）
0～149 と 1700～1799 に１本づつ、800～1100 に２本ある。
横黒縞（縦幅が非常に大きいので「縞」と呼びます）
150619[a,b] でのみ、投影画像のシーケンスすべてのスライス 330～800（総スライス数は1404）が黒くなっている。ただし、その黒さは 150619a の方が「濃い」。
(8) まず縦黒線の原因を知るため、それが発生しているサンプル回転角の投影画像の FP の画像を抽出しました。ただし、FZP-CT の投影画像は巨大なので動画にはせず、それぞれのサンプル回転角の画像
ファイルを計算機 vrm の以下のディレクトリにコピーしておきました：
vrm:/media/disk/bl47xu/150616_sfp/
150616e_fp/*.tif と 150619[a,b]_fp/*
それぞれの測定の FP の画像（"*" は前記の範囲のシーケンス番号）150616e_sfp.tif と 150619[a,b]_sfp.tif
それぞれの測定の SFP の画像（sfp.pdf のものの元画像）
注
このディレクトリには以下の画像もコピーしてあります。
150616e_p_0.tif と 150619[a,b]_p_0.tif
それぞれの測定の投影画像の最初のもの（０度の投影画像）
150616e_sp.tif と 150619[a,b]_sp.tif
それぞれの測定の SP の画像（sp.pdf のものの元画像）
これらの FP の画像を観察すればわかるように、多面体状のサンプルの大きな表面が X 線の照射方向と平行になっている撮影方向付近において SFP の画像の黒い縦線が発生しています。今朝の三宅さん
あての E-mail で「サンプルによる屈折・回折は（屈折 CT では）問題にならない」と書きましたが、これは吸収 CT には当てはまりません。サンプルの面の方向は吸収 CT の再構成画像の画質を左右します。
FZP-CT と同じサンプルを撮影した SIXM の吸収 CT 用のデータの場合、ノイズレベルが高いのでそれが目立たないだけだと思われます。やはり、X 線 CT では円柱状のサンプルがベストですね。
(9) 次に 150619[a,b] の横黒縞についてですが、FP の画像の観察からはその原因を特定できませんでした。ただし、そのスライスの範囲は FZP の「ピンホール」の位置に一致しているので、それによる「遮光」

90 / 124が原因かもしれません。前記の投影画像 150619[a,b]_p_0.tif を御覧下さい。測定 150619[a,b] の視野の四隅には X 線が来ていません。実験の場でも言いましたが、このような X 線が来ていない、入射光強度 I0 = 0 かつ透過強度 I = 0 の部分に対して、ぼくのプログラムでは投影値０で埋めて画像再構成をするようになっています。問題はその部分の測定値に暗電流起源のノイズが含まていることで、そのため
I0 = I = 0 にならず、投影値が点々状（スペックル）ノイズ的な空間分布になっている可能性があります。
これは投影画像では目立ちませんが、空間周波数に比例した再構成フィルタとの畳み込み演算を行った結果の FP の画像では強調されるかもしれません。この時、X 線が来ていない部分を含むスライスのSFP は高い値になり、それよりも低い、正しい FP の値の部分（150619[a,b]_sfp.tif の横黒縞）は SFP の
画像上で相対的に低い輝度で表示されている可能性があります。これは SP と SFP の「絶対値」の比較によってわかることですが、...。
(10)以上のような SP と SFP の解析から得た結果をまとめてみました：[1] SIXM の屈折 CT の画像再構成で用いている Pfeiffer の方法は SP と SFP から見る限り問題なし（これは画像再構成でも問題なしですが）。
[2] サンプルの面の方向による X 線の屈折・回折、および RID の空間勾配の急変は SP や SFP に影響を与える。つまり、これらは再構成画像の画質に影響するので、サンプルは円柱状に整形しておくのがベスト。
[3] X 線 CT 測定ではサンプルが視野からはみ出さないように設定するのはもちろんのこと、X 線が来ていない部分も無くすようにすべき。
長い E-mail になりました。とりあえず以上です。

91 / 124

添付ファイル sfp.pdf

92 / 124

Date:
From:
To:
Cc:

Wed, 14 Oct 2015 12:21:23 +0900
Tsukasa NAKANO
Akira TSUCHIYAMAMIYAKE, Akihisa TAKEUCHI, Kentaro UESUGI, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: Re: filtered_projection（FP）の計算プログラム

つちやまさま、
みやけさま、
なかのです。
> > 添付の位相像で左から右にかけて暗くなるのは、RID の空間勾配や X 線の回折・屈折の> > 影響なんですか？
> > こうした影響の場合、このようにグラディエントに変わらない気がするんですけど・・
添付の位相像は再構成画像ですよね。RID の空間勾配∂δ/∂r の「エラー」はサンプルを回転したそれぞれの角度θで観測したデータのものであり、(r,θ) と再構成画像上の座標系 (x,y) の関係は
r = x * cosθ + y * sinθです。そして、この関係式をもとにしてδ(x,y) を再構成しているので、∂δ/∂r のエラーはδ(x,y) ではθ 方向に伸びる線状の異常になるはずです。位相像が左から右にかけて暗くなるのなら、再構成画像の上下がθ
= 0 の方向なので、投影画像の最初の一枚が全然ダメ（r 方向に系統的に変化している投影値？）ということですね。しかし、7/17 の E-mail で紹介したθ = 0 の投影画像
vrm:/media/disk/tsukasa/150602/150603b/rxp.tif
は異常な感じではありません。再構成画像の左右方向の RID の系統的な変動は実在する密度勾配ではないのですか？とり急ぎ、
On Tue, 13 Oct 2015 19:35:01 +0900 Akira TSUCHIYAMA wrote:> 三宅さん>> 2015/10/13 8:45、MIYAKE のメール：>> > 土`山様> > 僕に尋ねられても・・・・。
>>> >> 円柱の方が確かに良さそうですけども、似たような影響が出そうなものに、> >> 結晶質試料の場合、回折の効果もあります。
> >> 現状では四角くて試料のほとんどが結晶質の流体包有物のような> >> ものでも微妙なコントラストの識別が出来ているので、> >> 致命的と言うほどではないような気がします。
>>> > という、上杉さんの意見に従うと、ですね。
>>> > 添付の位相像で左から右にかけて暗くなるのは、RID の空間勾配や X 線の回折・屈折の> > 影響なんですか？
> > こうした影響の場合、このようにグラディエントに変わらない気がするんですけど・・>> うーーん、竹内さん、上杉さん、中野さん、どうなのでしょうか？
>> >> もっともすべてのサンプルを円筒状にするのは大変でしょうから（高山さんにもお願いするとして）、

93 / 124

> >> 10/14 のセミナー後にでも、サンプルの打合せできないでしょうか？
>>> > １７時から JEOL と打合せなので、それまでなら。
>> お願いします。とにかく、標準サンプルの種類だけでも決めましょう。
>> ---------------------------------------------------------------> >> 三宅さん
> >>> >> 2015/10/12 12:28、Miyake のメール：> >>> 土`山様> >>> 中野さんと上杉さんのメールを読んでの、僕の理解では、> >>>> >>>> 中身までなかなかついていけてないのですが、12 月の SPring-8 の SIXM のビームタイムの> >>>> 標準サンプルですが、FIB 加工は円筒形にした方がいいということですよね。
> >>>> >>> 円筒形がベストで有るけど、時間がなければ、致命的ではないので通常のままでよい。
> >>> >> 添付したパワポのファイルを見て下さい。V 貝の AV_out01（梅-in 150603b）の例です。
> >> サンプルの左側の空気の部分は右側の空気の部分より明らかに明るいです。また、> >> サンプル自身も左側の方が右側より明るいように見えます。これが「サンプルの面の> >> 方向による X 線の屈折・回折、および RID の空間勾の急変」に起因するものですよね。
> >>> >> サンプルの真ん中だけの領域（cut：赤色）、サンプルと空気を含む左側の領域（cut-l：空色）、> >> サンプルと空気を含む右側の領域（cut-r：緑色）について、それぞれの RID のヒストグラムを> >> みて下さい。cut-l の空気の部分は RID>0 ですが、cut-r の空気の部分は RID<0（ヒストグラムには> >> 現れていないが）となっています。また、サンプルの部分も、RID(cut-I) > RID(cut) > RID(cut-r) と> >> なっており、これは標準サンプルについては致命的なのではないでしょうか？
> >>> >> ---------------------------------------------------------------> >>>> 三宅さん：> >>>> 現実的に（時間的に）大丈夫でしょうか？
> >>>> >>> 何を、どれくらい、作る必要があるかによりますが、今からなら大丈夫では？
> >>> ただし、僕が教える時間がどれくらいあるのかは不明。
> >>> >> もっともすべてのサンプルを円筒状にするのは大変でしょうから（高山さんにもお願いする> >> として）、10/14 のセミナー後にでも、サンプルの打合せできないでしょうか？
> >>>> >>> ----- Original Message ----> >>>> 中野さん
> >>>> c.c.に特推の研究員で京大に 10 月から来られた高山さんを追加しました。
> >>>>> >>>> 中身までなかなかついていけてないのですが、12 月の SPring-8 の SIXM のビームタイムの> >>>> 標準サンプルですが、FIB 加工は円筒形にした方がいいということですよね。
> >>>>> >>>> 三宅さん：> >>>> 現実的に（時間的に）大丈夫でしょうか？

94 / 124

Date:
From:
To:
Cc:

Wed, 14 Oct 2015 12:44:33 +0900
Tsukasa NAKANO
Akira TSUCHIYAMAMIYAKE, Akihisa TAKEUCHI, Kentaro UESUGI, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: Re: filtered_projection（FP）の計算プログラム

つちやまさま、なかのです。これはサンプルの表面（再構成画像の右上と左下）に付着している白金？起源の偽像のようですね。とり急ぎ、On Wed, 14 Oct 2015 12:35:34 +0900 Akira TSUCHIYAMA wrote:> 中野さん、竹内さん>> 再構成画像の左右方向の RID の系統的な変動は実在する密度勾配ではないです> （別の場所ですが、SEM でみても、こんな密度勾配は存在しません）。
> 原因は僕にはわかりませんが、とにかく内部のコントラストが周りの空気の. （実在しない）コントラストに引っ張られています。
>> 11/27-12/2 の撮影では、同じサンプルで円筒と直方体の両方も試してみたいと思っています。
Date:
From:
To:
Cc:

Wed, 14 Oct 2015 13:08:22 +0900
Tsukasa NAKANO
Akira TSUCHIYAMAMIYAKE, Akihisa TAKEUCHI, Kentaro UESUGI, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: Re: filtered_projection（FP）の計算プログラム

つちやまさま、なかのです。すみません。先程の「白金？」は FIB で爆撃に使っているイオンの残骸のことです。RID の勾配をもとにした位相 CT ではこの手のものが偽像となります。吸収 CT の場合を考えてみて下さい。内部に非常に高い LAC を持つ点状のものがある時、再構成画像上でその点の「影」ができます。これはその点を光が十分通っていない（== その光路の投影値が異常である）ために生じます。一方、位相 CT では RID の勾配が大きな部分（== モノの表面）がエラーの原因になります。そして、RID の空間勾配から RID を線積分するので、この大きな空間勾配はそれを挟む広がりを持った RID の異常になります。つまり、吸収 CT の再構成画像では線状の偽像となるものが位相 CT の画像では幅広のバンド状の偽像になるようです。今回の位相 CT の再構成画像上の空間勾配はサンプル表面の FIB イオンによるものではないでしょうか？ これが正しいなら、位相 CT に使う FIB サンプルの表面は十分に「掃除」しておかねばなりません。とり急ぎ、

95 / 124

Date:
From:
To:
Cc:

Wed, 14 Oct 2015 10:56:56 +0900
Tsukasa NAKANO
Kentaro Uesugi, Akihisa TAKEUCHI"TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: xp2tiff

みなさま、
GSJ/AIST のなかのです。言うのを忘れていました。SIXM 用の his2xp や xp2tg などのプログラムが使っている投影画像のシーケンスが格納されているバイナリファイル XP.bin から指定したシーケンス番号（０～）のものだけを抜き出し、それぞれを個別の TIFF 形式ファイルに変換するプログラム xp2tiff を書きました。これは以前に紹介した「XP.bin の上の最初の１枚の画像だけを抜き出す a2tiff」の上位互換プログラムです。
(0) xp2tiff（と a2tiff）のソースコードや 64 bit Windows 用の実行ファイルを SIXM 用の２個の書庫ファイルの両方に入れておきました：
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(1) a2tiff の使い方は以下の通りです（前にも説明した気がしますが、念のため）。
起動法
a2tiffXP.bin

{base

step}

BPS TIFF

起動パラメータXP.bin投影画像（のシーケンス）が入っているファイルの名前。
base と step
変換の際に使用する投影値と TIFF の整数画素値の関係式
投影値 = base + step × 整数画素値の係数値。ただし、これらの指定を省略すると画像の上の投影値の最小値と最大値の範囲を正規化した整数画素値に変換する。
base = 投影値の最小値
step = （投影値の最大値 - 最小値) / 整数画素値の最大値
BPS
TIFF 形式画像の画素値のビット数（１～16）。
TIFF
TIFF 形式画像を格納するファイルの名前。
標準エラー出力にタブコード区切りで書き出される２行×２個の値
１行目：画像の横および縦画素数
２行目：画像上の投影値の最小値と最大値(2) a2tiff とは異なり、xp2tiff は指定した画像を TIFF に変換せずに、指定した画像それぞれの上の投影値の最小値と最大値だけを調べることもできます。
起動法
xp2tiff XP.bin
もしくはrangeList

96 / 124

xp2tiff

XP.binrangeList

{base

step}

BPS

TIFF_format

起動パラメータXP.bin投影画像のシーケンスが入っているバイナリファイルの名前。なお、このファイルの先頭行（最初の改行文字まで）に投影画像の横・縦画素数とシーケンスの総数がタブコード区
切りの ASCII 文字列として書き込まれており、それは UNIX の標準コマンド "head" などを使えば閲覧できる。
head -1 XP.bin
rangeList
XP.bin の上の抽出したい画像のシーケンス番号（０～）を以下のような書式で指定する。
"-" もしくは "0-"すべての画像を抜き出す。
"0"
シーケンス番号０の画像だけを抜き出す。
"0-99"
０～99 の 100 枚の画像を抜き出す。
"0,2,4"０、２および４の画像。
"0-2,4"０～２および４の画像。
"1-"
１以上の画像すべて。
base、step および BPS
a2tiff の起動パラメータの意味と同じ（省略時の動作も同じ）。
TIFF_format
抜き出したシーケンス番号ごとの画像を格納する TIFF ファイルの名前を "DIR/%04d.tif"のような C 言語の標準関数 sprintf() 用の書式で指定する。ただし、この例の "DIR/"
は事前に自分で作成しておいた TIFF ファイルを格納するためのディレクトリの名前で、また、int 型整数値の変換書式 "%04d" により各画像のシーケンス番号は左０埋めの４桁の文字列（番号１なら "0001"）になる。
標準出力に書き出される画像ごとの行上のタブコード区切りの３個の値[1] シーケンス番号[2-3] その画像上の投影値の最小値と最大値
とりあえず以上です。

97 / 124

Date:
From:
To:
Cc:

Fri, 16 Oct 2015 11:49:11 +0900
Tsukasa NAKANO
Kentaro Uesugi, Akihisa TAKEUCHI"TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: 位相 CT シミュレータ

みなさま、なかのです。位相（屈折）CT の画像再構成プログラム xp2tg_dpc で使っている Pfeiffer 論文（PRL、98、10815、2007）の再構成フィルタ（投影データとして与えられた RID の空間勾配の線積分を畳み込み演算の際に行う再構成フィルタ）の性質をイマイチ理解できなかったので、普通の吸収 CT 用の投影データの値を数値的に微分して位相 CT で撮影したものと概ね等価な微分位相の投影データを合成するプログラム
xp2dp を書いてみました。
(1) xp2dp のソースコードや 64 bit Windows 用実行ファイルなどは SIXM 用の２個の書庫ファイルの両方に入れてあります。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip(2) xp2dp の起動法は以下の通りです。
xp2dp

XP.bin{Dr}

DP.binXP.binプログラム his2xp や hp2xp などで作成した吸収 CT 用の投影画像のバイナリデータファイルの名前。"-" を指定すると標準入力からデータを読む（ただし、これは Windows では NG）。
Dr
検出器の間隔。指定を省略すると値１だと見なされる。
DP.bin検出器 i の投影値 P(i) から「中心差分」で微分したその空間勾配 D(i) のバイナリデータを書き込むファイル名。"-" を指定すると標準出力にデータを書き込む（Windows では NG）。
注
上記の中心差分はプログラム xp2sp_dpc で行った投影値の「前方差分と後方差分の平均」と一致する微分の差分近似であり、また、これは「線形補間した測定値」を使った差分近似とも
言えます。
D(i)
= { P(i+1) - P(i-1) } / (2 × Dr)← 中心差分の定義= [ { P(i+1) - P(i) } / Dr + { P(i) - P(i-1) } / Dr ] / 2← 前方差分と後方差分の平均= [ { P(i+1) + P(i) } / 2 - { P(i) + P(i-1) } /2 ] / Dr
← 線形補間した測定値を使った差分(3) xp2dp で得た投影値の空間勾配のデータから xp2tg_dpc で再構成した画像と、xp2dp に入力した投影データから xp2tg で再構成した画像を比較しました。ただし、以下の測定の投影データを使いました。
On Thu, 01 Oct 2015 16:26:24 +0900 Tsukasa NAKANO wrote:98 / 124

>
>
>
>
>
>
>
>
>

020607a
BL47XU で撮影した金沢大学・森下くんの simplectite の CT 画像
http://www-bl20.spring8.or.jp/~sp8ct/hd2/morishita_symplectite/020607a/byte_222.gifから X 線 CT シミュレータで合成した、SPring-8 の通常の測定のものと同形式になっている測定データ。これは合成したデータなので ...
040711j
BL20B2 で撮影した直径が 1 mm 以下の均質な球形ビーズの CT 画像
http://www-bl20.spring8.or.jp/~sp8ct/tmp/040711j.gifの測定データ。上記のような良好な画質の画像を再構成できている。
具体的には、以下のような端末入力を行いました。
# 投影データ xp.bin の作成（ただし、１枚のスライス用だけを作成した）。
dit_tbl 040711j/raw | hp2xp 359 359 xp.bin# そのデータから普通に画像再構成
xp2tg xp.bin 0.000583 498.5 0 0 0.02930016 8 .
mv 0.tif 040711j_359.tif← カレントディレクトリ "." に作成された１枚のスライス画像を改名# 投影値の空間勾配のデータを作成し、それから画像再構成# Windows の場合は作業用ファイル（dp.bin）を使う必要がある。
xp2dp xp.bin Dr1 dp.bin
xp2tg_dpc dp.bin Dr2 498.5 0 0 0.02930016 8 .
rm dp.bin# Linux なら空間勾配のバイナリデータをパイプに流すことができる。
xp2dp xp.bin Dr1 - | ← 改行しない
xp2tg_dpc - Dr2 498.5 0 0 0.02930016 8 .
# いずれの場合も最後にスライス画像のファイル名を改名した。
mv 0.tif 040711j_359_dp.tif
ただし、xp2dp を使った処理では Dr1 と Dr2 のいずれか一方を１、他方の値を実際の検出器間隔（上の例では 0.000583 cm）にして下さい（その理由の説明は長くなるので省略します）。

(4) 以上のようにして作成した２測定の２枚ずつの再構成画像をこの E-mail に添付した PDF ファイル
dp.pdf に並べました。左が投影データから直接再構成したもの、右が xp2dp で得た微分位相のデータから再構成したスライス画像です。それぞれのスライスに出現した CT 値の値域を各画像の上部に記しておきました。左に比べて右の画像の CT 値の値域の幅が広いのは xp2dp による微分や xp2tg_dpc によるその回復処理に伴う計算誤差のためだと思われます。しかし、見た目には、左右の画像は概ね同じに見えます（CT 値と表示輝度の対応関係は左右の画像で同じにしてあります）。つまり、位相 CT シミュ
レータ（xp2dp）はそれなりに使えそうです。
とりあえず以上です。

99 / 124

添付ファイル dp.pdf

100 / 124

Date:
From:
To:
Cc:

Fri, 16 Oct 2015 17:41:18 +0900
Tsukasa NAKANO
Kentaro UESUGI
Akihisa TAKEUCHI, "TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: Re: 位相 CT シミュレータ

うえすぎさま、なかのです。E-mail ありがとうございます。
> xp2tg_dpc で再構成した場合、空気のあたりはきちんとゼロになるのでしたっけ？
> （Pfeiffer method で積分時の定数の取り扱いがどうなっているか忘れてしまいました）(1) 前に紹介したファイル pfeiffer.txt に書きましたが、Pfeiffer の方法では投影値 p(r) の Fourier（以下では"F" と略す）変換 P(k) を部分積分し、P(k) = ∫ p(r)・exp(-2・π・i・k・r) dr← F 変換の定義式= [ p(r)・exp(-2・π・i・k・r) / (-2・π・i・k)]_r=±∞- ∫ (∂p/∂r)・exp(-2・π・i・k・r) / (-2・π・i・k) drr =±∞ で p(r) = ０と仮定することにより、P(k) を「投影値の微分∂p/∂r の F 変換」に「フィルタ1/(2・π・i・k)」 を掛けた式に書き直します：
P(k) = ∫ (∂p/∂r)・exp(-2・π・i・k・r) dr / (2・π・i・k)
投影値 p(r) を使った画像再構成では「その F 変換 P(k) に対して周波数空間において再構成フィルタ
G(k) = |k| を掛けること」に相当する畳み込み演算を行うわけですが、上式より、∂p/∂r を用いた画像
再構成では
H(k) = G(k) / (2・π・i・k) = - i × sgn(k) / (2・π)ただし、sgn(k) は k の値の符号に応じて±1 を返す符号関数
を再構成フィルタとして使えば良いことになります（∂p/∂r に対する畳み込み演算のために H(k) の逆
F 変換をする必要がありますが、それは枝葉の話）。
(2) 以上のような Pfeiffer の方法で仮定しているのは「r = ±∞ で p(r) = ０」だけです。実際の画像再構成
処理では「r = ±∞」は「検出器列の両端」であり、そこでの p(r) =０の成否が Pfeiffer の再構成フィルタを用いた場合に問題になります。ただ、これは普通の画像再構成でも仮定していることなので、２つの手法で仮定していることの違いによる再構成画像の違いはないはずです（理論的には）。
> ところで> 位相 CT シミュレータ> は> 微分位相 CT シミュレータ> の方が名称としては正確ですかね。
(3) どちらの名称でもかまいません。ただ、古典物理的な観点から言うと、屈折角のような微分位相に関連した値以外の観測データを使う位相 CT はありえないのではないでしょうか？ つまり、Fermat の原理（光線は屈折率の線積分が最小となる経路を伝播する）に対する Euler 方程式（光線方程式）をもとにした位相 CT の観測データは屈折角を介した「微分位相の経路積分の値」です。それゆえ、位相 CT ==微分位相を観測値として使う CT なのでは？（Fermat の原理以外の物理的原理に基づく位相 CT はあ
るのか？）
とり急ぎ、

101 / 124

On Fri, 16 Oct 2015 14:12:10 +0900 Kentaro UESUGI wrote:> 中野さん>> 上杉です>> xp2tg_dpc で再構成した場合、空気のあたりはきちんとゼロになるのでしたっけ？
> （Pfeiffer method で積分時の定数の取り扱いがどうなっているか忘れてしまいました）>> ところで> 位相 CT シミュレータ> は> 微分位相 CT シミュレータ> の方が名称としては正確ですかね。
> xp2dp ですし
Date:
From:
To:
Cc:

Fri, 16 Oct 2015 18:14:40 +0900
Tsukasa NAKANO
Kentaro UESUGI
Akihisa TAKEUCHI, "TSUCHIYAMA, Akira", MIYAKE, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA
Subject: Re: 位相 CT シミュレータ

うえすぎさま、なかのです。先程の E-mail の (2) に書いた「Pfeiffer の方法による画像再構成で仮定していることは普通の画像再構成で仮定していることと同じ」は間違いでした。Pfeiffer の方法では「投影値の微分∂p/∂r の F変換」が可能だと仮定しています。つまり、r = ±∞（検出器の両端）で∂p/∂r = ０でないとダメです。それゆえ、Pfeiffer の方法を使った位相 CT の不完全 CT はヒドイことになりそうです。特に、スリットのため X 線の来ていない部分がある場合、そこでは暗電流起源の大きな∂p/∂r の値が発生する可能性がある。
とり急ぎ、

102 / 124

Date:
From:
To:
Cc:

Mon, 19 Oct 2015 20:04:16 +0900
Tsukasa NAKANO
Akihisa TAKEUCHI, Kentaro Uesugi
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: SIXM の屈折プロファイルの重心

たけうちさま、
うえすぎさま、
GSJ/AIST のなかのです。この E-mail に添付した PDF ファイル ri+rr+rp.pdf に示したような 8/19 の竹内さんからの E-mail で答えていただいた SIXM の refraction profile（SIXM のそれぞれのスキャンで測定した
暗電流補正済みの X 線強度分布）の「重心」に関する話です。
(1) 以前に質問した refraction profile の「左右非対称性」が「検出器への X 線の斜め入射」によって生じたと考え、SIXM の FZP による幾何学的な拡大を模した fan beam を使ったシステムの検出器で観測される
X 線強度プロファイル I(r) とその重心の位置 r0 の計算式を導出してみました。この E-mail に添付した
PDF ファイル fb.pdf を御覧下さい。そこに示したように、|r| << L（ただし、L はサンプルと検出器の間の距離）で、かつ、fan beam の X 線強度 J1(θ) が方位角θによらず一定値なら以下が成り立ちます。
[1] プロファイル I(r) は r に関して上に凸の放物線状になる。
[2] その重心 rc は X 線の検出器上の測定範囲のちょうど中央になる。
(2) サンプル内部で屈折する X 線はその過程でサンプルに吸収されます。その影響を知るため、この
E-mail に添付した PDF ファイル mb.pdf に示したような micro beam による測定を考えてみました。サンプル内部での X 線の r 方向への屈折は RID（δ）の空間勾配 －∂δ/∂r によって生じます。そして、RIDと LAC（μ）は概ね比例するので、∂δ/∂r < 0 なら∂μ/∂r < 0 になっているはずです。ここではそのような∂δ/∂r < 0 に対応する LAC の空間分布として r の一次式μ(r) = b－a・r を考えました（ただし、係数 a は LAC の r 方向の空間勾配の絶対値；a = |∂μ/∂r| ）。この場合、サンプルに入射した micro
beam 内部の r 方向の X 線強度分布が一様なら、サンプル透過後の強度は exp(a・S・r) に比例した値になります（ただし、S は X 線入射方向のサンプルの「厚み」）。そして、幅が R の micro beam に対して
積 a・S・R が十分小さい値なら以下が成り立ちます。
[3] 透過 X 線強度の空間分布の重心の位置（c）は屈折した micro beam のちょうど中央になる（micro beam の幅を R とすると、c = R / 2）。
(3) 以上の議論は現実的でしょうか？ また、その場合、竹内さんが SIXM の屈折角の推定に「屈折プロファ
イルの重心の値の差を使うこと」を推奨しているのは、上記の[2]と[3]のためでしょうか？ とり急ぎ、

103 / 124

添付ファイル ri+rr+rp.pdf

104 / 124

添付ファイル fb.pdf

105 / 124

添付ファイル mb.pdf

106 / 124

Date:
From:
To:
Cc:

Tue, 20 Oct 2015 09:45:17 +0900
Tsukasa NAKANO
Akihisa TAKEUCHI, Kentaro Uesugi
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira" <atsuchi@kueps.kyoto-u.ac.jp>
Subject: Re: SIXM の屈折プロファイルの重心

たけうちさま、
うえすぎさま、なかのです。昨日の E-mail の(2) に書いた「SIXM の吸収を考慮した屈折プロファイル」の議論に間違いを見つけました。以下の部分です。
> サンプル内部での X 線の r 方向へ> の屈折は RID（δ）の空間勾配 －∂δ/∂r によって生じます。そして、RID> と LAC（μ）は概ね比例するので、∂δ/∂r < 0 なら ∂μ/∂r < 0 になって> いるはずです。ここではそのような ∂δ/∂r < 0 に対応する LAC の空間分布> として r の一次式 μ(r) = b － a・r を考えました（ただし、係数 a は LAC の> r 方向の空間勾配の絶対値；a = |∂μ/∂r|）。
昨日の E-mail（に添付した PDF ファイル mb.pdf）では RID の負の空間勾配∂δ/∂r < 0 に対して LAC が負の空間勾配∂μ/∂r < 0 になっていることを仮定していましたが、重心の位置の「第一近似」の計算ではこれは不要でした。つまり、mb.pdf に記した重心 c の計算式
c = {1 / (a・S) }・F(a・S・R)
ただし、
F(t) = { (t - 1)・exp(t) + 1 } / { exp(t) - 1 }～ t/2← F(t) の「第一近似」の計算式
において分母と分子にある a = |∂μ/∂r| は打ち消しあいます。それゆえ、積 a・S・R（LAC の空間勾配とmicro beam の透過面積 S・R の積）が十分小さい値なら、∂μ/∂r の符号に関わらず、屈折プロファイルの重心の位置は幅 R の micro beam のちょうど中央に来ます（c = R/2）。とり急ぎ、

107 / 124

Date:
From:
To:
Cc:

Tue, 20 Oct 2015 10:26:42 +0900
Tsukasa NAKANO"TSUCHIYAMA, Akira"Akihisa TAKEUCHI, Kentaro Uesugi, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA, MIYAKE
Subject: Re: SIXM の屈折プロファイルの重心
つちやまさま、なかのです。先の E-mail に書いたように、LAC（RID ではありません）の空間勾配が大きい部分があるサン
プルでは、その部分に屈折角の測定エラーが発生し、それを使った位相 CT の再構成画像がおかしなことになる可能性があります。LAC の空間勾配が大きいのは粒子の表面（特に、高い LAC を持つ粒子の表面）ですから、結局、SIXM の位相 CT は石英の結晶のように LAC が低く、その空間分布が一様なサンプルの内部の微細構造の観察にしか使えないのかもしれません。とり急ぎ、
Date:
From:
To:
Cc:

Tue, 20 Oct 2015 15:19:54 +0900
Tsukasa NAKANO"TSUCHIYAMA, Akira"Akihisa TAKEUCHI, Kentaro Uesugi, Aki Takigawa, Junya Matsuno,
中村隆太, Takashi Sakurama, Yohei IGAMI, Akiko TAKAYAMA, MIYAKE
Subject: Re: SIXM の屈折プロファイルの重心
つちやまさま、なかのです。昼前の E-mail に書いた「SIXM で屈折角の測定エラーが発生する LAC の空間勾配（a）の大きさ」についてですが、その量的な見積もりができることに気づきました。
(1) 昨日の E-mail に書いたように、竹内さん推奨の屈折プロファイルの重心演算で SIXM の屈折角を推定するためには a・S・R << 1 が満たされている必要がある。ただし、S はサンプルの奥行き（～サンプルの径）、R はマイクロビームの幅。
(2) Dr を SIXM の吸収および位相の投影画像（の横方向）の画素の辺長とすると、サンプルの表面おける
LAC（μ）の空間勾配は概ね a = μ / Dr。そして、SIXM のスキャン幅はマイクロビームの幅と同程度の値（SIXM で「菊」のスキャンをした場合に相当？）なので、結局、a = μ/ R。
(3) 以上から、μ・S << １ なら SIXM の屈折角の推定に問題なし。より具体的には、サンプルの径（S）が100μm = 0.01 cm だと思えば、μ << 100（1/cm)。
(4) SIXM 測定で使った X 線エネルギー8 keV における LAC の理論値は以下の通り。
echo 8 | lac_cm density chemical_formulamaterial
density (g/cc)
chemical formulawater
1
H:2 O:1
quartz2.7
Si:1 O:2
calcite2.715
Ca:1 C:1 O:3
dolomite
2.86
Ca:1 Mg:1 C:2 O:6
forsterite
3.226
Mg:2 Si:1 O:4
fayalite
4.4048
Fe:2 Si:1 O:4
enstatite
3.21
Mg:1 Si:1 O:3
magnetite
5.2
Fe:3 O:4

LAC (1/cm)at 8 keV
10.372414
98.353975
204.278353
141.615899
103.982481
793.180449
107.496333
1166.599791と言うわけで、SIXM の正しい屈折角の推定は石英でも厳しいですね。とり急ぎ、

108 / 124

Date:
From:
To:
Cc:

Tue, 20 Oct 2015 18:40:22 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa"
上杉 健太朗, MIYAKE, Aki Takigawa, Junya Matsun, 中村隆太, "Takashi Sakurama",
Yohei IGAMI, Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: Re: SIXM の屈折プロファイルの重心

たけうちさま、なかのです。E-mail ありがとうございます。
(1) SIXM の吸収を伴う屈折に関して素朴な設定の計算ですみません。何せ、最初に「台地と山脈」のかたちをした屈折プロファイルを見てしまったので、マイクロビームは有限幅、かつ、X 線強度が一様ではないものだと思い込んでいました。
(2) 竹内さんが言われる通り「積 a・S・R が小さい値になる」は「マイクロビームの強度分布をデルタ関数と見なせる」に対応しているようですね。なお、ついでに言っておくと、積 a・S・R は無次元の値です。
とり急ぎ、
On Tue, 20 Oct 2015 16:06:55 +0900 "TAKEUCHI, Akihisa" wrote:> 中野様>> 竹内です。
> すみません、中野さんのメールを自分なりに理解し、考えをまとめてメールにするのに> 時間がかかり返答が遅くなりました。
>> 重心演算を RID 導出の根拠としている理由については、添付ファイルを参照してください。
> 自分用の備忘録的なものなので色々関係ないこともかいてあったり、読みづらいかもしれないですが> 最初の 2 ページが RID に重心演算を使う根拠に関する内容になっています。
> 一応下に要約も書きました。
>> もっとかいつまむと、> 検出器で得られる信号は、試料を透過したｘ線の複素振幅のフーリエ変換の絶対値自乗である> という考えから導かれてます。
> ひとつだけ近似があり、試料を透過するｘ線の幅が十分に小さい、という近似をしています。
> 試料を照射するビームの強度分布（点像分布関数）P(x)と試料の透過率分布 T(x)に対して、> 測定で得られる「透過像」P＊T が（＊はコンボリューション）、P＊T~ T が近似的に成り立つ、> つまり P(ｘ）をデルタ関数として扱えるという近似を使っており、これが成り立つ場合に、> 重心演算から RID が計算できます。
> なのでこの近似が成り立たない場合は、重心演算では正確な RID が得られないことになります。
> 「近似が成り立たない場合」に関しては実はまだちゃんと考えたことがありません。
>> 中野さんの計算のアプローチは、自分でやってみたことはなかったのですが、> 中野さんの提案されている「積 aSR が十分小さい」が、ひょっとしたら上記の> 「P をデルタ関数近似できる」と同義として考えられるのかもしれません。
>> うまく説明できてないようでしたらすみません。

109 / 124

Date:
From:
To:
Cc:

Wed, 21 Oct 2015 14:51:24 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", 上杉 健太朗
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, "Takashi Sakurama", Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: Martiono 論文

たけうちさま、
うえすぎさま、なかのです。そう言えば思い出しました。竹内さんが SIXM の位相画像の y or z 方向の「スムーシング」で使われている Martino 達の論文（この E-mail に添付しました）についてです。この論文にはイチャモンをつけたい細かいことが多々ありますが、ぼくが最も理解に苦しんでいるのは論文の最後の部分の式 (13) と(14) についてです。これらの式が言っている「トレードオフ」の意味がイマイチ理解できません。
詳しく言うと「位相画像の x 方向の微分値と微分位相の測定値の食い違いの自乗値 |∂φ/∂x - g|^2」と「位相画像の y 方向の微分値の自乗値 |∂φ/∂y|^2」を対等に扱うことが理解できない。これは、位相画像の y 方向の微分値が概ね ０、と仮定することになるのでは？
ところで、Fermat 原理から SIXM の屈折角を支配している「光線方程式」を導出するために「変分法」を復習していた時に、x 方向の微分位相の測定データから位相画像を推定できる、Martino 達のものとは別の手法を思いつきました。いわゆる「spline 関数」による推定法です。この E-mail に添付したファイル spline.pdf を御覧下さい。Martino 達の手法と併記するかたちで spline 関数による位相画像の推定法を記してあります。
ここで言っている spline 関数は「spline」の語源通りのもので、物理的には「薄板の曲げ」を支配する重調和（biharmonic）方程式を満たす関数です。与えられたデータを spline 関数で fitting すると、その曲面の「平均曲率の自乗和（total curvature）」は（同様な fitting で得た曲面の中で）最小になります。このE-mail に添付したファイル briggs.pdf にその証明が記されています（Briggs が 1974 年の論文で行った証明をぼくが別のやりかたで証明しました）。
１次元の重調和方程式は４階の常微分方程式 d^4 y / dx^4 = 0 です。そしてこの方程式を満足する１次元の spline 関数は３次までの多項式 y = a * x^3 + b * x^2 + c * x + d となります。一般に使われている spline 関数はデータの観測点を端点とする区間ごとに係数値 a～d を決めた、データの値を通る３次までの多項式ですが、この方法は２次元以上のデータに対しては使えません。その場合には spline
関数の原理に戻って考える必要があります。
このような spline 関数による位相画像の推定法には２つの利点があります。
[1] Martino 達の方法のようなトレードオフ・パラメータが不要。
[2] 人に説明する時に「spline 関数を使った」と言えば OK。
この手法をどうぞお試し下さい。
P.S.
実は、ぼくはこの新しい手法を実際に試していません。Martino 論文の手法を使っている竹内さんのプログラム HIDEKI なら、ちょっとした改造で spline 関数による新しい手法を試すことができると思います。

110 / 124

添付ファイル spline.pdf

111 / 124

添付ファイル briggs.pdf

112 / 124

Date:
From:
To:
Cc:

Mon, 23 Nov 2015 15:06:00 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", 上杉 健太朗
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, "Takashi Sakurama", Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: SIXM の位相回復プログラム

たけうちさま、
うえすぎさま、
GSJ/AIST のなかのです。10/21 の E-mail で紹介した Martino 論文の方法と、それに準じた spline 関数による SIXM の微分位相の画像の位相回復法に関する文書を修正しました。この E-mail に添付した
spline.pdf です。その２ページ目に記した数値計算法を用いて SIXM で測定した x 方向の微分位相の画像（のシーケンス）に対して位相回復を行う２個のプログラムを書きました：
dp2xp DP.bin { δx δy } λ XP.bin
Martino 論文の方法で位相回復を行うプログラム
dp2xp_spline DP.bin { δx δy } XP.bin
spline 関数を用いた手法で位相回復するプログラム
ただし、これらのプログラムの起動パラメータの意味は以下の通りです。
DP.binSIXM などで得た x 方向の微分位相画像（のシーケンス）のバイナリデータが入っているファイルの名前。
δx および δy
微分位相の画像の x および y 方向の画素の辺長。これらの値の指定を省略すると δx = 1 かつδy
= 1 と見なされる。
λMartino 論文の手法で使う trade-off parameter（TOP）。値０を指定すると単純な x 方向の線積分
で位相回復する。
XP.bin位相回復した結果の位相画像（のシーケンス）のバイナリデータを書き込むファイルの名前。
また、これらのプログラムが仮定している位相画像の境界条件（x 方向に point-refractive symmetry；PRS）に対応した x 方向の微分位相の画像を合成する、10/16 の E-mail で紹介した「位相 CT シミュレータ」用のプログラム xp2dp と概ね同じ処理を行うプログラム xp2dp_prs も書きました（xp2dp_prs の起動パラメータの意味は xp2dp のものとまったく同じです）：
xp2dp_prs

XP.bin{ Dx }

DP.bin

これら３個のプログラムのソースファイルや 64 bit Windows 用の実行ファイルを従来と同じ SIXM 用の２個の書庫ファイルの両方に入れておきました。
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip
なお、今回の位相回復プログラムでは Fourier 級数展開を行なっていますが、位相画像に対する幾何学的仮定を遵守するため FFT を使っていません。それを行列の乗算で実行しています。そのため、例えば縦横画素数がともに N の画像に対しては N^3 に比例した処理時間を要します（FFT なら (N * logN)^2～N^2）。

113 / 124

自作は面倒だったので、その処理時間軽減のためにマルチスレッド処理が可能な線形演算用ライブラリ
BLAS（Basic Liner Algebra Sub programs）を実装した以下の２種類のパッケージを利用しました（ただし、今回利用したのは BLAS の汎用行列乗算サブルーチン dgemm だけです）：
Intel C++ Compiler に付属の MKL（Math Kernel Library）の BLAS
何も指定しなくても計算機の CPU の物理コアの個数と同数のスレッドを用いてマルチスレッド処理を行ってくれる。また、使用するスレッド数を環境変数 MKL_NUM_THREADS で陽に指定することもできる。
OpenBLAS（http://www.openblas.net/)
Linux や Windows などで使用可能（ぼくのところでは MinGW64 用のものを利用）。環境変数
OPENBLAS_NUM_THREADS でスレッド数を指定しないとダメ。
そして、Intel MKL や OpenBLAS のライブラリのリンクが必要だったので、今回の３個のプログラムは前記の書庫ファイルに入れた以下の２個の「Makefile」のいずれかでコンパイルするようにしてあります：
Makefile.dp
Linux で Intel C++ Compiler と MKL を使うことを想定。変数 BLAS の設定を変更すれば Linux に付属の「reference BLAS」を使うことも可能。
Makefile.dp.mingw
Windows の MinGW64 で GNU C-compiler と OpenBLAS の利用を想定。OpenBLAS ライブラリを静的にリンクする設定にした（そのため、64 bit Windows 用の実行ファイルが巨大になりました）。
さて、プログラム xp2dp_prs で「Lenna 画像（画素数 512^2）」の微分位相の画像を合成し、それを dp2xp とdp2xp_spline で位相回復してみました。その結果はこの E-mail に添付した lenna.pdf に示した通りです。
lenna
もとの「Lenna 画像」
lenna_dp_prs
xp2dp_prs で合成した微分位相の画像
lenna_spline
spline 関数を用いた方法で位相回復した画像
lenna_λ（λ= 0、0.5 or 1）
TOP の値λを変えて Martino 論文の方法で位相回復した画像
注
微分位相の画像以外は位相値と表示輝度の対応関係を揃えてあります。また、画像の下の文字列は出現した位相値や微分位相の値の範囲です。
これら６枚の画像は 64 bit Windows 機上で以下のようにして作成することができます。どうぞお試し下さい。
[1] 前記の書庫ファイルからディレクトリ his/を解凍・展開[2] his/の下にこの E-mail に添付したテキストファイル lenna.bat.txt を lenna.bat に改名してコピー[3] his/の下にこの E-mail に添付した bzip2 で圧縮したバイナリファイル lenna.bin.bz2 を解凍したファ
イル lenna.bin をコピー
[4] コマンドプロンプトを開いて his/に移動し、"lenna.bat" を実行
注
バッチファイル lenna.bat は４スレッドのマルチスレッド処理で行列を乗算する設定にしてあります：
set OPENBLAS_NUM_THREADS=4
lenna.pdf を眺めればわかるように spline 関数による位相回復は全然ダメです。その原因は、...。また、もとの Lenna 画像はノイズ無しなので TOP の値λを変えた位相回復の良否は判断できませんが、少なくとも λ=0 の場合の単純な線積分による位相回復は dp2xp で概ね正しくできていることがわかりました。
とりあえず以上です。

114 / 124

添付ファイル spline.pdf

115 / 124

116 / 124

添付ファイル lenna.bat.txt@echo off
path=.;%path%set OPENBLAS_NUM_THREADS=4
if not exist lenna.bin (echo lenna.bin : file not found.
exit /b)echo lenna
a2tiff lenna.bin 0 1 8 lenna.tif
echo lenna_dp_prs
xp2dp_prs lenna.bin lenna_dp_prs.bin
a2tiff lenna_dp_prs.bin 8 lenna_dp_prs.tif
echo lenna_spline
dp2xp_spline lenna_dp_prs.bin lenna_spline.bin
a2tiff lenna_spline.bin 0 1 8 lenna_spline.tif
for %%L in (0 0.5 1) do (echo lenna_%%L
dp2xp lenna_dp_prs.bin %%L lenna_%%L.bin
a2tiff lenna_%%L.bin 0 1 8 lenna_%%L.tif)
添付ファイル lenna.bin.bz2
これは「Lenna 画像」のバイナリデータファイルなのでここに掲載しません。

117 / 124

添付ファイル lenna.pdf

118 / 124

Date:
From:
To:
Cc:

Sat, 05 Dec 2015 19:02:07 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", 上杉 健太朗
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, "Takashi Sakurama", Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: SIXM の位相回復プログラムの別バージョン

たけうちさま、
うえすぎさま、
GSJ/AIST のなかのです。11/23 の E-mail で紹介した Martino 論文の方法と、それに準じた spline 関数を用いた手法で SIXM の微分位相の画像（投影画像）の位相回復を行う dp2xp と dp2xp_spline を少しだけ書き換えたプログラムを書いてみました。この E-mail に添付したファイル spline.pdf を御覧下さい。その１～２
ページ目の内容は以前に差し上げたものと同じです。以前の版の位相回復プログラムでは spline.pdf の１
ページ目中程上に記した「Martino 論文の式 [10] の積分を x に関して部分積分した式」を使っていましたが、今回の版では微分位相 g の測定データを数値微分して h = ∂g/∂x の値を推定し、それを使って Martino
論文の式 [10] （を修正した式）を数値積分しました。この手法の詳細は spline.pdf の３ページ目に追記した通りです。
新しいプログラムの使用法は以前のものとまったく同じです。新旧のプログラムを試せるよう両者のソースファイルや Windows 用実行ファイルの名前を以下のようにしました。
部分積分（partial integration；PI）の式を使う以前のプログラム
ソースファイルの名前
dp2xp.c から dp2xp_pi.c に改名
それをコンパイルして得た Windows 用実行ファイルの名前
dp2xp_pi.exe
Martino 論文の手法で位相回復
dp2xp_pi_spline.exe
spline 関数によって位相回復
数値微分（numerical differentiation；ND）を行う新しいプログラム
ソースファイルの名前
dp2xp_nd.c
Windows 用実行ファイルの名前
dp2xp_nd.exe
Martino 論文
dp2xp_nd_spline.exe spline 関数
これらのファイルすべてを SIXM 用の書庫ファイルに入れておきました：
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.taz
http://www-bl20.spring8.or.jp/~sp8ct/tmp/his.zip
ただし、これらの書庫ファイルに入れた位相回復プログラムのコンパイル用スクリプト Makefile.dp とMakefile.dp.mingw は以前のものから変えていません。これらのスクリプトは単一のソースファイル dp2xp.cから２個の実行ファイル dp2xp と dp2xp_spline を生成するので、その実行に際してはソースファイルの名前の付け替えなどを行って下さい：
Linux 機の上での以前の版の位相回復プログラムのコンパイル
ln -s dp2xp_pi.c dp2xp.c
make -f Makefile.dp dp2xp dp2xp_spline
rm dp2xp.c

119 / 124

新しいプログラムのコンパイル
ln -s dp2xp_nd.c dp2xp.c
make -f Makefile.dp dp2xp
rm dp2xp.c

dp2xp_spline

また、事前に実行ファイル名の付け替えなどしておけば、11/23 の E-mail に添付したバッチファイル
lenna.bat を実行して新旧のプログラムの動作テストを行うこともできます：
64 bit Windows 機の上で以前のプログラムの動作をテスト
copy dp2xp_pi.exe dp2xp.exe
copy dp2xp_pi_spline.exe dp2xp_spline.exe
lenna.bat
新しいプログラムの動作テスト
copy dp2xp_nd.exe dp2xp.exe
copy dp2xp_nd_spline.exe dp2xp_spline.exe
lenna.bat
さて、Lenna 画像を用いた新しい位相回復プログラムの動作テストを行いました。この E-mail に添付した
lenna.pdf を御覧下さい。以前の結果の画像も並べて表示してあります。位相回復した画像 lenna_spline とlenna_[0,0.5,1] の下に記した出現した値の範囲が少しだけ違っていますが（新しいプログラムで得た値域の幅の方が狭い）、位相回復に失敗している lenna_spline のものを除けば新旧の画像の違いは見た目にはわ
かりません。手法 PI と ND の違いは計算誤差だけだと言えるでしょう。
とりあえず以上です。

120 / 124

添付ファイル spline.pdf（３ページ目のみ）

添付ファイル lenna.pdf
手法 ND で計算した画像は以前に手法 PI で得た、2015/11/23 の E-mail に添付した lenna.pdf のものと見分けがつかないので、ここには掲載しません。

121 / 124

Date
From:
To:
Cc:

Mon, 28 Dec 2015 18:23:33 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa", 上杉 健太朗
MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, "Takashi Sakurama", Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: Martino 論文の問題点

たけうちさま、
うえすぎさま、
GSJ/AIST のなかのです。これまで何度も E-mail した Martino 達の論文の手法についてですが、ものすご
いことに気づきました。端的に言うと、これは「タチ悪」の手法で、手法を支配する「トレードオフ・パラメータ」
を測定データから決めることができません（論文の構成を再度見直すと、Martino 達もこのことがわかっているようです）。この E-mail に添付した PDF ファイル spline.pdf を御覧下さい。面倒なので以前の PDF に追記しましたが、p.４～５だけで独立して読めるようになっています。その概要は以下の通りです。
[1] Martino 論文の手法では x 成分の微分位相の測定データ g(x,y) から以下の２つのことを仮定して
位相φ(x,y) を推定している：
仮定 (A)
推定すべき位相画像は式(1)の積分 I が停留値をとる場合の Euler 方程式である式(2)を満たす。ただし、これらの２つの式には未知のトレードオフ・パラメータλが含まれている。
仮定 (B)
φと∂g/∂x はどちらも x と y 方向のそれぞれに画像の幅の２倍の周期の奇関数である（この仮定により長方形の画像の４つの辺上でφ = 0 かつ∂g/∂x = 0 となる）。
[2] λの物理的な意味が不明なことを除けば、これらの仮定に特に問題はないように思える。問題は「λの値をどのようにして決めるか」で、これは最小自乗法の場合と同様に I(λ)が停留値（最小値）
をとるようなλの値を使えば良いように思える。
[3] 最小自乗法の場合と同様に考えると∂I/∂λ= 0 となるλの値で I(λ)は停留値をとるはずなので、∂I/∂λの計算式を導出してみた。
[4] しかしながら、最終的な式(13)に示されているように、∂I/∂λ= 0 となるλの有限の値は存在しない。つまり、測定データから I(λ)の値が最小になるλの値を決定することはできない。
[5] 以上から spline.pdf の最後に記した結論を得た：
Martino 達の論文（Opt. Lett., 38, 2013）に記されている手法は「閉じた手法」ではない。処理結果を支配するトレードオフ・パラメータの値を前提条件だけから決めることができない手法である。
と言う訳で、Martino 論文の手法を使うためには式(1)もしくは(2)の物理的な解釈を考えて、それに応じたト
レードオフ・パラメータの値を指定する必要があります。これは不可能に思えるので、今後ぼくは Martino 論
文の手法を使わないことにします。
なお、もしかしたら、この結論は式(1)から(2)を導出する際に用いた「変分法」では当然のこと（つまり、式(1)の積分 I を最小自乗法などと同様に取り扱えると考えることが間違い）だったのかもしれません。とり急ぎ、

122 / 124

添付ファイル spline.pdf（４～５ページのみ）

123 / 124

124 / 124

Date:
From:
To:
Cc:

Mon, 04 Jan 2016 13:37:00 +0900
Tsukasa NAKANO"TAKEUCHI, Akihisa"
上杉 健太朗, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama, Yohei IGAMI,
Akiko TAKAYAMA, "TSUCHIYAMA, Akira"Subject: Re: Martino 論文の問題点

たけうちさま、
GSJ/AIST のなかのです。E-mail ありがとうございます。
> 中野さんの勧める spline 法でそれが解決されているのであれば、誤解を招く名前の文書ファイル spline.pdf に書き足し・書き足ししたものを送り続けてすみません。昨年
11/23 の E-mail で紹介したように、spline 関数による位相回復は全然ダメです。それについては忘れてくだ
さい。とり急ぎ、On Mon, 4 Jan 2016 09:54:53 +0900 "TAKEUCHI, Akihisa" wrote:> 中野様、皆様>> 竹内です。あけましておめでとうございます。
>> 中野さんの考察拝見させていただきました。
> 返答が遅くなり申し訳ありません。
> Martino の論文すべてを理解しているわけでなく、> また、中野さんの考察もまだまだ僕の理解では追いつけておらず、> もっと勉強しないと分からないところだらけで申し訳ありません。
>> Martino の方法は、僕が Hideki で取り入れていた方法でも> トレードオフパラメータの値は主観的な判断で適当に決定しているだけなので> 中野さんの懸念される危うさは僕も認識しているところです。
> 中野さんの勧める spline 法でそれが解決されているのであれば、> もちろんそちらを利用すべきでしょう。
> 中野さんに Hideki での検証を昨年夏ごろに提案されてから未だに> ほとんど進んでいませんが、1-3 月中のビーム停止期間中にどうにか> やってみようと思っています。
> 取り急ぎ。