SIXM用の画 - X-ray Micro-Imaging at SPring-8

1 / 123
SIXM_mails
date
subject
4/28
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 / 123
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論文の問題点
3 / 123
Date:
From:
To:
Subject:
Tue, 28 Apr 2015 18:11:31 +0900
Tsukasa NAKANO
Kentaro Uesugi
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 / 123
(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 / 123
[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 / 123
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 / 123
使用例
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 / 123
省略します。
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 / 123
→
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 a
xp2tg a.bin 150e-7 200.5
0
屈折トモグラフィーの画像再構成
mkdir r
xp2tg_dpc r.bin 150e-7
200.5
-8
a
0
-8
> a.log
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 / 123
> センター値 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 / 123
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:
From:
To:
Cc:
Subject:
Fri, 19 Jun 2015 14:07:26 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
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 / 123
(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 / 123
# 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 / 123
Date:
From:
To:
Subject:
Tue, 23 Jun 2015 19:41:59 +0900
Tsukasa NAKANO
Kentaro UESUGI
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 / 123
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)
(2)
(3)
(4)
(5)
(6)
プログラム xp2tg（および xp2tg_dpc）のコードの修正
バイナリデータファイル XP.bin の（画像再構成以外への）利用法
his2xp と his2xp_top の処理内容の違いについて
xp2tg_dpc の起動パラメータ Dr に関する話
xp2tg のベンチマークテスト
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 / 123
１枚のスライスの画像再構成用の 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 / 123
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 / 123
さて、これらのすべて画像に「縦線」が入っています。つまり、スキャンごとの入射 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 / 123
# 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 に書いた通りです。
a2tiff
XP.bin
{base
step}
BPS XP.tif
20 / 123
起動パラメータ 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 / 123
これらに記したように、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 / δr
h(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) の畳み込み積分を行うことにより "filtered
projection" q(r) を計算します。
q(r) = ∫ p(ρ) ・[g,h] (r - ρ) dρ
22 / 123
ただし、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 として指定した値で filtered
projection を割るようになっています。つまり、吸収トモグラフィーの場合には 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 / 123
添付ファイル sixm.gif
添付ファイル sixm_ra.txt
SIXM で測定できる屈折角、Θr
SIXM の測定時の座標系 (s, r, z)
s
r
z
光源から投影面に向かう投影面に直交する座標軸
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 / 123
τと 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 / 123
p(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:
From:
To:
Cc:
Subject:
Tue, 07 Jul 2015 15:26:59 +0900
Tsukasa NAKANO
Kentaro UESUGI, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, Aki Takigawa, Junya Matsuno, 中村隆太, Takashi Sakurama
SIXM 用画像再構成プログラムに関する話の続きの続き
みなさま、
なかのです。6/25 と 7/2 の E-mail の続きです。
On Thu, 25 Jun 2015 17:19:47 +0900
Tsukasa NAKANO wrote:
26 / 123
> ... プログラム 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 / 123
上記の並びはプログラムの処理速度の降順になっていて（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 / 123
画像の画素数は常に従来の再構成画像の画素数以下の個数になっている）。さらに、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 / 123
添付ファイル 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 / 123
その他
シングルスレッドで処理
逆投影の計算法の違い（３種類）
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
040711j
size of projection images = 1000 * 720 pixels * 360 views
interval of detectors, dr = Dr = 0.000583 cm
31 / 123
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
020607a
size 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 / GPU
program
gsjgix / 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
040711j
37.000848
70.974262
130.093564
117.461060
36.611471
95.058632
145.413615
147.363683
38.605144
123.978912
181.122563
241.213953
020607a
16.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 / GPU
program
gsjgix / K20c
hp2tg_tm
/hp2tg_tint_OM
040711j
0.142750
0.273820
020607a
0.310532
0.216870
32 / 123
//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.453168
0.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
020607a
1e+06
log(pixels)
log(pixels)
hp2xp + xp2tg
hp2tg_tnai_OM
1e+06
100000
10000
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 / 123
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 / 123
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.bin axp_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 / 123
# 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 / 123
添付ファイル 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 / 123
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 と Ry
RMSD 計算の対象とする２画像の x および y 方向の位置ズレの相対的な範囲を指定するための
パラメータ。例えば Rx = 0.5 を指定すると、０度の画像に対して 180 度の画像を画像の横幅の半
分だけ左に移動した状態から右に半分移動した状態までのすべての配置に対する RMSD を計算
する。なお、xp2rc.csh と等価な処理を行う、HASPET 用の画像再構成プログラムのひとつである
C-shell script "center.csh" は Rx = 0.5、Ry = 0.25 を採用している。
RMSD.txt
38 / 123
位置ズレごとの 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 / 123
> 具体的には 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 / 123
添付ファイル 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 / 123
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 / 123
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 / 123
Date:
From:
To:
Cc:
Thu, 16 Jul 2015 11:05:16 +0900
Tsukasa NAKANO
MIYAKE, "中村隆太"
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 <[email protected]>
MIYAKE, "中村隆太"
Akira TSUCHIYAMA, "MATSUNO Junya", Aki Takigawa, Takashi Sakurama,
Kentaro Uesugi, "TAKEUCHI, Akihisa"
Subject: SP8_BL47XU_SIXM_150603c_a.HIS_訂正
44 / 123
みなさま、
なかのです。先程の 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 ‐> 1755732483 27783218835 150603b/a.HIS # [b] % cksum 150603c/a.HIS ‐> 1103610525 27772972435 150603c/a.HIS # [c] % cksum 150603d/a.HIS ‐> 3027004535 27783218832 150603d/a.HIS # [d] # [b],[c],[d] % head ‐1 150603b/a.log ‐> 601,451,100,1,0 % head ‐1 150603c/a.log ‐> 601,451,100,1,0 % head ‐1 150603d/a.log ‐> 601,451,100,1,0 % expr 601 ¥* 451 + 100 ‐> 271151 # [e] % check_his 150603b/a.HIS | cat ‐n | tail ‐1 ‐> 27783218835 # $1 == [b] 271151 128 400 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 128 400 2 # $1 == [e]
Date:
From:
Fri, 17 Jul 2015 16:34:28 +0900
Tsukasa NAKANO
45 / 123
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:
From:
To:
Cc:
Subject
Thu, 16 Jul 2015 11:20:22 +0900
Tsukasa NAKANO
Kentaro Uesugi, "TAKEUCHI, Akihisa"
Akira TSUCHIYAMA, MIYAKE, "MATSUNO Junya", Aki Takigawa, Takashi Sakurama, "中村隆太"
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 / 123
たけうちさま、
なかのです。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 が食い違って
います。
測定
150605d
150605f
150605i
150605j
a.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 / 123
たけうちさま、
なかのです。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 / 123
ディレクトリ /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 / 123
その際に先のファイル 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 / 123
(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 / 123
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 / 123
(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 / 123
(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
[email protected]:/media/disk/tsukasa/150602/ 150602d.org/, 150603j.org/ sample.xlsx log.pdf sixm.csv, sixm.ods 150603c/, 150605a/ a.HIS, a.log, a.par axp.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 / 123
150604b/, 150604c/, 150604d/, 150604e/, 150604f/, 150604g/, 150605b/, 150605c/, 150605d/, 150605e/, 150605f/, 150605h/, 150605i/, 150605j/ a.HIS, a.log, a.par axp.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.xp i0s.tbl his2xp_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 / 123
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 / 123
[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 ビット画素値
の値域は
150602d
1514.00 ～ 1750.07
150603j
1514.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 / 123
屈折プロファイルの「重心」の位置
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 / 123
れ青と赤色で示してあります。また、３個の 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 / 123
添付ファイル ri+rp.pdf
添付ファイル rr+sd.pdf
60 / 123
添付ファイル rr+lp.pdf
添付ファイル rr+gp.pdf
61 / 123
添付ファイル msd+gp.pdf
62 / 123
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 / 123
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
FZP の一次回折光を集光プローブとして試料をスキャンして、
プローブが広がった像を後ろのカメラで取り込んでいるので
この広がったプローブが「台地」のような像として検出されます。
ただ、FZP はレンズとして理想的な素子ではなく、集光する光以外にもいろんな光が出てきます。
この場合一番大きな影響があるのが非回折光、つまり FZP をそのまま透過してくる光です。
これをできるだけ止めるために、試料面直前に OSA と呼んでいるスリットをおいてます。
「山脈」に見えるのは、このスリットを抜けてきた FZP の非回折光です。
非回折光はほとんど広がらないのに対して、回折光は約 60 倍に広がるので、
その分ピーク強度としては弱くなります。
ですが、全体強度としては「台地」の方が圧倒的に大きく（だいたい 25:1）なので
「山脈」=ノイズ は、無視してそのまま一緒に計算しています。
もちろん、本来は、「山脈」は分離して計算すべきです。
ただ、光学的にセンタービームストップを置くなどして「山脈」を分離した実験も試みましたが、
実際のところほとんど変わらない
（むしろ位相像では若干感度が落ちる。重心演算上では重心位置をキャンセルする向きに働くので）
ので、「山脈」は気にしないで計算しているのが実情です。
実際、.his のデータでは「山脈」の強度が飽和しているものも多くありますが、それも気にしません。
それよりも「台地」の信号量を増やす方が結局 SN の向上につながるので。
「山脈」が出ない素子としては、ミラーや屈折レンズなどがⅩ線領域でもありますが、
ミラーは装置として安定性が問題＋スペックルノイズで感度が落ちる
レンズは 8keV 前後でいい素子がない
ということで FZP を使っています。勿論、いいミラーやレンズが出てくれば
そちらに乗り換えるのがベストかと思います。
[a][c]についてはこれで回答になったかと思います。
[b]については、勿論左右対称であるに越したことはありません。
データによって対称性が崩れてるのは、OSA、Diffuser、Δθなどのずれが考えられます。
Δθのズレ自体は、本来左右対称性には関係ないはずですが、
47XU での実験の場合、FE スリット開口を極端に狭めていてしかもその開口の形がいびつなため、
Δθがずれるとその影響が出てきてしまいます。
左右非対称であっても、測定の間に大きく変化なければ大した問題ではないと考えていますが
中野さんも指摘されている短波長の縞状のノイズはΔθが原因かと考えています。
Martino の方法は僕自身もあれが最適な方法とは思っていません。現状可能な妥協案と考えています。
Pfeiffer の方法ベースでノイズ除去も納得いく方法があればそれがベストと思います。
その点は中野さんと完全に同意です。
単に僕がいい方法を見つけられてないだけという。
以上、長くなりました。
64 / 123
竹内さんからの E-mail の添付ファイル
65 / 123
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 / 123
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 / 123
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.h
run.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 / 123
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・Σ_j
w(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 / 123
(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 coefficients
index of refraction
n = 1 - δ - i・β
= 1 - re / ( 2・π )・λ^2・Σ_j N(j)・f(j)
δ
β
re
λ
N(j)
f(j)
refractive index decrement (RID)
absorption index
classical electron radius = 2.817940289458e-15 [m]
wave length
number of element "j" per unit volume
atomic scattering factor of element "j"
---λ [m] = ( h・c / e ) / E [eV] = 1.239841576558181422983e-6 / E [eV]
h
c
e
E [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
ρ
w(j)
A(j)
Avogradro's number = 6.0221417930e+23 [1 / mol]
density of material (including element "j")
weight ratio of element "j"
atomic weight of element "j"
(1.7.1)
70 / 123
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
# Z 01 H 02 He 03 Li 04 Be 05 B 06 C 07 N 08 O 09 F 10 Ne 11 Na 12 Mg 13 Al 14 Si 15 P 16 S 17 Cl 18 Ar 19 K 20 Ca 21 Sc element name http://henke.lbl.gov/optical_constants/getdb2.html “echo 8 | rid_cm 1 Z:1” ↓ Hydrogen Helium Lithium Beryllium Boron Carbon Nitrogen Oxygen Fluorine Neon Sodium Magnesium Aluminum Silicon Phosphorus Sulfur Chlorine 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 / 123
22 Ti 23 V 24 Cr 25 Mn 26 Fe 27 Co 28 Ni 29 Cu 30 Zn 31 Ga 32 Ge 33 As 34 Se 35 Br 36 Kr 37 Rb 38 Sr 39 Y 40 Zr 41 Nb 42 Mo 43 Tc 44 Ru 45 Rh 46 Pd 47 Ag 48 Cd 49 In 50 Sn 51 Sb 52 Te 53 I 54 Xe 55 Cs 56 Ba 57 La 58 Ce 59 Pr 60 Nd 61 Pm 62 Sm 63 Eu 64 Gd 65 Tb 66 Dy 67 Ho 68 Er 69 Tm 70 Yb 71 Lu 72 Hf 73 Ta 74 W 75 Re 76 Os Titanium Vanadium Chromium Manganese Iron Cobalt Nickel Copper Zinc Gallium Germanium Arsenic Selenium Bromine Krypton Rubidium Strontium Yttrium Zirconium Niobium Molybdenum Technetium Ruthenium Rhodium Palladium Silver Cadmium Indium Tin Antimony Tellurium Iodine Xenon Cesium Barium Lanthanum Cerium Praseodymium Neodymium Promethium Samarium Europium Gadolinium Terbium Dysprosium Holmium Erbium Thulium 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 / 123
77 Ir 78 Pt 79 Au 80 Hg 81 Tl 82 Pb 83 Bi 84 Po 85 At 86 Rn 87 Fr 88 Ra 89 Ac 90 Th 91 Pa 92 U Iridium Platinum Gold Mercury Thallium Lead Bismuth Polonium Astatine Radon Francium Radium Actinium 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 formula RID at 8 keV RID / density material water 1 H:2 O:1 3.620619e‐06 3.620619e‐06 quartz 2.7 Si:1 O:2 8.853034e‐06 3.278901e‐06 calcite 2.715 Ca:1 C:1 O:3 8.898499e‐06 3.277532e‐06 Ca:1 Mg:1 C:2 O:6 9.350504e‐06 3.269407e‐06 dolomite 2.86 Mg:2 Si:1 O:4 1.053804e‐05 3.266596e‐06 forsterite 3.226 Fe:2 Si:1 O:4 1.347144e‐05 3.058355e‐06 fayalite 4.4048 Mg:1 Si:1 O:3 1.049759e‐05 3.270278e‐06 enstatite 3.21 magnetite 5.2 Fe:3 O:4 1.553156e‐05 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 quartz 2.7 SiO2 8.89687453E‐06 3.27693374E‐06 calcite 2.715 CaCO3 9.34872423E‐06 3.26878489E‐06 dolomite 2.86 CaMgC2O6 1.05358367E‐05 3.26591316E‐06 forsterite 3.226 Mg2SiO4 fayalite 4.4048 Fe2SiO4 1.34696238E‐05 3.05794242E‐06 1.04954597E‐05 3.26961344E‐06 enstatite 3.21 MgSiO3 magnetite 5.2 Fe3O4 1.55297585E‐05 2.98649184E‐06 73 / 123
Date:
From:
To:
Cc:
Wed, 30 Sep 2015 09:41:54 +0900
Tsukasa NAKANO
Akihisa TAKEUCHI
MIYAKE, 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 / 123
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
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 / 123
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>>> なので、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 / 123
>
>
>
>
>
>
>
>
>
>
>
>>>> は、分かりました。
>>>>
>>>>> そして、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 / 123
>
>
>
>
>
>
>
>
>
>
>
>
>
>
中野さんが指摘された事以外に、位相の場合は
計算エラーで生じるバックグラウンドの緩い傾斜なんかもあり、
それが土'山さんの云う「ふらつき」につながるのかな、と
考えてますが、ちょっとまだわかりません。
---------------------------------------------------------------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 / 123
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 / 123
xp2sp XP.bin {Dr} SP.bin
xp2sp_dpc XP.bin {Dr} SP.bin
XP.bin
his2xp、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 / 123
--- 以下の２個は 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 / 123
添付ファイル sp.pdf
82 / 123
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 / 123
(5) ７測定・13 セットの X 線 CT の投影画像の SP を再度計算してみました。この E-mail に添付したファイル
sp.pdf を御覧下さい。赤枠で囲んだ４個の SP の画像だけが 10/1 の E-mail に添付した sp.pdf に示し
たものと違っています。改造した xp2sp_dpc は微分位相のデータに含まれる誤差をそれぞれの投影値
に割り振っているので、それらから計算した SP の値域の幅が狭くなっています。
dataset
150602d_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 の画像上での横座標値（== 投影画像のシーケンスの
順番、０～）は測定ごとに概ね同じ値なので、これらの黒い縦線はデータ処理の過程で生じたものでは
ありません。
測定番号
150602d
150603j
黒い縦線の横座標値
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 / 123
[1] RID の空間勾配∇δは facet を横切る時に非常に大きな値になる。
[2] その facet が X 線の照射方向と平行になっている投影画像でベクトル∇δの絶対値 |∇δ| を
測定することになる。
[3] この |∇δ| は facet の位置で「デルタ関数」的に空間分布しているので、そのデータを数値的
に線積分しても facet の位置で「階段関数」的に変動する投影値の空間分布を正確に推定
することは不可能。
[4] この「デルタ関数」の「高さ」として真の値よりも低い値を測定することになるので、それを線積
分した「階段関数」の「段差」も真の値より小さくなる。すなわち、この撮影方向の画像の投影
値の推定値は（他の方向に比べて）全体的に低い値になり、SP の値も小さくなる（== SP の
画像上に黒い縦線が発生する）。
とりあえず以上です。
85 / 123
添付ファイル sp.pdf
86 / 123
Date:
From:
To:
Cc:
Thu, 08 Oct 2015 09:19:16 +0900
Tsukasa NAKANO
MIYAKE
"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 / 123
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 / 123
(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
or
export 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 / 123
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 / 123
が原因かもしれません。前記の投影画像 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 / 123
添付ファイル sfp.pdf
92 / 123
Date:
From:
To:
Cc:
Wed, 14 Oct 2015 12:21:23 +0900
Tsukasa NAKANO
Akira TSUCHIYAMA
MIYAKE, 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 / 123
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>> 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 / 123
Date:
From:
To:
Cc:
Wed, 14 Oct 2015 12:44:33 +0900
Tsukasa NAKANO
Akira TSUCHIYAMA
MIYAKE, 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 TSUCHIYAMA
MIYAKE, 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 / 123
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 の使い方は以下の通りです（前にも説明した気がしますが、念のため）。
起動法
a2tiff
XP.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 / 123
xp2tiff
XP.bin
rangeList
{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 / 123
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.bin
XP.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 / 123
>
>
>
>
>
>
>
>
>
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 / 123
添付ファイル dp.pdf
100 / 123
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) dr
r =±∞ で 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 / 123
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 / 123
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 / 123
添付ファイル ri+rr+rp.pdf
104 / 123
添付ファイル fb.pdf
105 / 123
添付ファイル mb.pdf
106 / 123
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" <[email protected]>
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 / 123
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_formula
material
density (g/cc)
chemical formula
water
1
H:2 O:1
quartz
2.7
Si:1 O:2
calcite
2.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 / 123
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 / 123
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 / 123
添付ファイル spline.pdf
111 / 123
添付ファイル briggs.pdf
112 / 123
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.bin
SIXM などで得た 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 / 123
自作は面倒だったので、その処理時間軽減のためにマルチスレッド処理が可能な線形演算用ライブラリ
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 / 123
添付ファイル spline.pdf
115 / 123
116 / 123
添付ファイル 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 / 123
添付ファイル lenna.pdf
118 / 123
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 / 123
新しいプログラムのコンパイル
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 / 123
添付ファイル spline.pdf（３ページ目のみ）
添付ファイル lenna.pdf
手法 ND で計算した画像は以前に手法 PI で得た、2015/11/23 の E-mail に添付した lenna.pdf のもの
と見分けがつかないので、ここには掲載しません。
121 / 123
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 / 123
添付ファイル spline.pdf（４～５ページのみ）
123 / 123