PythonによるX線回折顕微法（位相回復）

X線回折顕微法（位相回復）とは

全てを理解するのは大変ですので、ここでは簡単に説明します。
詳しく知りたい方は、下のサイトをご覧ください。
本記事は、こちらを参考にしております。

---

X線を試料に照射すると、回折したX線がスクリーンに現れます。
回折したX線は、複素数のデータを持っています。
しかし、スクリーンに現れている像（観測可能なもの）は、実数部分の強度のみで位相の部分は観測できません。
試料の電子密度分布（オリジナルのデータ）をフーリエ変換して絶対値の2乗を取ったものです。
ゆえに、位相のデータは消えてしまいます。
これでは、元に戻そうと逆フーリエ変換をしても位相のデータがないために元には戻りません。
これを元に戻す手法が、今回紹介するX線回折顕微法（位相回復）というものです。

回復の流れは下の概略図です。

逆空間の制約は、位相データが失われることを示しています。

実空間の制約は、最も肝である位相を回復させるためのアルゴリズムを適用することを示しています。
そのアルゴリズムは、Fienup により提案された HIO (hybrid input-output) アルゴリズムというものです。
下の式がそのアルゴリズムの内容です。

Sはサポートと呼ばれるもので、実空間においてデータが存在している領域を示します。
βの値は、電子密度を減らす度合いを決めるパラメータです。

上で説明したことを繰り返すと下のように回復していきます。

回復していく様子

元の画像

単に逆フーリエ変換して戻した画像

現時点ではどう回復していくのかわかりずらいかと思いますので、下に進んでプログラムを書きながら理解してもらえたらと思います。

プログラム

ライブラリの導入

3つのライブラリを使用します。

1. NumPy
2. ImageIO
3. math

それぞれについて簡単に説明します。

NumPyは、多次元配列を効率的に扱うライブラリです。 Pythonの標準ライブラリではないが、科学技術計算や機械学習など、ベクトルや行列の演算が多用される分野では、ほぼ必須のライブラリです。

Imageioは、アニメーション画像、ボリュームデータ、科学的フォーマットなど、さまざまな画像データの読み書きを簡単に行うためのPythonライブラリです。

mathは、数学計算用の変数や関数、定数が入っているライブラリです。

それでは、以下のコードを追加してください。

import numpy as np
import numpy.fft as fft
import imageio
import math

これでライブラリの導入は完了です。

画像の取り込み

位相回復に用いるオリジナルの白黒画像を用意します。
白黒である理由は、X線を当ててスクリーンに現れるものは電子密度の高低だけの情報であるからです。
理論上はカラー画像を用いてもできますが、実験で用いないというのと、少し複雑になるといいう理由でここでは取り扱いません。

白黒の場合は、1pxの情報が1つの値で表現される。
それに対しカラーの場合は、RGBの3色、すなわち3つの値で表現される。
つまり、白黒に比べてカラーは操作が３倍になってしまう。

サイズが大きすぎると、PCが重くなりますので、小さいものを推奨します。

ここでは、この画像を用います。150px✖150pxのPNG画像です。(.pngと.bmpの画像は動作確認しています。)

準備ができたら、こちらのコードを追加します。

data = imageio.imread("画像のパス")

ファイルのパスの部分には、画像のパスを入れてください。

print(data)で中身を見てみます。

---

dataの中身（一部）

[174, 174, 174, 174, 174, 174, 174, 174, 174, 174, 174, 174, 174,174, 174, 174, 174, 174, 174, 168, 144, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136, 136,136, 136, 136, 136, 136, 136, 136]

要素数が大きすぎて省略表現になる場合は、こちらをファイル上部に加えてください。すべて表示されるようになります。

np.set_printoptions(threshold=np.inf)

---

上のような大量の数字が格納されています。
これは、画像の１ピクセルごとの白黒の濃淡の度合いを示す値を格納しています。
白黒の濃淡を0~255までの256段階で表現しています。

print(data.shape)で、dataの形を見てみます。

(150, 150)と表示されます。
これは、行列の形が150×150であることを意味しており、150px✖150pxの画像であることと一致しています。

これで画像の取り込みについては完了です。

オーバーサンプリングを模擬したパッド画像の作成

オーバーサンプリングを模擬したパッド画像の作成を行います。
サンプル画像をフーリエ変換し、位相回復を始める元のデータを作成します。
下の画像の→の部分です。

次に、そのデータをオーバーサンプリングします。
フーリエ変換された逆空間のデータは、実空間に比べて1/2の大きさであるため、元に戻すためにデータ数を2倍以上に増やす必要があります。

以下のコードを追加してください。

pad_len = len(data)
padded = np.pad(data, ((pad_len, pad_len),(pad_len, pad_len)), 'constant', 
                constant_values=((0,0),(0,0)))

ft = fft.fft2(padded)

何をしているのか見てみます。
まず、以下のコードです。

pad_len = len(data)

print(pad_len)で中身を見てみます。
150という値が入っています。

len()関数は、引数に指定したオブジェクトの長さや要素の数を取得します。

よって、sourceの要素の長さをpad_lenに代入しているということです。

次に、以下のコードです。

padded = np.pad(data, ((pad_len, pad_len),(pad_len, pad_len)), 'constant', 
                constant_values=((0,0),(0,0)))

print(padded)で中身を見てみます。
paddedの中身は、中央付近にdataの行列があり、その回りが0で囲まれています。
print(padded.shape)で行列の形は、450×450であるとわかります。

つまり、縦、横3倍ずつののオーバーサンプリングになっています。

np.pad()関数は、第一引数に置いたオブジェクトに対し、その外側に第二引数で指定した数だけ、第三引数で指定したモードで、第四引数で指定した値を付け加えるというものです。
第二引数の((a,b),(c,d))は、aが上、bが下、cが左、dが右にどれだけ加えるかを決めます。第四引数においては、加える値を決めます。

よって、ここではdataの行列の上下左右にpad_len(=150)ずつ0を加えています。

以下に分かりやすい例を示します。
値を変えて理解を深めてみてください。

x=np.arange(1, 3*3 + 1).reshape(3, 3)
padded_ex = np.pad(x, ((2,1),(1,2)), 'constant', 
                constant_values=((1,0),(2,0)))
print(x)
print(padded_ex)

---

xの中身

[1 2 3]
[4 5 6]
[7 8 9]

padded_exの中身

[2 1 1 1 0 0]
[2 1 1 1 0 0]
[2 1 2 3 0 0]
[2 4 5 6 0 0]
[2 7 8 9 0 0]
[2 0 0 0 0 0]

---

最後に、以下のコードです。

ft = fft.fft2(padded)

print(ft)で中身を見ます。

---

[ 3.84266800e+06 +0.j -3.17392520e+06-186627.01427012j
1.57691612e+06+252617.37925885j ... 1.64735118e+04+169556.99706177j
1.57691612e+06-252617.37925885j -3.17392520e+06+186627.01427012j]
[-3.17531855e+06-248438.65644106j 2.61590728e+06+362166.12648323j
-1.29399605e+06-318845.10560134j ... -2.43796108e+03-127288.30912719j
-1.31543615e+06 +97850.85466241j 2.63063099e+06 +53638.44218094j]
[ 1.58118214e+06+341171.48173816j -1.30082797e+06-359974.1270297j
6.45505335e+05+247508.34246826j ... -7.83254508e+03 +56634.46219532j
6.64698324e+05 +44901.31504278j -1.31528247e+06-207599.70418091j]

---

paddedとは、明らかに違う値が入っています。
print(ft.shape)で形を見ると、450×450でpaddedと違いはないです。

fft.fft2(x)関数は、xの値をフーリエ変換します。
よって、ここではpadded内のすべての値をフーリエ変換しています。

以上で、オーバーサンプリングを模擬したパッド画像の作成が完了です。

回折パターンのシミュレート

フーリエ変換した画像データの絶対値を取り、下の画像の|F(K)|を作成します。

以下のコードを追加してください。

g = np.abs(ft)

np.abs(x)関数は、xの絶対値を返します。

次に、以下のコードを追加してください。

padded_len = len(padded)

後で使うので、paddedの長さを取得しておきます。

以上で、回折パターンのシミュレートは完了です。

画像の位置とパディングの位置を把握する（サポートとその外の判定）

サンプル画像のある所（サポート）を1、ない所を0として、サポート領域か、その外かを判定するためのマスクを作成します。

以下のコードを追加してください。

mask = np.ones((pad_len+2,pad_len+2))

152×152の行列で、要素がすべて1になっています。

np.ones()関数は、要素がすべて1の行列を作成する関数です。

次に、以下のコードを追加してください。

mask = np.pad(mask, ((pad_len-1, pad_len-1),(pad_len-1, pad_len-1)), 'constant', 
                constant_values=((0,0),(0,0)))

要素がすべて１の152×152の行列の上下左右を0でパディングして、450×450の行列を作成しています。

以上で、画像の位置とパディングの位置の把握（サポートとその外の判定）は完了です。

ランダムな位相情報を使用した初期推測

位相回復を始めるときの一番最初の位相はランダムに与える必要があります。

以下のコードでランダムな位相を作成します。

phase = g * np.exp(1j * np.random.rand(l,l) * 2 * math.pi)

|G(K)|にランダム位相e^(iφK)をかけて、G(K)を作成しています。

以上で、ランダムな位相情報を使用した初期推測は完了です。

反復回数の設定

下の画像のサイクルを繰り返すことでより元の画像の復元ができます。

ここでは、は800回繰り返します。
以下のコードを追加してください。

n = 801

以上で、反復回数の設定は完了です。

回復する時の係数の設定

下の画像のβの値を設定します。

理論的には、1.0が最適だとわかっています。
ここでは、特に意味はないですが0.9に設定してみます。

以下のコードを追加してください。

b = 0.9

以上で、回復する時の係数の設定は完了です。

回復

準備は整ったので、下の画像のサイクルを作成します。

以下のコードを追加してください。

#前回の結果
pre = None

for s in range(0,n):
    #フーリエドメイン制約を適用する
     reci = g * np.exp(1j * np.angle(phase)) 
    
    inv = fft.ifft2(reci)
    inv = np.real(inv)
    if pre is None:
        pre = inv
    
    #実空間制約を適用する
    temp = inv
    for i in range(0,l):
        for j in range(0,l):
            #画像領域が正であること
            if inv[i,j] < 0 and mask[i,j] == 1:
                inv[i,j] = pre[i,j] - b*inv[i,j]
            #サポート領域の強度をゼロに近づける
            if mask[i,j] == 0:
                inv[i,j] = pre[i,j] - b*inv[i,j]
    
    pre = temp
    
    phase = fft.fft2(inv)

    #進行状況を画像で保存する
    if s % 10 == 0:
        imageio.imwrite("保存先のパス/progress" 
                    + str(s) + "任意の画像の拡張子", pre)
        print(s)

何をしているのか見ていきます。

まず、以下のコードです。

pre = None

preは前回の回復結果です。
Noneで初期化しておきます。

次に、以下のコードです。

for s in range(0,n):

サイクル全体の繰り返しを行っています。

次に、以下のコードです。

reci = g * np.exp(1j * np.angle(phase))

下の画像のピンクの部分の操作です。

np.angle(z)関数は、複素数zの偏角を返します。

次に、以下のコードです。

inv = fft.ifft2(reci)
inv = np.real(inv)
if pre is None:
        pre = inv

fft.ifft2(x)関数は、xの逆フーリエ変換した値を返します。

np.real(z)関数は、複素数zの実数部の値を返します。

if文の所は、1回目のサイクルの時は1個前の結果が存在しないため、補完しています。

次に、以下のコードです。

temp = inv
for i in range(0,l):
    for j in range(0,l):
        if inv[i,j] < 0 and mask[i,j] == 1:
            inv[i,j] = pre[i,j] - b*inv[i,j]

        if mask[i,j] == 0:
            inv[i,j] = pre[i,j] - b*inv[i,j]

下の画像の緑の部分の操作です。

具体的には、下の画像の計算を行っています。
上側は、temp = invです。
それ以外が下側の計算です。

preが$ρ_n$、invが$ρ'_n$です。

次に、以下のコードです。

pre = temp
    
phase = fft.fft2(inv)

次のサイクルへの移行部分です。
preに関してはそのまま次のサイクルに使用したいため、フーリエ変換する必要はないです。
invはフーリエ変換します。

最後に、以下のコードです。

    if s % 10 == 0:
        imageio.imwrite("保存先のパス/progress" 
                    + str(s) + "任意の画像の拡張子", pre)
        print(s)

回復の途中経過の画像preを書き出しています。

保存先のパスの所に保存したいフォルダのパスを入力してください。

任意の画像の拡張子の所に.bmpや.pngなどを入力してください。

以上で回復の操作は完了です。

ソースコード全体

import numpy as np
import numpy.fft as fft
import imageio
import math

#ここから画像の取り込み

data = imageio.imread("画像のパス")

#ここからオーバーサンプリングを模擬したパッド画像の作成

pad_len = len(data)
padded = np.pad(data, ((pad_len, pad_len),(pad_len, pad_len)), 'constant', 
                constant_values=((0,0),(0,0)))

ft = fft.fft2(padded)

#ここから回折パターンのシミュレート

g = np.abs(ft)
padded_len = len(padded)

#ここから画像の位置とパディングの位置を把握する（サポートとその外の判定）

mask = np.ones((pad_len+2,pad_len+2))
mask = np.pad(mask, ((pad_len-1, pad_len-1),(pad_len-1, pad_len-1)), 'constant', 
                constant_values=((0,0),(0,0)))

#ここからランダムな位相情報を使用した初期推測

phase = g * np.exp(1j * np.random.rand(padded_len,padded_len) * 2 * math.pi)

#ここから反復回数の設定

n = 801

#ここから回復する時の係数の設定

b = 0.9

#ここから回復

#前回の結果
pre = None

for s in range(0,n):
    #フーリエドメイン制約を適用する
    reci = g * np.exp(1j * np.angle(phase)) 
    
    inv = fft.ifft2(reci)
    inv = np.real(inv)
    if pre is None:
        pre = inv
    
    #実空間制約を適用する
    temp = inv
    for i in range(0,padded_len):
        for j in range(0,padded_len):
            #画像領域が正であること
            if inv[i,j] < 0 and mask[i,j] == 1:
                inv[i,j] = pre[i,j] - b*inv[i,j]
            #サポート領域の強度をゼロに近づける
            if mask[i,j] == 0:
                inv[i,j] = pre[i,j] - b*inv[i,j]
    
    pre = temp
    
    phase = fft.fft2(inv)

    #進行状況を画像で保存する
    if s % 10 == 0:
        imageio.imwrite("保存先のパス/progress"
                    + str(s) + "任意の拡張子", pre)
        print(s)

Gifの作成

せっかくなので、回復されていく様をGifにします。

コードファイルと同じフォルダ内に画像を設置して、以下のコードを実行します。

from PIL import Image
#画像を入れるリスト
pictures=[]
#画像を入れる
for i in range(79):
    pic_name='progress'+str((i+1)*10)+ '画像の拡張子'
    img = Image.open(pic_name)
    pictures.append(img)
#gifアニメを出力する
pictures[0].save('anime.gif',save_all=True, append_images=pictures[1:],
optimize=True, duration=500, loop=0)

'画像の拡張子'には、ご自分の作成した画像の拡張子を入力してください。

おわりに

無事に画像は回復されましたか？
反復数を増やすなど色々つついて遊んでみてください。

コード自体は短いのですが、原理を理解するのは骨が折れます。

生物や医療の分野でホットな技術のようです。

参考

1. http://www.jssrr.jp/journal/pdf/19/p003.pdf

2. J. R. Fienup: Appl. Opt. 21, 2758 (1982).