高精度輻射輸送計算と機械学習との融合による近赤外光トモグラフィ-現状とこれから筑波⼤学

計算科学研究センター
⽮島 秀伸
グループメンバー: 梅村雅之、安部牧⼈、⾼⽔裕⼀、佐藤⼤樹
星詳⼦（浜松医科⼤） 、三村徹也（浜松フォトニクス）
計算メディカルサイエンスワークショップ 2021/9/6

1Hillman (2007)より

生体光イメージングとは？
波⻑700-1000ナノメートルの近⾚外線を⼈体へ照射し、透過もしくは戻ってきた光を検出して⽣体内を診断する⼿法
利点
ゼロ被曝
⾮侵襲
ベッドサイドモニタリング

Fig. 1 Intrinsic chromophores and structures of the brain. !a" For small animal imaging, it is possible to use visible light to look at oxygenationdependent hemoglobin absorption !as well as exogenous absorbing and fluorescing dyes". The cortex can be exposed and imaged at highresolution, or imaged through the intact skull and scalp at lower resolution. !b" For human brain imaging, near-infrared light will penetrate more
readily through the scalp and skull to sample the brain thanks to lower scatter and absorption. The human brain has a more complex convoluted
structure compared to lower mammals. Spectra: major chromophores in brain are oxy- and deoxyhemoglobin and water !lipid is omitted here". For
generation of these spectra: water content is assumed to be 90% !quoted value in newborn brain, 71 to 85% in adult brain34". Hemoglobin
absorption is shown assuming 2-mM concentration of hemoglobin in blood and 3% content of blood in tissue !60 mM in tissue34" calculated fromspectra.257 HbR spectrum represents 3% fully deoxygenated blood in tissue. Scatter spectrum is approximate using: !!s = A"−b, where A = 1.14
# 10−7 m and b = 1.3.228

安⼼・安全
新⽣児・幼児へも適⽤可能

ただし、画像再構成には光の伝播について理論モデルが必要
⼤規模な並列計算により可能

may be removed to provide a direct view of an area of thebrain’s surface. In mice, the brain can often be visualized
without even thinning the skull. Methods have been carefully
developed to ensure that normal brain function is not affectedby these approaches. Most often, the animal will remain anesthetized during measurements.23 In some cases, however, apermanent window can be surgically implanted and the animals can survive for over a year while an area of their brain isaccessible for repeated direct optical imaging.14,24
Imaging of pathologies requires only that the animal beanesthetized. Images during acute interventions, such as adJournal of Biomedical Optics

ministration of drugs or introduction of a stroke, must be carefully synchronized with image acquisition. While the brain’sfunctional response to stimulus is being studied, a carefully
controlled stimulus is presented during imaging. Such stimuli
may include visual !e.g., showing patterns or pictures", somatosensory !e.g., electrical or tactile stimulation of the whiskers
or forepaw", or auditory stimuli !different frequencies and amplitudes". Stimuli may be repeated in regular blocks, or presented in a randomized pattern.18,25 Data are usually shown as
the average response to multiple stimulus repetitions to obtain

051402-3

2September/October 2007!

Vol. 12!5"

⽣体光イメージング
プロセス

近⾚外光計測本講演
1.TRINITY開発
図1 生物学的酸素濃度指示物質
2.ファントム計算
cyt. ox.の酸素への親和性はHbやMbより高く、正常細胞におにあるが、がん組織で報告されているような低酸素濃度では還3.Wavelet法/GPU順問題解析て最も感度が高い。酸素濃度の定量的評価には酸化型cyt. ox.濃診断
近赤外光を用いる以外に計測手段はなく、さらに近赤外光計測
輻射輸送計算画像再構成
可能なのは、近赤外光トモグラフィ（near-infrared
optical tomog
[A]

[B]

気管

[C]

画像再構成領域
逆問題解析
機械学習 図2 甲状腺NIROT. A,ヒト頸部計測風景;
B,水平断核磁気共鳴（Mimura
et al. (2021)

数(P4.機械学習(LSTM)a)画像（801 nm）.MR画像の甲状腺位置よりやや前よりでは成されている（実線枠内）; C, 組織酸素飽和度(StO2)の画像 (5.光⾳響トモグラフィ

3
NIROTは、図2Aのように複数の照射－検出ファイバを用い解析することによって光学特性値（吸収係数・散乱係数）をもまる光吸収物質の濃度などの空間分布を画像化する技術である

るX線とは異なり、多重散乱によってさまざまな方向に広がっ

研究の構想（生体光イメージング技術の開拓）
Yajima et al.

①⾰新的技術開発
⼤規模シミュレーション(⾼精度輻射輸送計算、ビッグデータ)②巨⼤データバンクの創設
国際的に開かれた医学・天文学の構造解析拠点へ

edshift at z ∼ 14, 10.4, 8.5, and 6.2, respectively. The box size is 1 Mpc in the comovingof erg s−1 cm−2 arcsec−2.

最先端技術(Wavelet法、 GPU/FPGA
機械学習/深層学習)
研究環境(スーパーコンピュータCygnus、専有PCクラスタ)

国際ひかりイメージング
データバンク
医学計測
天文観測

医療診断
構造解析

luding, in clockwise direction, star formation rate, emergent Lyα luminosity, equivalentcircles) and UV continuum (1300 - 1600) Å, open circles). Open triangles and squaresof the surface brightness with 10−18 and 10−19 erg s−1 cm−2 arcsec−2, respectively.

54

保存する計算データ規模
輻射強度 I: 各場所、各時刻で⾓度bin数分の情報
空間グリッド数: (100)^3、⾓度グリッド数: 3000、スナップショット数: 8000
３次元情報をフルに全て保存する場合
8 byte x 106 x 3000 x 8000 = 192 TB / 1パラメータ
表⾯情報(2次元)のみ、時間保存も1/10にする
8 byte x 104 x 3000 x 800 = 192 GB / 1パラメータ=> 1PBストレージで~5000パラメータ保存可能
JST創発的研究⽀援事業 (PI:⽮島秀伸) の予算にて購⼊予定「宇宙物理輻射輸送計算で拓く新しい⽣体医⽤光学」(2021年度-2027年度)
5大規模輻射輸送シミュレーション(Yajima et al. 2021, submitted to JQSRT)

新汎⽤型輻射輸送計算コード(100%オリジナルで開発)

TRINITY(Yajima, Abe, Umemura, Takamizu, Hoshi, submitted to JQSRT)(Time-dependent Radiation Transfer in Near-infrared Tomography)

Short輻射輸送⽅程式

Ij+1*⾼速⾼精度な光線デザイン*⼤規模並列計算が可能（MPI+OpenMPのハイブリッド並列）

⽣体を1 mmの空間分解、1 ピコ秒の時間分解で計算可能に

TRINITY

Ij-1

Ij

Ij

Ij-1

Ij-1

6

既存の⽅法との⽐較
真空に斜めビーム照射
経路上(y=2cm)に吸収体を置く
Short法(有限要素法)では数値的な拡散により、ビームが広がる
その結果吸収体を⼈⼯的に回り込
んでいる7

散乱がある場合
散乱係数 1.0 cm-1 を⼀様にセット
経路上に⼤きさをもつ吸収体
⼈⼯的な拡散よりも、物理的な拡散が⽀配的になり、スキーム間の差が⼩さくなっている
吸収体背⾯での光の散乱の回り込
みはTRINITYの法が⼩さい
8Figure 7: Left panel: Schematic view of a pulse experiment for phantom. S1-S8 show positions
計測平面：ファントム底面から、35mm, 45ofpulse
injection.
D1-D8
represent
thepositionsof detectors. Right panel: radiation fieldmm, 55 mm
obtained from the radiative transfer simulation. The color shows the photon energy density
in log-scale.

光パルス数値実験

測システム(TRS-80)による計測に用いたファントム
⽣体模擬物質（ポリウレタン製）にトム計測計測に用いたファントム
1. ジュラコン（ポリオキシメチレン）
対するパルス照射実験との⽐較
40 X 40 X 180 mm

Properties of phantom

1. ジュラコン（ポリオキシメチレン）µa = 0.019/mm, µ40s’ X= 400.86/mm(800 nm), g = unknown, n =1.48
X 180 mm
µa = 0.019/mm, µs’ = 0.86/mm (800 nm), g = unknown, n =1.48
Table 1:

Box size

4.0 cmµa0.22 cm 1µs

22.6 cm 1

g0.62

n1.51

Properties of phantom made by polyurethane: µa and µs are absorption andscattering coefficients, g is an anisotropy factor, and n is the refractive index.

18S, 照射；D, 受光、最短SD間隔, 8 mm

Figure 6: Left panel: Schematic
view of a pulse experiment for phantom. S1-S8 show positions
S, 照射；D,計測平面：ファントム底面から、35
受光、最短SD間隔,8 mm
mm, 45Figure 6: Leftofpanel:
SchematicD1-D8
view representof a pulsethe
experiment
fordetectors.
phantom.
S1-S8show
positionspulseinjection.
positionsofRightpanel:
radiation
fieldmm, 55 mm
計測平面：ファントム底面から、35mm, 45
obtainedfrom represent
the radiativetransfersimulation.
The colorRightshows panel:the photonenergy density
pulse
injection.
D1-D8
thepositionsof detectors.
radiationfield
2018/07/20
mm,of55mm
log-scale.
obtained frominthe
radiative transfer simulation. The color shows the photon energy density

in log-scale.

near the injection point. Because of high scattering efficiency, the trajectories

9

Table 1:µs

22.6 cm 1

g

0.62

n1.51

計測に用いたファントム
1. ジュラコン（ポリオキシメチレン）
40 X 40 X 180 mmµa = 0.019/mm, µs’ = 0.86/mm (800 nm), g = unknown, n =1.48

Properties of phantom made by polyurethane: µa and µs are absorption andscattering efficiencies, g is a anisotropy factor, and n is the refractive index.

光パルス数値実験

S, 照射；D, 受光、最短SD間隔, 8 mm
Figure 6: Left panel: Schematic view of a pulse experim
計測平面：ファントム底面から、35mm, 45
pulse
injection. D1-D8 represent the positions of dmm,of55mm

obtained from the radiative transfer simulation. The cin log-scale.

計測データ

near the injection point. Because of high scatt

数値計算

of photons are bent typically with the mean fre

If there is no absorption, the typical number(⌧ (1205g))2 ⇠ 3.2 ⇥ 103 . Thus, the box is filladdition, even after the stop of injection, anare trapped due to multiple scattering processe
crossing time, tcross = 4.0/c/n = 8.8 ⇥ 10 2panel. Thereafter, most photons gradually dis
absorption or escape.
210

Figure 8 presents the comparison of simulati

data. We find that the shape of the output signa

tal data. For the detectors further from the injeat later time. In addition, the full width at hawith the distance from S2. For example, FW
215

0.8 ns for D1. This means that there is wider

各場所での検出光の波形がFigure 8: Time-resolved light signal. Each panel shows the light signal taken from di↵erent
計測データと⼀致
detectors as shown Figure 6. Black solid lines represent the injected pulse. Thick gray and(i.e., traveling time) in the case of D5.

Figure 7:xy two-dimensional map of photon energy density at z = 2.0 cm in log-scale.

The values are normalized by the energy peak at t = 1.4 ns.

1510

red curves show the experimental data and the simulation results.
163.2. Phantom test with an absorber脳・甲状腺計算に向けて(星教授より提供)

MRIより構造抽出

甲状腺 吸収率 x 3

テスト計算結果甲状腺 吸収率 x 30シグナル

甲状腺なし
甲状腺 吸収率 x 3
甲状腺 吸収率 x 10
甲状腺 吸収率 x 30
時間[ns]

11計算高速化 (Wavelet法+GPU化) Abe et al. in prep.
Wavelet法
⾓度分解能を⾃動的かつ局所的に最適化（計算メモリを激減、計算の⾼速化）

0.00012Full-resolution
wavelet

D10.0001
8e-05
6e-05

D5

4e-05
2e-05
0

1

20.9(%) (%)(𝑛) = 𝜆$ 𝜑$ (𝑛) +"(&) (&)∑ ∑𝛾$,( 𝜓$,( (𝑛)&'%(

精度を保ったまま、計算メモリを⼤幅に減らす事に成功

33.5
time [ns]

44.5

5

4.5

5

Full-resolution runに対するメモリ使⽤率0.8
memory utilization𝐼(")

2.5

0.7
0.6
0.50.4
0.3
0.2
0.1
02

2.5

3

3.5
time [ns]4

計算高速化 (Wavelet法+GPU化) Abe et al. in prep.
• GPU版TRINITYコードも開発済み• CPUコードに⽐べて~6-7倍程度の⾼速化に成功• ボトルネックは通信部分 (40%程度を占める)

1 GPU/CPUあたりの計算規模を323 meshに固定、300 stepの計算時間を⽐較

要求メモリ 約5分の1
計算時間はGPUとwaveletを組み合わせれば約10倍程度早くなる→ 10000パラメータのデータ構築へ

CPU (@TRINITY):Intel Xeon Gold 6140
(18 core)機械学習による画像再構成

元の時系列データ

レーション(Takamizu et al. in prep.)
Google TensorFlowを⽤いて機械学習の光伝播シミュレーション
2次元検出モデル
モデルを作成（LSTM法）

機械学習で用いるデータ：8次元の時系列データ

吸収係数はファントムと同じ、直径を半分にした吸収体をファントムに吸収体をあらゆる領域内に敷き詰める

A:吸収率=𝐼𝑎𝑏𝑠
𝐼𝑛𝑜𝑎𝑏𝑠5mm ー＞2.5mm

場所において学習−1

８つの検出器での時系列データを分類する＝教師あり機械学習
検出位置：28か所

ガンの大きさ、形状、個数など、より任意性の高い設定における検出を目指す

検出位置：168か所

14

をる

異なる大きさの５０点未満非表示
２次元でのモデル化＝吸収体の2.5倍＝吸収体の２倍＝1.5倍

キレイに再現されるよ複数かつ異なる異常部位に関しても⾼精度に任意の形の複数体のガ判定出来るようになった
真のガンの位置がかる
15３次元でのモデル化１面モデル

１面モデル

大きい1体も位置がずれ柱方向によく反応しているが、大きい1体も位置がずれている
位置がずれている

３次元検出モデル

Ｚ軸
柱方向によく反応しているが、位置がずれている

個の吸収体㎜の球体）

X軸

１２グループ各グループ 18個の球
＝2段階

Y軸

216箇所の吸収パターンで学習データ作成

12分類×18分類

ある程度正解しているが、形や位置はズレている
現在精度向上のために、検出⾯などいろいろと⼯夫中16

光音響トモグラフィ(PAT:Photoacoustic Tomography)
近⾚外光と⾳波を組み合わせた診断法

近⾚外光

佐藤⼤樹さん(M2)
修⼠論⽂予定光⾳響波を⽤いた⼈の⼿の⾎管の画像
Canon(http://www.canon.com/technology/future/index.html)

⾳波＊光イメージングのみよりもさらに深部を調べられる＊⾎管・⾎流の画像化に強い（脳虚⾎、薬剤輸送）＊条件によっては細胞レベルの可視化

⾚⾎球のイメージング
Shintate et al. (2020)

基礎⽅程式
照射（近⾚外光）

検出（光⾳響波）

輻射輸送⽅程式

光⾳響波動⽅程式

CIP法を使⽤してコードを作成

⽬的：これら輻射輸送、光⾳響波の式を正確に解いたPATモデルを作成する。脳内の虚⾎・出⾎の診断を⽬指す。

ヒト頭部を想定した光⾳響波伝搬シミュレーション
⼆次元空間( 4cm × 4cm)
頭蓋⾻

背景組織• 光学特性値
吸収係数(cm-1)

散乱係数(cm-1)

背景組織

080.0

⾎管

4.758675.0

頭蓋⾻

086.0

⾳速
(cm/μs)密度(g/cm2)

背景組織0.15

1.0

⾎管0.151.0

頭蓋⾻

0.45

1.0

1.0cm

近⾚外光照射位置

2.0cm

2.2cm

⾎管

• ⾳響特性値

1.5cm

ヒト頭部光音響波シミュレーション
3"all_35-25.dat" using 1:2:4

輻射輸送計算の結果
光エネルギーの分布

3.5

102-1

1.5-2

1-3

0.5-4

0-5
0.5

1

1.5

2

2.5

33.54)

2.5

0log10(32頭蓋⾻

4

⾎管による光吸収音響シグナル2.5"pa_2d_cip14_dec_4.dat"

妥当な時間で⾳波が表⾯に到着頭蓋⾻での反射により複数ピークが⾒られる

直接届いた光⾳響波

21.51

多重反射された光⾳響波

0.50-0.5051015202530354045光⾳響波に対する機械学習
ヒト頭部を想定して、光源を1箇所、⾎管を1本配置して光の伝搬・吸収のシミュレーション + 光⾳響波伝搬のシミュレーション

シミュレーションから得られた9箇所の検出位置での光⾳響波の波形を使⽤して機械学習を⾏い、⾎管が複数本ある場合の位置や⼤きさを判断する。

• 光源位置と光⾳響波検出位置は⽚側のみ
光源位置：4箇所
光⾳響波検出位置：9箇所• ⾎管の配置パターン数は
半径：0.25cm → 28パターン
半径：0.5cm → 28パターン
半径：0.75cm → 15パターン

機械学習による⾎管の同定判定された位置に点数を付与し、各光源位置（4箇所）での判定結果を重ね合わせる。

⽩丸が⾎管の正しい位置。
⾎管が⼆本や三本ある時でも近しい位置もしくは、正確な位置を判定できている。
⼀⽅で、まだまだ近しい位置すら判定できていない場合もある。
様々な場合で正確な判定ができるように、精度を向上させていく必要がある。

まとめ＊輻射輸送計算コードTRINITYを開発し、ファントムで検証した＊ Wavelet法+GPU化によって10倍程度の⾼速化、低メモリ化＊光⾳響波の計算コードを開発した＊機械学習で簡単な状況では95%以上の精度で判定可能になった
今後＊甲状腺・脳などの複雑な部位への適⽤＊光イメージング・光⾳響のお互いの⻑所・短所、組み合わせの効果について定量評価＊FPGA化＊ビッグデータ処理について
24