△ 15」/「 一
V) xelX 一
印Vbxel
̲y
ハ
N
A
4ε1'α ツ
4
一
・PnextjPc。rrent ε〃 α.κ
一
■ ■ ■■1
■
γc砂=
, γc∫.κ
1
(》 《 ,
殉 κε1、P
一
殉 κθ13舵 κrκ,ight‑Pκ ll
レわXε18∫zεツ ニ〜そ)xel ̲y
一
図7‑53D‑DDAに よ る ボ ク セ ル 間 の 進 行 距 離 算 出 の 概 要(二 次 元 表 示)
54
7.4散 乱 角 決 定 の た め の ク ラ イ ン ー仁 科 の 公 式 に 即 し た 二 次 元 テ ー ブ ル の 作 成 従 来 のMCシ ミ ュ レ ー シ ョ ン の コ ン プ ト ン 散 乱 処 理 に は,散 乱 角 の 決 定 の た め に 合 成 棄 却 法 が 用 い ら れ る22)。 合 成 棄 却 法 を 用 い た コ ン プ トン 散 乱 処 理 で は,条 件 分 岐 や 複 数 回 の 乱 数 生 成 を 含 み,散 乱 角 度 が 採 用 さ れ る ま で 反 復 処 理 が 行 わ れ る 。 よ っ て,GPUコ ン ピ ュ ー テ ィ ン グ で はwarpdivergenceが 生 じや す い 。 そ こ で 開 発 コ ー ドで は,入 射 光 子 エ ネ ル ギ ー に 対 す る 散 乱 角 の 確 率 密 度 分 布 を 累 積 確 率 関 数 に 変 換 し た 二 次 元 テ ー ブ ル を 作 成 し,一 様 乱 数[0,1]
の 値 に 対 応 し た コ ン フ.ト ン 散 乱 角 を サ ン プ ル す る ア ル ゴ リ ズ ム で 効 率 化 を 図 る 。2章 に 記 述 し た よ う に,MV‑X線 領 域 に お け るIncoherentscatteringfunctionS(x,紛 の 補 正 の 効 果 は 非 常 に 小 さ い た め,コ ン プ トン 散 乱 角 の 二 次 元 テ ー ブ ル は ク ラ イ ン 仁 科 の 式 に よ る 微 分 断 面 積 の み を 用 い て 作 成 し た16)。
ク ラ イ ン ー仁 科 の 式 に よ る 単 位 立 体 角 に お け る 電 子 あ た り の 散 乱 断 面 積 は,光 子 エ ネ ル ギ ー をE,散 乱 角 を θ,古 典 電 子 半 径 をro=2.818×10'13cmと し た 場 合,次 式 で 表 さ れ る 。
dσ 宿 d9=7
1 1一 ト α(1‑cosθ)
2
、+,。s2θ+α2(1‑… θ)21
+α(1‑c・sθ)
(7.8)
こ こ で α は 一E/0.511を 表 す 。 立 体 角 と散 乱 角 と の 間 に はd9ニ2πsinθdθ の 関 係 が あ る た め,式 (7.8)は 以 下 の よ う に 変 形 で き る 。
dσ=γ3πsinθ 1
1一 ト α(1‑cosθ) 2
、+。 。s2θ+α2(1‑… θ)21
十 α(1‑COSθ)
dθ
(7.9)
よ っ て,(θ α≦ θ ≦ θう)の 条 件 を 満 た す 球 帯 断 面 積 σ(θ)は次 式 よ り 求 ま る 。
σ(θ)一 倉 血 θ
α
1 1十 α(1‑cosθ)
2
、+,。s2θ+α2(1‑… θ)2
1+α(1‑c・sθ)
dθ
(7.10)
算 出 した 球 帯 断 面 積 σ(θ)を正 規 化 し,コ ン プ トン散 乱 角 θの 確 率 密 度 分 布 を 取 得 した 。図7‑6 に 光 子 エ ネ ル ギ ー 対 す る コ ン プ トン 散 乱 角 の確 率 分 布 を数 例 示 す 。 そ の 後,図7‑7に 示 す よ
うに確 率 密 度 分 布 を累 積 確 率 分 布 に 変 換 した 。 一 様 乱 数[0,1]と 正 規 化 され た 確 率 を 関連 づ け る こ と に よ り,一 様 乱 数[0,1]の 値 に対 応 した 累 積 確 率 分 布 の 散 乱 角 の値 を 参 照 す る よ うに し た 。
(xlO‑3) 2.0
1.5
m
ヨ O σ3 ρ o 乙
0.5
0.OlMeV O.lMeV
‑・ ・一 一1 .OMeV
‑一・一一2 .5MeV 6.OMeV
0.0
0306090120150180
scatterangle[deg]
図7‑6光 子 エ ネ ル ギ ー に よ る コ ン プ ト ン 散 乱 角 の 確 率 分 布
180
150
0 2 1 90 60
[bD Φ ℃ ﹂ ￡ bD 自 憩 § の
30
0.OlMeV O.lMeV
… ・‑1
.OMeV
‑一 一一 一2 .5MeV
0
0.00.20.40.60.81.O
probability
図7‑7光 子 エ ネ ル ギ ー に よ る コ ン プ ト ン 散 乱 角 の 累 積 確 率 分 布
56
7.5GPUに よ る 乱 数 生 成
乱 数 と は,生 成 さ れ た 数 か ら次 の 数 が 任 意 の 時 点 ま で 予 測 不 可 能 で あ り,さ ら に 規 則 性 が な い た め 状 況 の 再 現 が 不 可 能 な 数 列 の こ と を さ す 。 よ っ て コ ン ピ ュ ー タ 上 で 生 成 さ れ た 乱 数 は,数 値 の 予 測 が 困 難 か つ 数 値 の 偏 りが 小 さ い 必 要 が あ る 。 コ ン ピ ュ ー タ に よ り生 成 さ れ た 乱 数 は 擬 似 乱 数 と 表 さ れ,RandomNumberGenerator(RNG)に よ り生 成 さ れ る 。RNGの 良 し 悪 し は,周 期 の 長 さ,分 布 の 偏 り,メ モ リ使 用 量 お よ び 速 度 で 決 ま る 。 全 て の 擬 似 乱 数 に は 周 期 が あ り1周 期 が 終 了 す る と ま た 同 じ 数 列 を 繰 り 返 し 生 成 す る 。MC計 算 は 多 く の 乱 数 の 生 成 を 必 要 と す る ア ル ゴ リ ズ ム の た め,周 期 が 短 い 疑 似 乱 数 で は 計 算 結 果 が 偏 り,正 確 な 処 理 が 不 可 能 と な る 。
MC計 算 に 一 般 的 に 用 い ら れ る 高 精 度RNGと し てMersienneTwister(MT)ア ル ゴ リ ズ ム が あ る28)。MTは 最 大2216091‑1の 非 常 に 長 い 周 期 を も ち,均 等 な 分 布 を 可 能 と す る 。 ま た, 乗 算 ・除 算 な ど の 処 理 の 遅 い 命 令 を 避 け る ア ル ゴ リズ ム と な っ て お り,CPUで は 比 較 的 高 速 で あ る 。 た だ し,2492バ イ トの 状 態 ベ ク トル を 使 用 す る た め,GPUへ 実 装 す る た め に は 容 量 の 大 き いglobalmemoryへ 格 納 す る 必 要 性 が あ る 。 そ こ で,GPU向 け に 改 善 さ れ たMTア ル ゴ リ ズ がSaitoら に よ り 開 発 さ れ て い る29)。 こ れ は,CUDAの 乱 数 ラ イ ブ ラ リ(CURAND) に 実 装 さ れ て お り,51200threadで 制 限 が か か る よ う に な っ て い る 。
ま た,Xorshiftア ル ゴ リ ズ ム は,Marsagliaら に よ っ て2008年 に 導 入 さ れ た30).こ の ア ル ゴ リ ズ ム は,周 期 はMTよ り短 い が 状 態 ベ ク トル が 小 さ く 処 理 がMTよ り速 い 。CURANDラ イ ブ ラ リ で は,わ ず か な 修 正 が 加 え られ てXOR‑shiftwithWeylsequence(XORWOW)ア ル ゴ リ ズ ム と し て 実 装 さ れ て い る 。 開 発 コ ー ドで は,CUDAに 実 装 さ れ たRNGの な か で,乱 数 生 成 が 高 速 か つMTの 次 に 周 期 が 長 いXORWOWを 用 い る こ と に し た 。
7.6GPU光 子 輸 送MonteCarloシ ミ ュ レ ー シ ョ ン ア ル ゴ リ ズ ム
図7‑8にGPUMCコ ー ドの フ ロ ー チ ャ ー トを 示 す 。GPUMCシ ミ ュ レー シ ョ ン で は,並 列 数 を 設 定 し,CPUメ モ リ か らGPUメ モ リ ヘ デ ー タ を 転 送 し た 後,ホ ス トか ら デ バ イ ス を 呼 び 出 す こ と に よ りkemelに お け る 並 列 計 算 を 施 行 す る 。 デ バ イ ス に お け る 光 子 輸 送 の 並 列 処 理 に よ り,設 定 し た 全 ヒ ス ト リ ー 数 が 終 了 し た 後,GPUメ モ リか らCPUメ モ リヘ サ ン プ ル デ ー タ を 転 送 し,出 力 す る 。
7.6.lMonteCarlo計 算 前 の 処 理
GPUMCシ ミ ュ レ ー シ ョ ン で は,光 子 輸 送 計 算 に 入 る 前 にCPUメ モ リ 上 の デ ー タ をGPU メ モ リへ 転 送 さ せ る 必 要 が あ る 。 転 送 す る デ ー タ を 以 下 に 示 す 。
●
●
●
●
●
シ ミ ュ レ ー シ ョ ン に 用 い るCT画 像(3Dボ ク セ ル 画 像)
3D‑CT画 像 を マ ク ロ セ ル 法 に よ り 変 換 し た 画 像(3Dマ ク ロ セ ル 画 像) 光 子 エ ネ ル ギ ー に 対 す る 断 面 積 デ ー タ
(相 互 作 用 の 種 類:光 電 効 果,コ ン プ ト ン 散 乱,電 子 対 生 成,全 断 面 積) 光 子 エ ネ ル ギ ー に 対 応 し た コ ン プ ト ン 散 乱 角 の 累 積 確 率 分 布
(cumulativeprobabilitydistributionofComptonscattering;CPDc。mpt。n)
線 源 か ら 射 出 さ れ る 一 次 光 子 の エ ネ ル ギ ー ス ペ ク トル の 累 積 確 率 分 布 (cumulativeprobabilitydistributionofenergyspectrum;CPDprimary)
3Dボ ク セ ル 画 像 と3Dマ ク ロ セ ル 画 像 は 容 量 が 大 き い た めglobalmemoryに 格 納 す る 。 断 面 積 デ ー タ,散 乱 角2Dテ ー ブ ル お よ び 光 子 の エ ネ ル ギ ー ス ペ ク トル は,光 子 の エ ネ ル ギ ー に 対 応 し た 各 々 の 離 散 デ ー タ をtexturememoryに 格 納 し,texturememory特 有 のbuild‑in interpolationを 用 い る こ と に よ り線 形 補 間 し た デ ー タ を 効 率 良 く 取 得 可 能 に し た 。 そ の 他, GPUメ モ リ に あ ら か じ め 確 保 す る 項 目 を 記 す 。
・ 各threadに お い て 追 跡 す る 光 子 の 情 報(光 子 エ ネ ル ギ ー,光 子 の3D座 標,単 位 ベ ク ト ル,光 子 の 状 態 を 表 す フ ラ グ)を 格 納 す る た め の メ モ リ領 域
・ 各threadに お い て 用 い る 乱 数 の 値 を 格 納 す る た め の メ モ リ領 域
・ 画 像 処 理 の た め に 保 存 す る 情 報(EPIDに 入 射 し た 一 次 光 子,散 乱 光 子,再 投 影 光 子 の エ ネ ル ギ ー とEPID上 で の 二 次 元 座 標)を 格 納 す る た め の メ モ リ領 域
上 記 の 三 項 目 は,thread数 だ け 確 保 す る た め 大 容 量 を 必 要 と す る 。 こ れ ら の デ ー タ の た め の メ モ リ領 域 はglobalmemory上 に 確 保 す る が,3章 に 記 載 し た よ う にglobalmemoryは レ イ テ ン シ が 大 き い 。 そ こ で,本 研 究 で はCUDAに 組 み 込 ま れ たC++ベ ー ス のStandardTemplateLibrary
(STL)で あ るthrustを 用 い る31)。thread数 に 依 存 す る 膨 大 な デ ー タ は,globalmemory上 でthrust devicevectorと し て 処 理 す る こ と に よ り,メ モ リへ の 読 み 書 き だ け で な く,ソ ー トや 削 除 な ど の デ ー タ 整 理 を 効 率 良 く 行 う こ と が 可 能 と な る 。 ま た,kemelに お け る 並 列 処 理 で は,デ ー タ をglobalmemoryか らregisterに 格 納 す る こ と に よ り計 算 時 の デ ー タ ア ク セ ス を 効 率 化 さ せ た 。
58
HOST
CPUコ ン ピ ュ ー テ ィ ン グSTART
各 種 デ ー タ をCPUメ モ リ へ 格 納 シ ミ ュ レ ー シ ョ ン に 用 い るCT画 像(3Dボ ク セ ル 画 像)
・ 光 子 エ ネ ル ギ ー に 対 す る 断 面 積 デ ー タ
(相 互 作 用 の 種 類:光 電 効 果,コ ン プ ト ン 散 乱,電 子 対 生 成,全 断 面 積) 光 子 エ ネ ル ギ ー に 対 応 し た コ ン プ ト ン 散 乱 角 の 累 積 確 率 分
(cumulativeprobabilitydistributionofComptonscattering;CPDc。mpt。n)
・ 線 源 か ら 射 出 さ れ る 光 子 の エ ネ ル ギ ー ス ペ ク トル の 累 積 確 率 分 布 (cumulativeprobabilitydistributionofenergyspeCtrum;CPDprimary)
マ ク ロセ ル 法 に よ りCT画 像 か らマ ク ロセ ル 画 像 を 生 成(3Dマ ク ロ セ ル 画 像)
光 子 輸 送 計 算 時 に 用 い るGPUメ モ リ の 確 保 光 子 情 報 用
(光 子 エ ネ ル ギ ー,光 子 の3D座 標,ベ ク トル,光 子 の 状 態 を 表 す フ ラ グ) 生 成 した 乱 数 の 格 納 用
・ 光 子 サ ン プ リ ン グ デ ー タ の 格 納 用
(EP■)に 入 射 し た 一 次 光 子 、 散 乱 光 子 、 再 投 影 光 子 の エ ネ ル ギ ー お よ び EPII)入 射 座 標,(hVp,x,y),(hv、,x,y)(hv,,x,y))
GPUメ モ リ容 量 に応 じた 並 列 蜘 の設 定
lthreadあ た り の ヒ ス ト リ ・一・lkN,,,eadの 設 定
全 ヒ ス ト リ ー 数1隣 。t、1の決 定Nt。t、1ニNth,e、d×M
回
図7‑8(a)GPUMCコ ー ド の フ ロL‑一一・チ ャ ー ト
前項続き
!
DEVICE
GPUコ ン ピ ュ ー テ ィ ン グ 、
種 々 の デ ー タ をGPUメ モ リ へ 転 送(CPUmemorytoGPUmemory)
, 冒
Firstkemelの 実 行 命 令 , Firstkemelの 並 列 計 算 実 行
■
「F
一
Secondkemelの 実 行 命 令 , Secondkemelの 並 列 計 算 実 行
層 u
Thirdkemelの 実 行 命 令 , Thirdkemelの 並 列 計 算 実 行
ヒ ス ト リ ー 数 し な い く妬ね1の追 跡 終 了
す る
層
光 子 サ ン プ リ ン グ デ ー タ をCPUメ モ リ へ 転 送 (GPUmemorytoCPUmemory)
ー︑
サ ン プ リン グデ ー タの保 存
ノ
サ ン プ リング デ ー タ格 納 GPUメ モ リの再 初 期 化
図7‑8(b)GPUMCコ ー ド の フ ロ ー チ ャ ー ト
60
7.6.2ヒ ス ト リ ー 数 の 設 定
GPUMCシ ミ ュ レ ー シ ョ ン で は,lthreadあ た り1光 子 の 追 跡 が 行 わ れ る 。 よ っ て 並 列 に 輸 送 さ れ る 光 子 の 数 は,「lblockあ た り のthread数 」 と 「block数 」の 積 と な る 。 並 列 数M‑thread 数 ×block数 と お く と,並 列 数MはGPUの メ モ リ 容 量 に よ っ て 上 限 が 与 え ら れ る 。 全 ヒ ス ト
リ ー 数Nt。t、1が 並 列 数Mよ り 小 さ い 時,lthreadあ た り の ヒ ス ト リ ー 数 は1と な る 。 一 方,Nt。t、1 がMよ り 大 き い 時,lthreadあ た り 射 出 す る 一 次 光 子lftN,,,eadは 次 式 よ り 算 出 さ れ る 。
Nt。tal Nthread=
M
(7.ll)
算 出 さ れ たNth,eadの 小 数 点 以 下 を 切 り 捨 て た 場 合,設 定 し た ヒ ス ト リ ー 数 に 対 し て 統 計 不 足 と な る 。 し か し,Nth,eadは 整 数 値 で 定 義 し な け れ ば い け な い た め,Nth,eadは 整 数 値 で 算 出 し た の ち1を 加 算 す る こ と に よ りNt。t、1以上 の 十 分 な 光 子 数 を 射 出 し た 。 よ っ て,Nt。t。1は 「lthread あ た り の ヒ ス ト リ ー 数 」,「lblockあ た り のthread数 」 お よ び 「block数 」 の 積 と な る 。 開 発 コ ー ドで は 全threadの 射 出 数 をNth,eadに 設 定 す る こ と に よ り,全threadを 可 能 な 限 り ア ク
テ ィ ブ な 状 態 に 保 つ よ う に し た 。
7.6.3kernel関 数 の 反 復 計 算 に よ る 光 子 輸 送
開 発 し た コ ー ドで は,デ バ イ ス に お け る 並 列 処 理 の た め のkernel関 数 を3つ に 分 け,kemel を 反 復 す る こ と に よ り,マ ク ロ セ ル 内 の 光 子 輸 送,相 互 作 用 判 定 お よ び 光 子 サ ン プ リ ン グ を 施 行 す る 。
(a)Firstkemel
図7‑9に,Firstkernel関 数 の フ ロ ー チ ャ ー ト を 示 す 。Firstkernelは,一 次 光 子 お よ び 散 乱 光 子 の エ ネ ル ギ ー と 単 位 ベ ク トル の 決 定 の た め に 実 行 さ れ る 。 一 次 光 子 お よ び 散 乱 光 子 の 進 行 方 向 の 算 出 は,一 部 に 同 じ 計 算 過 程 が 含 ま れ る た め,同 一 関 数 内 で 処 理 す る こ と に
よ っ てwarpdivergenceの 低 減 を 図 る 。
ま ず,一 次 光 子 に 関 す る 処 理 に つ い て 記 述 す る 。 図7‑10に 一 次 光 子 の 極 角 θ お よ び 方 位 角 レ 算 出 の 概 念 図 を 示 す 。 従 来 のMC計 算 で は,(θ,V)は 一 様 乱 数 に よ り ラ ン ダ ム に 決 定 さ れ る 。 一 方,本 研 究 に お け るGPUMCコ ー ド で は,処 理 の 効 率 化 を 図 る た め,thread
とblockの イ ン デ ッ ク ス(threadld,blockldx)を ピ ク セ ル 座 標(Up,Vp)に 対 応 さ せ て,全 ピ ク セ ル に 一 定 量 の 光 子 を 同 時 に 射 出 さ せ る 手 法 を と る 。threadの 次 元 は 一 次 元 に 設 定 し, threadの イ ン デ ッ ク スthreadldx.xをEPIDの σ 方 向 の ピ ク セ ル と 対 応 さ せ る 。 ま た,block 次 元 は 二 次 元(blockldx.x,blockldx.y)に 設 定 し,blockの κ 方 向 の イ ン デ ッ ク スblockldx.xを EPIDのV方 向 の ピ ク セ ル に 対 応 さ せ る 。 こ れ に よ り,EPID平 面 の ピ ク セ ル 数 が σ ×V(1≦
Up≦U,1≦ γp≦V)の 時,ピ ク セ ル 座 標(Up,Vp)は(Up,γp)=(threadldx.x,blockldx.x)の 関 係 と な る 。 本 研 究 で 用 い たGPUのlblockあ た り のthreadldx.xの 最 大 は1024で あ り,EPIDの σ 方 向 の ピ ク セ ル 数 の 二 倍 で あ る 。 よ っ て,blockldx.κ も 同 様 にV方 向 の ピ ク セ ル 数 の 二 倍 の768 に 設 定 し,threadldx.xお よ びblockldx.κ を2で 除 算 し て 小 数 点 以 下 を 切 り 捨 て る こ と で, 各 ピ ク セ ル(Up,Vp)へ 光 子 を4個 同 時 に 射 出 で き る よ う に し た 。 ま た,blockldx.yの 最 大 値 を メ モ リ の 容 量 が 許 す 限 り 大 き く し,一 度 に(Up,Vp)へ 向 け て 射 出 さ れ る 光 子 数 を さ ら に 増 加 さ せ た 。 し た が っ て,各 ピ ク セ ル(Up,Vp)に 同 時 に 射 出 さ れ る 光 子 数 は4×blockldx.yに な る 。 こ れ に よ り,GPUMCシ ミ ュ レ ー シ ョ ン 全 体 の 並 列 数Mは,threadldx.κ ×blockldx.x×blockldx.y と な る 。
62
Firstkemelの 並 列 計 算 実 行
1
一 様 乱 数 の 生 成 ア lr2[0,1]
↓
フ ラ グ が 一 次 光 子 射 出
コ ン プ トン散 乱
threadldx.xお よ びblockldx.xに 対 応 し た ピ ク セ ル 座 標(Up,Vp)を 距 離 単 位(u,v)に 変 換 し 、 線 源 一EPID間 距 離1か ら一 次 光 子 の 極 角 θお よ び 方 位 角 レを 算 出
θ一一(翠)ψ 一t・n‑・(老)
一 様 乱 数r
lと エ ネ ル ギ ー ス ペ ク トル のCPDprim、ryに よ り、 各threadに お け る 一 次 光 子 エ ネ ル ギ ーEを サ ン プ リ ン グ
一 様 乱 数r
lと 散 乱 前 の 光 子 エ ネ ル ギ ・一一Eに よ っ て コ ン プ ト ン 散 乱 角 の 累 積 確 率 分 布CPD、 。npt。.
か ら散 乱 極 角 θを サ ン プ リ ン グ
散 乱 光 子 エ ネ ル ギ ・一一・E'の 算 出 E
1+α(1‑cosθ)
暫 定 的 相 互 作 用 半1」定 の た め の 値 算 出 一109(1‑r2)
Cutoff エ ネ ル ギ ー
以 上
未満
散 乱 方 位 角 レの 算 出 ψ=2πr2
極 角 θお よ び 方 位 角 レよ り光 子 の 単 位 ベ ク トル(κ。,y。,z、)を算 出
Secondkemelの 並 列 計 算 へ
Firstkemelの 1.一 次 光 子 光 子 射 出 へ
図7‑9Firstkemel関 数 の フ ロL‑一一・チ ャ ー ト