はじめに

日本に画像再構成コミュニティを作りたい！

かつて、日本には画像再構成の研究者や画像再構成を学ぶ学生が結構いました。今も一定数はいますが、活発に研究にしていた方が年齢的に第一線を退きつつあり、研究は緩やかに衰退しているように感じます。

対照的に、欧米や中国のコミュニティは非常に活発です。それはもう凄いです。なぜなら、研究者の母数も違いますが、研究者は研究に専念出来ること、画像再構成を学んで働けるポスト・企業が多いことが大きく影響していると感じます。

あとは、画像再構成の方法論に関する取っつきやすい和書が少なく、方法論を体系的にまとめた和書やWeb情報がありません。そして、画像再構成に興味のある人が質問・議論できるコミュニティがないことも、研究の緩やかな衰退を招いているように思います。深層学習コミュニティが羨ましい！

画像再構成の面白さを感じてもらえるように、画像再構成の「古典」から「最先端」までの方法論を、実装の話も踏まえつつ「体系的に」まとめていきます。

コメント・質問歓迎です！なお、私は「数学・物理に面白みを感じるが、決して得意ではない」という人種なので、数学・物理の厳密な観点では間違ったことも書くと思います。見つけたら、コメントください。質問にびしっと答えられないこともあると思います。その際は一緒に考えてもらえると嬉しいです。

なお、全ての記事に共通で、事実やコンセンサスが取れている内容を私見と区別するために、私見は「・・・と思います」「・・・と感じます」と書きます。

画像再構成とは

画像再構成(image reconstruction)とは、人間の目やセンサーによって直接的に計測できない物理量の空間的・時間的な分布(すなわち画像)$\boldsymbol{x}$を、その物理量に関係した別の物理量$\boldsymbol{y}$を頼りに「数学の力を借りて」推定することです。以降、$\boldsymbol{y}$を観測データと呼びます。

観測データ$\boldsymbol{y}$から画像$\boldsymbol{x}$を数学的に推定するのですから、$\boldsymbol{y}$の観測過程を$\boldsymbol{y}=H(\boldsymbol{x})$として数学的に表現できないと画像再構成は出来ません。ここで、$H(*)$は観測過程を表す関数です。以降、$H(*)$をシステム関数と呼びます。

一般に、何らかの入出力システムがあるとき、入力$\boldsymbol{x}$から出力$\boldsymbol{y}$を求める問題を「順問題(forward problem)」、反対に出力$\boldsymbol{y}$から入力$\boldsymbol{x}$を求める問題を「逆問題(inverse problem)」と言います。画像再構成は典型的な逆問題です。

画像再構成を必要とする装置は色々ありますが(別記事でまとめます)、我々の最も身近なところでいくと、(それでも遠いかもしれませんが)Ｘ線ＣＴだと思います。ＣＴとはcomputed tomographyの略で、日本語では「計算機トモグラフィー」や「計算機断層撮影」と訳されます。Ｘ線ＣＴは関心対象物に多数の方向からＸ線を照射して、多数のＸ線の影絵(物体を透過したＸ線の強度画像)を測定します。手を光にかざして、机や壁に影ができることと基本的に同じです。Ｘ線ＣＴでは、観測データである多数の影絵を数学的に巧妙に処理して組み合わせることで、物体の内部を輪切り画像として画像再構成します。Ｘ線によるイメージング装置としては、健康診断でも登場するレントゲン撮影装置の方が馴染みがあると思いますが、レントゲン画像はＸ線の影絵そのものであり、再構成された画像ではありません。

少し脱線しますが、放射線技師にスポットを当てたテレビドラマ「ラジエーションハウス」が放送されていましたが、その中で何度か「放射線技師は医師にオーダーされた写真だけ撮ってればいいんだよ！」的なセリフがありましたが、「写真」という表現には違和感があります。なぜなら、劇中でも何度か登場したレントゲン画像(Ｘ線静止画像)やＸ線動画像は２次元的な画像センサで撮影されたＸ線の影絵なので写真と呼んでもよいのですが、Ｘ線ＣＴ画像やＭＲＩ画像は観測データ(写真的なものと言えなくはない)を再構成した３次元画像なので、２次元の写真とは似て非なるものなんですよね。

画像再構成の数学的枠組み

画像再構成には、様々な数学的枠組み(アプローチ)があります。その中で基本的といえる枠組みは、システム関数$H(*)$が既知であることを前提に、$\boldsymbol{y}=H(\boldsymbol{x})$というシステム方程式を$\boldsymbol{x}$について解くこと、そして、前もって導き出したシステム関数$H(*)$の逆関数$H^{-1}(*)$を$\boldsymbol{y}$に作用させること、の２つです。「逆関数」とさらっと書きましたが、解析学を中心とした数学に関する深い知識と数学的センスがないと逆関数は導出できません(私にはできません)。

さて、観測には誤差がつきものです。誤差を含まない観測データが得られることはまずありません。誤差はランダムな要素(統計誤差)と非ランダムな要素(系統誤差)に分解できます。ランダムは確率論的、非ランダムは決定論的ともいいます。何回も同じ観測を行い、観測データを平均化すれば、再現性のない統計誤差は小さくなりますが、再現性のある系統誤差は変わりません。

したがって、現実的なシステム方程式は$\boldsymbol{y}=H(\boldsymbol{x})$ではなく、$\boldsymbol{y}=H(\boldsymbol{x})+\boldsymbol{b}+\boldsymbol{e}$です。$\boldsymbol{b}$が系統誤差、$\boldsymbol{e}$が統計誤差です。この場合でも「システム方程式を解く」、「逆関数を作用させる」が画像再構成の基本ですが、前もって系統誤差$\boldsymbol{b}$を求めておく必要があります。場合によっては、画像再構成よりも系統誤差推定の方が難しいということもあるかもしれません。

一方、統計誤差$\boldsymbol{e}$はランダムなので、どうあがいても求められません。具体的な値は神のみぞ知るです。しかし、我々人間も、観測系の物理に関する知識を背景に、統計誤差$\boldsymbol{e}$(もしくは観測データ$\boldsymbol{y}$)が従う確率分布の(近似的な)形状を知っている場合があります。

その場合、誤差の確率分布に基づいて、観測データ$\boldsymbol{y}$が得られる確率(生起確率)を画像$\boldsymbol{x}$の関数 $p(\boldsymbol{y}|\boldsymbol{x})$として定式化できれば、生起確率を最大にする$\boldsymbol{x}$が統計的な意味で最も確からしい解といえます。簡単に言えば、観測データに対する「当てはまり」が最良となる$\boldsymbol{x}$が尤もらしいと考えます。観測データの生起確率は、統計学で「尤度」と呼ばれます。

この「尤度最大化問題を解く」という枠組みは、「システム方程式を解く」、「逆関数を作用させる」に次ぐ第３の枠組みです。尤度最大化問題は一般に「最尤推定(maximum likelihood estimation)」と呼ばれ、その解は最尤推定解(maximum likelihood estimator)と呼ばれます。それぞれの枠組みについては、別記事で掘り下げます。

なお、観測データの統計誤差はデノイジング処理により低減することが可能ですが、副作用として系統誤差の増大を伴います。また、観測データには手を加えずに、統計誤差が伝搬した推定結果の方をデノイジングすることも可能ですが、これも系統誤差の増大を伴います。つまり、系統誤差と統計誤差は本質的にトレードオフ(反比例)の関係です。

誤差に頑健な画像再構成

系統誤差と統計誤差のトレードオフの改善方法は大別すると３つあります。

ハードウェアもソフトウェアも変えずにデータ収集の条件を変えて観測データの情報をリッチにする(例:データ収集時間を長くする)
ハードウェアを改善して観測データの情報をリッチにする(例:センサーの高感度化・高解像度化)
ソフトウェアを改善して誤差に対する頑健性を高める

画像再構成を研究する最大の動機は「３」です。追加の情報がないと打つ手なしなので、頑健性を高めるためには、観測データ以外の情報、つまるところ、観測対象自体に関する先見情報が必要です。何らかのコスト(ヒト・モノ・カネ)をかけて、先見情報を取得し、そして画像再構成の計算フローに組み込まなければなりません。

少し脱線しますが、深層学習による画像認識精度の劇的な改善は、まさに先見情報の力ではないでしょうか？ヒト・モノ・カネをかけて大量の学習データ(＝先見情報)を集めて、ヒト・モノ・カネをかけて推論モデルを構築しています。

先見情報を用いた画像再構成の代表的な枠組みは「事後確率最大化問題を解く」です。これが画像再構成の第４の枠組みです。具体的には、先見情報を確率分布の形で(作為的かつ便宜的に)表現した事前確率$p(\boldsymbol{x})$を導出し、観測データの尤度と観測対象の事前確率の積$p(\boldsymbol{y}|\boldsymbol{x})p(\boldsymbol{x})$を最大化する$\boldsymbol{x}$を解とします。事後確率最大化問題は一般に「最大事後確率推定(maximum a posterior estimation)」と呼ばれます。尤度最大化も事後確率最大化も何らかの関数を最大化するだけなので、数値計算の観点では大きな差はありません。

尤度最大化では観測データへの当てはまりをひたすら追求するのに対して、事後確率最大化では、観測データへの当てはまりと先見情報への当てはまりに「バランス」を要求します。事後確率最大化問題には、そのバランス(先見情報をどの程度加味するか)を調整するパラメータ(ハイパーパラメータ)があり、パラメータを０にすると、尤度最大化問題と一致します。よって、尤度最大化問題は事後確率最大化問題の特殊形です。

Ｘ線ＣＴやＭＲＩ(magnetic resonance imaging)において実用化されている圧縮センシング再構成は、スパース性を先見情報に用いた画像再構成の典型例です。

さて、正しい先見情報であれば系統誤差と統計誤差のトレードオフは改善されますが、間違っていれば反対に改悪となるので、先見情報の選択は非常に重要です。どんな先見情報が有効かは完全にケースバイケースで、汎用的な先見情報というものは残念ながら存在しません。つまり、観測対象に特化した事前確率を自ら設計する、もしくは既存の事前確率から良さげなものを選択することになります。

しかし、これはおかしな話ですよね。画像再構成によって観測対象の性質を知りたいのに、前もって観測対象の性質を知っていなければ事前確率を設計・選択できません。これ「ニワトリが先か、タマゴが先か」の状態です。なので、真に正しい先見情報(＝真理)は基本的に知り得ません。したがって、真に正しい事前確率は基本的に定義できません。我々ができることは、「それなりに正しいと期待できる」先見情報を見出だして事前確率を定義し(ハンドメイドし)、それを信じることだけです。なので、前の文章でわざわざ「作為的かつ便宜的に」と括弧書きしました。

事前確率を信じた結果、所望の結果が得られなければ、事前確率が不適切であったと、事後的に結論します。つまり、出たとこ勝負です。この意味で、事後確率最大化問題を解くことによる画像再構成は極めて帰納的な(結果を出発点とするトップダウンの)アプローチになります。他の３つの枠組みは、基本的に演繹的な(観測原理を出発点とするボトムアップの)アプローチです。

学習型画像再構成

画像認識を震源地とする深層学習の爆発的な拡がりは、画像再構成の分野にも波及しており、観測データを画像に変換する画像再構成関数、つまりシステム関数の逆関数ライクな関数を、観測過程に関する数学的・物理的な事前知識を全く入れずに、深層ニューラルネットワーク(ＤＮＮ)で強引に表現し、完全にデータドリブンで学習するという、古典的・今日的な画像再構成の方法論とはかけ離れた、超ラジカルな方法も提案されています。これが「再構成モデルを学習し適用する」という画像再構成の第５の枠組みです。

また、画像再構成を完全にＤＮＮ化するのではなく、画像再構成に関係する定型的な計算をモジュール(重み固定の層)として使ったＤＮＮ再構成や、これまで述べてきた４つの画像再構成の枠組みの中に、(事前学習した)ＤＮＮを組み込んだＤＮＮ応用再構成も提案されています。今後の展開が楽しみです。

しかし、忘れてはならないことは、学習データとなる画像は、基本的に学習に頼らない「非学習型画像再構成」の結果なので、非学習型画像再構成にも更なる進化が求められます。ただ、計算ハードウェアの進化とともに、データ生成やシミュレーションといったソフトウェア技術が高度に進化すると、学習データは「集める」から「作る」にパラダイムがシフトする可能性が高く、そうなると非学習型再構成が駆逐される可能性がもゼロではありません。

Questions

さて、ここまではどんな画像再構成問題にも当てはまるように、あえて具体例を出さず、抽象的に説明してきましたが、抽象的であったが故に、所々で「？」が頭に浮かんだのではないでしょうか。

Q.システム関数はどうやって知るの？
Q.系統誤差はどうやって知るの？
Q.システム関数は非線形でもいいの？
Q.生起確率はどうやって知るの？
Q.先見情報をどうやって事前確率として表現するの？

Q.システム方程式の解は常に存在するの？それは一意(１つ)なの？
Q.システム関数の逆関数は常に存在するの？それは一意なの？
Q.生起確率の最大化解は常に存在するの？それは一意なの？
Q.事後確率の最大化解は常に存在するの？それは一意なの？

Q.システム方程式はどうやって解くの？
Q.逆関数はどうやって求めるの？
Q.尤度最大化問題はどうやって解くの？
Q.事後確率最大化問題はどうやって解くの？

Q.推定結果の系統誤差や統計誤差は評価できるの？
Q.学習型画像再構成って信頼できるの？

答えは・・・

全て「ケースバイケース」

です。答えになっておらず、すいません！しかし、それが事実です。

画像再構成と一言で言っても、その原理や方法は基本的に装置毎に異なります。したがって、画像再構成マスターへの道のりは長く険しいのですが(ちなみに私はマスターにはほど遠いです)、「とある特徴」を共有する装置群では、画像再構成の方法論が共通しており、その方法論を学べば、幅広く応用が出来ます。その装置群とは、直接的・間接的に関わらず「物理量の積分値を計測する」という「線形」な観測過程に基づく装置です。

ちなみに、最後の疑問に関しては「学習データの質や多様性が確保されれば」というのがその場しのぎの政治家的な回答ですが、質や多様性はとても抽象的な表現で、数値的にも評価出来ないと思うので、結局はデプロイ先での性能が全てなんだと思います。

投影からの画像再構成

物理量の積分値を計測する装置の場合、誤差のない観測データ$H(\boldsymbol{x})$は$T(P\boldsymbol{x})$と表現できます。ここで、$P$は行数が積分データの次元$M$、列数が画像$\boldsymbol{x}$の次元(画素数)$N$である$M×N$行列であり、$P\boldsymbol{x}$が積分データを表します。また、$T(*)$は積分データに対する変換関数を表します。この場合、システム方程式は$\boldsymbol{y}=T(P\boldsymbol{x})+\boldsymbol{b}+\boldsymbol{e}$となります。

さらに、変換関数$T(*)$の逆変換$T^{-1}(*)$が既知であれば、画像$\boldsymbol{x}$に関する線形方程式$T^{-1}(\boldsymbol{y}-\boldsymbol{b}-\boldsymbol{e})=P\boldsymbol{x}$が得られます。このように線形方程式に帰着できれば、方程式を数値的に解いて画像再構成することのイメージができるのではないでしょうか。

物理量の積分計測に基づく装置の中でも、物理量の「線積分値」を計測する装置は重要です。なぜ重要かというと、Ｘ線ＣＴをはじめとする断層撮影装置がそのクラスに属するためです。物理量の線積分値はしばしば「投影(projection)」と呼ばれ、Ｘ線ＣＴをはじめとする断層撮影装置は「投影からの画像再構成」と呼ばれる逆問題を解くことで画像再構成します。

今後の記事では「投影からの画像再構成」を中心に、その方法論を体系的にまとめていくつもりです。

おわり