GPU を用いた波動光学計算用ライブラリ GWO ライブラリの開発とその応用 Numerical calculation library for wave optics, the GWO library, using the graphic processing unit and its applications 下馬場朋禄*1，三浦潤也*1，佐藤芳邦*1，阿部幸男*2，白木厚司*3，市橋保之*2， 中山弘敬*2，増田信之*2，杉江崇繁*4，高田直樹*5，伊藤智義*2 Tomoyoshi Shimobaba, Junya Miura, Yoshikuni Sato, Yukio Abe, Atsushi Shiraki, Yasuyuki Ichihashi, Hirotaka Nakayama, Nobuyuki Masuda, Takashige Sugie, Naoki Takada, Tomoyoshi Ito *1 山形大，*2 千葉大，*3 木更津高専，*4 東京工科大，*5 湘北短期大 Yamagata University, Chiba University, Kisarazu National College of Technology, Tokyo University OF Technology, Shohoku College ABSTRACT In optics, several diffractions are used for calculating scalar light propagation. The calculation result provides us with the generation of a computer-generated-hologram (CGH), the numerical reconstruction image from a hologram, the optical characteristics of an optical device, and so forth. The acceleration of the calculation commonly uses the fast Fourier transform; however, in order to compute a CGH and real-time reconstruction from holograms, recent computers do not have sufficient computational power. In this paper, we develop a numerical calculation library for the diffraction integrals using the graphic processing unit (GPU), the GWO library, and report the performance of the GWO library. The GPU chip allows us to use a highly parallel processor. The maximum computational speed of the GWO library is about 20 times faster than a personal computer. Keywords: CGH，GPU，GPGPU，ホログラフィ，回折計算，計算機合成ホログラム １．はじめに 各 種 光 学 素 子 ， 回 折 光 学 素 子 ， 光 MEMS （Micro Electro Mechanical Systems）といった 光 学 部 品 の 光 学 特 性 解 析 [1][2] や ， CGH （Computer Generated Hologram）の高速計算 [3][4]，CCD カメラなどで撮影したホログラム からの像再生（ディジタルホログラフィ）[5]， 電子ホログラフィックディスプレイ設計用 CAD ツールにおける再生像計算[6]など，回折 計算は光学の幅広い分野で使用されている． 下馬場 朋禄 [email protected] 山形大学大学院理工学研究科情報科学専攻 〒992-8510 山形県米沢市城南 4-3-16 TEL 0238-26-3729 この回折計算をプログラムに容易に実装する ために，各種の回折計算や付加機能などを実 装したC++クラスライブラリLightWave があ る[7]． 回折計算は FFT（Fast Fourier Transform）に よる高速計算アルゴリズムが知られているが， 現在の高性能な CPU を使用したとしても，リ アルタイムでの CGH 計算やディジタルホロ グラフィを行うことは困難である． 一方，GPU（Graphics Processing Unit）は，3 次元グラフィックス処理に特化した専用ハー ドウェアであったが，近年では，浮動小数点 演算を行う小型プロセッサをそのチップ内に 多数搭載し，プログラム可能な高並列 1 チッププロセッサへ進化している．この特徴 に着目し GPU 上で汎用的な数値計算を行う GPGPU（General-Purpose GPU）（GPU コン ピューティングとも呼ばれる）は，最近のハー ドウェアによる高速数値計算の一手法として 注目を集めている． そこで本論文では，GPU 上に各種の回折計 算を実装することで，GPU に関する知識の無 いユーザに対しても手軽に GPU パワーを享 受できる回折計算用ライブラリの開発とその 性能について述べる．回折計算は用途に応じ て各種の計算方法が提案されているが，現在 の GWO ラ イ ブ ラ リ は ， Fraunhofer 回 折 ， Fresnel 回折（畳み込み表現とフーリエ変換表 現の 2 種類），角スペクトル法，Shifted-Fresnel 回折の計 5 種類の回折計算をサポートする． GWO ラ イ ブ ラ リ は ， Windows 上 の DLL （DynamicLink Library）として実装されており， 既存のプログラミング環境（例えば，Microsoft 社の Visual Studio 環境など）との併用も可能 となっている． ２．回折計算 2.1 角 ス ペ ク ト ル 法 ・ Fresnel 回 折 ・ Fraunhofer 回折 角スペクトル法については次式を計算すれ ばよい． u (m2 , n2 ) = FFT −1[FFT[a (m1 , n1 )]× exp(iz k 2 − 4π 2 (( sΔf x ) 2 + (tΔf y ) 2 ) ] フレネル回折積分は，その表現方法によって ２種類あり，畳み込み表現は， [ FFT exp(iπ ((m1Δx) 2 + (n1Δy ) 2 ) /(λz ) y2 u x 2, y2 Aperture a x1, y1 y1 x2 w avelength z x1 Fig.1 Arrangement of aperture and observation planes となる．また，フーリエ変換表現は， u (m2 , n2 ) = C2 FFT[a(m1 , n1 ) × exp(iπ ((m1Δx) 2 + (n1Δy ) 2 ) /(λz )] となる．フラウンホーファ回折は， 図 1 に回折計算における開口面と観察面の 配置図を示す．回折積分の数値計算を行う場 合は，各変数を離散化し，フーリエ変換に FFT を用いる． まず，開口面座標を ( x1 , y1 ) = ( m1Δx, n1Δy ) ， 観察面座標を ( x2 , y2 ) = ( m2 Δx, n1Δy ) と離散化 する． Δx, Δy は各座標軸のサンプリング間隔 を表す．複素振幅分布（開口関数）を a ( m1 , n1 ) ， 開口面によって回折された光を観察する面を u (m2 , n2 ) とする．開口面と観察面の距離を z ， 光の波長を λ ，波数を k = 2π / λ とする． また FFT と逆 FFT の演算子をそれぞれ FFT, FFT −1 と し ， 周 波 数 領 域 で の 座 標 を ( f x , f y ) = ( sΔf x , tΔf y ) とする．Δf x , Δf y は周波 数領域でのサンプリング間隔を表す． N x , N y は開口面と観察面の x 軸， y 軸方向の画素数 を表す． u (m2 , n2 ) = C1 FFT −1[FFT[a(m1 , n1 )]× Observation Plane ] u (m2 , n2 ) = C3 FFT[a(m1 , n1 )] となる． C1～C3 は複素数の定数で，詳細は紙 面の都合上，割愛する．各回折積分を比較す ると，フーリエ変換回数が異なり，また，観 察面の座標変換を伴うものがあることに注意 する必要がある．各回折計算の特徴を表 1 に まとめる．詳細は文献[8]などを参照されたい． Table.1 Characteristic of diffraction calculations Number Type Memory Angular Spectrum 2Nx×2Ny 2 Fresnel(Convolution) 2Nx×2Ny 3 same Fresnel(Fourier) Nx×Ny 1 1/λz Fraunhofer Nx×Ny 1 1/λz Shifted-Fresnel 2Nx×2Ny 3 arbitary of FFT Sampling same 2.2 Shifted-Fresnel 回折 2.1 節で述べた回折計算は，開口面と観察面 のサンプリング間隔を自由に設定することは できない．表 1 のように開口面と観察面のサ ンプリング間隔は同一か，もしくは，観察面 のサンプリング間隔が 1 /(λz ) 倍されたものと なる． 近年，この制約に縛られないShifted-Fresnel 回折と呼ばれる新しい回折計算手法が提案さ れた[9]．この手法は，Fresnel回折から導出さ れているため，距離計算に近軸近似を用いて いるが，開口面と観察面のサンプリング間隔 を自由に設定できるため新しい応用が期待で きる．Shifted-Fresnel回折は次式を計算すれば よい． u (m2 , n2 ) = Cs FFT −1[ FFT [ A(m1 , n1 )] × FFT [ B(m1 , n1 )]] ここで， C s ， A( m1 , n1 ) ， B ( m1 , n1 ) は， Cs = π 1 exp(ikz ) exp(i ( x22 + y22 )) × iλ z λz 2π exp(−i ( x01m2 Δx2 + y01n2 Δy2 )) λz π A( m1 , n1 ) = a ( m1 , n1 ) exp(i ( x12 + y12 )) × λz π 2 exp(i ( x2 + y22 )) × λz 2π exp( −i ( x x + y y )) × λz 1 o 2 1 o 2 exp( −iπ ( s x m12 + s y n12 )) B (m1 , n1 ) = exp(iπ ( s x m12 + s y n12 )) ここで，Δx1 , Δy1 は開口面のサンプリング間 隔， Δx2 , Δy2 は観察面のサンプリング間隔， ( xo1 , yo1 ) と ( xo 2 , yo 2 ) は，原点からの開口関数 と観察面のオフセット座標を表す．また s x = Δx1Δx2 /(λz ) ， s y = Δy1Δy2 /(λz ) ， x1 = m1Δx1 + xo1 ， y1 = n1Δy1 + yo1 ， y2 = m2 Δx2 + xo 2 ， y2 = n2 Δy2 + yo 2 と定義し た． ３．GPU の概要 GWO ラ イ ブ ラ リ は ， NVIDIA 社 の CUDA[10][11][12]に対応したGPUチップ上で 動作する．CUDAは， C言語ライクな言語で GPUのプログラミングが行える．CUDA対応の GPUの構造を図 2 に示す． ホストコンピュータは，CUDA 対応の GPU を高性能なコプロセッサとして扱うことがで きる．ホストコンピュータと GPU ボード間は PCI-Express バス経由で接続され，DMA（Direct Memory Access）操作によりデータ転送がボト ルネックにならないよう工夫されている．計 算対象となるデータは，GPU ボード上のデバ イスメモリへ格納する． Fig.2 Structure of GPU for CUDA GPU が計算を行う際は，その処理をスレッ ドと呼ばれる処理単位に分割を行う．図中の マルチプロセッサは 8 個のストリームプロ セッサ（SP）で構成されており，同じマルチ プロセッサに属する SP は同一の命令で異なる デ ー タ を 処 理 す る SIMD（ Single Instruction Multiple Data）的な動作をする． 分割されたスレッドの各マルチプロセッ サ・SP への割り当てはハードウェアで自動的 に行われる．GPU がデバイスメモリのデータ を読み込むには，GPU のクロックで数 100 ク ロックのレイテンシを要する．スレッド化に より，あるスレッドがこのレイテンシを待つ 間に，他のスレッドが SP を使えるように切り 替えが行われ，より効率的な処理が行える工 夫が施されている． このように，GPU-メモリ間の広帯域データ 転送能力と，複数個の SP が多数のスレッドを 並列処理することにより，GPU は計算の高速 化を図っている．本論文で使用した GPU は Geforece8800GTS で，12 個のマルチプロセッ サと 96 個の SP が搭載されている． ４．GWO ライブラリ 現バージョンのGWOライブラリは表 1 に示 す 5 つの回折計算をサポートしている．回折 計算はFFTを多用するが，GWOライブラリ内 部ではGPU上でFFTを実行できるCUFFTライ ブラリを使用している[13]． 図 3 に GWO ライブラリを使用した Fresnel 回折（畳み込み表現）のソースコード例を示 す．gwoInit 関数は初期化関数で，この関数の 引数として以下のパラメータを指定すると回 折積分の種類を指定できる． ・GWO_FRESNEL_CONV Fresnel 回折（畳み込み表現） ・GWO_FRESNEL_FOURIER Fresnel 回折（フーリエ変換表現） ・GWO_FRAUNHOFER Fraunhofer 回折 ・GWO_ANGULAR 角スペクトル法 ・GWO_SHIFTED_FRESNEL Shifted-Fresnel 回折 gwoSetPitch と gwoSetWaveLength はそれぞ れ，開口面のサンプリング間隔と入射光の波 長を設定する．gwoSendData で開口データを GPU へ送信する．gwoCalc は gwoInit で指定し た回折計算方法に基づいて計算を行う．計算 終了後に，で指定した回折計算方法に基づい て計算を行う．計算終了後に，gwoRecieveData で GPU から計算結果を回収する．このように， 数ステップの手続きのみで，GPU の演算性能 を享受できるライブラリとなっている． 表 2 に各回折計算（1024×1024 サイズ）を CPU 単体で実行させた場合と GWO ライブラ リで実行させた場合の計算速度の比較を示す． Table.2 Calculation time Type CPU only GWO Angular Spectrum 1873(ms) 130(ms) Fresnel (Convolution) 2995(ms) 153(ms) Fresnel (Fourier) 97(ms) 27(ms) Fraunhofer 95(ms) 25(ms) Shifted-Fresnel 3030(ms) 182(ms) GWOライブラリの利用して，ディジタルホ ログラフィック顕微鏡で撮影したホログラム からの像再生や[14]，CGHのリアルタイム生成 などを行っている． まとめ GPU に関する知識の無いユーザに対して も手軽にGPU パワーを享受できる回折計算 用ライブラリの開発を行った．本ライブラリ はSourceForge.jp上で近日配布する予定である [15]． 参考文献 [1] T.P.Kurzweg, S.P.Levitan, P.J.Marchand, J.A.Martinez, K.R.Prough, D.M.Chiarulli, “A CAD Tool for Optical MEMS”, Proc.36th ACM/IEEE conf. on Design automation, 879-884 (1999) gwoInit(GWO_ FRES NEL_ CONV,WIDTH,HEIGHT);//初期化 gwoS etPitch(10.0e-6,10.0e-6);//ピッチ設定 gwoS etWaveLength(633.0e-9);//波長設定 gwoS endData(ape);//GPUに開口面データを送る gwoCalc(0.2);//フレネル回折計算（0.2m伝播） gwoIntensity();//回折パターンの強度計算 gwoRecieveResult(ape);//GPUから強度データ取得 Fig3. Example of source code using the GWO library [2] T.P.Kurzweg, S.P.Levitan, J.A.Martinez, M.Kahrs, D.M.Chiarulli, "An Efficient Optical Propagation Technique for Optical MEM Simulation", Fifth International Conference on Modeling and Simulation of Microsystems (MSM2002), 352-355 (2002) [3] K.Matsushima,"Computer-generated holograms for three-dimensional surface objects with shade and texture", Appl. Opt., 44, 22, 4607-4614 (2005) [4] 坂本雄児，長尾智大，“パッチモデル用いた 計算機合成フーリエ変換ホログラムの高速 計算法”，電子情報通信学会誌，J85-C，3， 150-157 (2002) [5] T.C.Poon (ed.), ”Digital Holography and Three Dimensional Display – Principles and Applications”, Splinger (2005) [6] 三浦潤也，下馬場朋禄，” 電子ホログラ フィックディスプレイの構築支援シミュ レータの開発”，3 次元画像コンファレンス 2008 講演論文集 (2008) [7] 松島恭治，“オブジェクト指向波動光学シ ミュレーションライブラリ LightWave”， Optics Japan 2004 講演予稿集, 304 – 305 (2004) [8] Okan K. Ersoy : ” Diffraction, Fourier Optics And Imaging”, Wiley-Interscience (2006) [9] R.P.Muffoletto, J.M.Tyler and J.E.Tohline, “Shifted Fresnel diffraction for computational holography” Opt. Express, vol.15 No.9 5631-5640 (2007) [10] NVIDIA, "NVIDIA CUDA Compute Unified Device Architecture Programming Guide Version 1.1”, nVIDIA, (2007) [11] 下馬場朋禄，伊藤智義，”CUDA 技術を利用 した GPU コンピューティングの実際（前 編）”，インターフェース 6 月号，CQ 出版 社，2008 年 6 月号，pp.144-153 (2008) [12] 下馬場朋禄，伊藤智義，”CUDA 技術を利用 した GPU コンピューティングの実際（後 編）”，インターフェース 7 月号，CQ 出版 社，2008 年 8 月号 (in print) [13] NVIDIA, "CUDA FFT Library Version 1.1 Reference Documentation”, NVIDIA, (2007) [14] 佐藤芳邦，下馬場朋禄，三浦潤也，伊藤智 義，“リアルタイム・ディジタルホログラ フィック顕微鏡” 3 次元画像コンファレン ス 2008 講演論文集 (2008) [15] https://sourceforge.net/projects/thegwolibrary/
© Copyright 2024 ExpyDoc
