コンピュータホログラフィによる光伝搬動画像記録・再生法の高速シミュレーション
角江 崇
Takashi KAKUE
千葉大学 大学院工学研究院〒263-8522 千葉県千葉市稲毛区弥生町 1-33
Graduate School of Engineering, Chiba University
1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan

1. はじめにな計算機リソースかつ計算時間で処理可能な光コンピュータホログラフィでは，計測対象とな伝搬動画像記録・再生法の計算機シミュレーショる光の情報をホログラムと呼ばれる 2 次元の干渉ンモデルを新たに開発することである．本稿では，縞画像としてイメージセンサで記録する．記録し高速フーリエ変換（FFT）に基づく高速化手法[6]たホログラムに対して計算機上で再生処理計算について紹介する．

を行い，計測対象の情報を有した光の強度分布と位相分布，つまり 3 次元像を数値的に再生できる．
この特長を活かして，透明ガス流や超音波場の可視化へと適用した例[1-3]が報告されている．コン
ピュータホログラフィの応用事例のうち，光が伝搬する様子を記録，再生可能な技術（以下，光伝搬動画像記録・再生法と呼ぶ）[4,5]は，光源に超短光パルスレーザを用いる．時間分解記録手法な
どのような他技術では成し得ない時空間的に連続な 3 次元動画像，つまり 4 次元（空間 3 次元＋時間 1 次元）像を単一露光で取得できる大きな利点がある．その一方で，得られる 4 次元像には時空間歪みが生じることが明らかにされていた[4]．
光が伝搬する様子を精確に再現するにはこの時空間歪みの補償が必要であり，そのためには光伝搬動画像記録・再生法の計算機シミュレーションが有用であると考えられる．しかし光の 4 次元的な振る舞いを詳細に解析するためには，膨大な計算機リソースと長大な時間，双方を必要とする．
このことから，光伝搬動画像記録・再生法における 4 次元像の時空間歪み補償法を確立するには至っていなかった．
本研究課題では上述の問題に焦点を当てた．すなわち，時空間歪み補償法の確立に向け，現実的2. 光伝搬動画像記録・再生法
2.1 記録の原理
図 1 に光伝搬動画像記録・再生法における記録の原理図を示す．光源となる超短光パルスレーザから出射された光パルスはビーム分割素子で二分される．それぞれ，凹レンズと凸レンズを組み合わせたビーム径拡大素子によって広げられた
後，一方は拡散板，他方はホログラム記録材料へと導入される．拡散板に入射した光パルスは散乱され，様々な方向へと広がっていく（図 1(a)）．この光パルスを物体光パルスと呼ぶ．ホログラム記
録材料へと入射される光パルスは参照光パルスと呼ばれる．図 1(a)の瞬間においては，参照光パルスはホログラム記録材料へは到達していない．
次に，図 1(a)で発生した物体光パルスが，ホログ
ラ ム記 録材料 へと 到達す る瞬 間を考 える （図1(b)）．物体光パルスと同時に参照光パルスが到達する場合，両者は互いに干渉するため，干渉縞が形成される．この干渉縞がホログラムとして，ホログラム記録材料へ記録される．同時に到達しない場合には干渉縞は形成されない．一方，拡散板を照明している光パルスは移動し，図 1(a)の瞬間とは別の位置から新たに物体光パルスが発生

BSM

超短光パルスレーザ

NL CL
M

BSM

M

光パルス

NL CL

拡散板

M

超短光パルスレーザ

NL CL

M
NL CL

拡散板

M

超短光パルスレーザ

NL CL

光パルス

NL CL

物体光パルス
参照光パルス

BSM

拡散板

M

物体光パルス

物体光パルス

参照光パルス

参照光パルスホログラム記録材料(a)

ホログラム記録材料

ホログラム記録材料

(b)

(c)

図 1 光伝搬動画像記録・再生法の記録の原理[6]．(a)から(c)の順に時間が経過．図中の略称はそれぞれ，BS: ビーム分割素子，M: ミラー，NL: 凹レンズ，CL: 凸レンズを表す．
する．この物体光パルスは，図 1(c)の瞬間において，ホログラム記録材料へと到達する．このとき物体光パルス

参照光パルスに目を移すと，図 1(b)の瞬間とは別の位置においてホログラム記録材料を照明して光パルス

参照光パルスいる．図 1(b)と同様に，物体光パルスと参照光パ
拡散板

ルスが同時にホログラム記録材料を照明する場

ホログラム
記録材料
(a)

合，干渉縞が形成され，記録される．つまり光伝光パルスの断面

搬動画像記録・再生法においては，物体光パルス

ホログラム
記録材料の有する時間情報が，ホログラム記録材料において空間情報に変換されて記録される．この性質により，複数回の露光を必要とせず，単一露光によ物体光パルスって 4 次元像を取得できる．
2.2 シミュレーションモデル

参照光パルスの断面
拡散板(b)

図 2 に，光伝搬動画像記録・再生法のシミュレーションモデルを示す．計算処理およびモデルの

図 2

光伝搬動画像記録・再生法のシミュレーシ

簡単化のため，以下のような仮定を導入した．拡ョンモデル[6]．(a) 上面図，(b) 俯瞰図．

散板とホログラム記録材料とは平行であり，その間の距離を𝑑とする．拡散板を照明する光パルスの点P𝐻 (𝑥𝐻 , 𝑦𝐻 , 𝑧𝐻 )に記録される条件は以下のようは，𝑧軸に対して𝑥方向に𝜃𝑂 の角度で入射させる．に与えられる[6]．

光パルスの強度分布は矩形波状であり，光パルスと拡散板との断面像は直線状になる．この直線状の断面から物体光パルスが発生する．散乱の指向性は考慮せず，散乱によって物体光パルスの位相がランダムに乱されたうえで，物体光パルスは球面波状に一様に広がっていくと考える．物体光パ
ルスがホログラム記録材料へと到達する時刻は|𝑥𝐻 sin 𝜃𝑅 − 𝑥𝐷 sin 𝜃𝑂 + 𝑑 − 𝑟𝐷𝐻 | < 𝑐𝑊

(1)

𝑟𝐷𝐻 = √(𝑥𝐻 − 𝑥𝐷 )2 + (𝑦𝐻 − 𝑦𝐷 )2 + 𝑑 2(2)

ここで，𝑐は光速，𝑊は光パルスの時間幅を表す．
(1)式を満たす光のみをホログラム記録材料上で重ね合わせることにより，光伝搬動画像記録・再生法で記録されるホログラムを数値計算できる．

その経路，すなわち光路長に依存して異なる．ホ3. 提案手法

ログラム記録材料へ角度𝜃𝑅 で入射する参照光パ

3.1 PSF を利用した光の回折の計算

ルスを考えた場合，拡散板上の点P𝐷 (𝑥𝐷 , 𝑦𝐷 , 0)から光の伝搬および重ね合わせ，すなわち光の回折発生した物体光パルスが，ホログラム記録材料上は，時間項を考慮しない場合においては以下のように記述できる．
𝑢2 (𝑥2 , 𝑦2 , 𝑑) =×

いては(8)式ではなく(3)式を直接的に計算してい1∬ 𝑢1 (𝑥1 , 𝑦1 , 0)𝑖𝜆たため，計算に莫大な時間を要していた．そこで(3)𝑑2𝜋2 exp (𝑖 𝜆 𝑟21 ) 𝑑𝑥1 𝑑𝑦1𝑟21𝑟21 = √(𝑥2 − 𝑥1 )2 + (𝑦2 − 𝑦1 )2 + 𝑑 2

新たに，光路長差を考慮して PSF の分布を変調する手法を提案する．まず，以下のような関数ℎを新たに定義する．

(4)

𝑢1 と𝑢2 はそれぞれ伝搬元と伝搬先における光のℎ(𝑥1 , 𝑦1 ; 𝑥2 , 𝑦2 ) =𝑑2𝜋2 exp (𝑖 𝜆 𝑟21 )𝑖𝜆𝑟21(9)複素振幅分布を表す．𝑖は虚数単位，𝜆は光の波長

ここで，𝑟21 は(4)式で与える．(9)式は，伝搬元のである．伝搬元と伝搬先の面は平行であり，その面上にある任意の点(𝑥1 , 𝑦1 , 0)からの光が，伝搬先距離は𝑑である．(3)式は次のように，畳み込みをの面上にある任意の点(𝑥2 , 𝑦2 , 𝑑)にどのような分利用して書き換えることができる．

布を形成するのかを表している．(9)式と(6)式を𝑢2 (𝑥, 𝑦, 𝑑) = 𝑢1 (𝑥, 𝑦, 0) ∗ ℎ𝑅𝑆 (𝑥, 𝑦)(5)

用いると以下の関係が得られる．
ℎ(𝑥0 , 𝑦0 ; 𝑥0 + 𝑥, 𝑦0 + 𝑦) = ℎ𝑅𝑆 (𝑥, 𝑦)𝑑2𝜋ℎ𝑅𝑆 (𝑥, 𝑦) =exp (𝑖𝑟)
2𝑖𝜆𝑟
𝜆(6)

𝑟 = √𝑥 2 + 𝑦 2 + 𝑑 2(7)∗は畳み込み演算子を表す．ℎ𝑅𝑆 (𝑥, 𝑦)は点拡がり関数（PSF）と呼ばれ，光がどのように広がっていくのかを表すものである．(5)式から，PSF の分布が不変であれば，点拡がり関数との畳み込みによ(10)(10)式は，点(𝑥0 , 𝑦0 , 0)からの光が点(𝑥0 + 𝑦, 𝑦0 +𝑦, 𝑑)に形成する分布が PSF に相当することを示している．つまり，光の発生点が(𝑥0 , 𝑦0 )だけ変化しても光の分布は(𝑥0 , 𝑦0 )だけシフトするのみで，光の分布自体は変化しないことを意味する．
次に，ℎに対して以下の制約を与える．って光の回折を記述できることが分かる．畳み込|𝑥 sin 𝜃𝑅 − 𝑥0 sin 𝜃𝑂 + 𝑑 − 𝑟0 | < 𝑐𝑊

(11)

みはフーリエ変換により計算できるので，(5)式は𝑟0 = √(𝑥 − 𝑥0 )2 + (𝑦 − 𝑦0 )2 + 𝑑 2(12)

さらに以下のように表せる．
𝑢2 (𝑥, 𝑦, 𝑑)= ℱ −1 [ℱ[𝑢1 (𝑥, 𝑦, 0)]ℱ[ℎ𝑅𝑆 (𝑥, 𝑦)]](11)式は，光伝搬動画像記録・再生法において物(8)体光パルスがホログラム記録材料へと記録される た め の 条件 で ある (1)式 と 同 様 であ る ．𝑑 =ℱ[ ]およびℱ −1 [ ]はそれぞれフーリエ変換，逆フー270 mm，𝜆 = 522 nm，𝑝 = 2.2 μmを与え，計算対リエ変換の演算子を表す．計算機においては，フ象とする領域を 4,096×4,096 画素で表した場合にーリエ変換は FFT を利用して高速に計算可能でおけるℎの分布を図 3 に示す．各画像中に示したある．したがって，(8)式に基づいて光の回折を計

座標は(𝑥0 , 𝑦0 )を表す．前述したように，ℎの分布算すれば，(3)式を直接的に計算するよりも高速には(𝑥0 , 𝑦0 )に応じてシフトしているだけであるこ

結果を取得できる．とが分かる．一方，図 3(b)は𝜃𝑅 = 2°，𝜃𝑂 = 3°と3.2 光路長差を考慮した PSF による光伝搬動画し，(11)式を用いて光路長を考慮した場合におけ像記録・再生法のシミュレーションるℎの分布を表す．(11)式を満たさない領域につ光伝搬動画像記録・再生法においては，参照光いては黒色で表現した．図 3(b)では，(𝑥0 , 𝑦0 )が変と物体光が同時にホログラム記録材料へと到達化するとその位置がシフトするだけでなく，分布した場合にのみ干渉が生じるため，光の伝搬にお中のパターンも変化している様子が分かる．このいて時間項，すなわち光路長差を考慮する必要がことは，図 3(b)で与えられる分布を(5)式におけある．したがって 3.1 節に示した計算をそのままる PSF として利用できないことを表している．そ適用できない．この点を考慮し，従来手法[5]におこで今回は，(11)式に対してさらに𝜃𝑅 = 𝜃𝑂 という

(a)

(b)

(c)

図 3 ℎ(𝑥0 , 𝑦0 ; 𝑥0 + 𝑥, 𝑦0 + 𝑦)の分布[6]．(a) 光路長差を考慮しない場合，(b) 光路長差を考慮した場合（𝜃𝑅 = 2°，𝜃𝑂 = 3°），(c) 光路長差を考慮した場合（𝜃𝑅 = 𝜃𝑂 = 2°）．

た場合におけるℎの分布を示す．図 3(c)の各画像
は，図 3(a)と同様に，(𝑥0, 𝑦0 )に応じて分布中のパターンがシフトしているだけであり，分布そのものは不変であることが分かる．提案手法においては，図 3(c)の(𝑥0 , 𝑦0 ) = (0, 0)における分布を PSFとして利用する．すなわち，提案手法におけるPSF は以下のように表現できる．

1.0E+05
従来手法

104
1.0E+04
1.0E+03

計算時間 [秒]条件を加えた．図 3(c)に，𝜃𝑅 = 𝜃𝑂 ≡ Θ = 2°とし

102
1.0E+02
1.0E+01

2.1×105100
1.0E+00

1.0E-01

高速化比

5.5×105

7.3×104
2.1×10410-2
1.0E-02

提案手法

1.0E-0310-4
1.0E-04
256 5121,0241280

2,0482304

3328

4,0964352

N [画素]

図 4 計算時間の計測結果[6]．

ℎ𝐿𝐼𝐹 (𝑥, 𝑦)𝑑2𝜋exp (𝑖𝑟) (𝑠𝑎𝑡𝑖𝑠𝑓𝑦 (14))
2= {𝑖𝜆𝑟𝜆(𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒)0(13)|𝑥 sin Θ + 𝑑 − 𝑟| < 𝑐𝑊

(14)で必要となる FFT の計算には，CUDA で提供される cuFFT ライブラリを用いた．図 4 に計算領(13)式で与えられる PSF を用いて(8)式に示した
畳み込みを計算することで，光伝搬動画像記録・再生法においてホログラム記録材料上に形成される干渉縞の分布を得ることができる．

4. 提案手法の評価域の画素数と計算時間との関係を示す．ここでは，(3)式を直接的に計算する手法を従来手法とした．
計算領域の画素数を𝑁 × 𝑁画素としたとき，𝑁が512，1,024，2,048，4,096 の場合における結果を示した．図 4 の結果はいずれも，計算時間を 10 回計測して得られた平均値である．従来手法においては，計算時間がそれぞれ1.2 × 101 ，1.5 × 102 ，提案手法の有効性を示すために計算時間を評

1.7 × 103，1.9 × 104秒であった．一方，提案手法

価した．計算機環境は以下のとおりである．OS とにおいてはそれぞれ5.9 × 10−4，2.0 × 10−3，7.8 ×

して Windows 10 Enterprise，CPU として Intel10−3，3.4 × 10−2秒であり，最大で約 56 万倍の高Core-i9 10900K，GPU として GeForce RTX 3090，速化を達成できた．次に，従来手法を基準としたメモリ 64GB を搭載した PC を用いた．プログラ

場合に，提案手法で得られる再生像の画質を二乗

ム 開 発 環 境 と し て Microsoft Visual Studio

平均平方根誤差（RMSE）により評価した．RMSE

Community 2019 と CUDA 10.2 を用いた．(8)式は，得られた再生像を 256 階調化した後で計算した．𝑁が 512，1,024，2,048，4,096 の場合において，それぞれ 38.0，25.6，19.8，15.5 という結果が得られた．光伝搬動画像記録・再生法に限らずBS

超短光パルスレーザ

M

拡散板の開口パターン

コンピュータホログラフィにおいては，散乱光をBE記録する場合には再生像に粒上のノイズ（スペックルノイズ）が重畳し，画質が低下することが知BEられている．この点を考慮すると，従来手法と提拡散板案手法とで再生像画質には大きな差がないことを確認できた．

MBC

次に，提案手法による計算機シミュレーションで得られた結果を，光学実験で得られた結果と定性的に比較した．図 5(a)に示すような光学系を構

イメージセンサ(a)

築してホログラムを取得した．中心波長 522nm，時間幅 178fs のフェムト秒光パルスレーザを光源として用いた．光パルスが伝搬する様子を分かりやすくするため，拡散板には解像力テストターゲ

時間ットのパターンを開口として設定した．画素数(b)2,592 × 1,944画素，画素ピッチ2.2 μmのイメージ

図 5 光学実験の概略[6]．(a) 光学系の模式図，センサを用い，𝑑 = 270 mm，Θ = 2°の条件でホロ(b) 期待される再生像の模式図．

グラムを記録した．図 5(b)には，図 5(a)の光学系

図中の略称号はそれぞれ，BS: ビーム分割素子，で記録したホログラムから得られる 4 次元像の模

M: ミラー，BE: ビーム径拡大素子，BC: ビーム式図を示す．拡散板を用いる場合，光パルスと拡結合素子を表す．

散板との断面像が記録されるため，実際には 3 次元（空間 2 次元＋時間 1 次元）像となる．断面像覧可能であるので，よければそちらも参照されたが左から右へと伝搬していく様子が観察できる．

い．いずれの場合においても図 5(b)で期待された

また，再生像中に現れる解像力テストターゲット

ように，弓状に歪んだ光パルスの断面像が左からのパターンの変化する様子から，光パルスが伝搬右へと伝搬していく様子を観察できた．両者の結する様子を容易に観察できると期待される．ここ

果は定性的によく一致しており，提案手法の有効で，断面像は直線状であるが，1 章でも述べたよ性をじゅうぶん示しているといえる．一方で，光うに得られる 4 次元像には時空間歪みが生じるた学実験から得られた再生像はエッジ部分がぼけめ，図 5(b)に示したように実際の像は弓状に歪む．ているのに対し，計算機シミュレーションで得ら記録したホログラムに対する再生処理についてれた結果はエッジ部分がはっきりしていること

は，紙面の分量の都合からここでは割愛する．詳が確認できた．この差は，計算機シミュレーショ

細は参考文献[6]を参照されたい．図 6 に得られンにおける処理を簡単化する目的で導入していた再生動画像から抜き出した静止画の一部を示た仮定によって生じたものと考えられる．計算機

す．(a)，(b)はそれぞれ光学実験，計算機シミュレ

シミュレーションにおいては，光パルスは矩形波ーションで得られたホログラムからの再生結果状の強度分布を有するという仮定を導入していである．画像間の時間差は 8.2fs に相当する．実た．一方で，実際の光パルスはガウス分布状の強際の再生動画像は参考文献[6]中のリンクから閲度分布を有することが知られている．この違いが，

5 mm時間(a)

5 mm時間(b)

図 6 得られた再生動画像から抜き出した 6 コマの静止画[6]．(a) 光学実験結果，(b) 計算機シミュレーション結果．画像間の時間間隔は 8.2fs．実際の再生動画像は参考文献[6]中のリンクから閲覧可．
再生像におけるエッジのぼけ具合の差に影響を与えたと予想できる．

5. まとめコンピュータホログラフィの応用事例の一つである光伝搬動画像記録・再生法の計算機シミュレーションを，FFT に基づいて高速に処理可能な手法を提案，実証した．計算時間の評価結果から，提案手法は従来よりも約 56 万倍高速化可能であることが示された．RMSE を利用した画質評価に謝辞
本研究は，公益財団法人高柳健次郎財団の助成によるものである．また，本研究の一部は JSPS 科研費（20H05887，20J23542），公益財団法人島津科学技術振興財団の助成を受けたものである．ここに記し，関係各位に謝意を表す．
共同研究者である千葉大学の伊藤智義教授，下馬場朋禄教授，京都工芸繊維大学の粟辻安浩教授，井上智好氏にこの場を借りて感謝申し上げる．

より，提案手法が従来手法と同等の再生像画質を参考文献得られることを確認した．光学実験によって得ら[1] T. Kakue, et al., Opt. Lett. 36, 4131-4133

れたホログラムからの 4 次元像と定性的に比較し(2011).

た結果，提案手法の有効性を示すことに成功した．

[2] T. Kakue, et al., Sci. Rep. 7, 10413 (2017).

今後は，提案手法の拡張性および汎用性を高める[3] Y. Takase, et al., Appl. Opt. 60, A179-A187

ために，本研究において導入していたいくつかの仮定を満たさない場合，例えば𝜃𝑅 ≠ 𝜃𝑂 の条件においても適用可能な，新たな計算手法およびシミュレーションモデルを開発することを目指す．将来的には，4 次元像に生じる時空間歪みを補正し，光パルスの振る舞いを精緻に捉えることが可能な未踏技術の確立へと挑戦する予定である．

(2021).
[4] T. Kakue, et al., IEEE J. Sel. Top. Quantum
Electron. 18, 479-485 (2012).
[5] T. Kakue, et al., OSA Continuum 4, 437-454(2021).
[6] T. Kakue, et al., JOSA A 39, A7-A14 (2021).