はじめに
以前浜松ホトニクス社のセミナー動画でLCOS-SLMの解説動画を見ることがありました。
LCOS-SLMは空間位相変調器のことで、液晶を使って位相を制御することで任意のビームを形成する素子です。色々設定した後にONすると複雑なビーム形状を整形する様は魔法のようでした。
ここで思いました。“位相制御するのは分かるけど、どうやって必要な位相求めているんだ？”
動画では勿論アルゴリズム詳細は説明がありません。調べてみると、Gerchberg–Saxton法(以下、GS法)で位相を求められるみたいなので、Pythonで実装してみます。
(LCOS-SLMがGS法使っているかは分かりません。)
私が見たセミナー動画は一般向け公開されているものではなかったのですが、浜松ホトニクス社から原理の解説動画がyoutubeで複数公開されています。
製品ラインナップも色々ありますね。
高パワー向けもありますが、レーザ加工を意識していて確かに相性が良さそうです。
https://www.hamamatsu.com/jp/ja/product/optical-components/lcos-slm.html
GS法の詳細について
GS法は、入力面（SLM）と出力面（目標）の間を伝播計算で往復し、各面で振幅を強制的に書き換えることで、最適な位相パターンへ収束させるアルゴリズムです。
入力面と出力面の光分布を考えます。
出力面で所望のビームプロファイルを得たいと思っています。(例えば出力面で特定の文字列を形成したい等)
出力面での振幅は既知(ターゲット、得たい分布)、位相は不明です。
入力面は振幅と位相は既知です。(振幅はガウシアンビーム、平行光なら位相は均一とする等)
上記の状況で以下プロセスを繰り返し(約100回)行います。
1. 入力面から出力面に伝播計算を行う
今回は単純にフーリエ変換を行います。(フラウンフォーファー回折)
2. 出力面での更新
得られた振幅はターゲットの振幅と全く異なった分布です。そのため、振幅をターゲットのものに変更します。ただし、位相はそのまま維持します。
3. 入力面へ逆伝播
逆フーリエ変換により、出力面から入力面の分布を求めます。
4. 入力面での更新
2.と同様に、位相は維持したまま、入力の振幅に戻します。
PythonによるGS法の実装
GS法をPythonで実装してみます。
入力面はガウシアンビーム、出力面は”YOHAKU”を形成することとします。
入力面の位相はランダムとします。この方が収束性が良いためです。
収束の様子を確認するのにMSE(Mean Square Error)を記録しておきます。MSEは各ステップでのターゲットと計算結果の差分の2乗の平均したものでターゲットからのずれを表します。機械学習などでよく使われますね。
ターゲットとなる”YOHAKU”の画像を用意しました。まずはこれを読み込みます。
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
# 白黒画像を読み込む
img = Image.open("YOHAKU.png").convert("L")
# サイズを512×512に
img = img.resize((512, 512))
# numpy配列へ
target = np.array(img)
# 白文字を1、背景を0
target = target / 255.0
# 振幅
target = np.sqrt(target)ターゲットが決まったので、以下でGS法を実施していきます。
反復回数は200回とします。
# =========================
# パラメータ
# =========================
N = 512
n_iter = 200
# =========================
# 入力面の振幅（ガウシアンビーム）
# =========================
x = np.linspace(-1, 1, N)
X, Y = np.meshgrid(x, x)
r2 = X**2 + Y**2
beam_radius = 0.3
A_in = np.exp(-r2 / beam_radius**2)
# =========================
# 出力面の目標強度
# =========================
A_target = target
# =========================
# 初期位相
# =========================
phase = np.random.rand(N, N) * 2 * np.pi
# =========================
# GS 反復
# =========================
errors = []
for i in range(n_iter):
# 入力面の複素振幅
U_in = A_in * np.exp(1j * phase)
# 出力面へ伝搬（FFT）
U_out = np.fft.fftshift(np.fft.fft2(U_in))
# 出力面の位相を保存
phase_out = np.angle(U_out)
# 出力面で振幅を目標値に置換
U_out_new = A_target * np.exp(1j * phase_out)
# 入力面へ戻す（逆FFT）
U_in_back = np.fft.ifft2(np.fft.ifftshift(U_out_new))
# 入力面の位相を更新
phase = np.angle(U_in_back)
# 誤差を記録（表示用）
recon = np.abs(U_out)**2
recon /= recon.max() + 1e-12
err = np.mean((recon - target)**2)
errors.append(err)
# =========================
# 最終結果
# =========================
U_final = A_in * np.exp(1j * phase)
U_out_final = np.fft.fftshift(np.fft.fft2(U_final))
I_final = np.abs(U_out_final)**2
I_final /= I_final.max() + 1e-12
# =========================
# 表示
# =========================
fig, ax = plt.subplots(1, 4, figsize=(14, 4))
ax[0].imshow(target, cmap='gray')
ax[0].set_title('Target intensity')
ax[0].axis('off')
ax[1].imshow(I_final, cmap='gray')
ax[1].set_title('Reconstructed intensity')
ax[1].axis('off')
ax[2].imshow(phase, cmap='hsv')
ax[2].set_title('Phase hologram')
ax[2].axis('off')
ax[3].plot(errors)
ax[3].set_title('Error')
ax[3].set_xlabel('Iteration')
ax[3].set_ylabel('MSE')
plt.tight_layout()
plt.show()計算結果は下図です。
左からターゲット、計算した強度分布、計算した位相分布、MSE です。
振幅・輝度
計算した強度分布は綺麗にパターンが形成出来ていますが、明るい領域では色が均一になっていないことが分かります。試しに中央で輝度を切り出してみます。
輝度は明るい領域では一定であることが理想ですが、かなりばらつきが大きいことが分かります。
これはGS法のよく知られた特性であるようで、改善したWGS法(Weighted Gerchberg–Saxton)というものがあるようです。
位相
位相マップは一目見てもランダムのようにしか見えません。輝度と同じく中央ラインで切り出してみます。
やはりランダムにしか見えず、これを人間が解釈するのは難しそうです。
最初に位相はランダムとした方が収束が早いのが何となく理解できますね。
MSE
MSEは反復回数100回程度で、最小値を取っているように見えます。コードでは200回としましたが、反復回数100回と200回で差は分かりませんでした。100回程度で収束したと思って良さそうです。
最初はより単純な十字をターゲットにしていたのですが、それだと5回程度で収束しました。当たり前ですが、ターゲットの複雑さで必要なループ回数が変わるようです。
GS法だとMSEエラーは進捗状況の可視化指標にすぎませんが、ニューラルネットワークの学習ではMSE(若しくは損失関数)を直接下げるように計算を行います。この手法も最適な位相を求める強力な手法です。
次回記事
作成した画像はボツボツしていたので、どこまで綺麗にできるか挑戦しました。
今回は伝播計算としてFFTを行いましたが、これは無限遠への伝播を意味しています。
実際には、レンズで集光(フラウンフォーファー回折)するか、有限距離の伝播を考える必要があります。今回の計算にてFFTを角スペクトル法(ASM)に置き換えてみるのも面白そうなので、その内やれればいいなと思います。
角スペクトルについてはこの記事を参照ください。
レンズを使ったフラウンフォーファー回折を考える際はこちらの記事シリーズを参照ください。
まとめ
本記事ではGS法の概要を説明し、Pythonを用いて”YOHAKU”の文字を生成するための位相分布を計算しました。
- GS法では、生成されるパターンの輝度に大きなばらつきが見られました。
- 求められた位相分布はランダムに近く、人間が直感的に理解することは困難です。
- MSEエラーは100回程度で最小値をとるが、200回時と得られる画像との差は分かりませんでした。
本記事で扱ったGS法は、位相回復（Phase Retrieval)と呼ばれる問題を解く代表的なアルゴリズムの一つです。強度しか観測できない状況から位相を推定するという点は非常に興味深く、光学だけでなくX線回折や電子顕微鏡など、さまざまな分野で重要なテーマとなっています。
また、「欲しい強度分布から位相分布を求める」という意味では、GS法は逆問題を解く手法とも考えられます。逆問題は一般に解くことが難しい一方で、工学や科学の多くの課題は最終的に逆問題へ帰着します。
GS法自体は1970年代に提案された古典的なアルゴリズムですが、その後もWeighted Gerchberg–Saxton法や勾配法、深層学習を用いた位相設計など、さまざまな発展手法が提案されています。シンプルでありながら現在の計算光学にもつながる、非常に面白いテーマだと感じました。
参考文献
フーリエ光学を日本語で学びたい方におすすめの一冊です。 回折やフーリエ変換、ホログラムも分かりやすく解説しており、GS法の背景となる光学理論を体系的に理解できます。
GS法を理解するうえで土台となるフーリエ光学を学べる定番の名著です。フーリエ変換、回折、レンズのフーリエ変換、空間周波数など、GS法の背景となる理論が体系的に解説されています。「なぜFFTで伝搬を計算できるのか」を深く理解したい方におすすめの一冊です。
コメント