が挙げられる．これらの問題が解決されれば，医療用 CT では低被曝化・測定時間短縮・高画質化，また電子線トモグラ
フィーでは測定時間短縮・試料損傷低減などの様々な長所が圧縮センシングを用いたトモグラフィー画像再構成
Tomographic Image Reconstruction Using
Compressed Sensing

生まれる．そこで，歴史的にこれらの問題を解決するために膨大な数の研究が行われたが，2000 年代中盤までは解決が困難であると考えられてきた 1,2）．
これに対して，2006 年に Donoho や Candes らにより圧縮センシング（Compressed Sensing, CS）と呼ばれる不足した測定データや SN 比が低い測定データから高精度または厳密な信号（画像）復元を行う逆問題の新解法が考案され，上述

工藤 博幸ａ，ｂ，董建ｃ，加茂 勝己ｄ，ｄ
堀井 則孝 ，古河 弘光ｄ，波多聰ｅ，
村山 光宏ｆ，佐藤 和久ｇ，宮崎 伸介ｈ

では，Donoho や Candes らの論文が発表されてから直ちに

Hiroyuki Kudo, Dong Jian, Katsumi Kamo, Noritaka Horii,

線量を落として低被曝で撮影を行う低線量 CT の画像再構成

Hiromitsu Furukawa, Satoshi Hata, Mitsuhiro Murayama,

法として実用化されるに至った．また，電子線トモグラフィー

Kazuhisa Sato and Shinsuke Miyazakiの分野では 2010 年以降 Midgley らや Murayama らがいち早

ａｂ

筑波大学システム情報系情報工学域

JST-ERATO，百生量子ビーム位相イメージングプロジェクトｃ

筑波大学システム情報系

ｄ（株）システムインフロンティアｅ

九州大学大学院総合理工学研究院

ｆｇ

バージニア工科大学材料工学科

大阪大学超高圧電子顕微鏡センターｈ（株）メルビル

要旨

医療用 CT や電子線トモグラフィーの分野において，少数方向の投影データ・角度欠損がある投影データ・SN比が低い投影データから高画質の画像を再構成する手法として注目されている圧縮センシング（CS）を用いた画像再構成法に関して，数学的原理を中心として平易に解説する．また，最近筆者らが開発した電子線トモグラ
フィー用途の CS を用いた画像再構成アルゴリズム ISER（Iterative SEries Reduction）を紹介する．

キーワード：コンピュータトモグラフィー（CT），電子線トモグ
ラフィー，画像再構成，圧縮センシング，トータルバリエーションの問題が解決できる見通しが出てきた 3,4）．医療用 CT の分野
膨大な数の CS 画像再構成に関する研究が行われ，例えば X

く CS の適用を提唱する論文を発表して，その後研究のアクティビティやユーザの期待は急速に増大している 5～8）．特に，材料系電子顕微鏡分野を牽引する Midgley らは，画像再構成法の研究を SIRT 法から開始して，DART（Discrete Algebraic
Reconstruction Technique）法 9），そしてその後 CS に行き着いた研究歴をもち，成功を収めている 5,6）．更に，CS の数学的理論に関しても体系化や性能向上が進み，幾つかの書籍も出版されている 10～12）．
本論文では，CS に馴染みのない電子線トモグラフィー分野の読者を対象として，CS の数学的原理を中心に平易に解説し，更に筆者らが電子線トモグラフィーへの応用を目的として開発した CS 画像再構成アルゴリズム ISER（IterativeSEries Reduction）を紹介する．
2.

圧縮センシング（Compressed Sensing, CS）の数理

2.1

CS の原理まず，CS を適用する典型的な例である少数方向や角度欠損がある投影データからの CT 画像再構成を想定して，CS の原理を述べる．画像の画素値を一列に並べた J 次元ベクトルをx，測定した投影データ（電子線トモグラフィーの場合は連続傾斜像データ）を一列に並べた I 次元ベクトルを y，x と1.

はじめにy を関係づける I × J システム行列を A = {aij} で表す注 1）．このとき，画像再構成は線形方程式 Ax = y を解き y から x を

筆者は，CT（コンピュータ・トモグラフィー）において求める逆問題であるが，少数方向や角度欠損がある投影デー測定した投影データから断層像を生成する画像再構成法を専タしか測定できない場合には I << J となり，Ax = y の解は門とする研究者である．CT 技術における長年の未解決問題無数に存在するため再構成画像 x は一意に定まらない．このとして，1）投影方向数が制限された少数方向投影データかような状況において高品質な（または厳密な）画像の復元をらの画像再構成，2）角度欠損がある不完全投影データから可能にする逆問題解法が CS である．もちろん，投影データの画像再構成，3）線量を落として（または露光時間を短くy の情報だけでは Ax = y の解は一意に定まらず，画像 x にして）測定した SN 比が低い投影データからの画像再構成，関する何らかの先験情報が必要となる．先験情報を用いて逆ａ〒 305–8573 茨城県つくば市天王台 1–1–1
E-mail: kudo@cs.tsukuba.ac.jp
2015 年 12 月 14 日受付
48

問題の解を一意に定まるようにすることを正則化（regularization）と呼び，CS では対象画像が『疎（スパース，sparse）』である性質を正則化に利用する．例として，楽器音の連続波
顕微鏡 Vol. 51, No. 1（2016）【著作権者：公益社団法人 日本顕微鏡学会】

勾配変換に対応する W の例を示す．

形を測定標本値の集合から再生する問題を考える．正しく波形を再生するには，通常シャノン標本化定理に基づき細かく次に，具体的な画像再構成の手法について述べる．CS では，標本化したデータが必要である．しかし，図 1（a）に示すよa）Wx のどの成分が零に近く無視できるかを同定すること，うに，楽器音は通常フーリエ変換した空間で少ない個数の周b）Wx の零成分を無視して x を復元すること，の 2 つを行波数の正弦波の和で表現される性質を利用することで，未知う必要がある．a）と b）を同時に単一の数学的枠組みで実パラメータの個数を大幅に削減して粗く標本化したデータか現するため，ベクトル x の L1 ノルムが疎な性質を上手く評らでも正しい再生が可能になる．この例から分かるように，価できることに注目した，以下に述べる巧妙な解法が用いら疎であるとは，x にスパーシファイ（sparsify）変換と呼ばれる．ここで，パラメータ p ≧ 0 に対して，Lp ノルムとはれる線形変換 W を施したベクトル Wx の多くの成分が零ま次式で定義されるベクトルの大きさやベクトル間の距離を評たは零に近い無視できるほど小さい値をとる性質である．無価する量である（p をノルムの次数と呼ぶ）．

視できる成分を零とみなせば Ax = y における未知数の数を1/ p⎧⎪⎛ K⎞⎪⎜p⎟⎪⎟t⎜⎪ ∑ k ⎟⎟ (p > 0のとき)⎠⎪⎪⎝⎜ k=1t p =⎨1/ ε⎪⎪Kε⎞⎟⎟ (p = 0のとき)⎪⎪ lim ⎛⎜⎜tk∑⎟⎪⎪ε →+0 ⎜⎝⎟⎠
k=1⎪⎩

J から少数（上手くいけば I 以下）に削減することができ，
Ax = y の解を一意に定める（または一意に定まりやすくする）ことが可能になる．具体的なスパーシファイ変換 W としては，望遠鏡で撮影した天体の画像は背景が真っ暗で非零（2）の画素数が少ないので恒等変換，楽器音は少数個の周波数の正弦波の重ね合わせで表現できるのでフーリエ変換，自然画

CS の解法としては，『合成（Synthesis）型』と呼ばれる手像の場合はエッジ部を除いた同一領域内でほぼ一様な濃度値法と『分析（Analysis）型』と呼ばれる手法の 2 種類が存在をもつのでエッジ検出を行う濃度勾配変換やウェーブレットすることが知られている 13）．どちらを利用するかは応用分野

変換など，扱う信号（画像）の種類により様々なものが考えにより趣向が分かれるが，CT 画像再構成では後述するトーられる．以降では，スパーシファイ変換を K × J の行列 W（通タルバリエーションに代表される分析型がよく用いられるた

常は K ≧ J）で表し次のように表現する．

め，以降では分析型の解法について述べる．分析型の解法でt = Wx

は，画像 x を未知数とする以下の最適化問題を解くことに（1）

ただし，t = (t1, t2,…, tK)T は変換係数を一列に並べた K 次元ベクトルである（T は転置を表す）．CS では，スパーシファより，画像再構成を行う．
［測定データ y に雑音がない場合］(イ変換 W の選択が性能に大きな影響を与える．疎なベクト

minimize W x p)psubject to Ax = y（3）ルに変換する効果が大きければ，未知数の数をより少数に削減でき CS の性能は向上する．具体的な W のイメージをつ［測定データ y に雑音がある場合］

かんでもらい以降の説明を容易にするため，図 1（b）に後述
するトータルバリエーション正則化で用いられる画像の濃度(

minimize β W x p) + ( Ax − y )
p2

（4）

2

ただし，式（3），（4）におけるノルムの次数 p の選択に関しては，以下に述べる理由により CS では p = 1 が好んで用いられる．ベクトル t の次元が 2 の場合（t = (t1, t2)T）における p = 0, 1, 2 に対応する L0, L1, L2 ノルムの等値線を図 2 に示す．この図から分かるように，L0, L1 ノルムは L2 ノルムと比較して座標軸から離れた疎でないベクトルを過度に大きく評価する性質をもつ．よって，L0, L1 ノルムを用いることで，座標軸に近い疎な解を上手く抽出することができる．また，本来ベクトルの非零成分の個数を表すのは L0 ノルムであり，これを用いたいが非凸関数で厳密な最小化が困難であり，ほとんど全ての CS では L0 ノルムを近似する凸関数の L1 ノル（4）の定式化は，Wx が疎になムが用いられる注 2）．式（3），るように罰金（ペナルティ）を課して Ax = y を解くことに相当する．
また，式（4）において β > 0 は正則化の強さ（疎な度合い）を調節するハイパーパラメータであり，β が大きいとアーティファクトや雑音の低減効果は大きくなるが画像中の小さ図1

楽器音信号と画像に対するスパーシファイ変換 W の例．

最近の研究と技術

圧縮センシングを用いたトモグラフィー画像再構成な構造物が消滅しやすくなり，逆に β が小さいとアーティ
49【著作権者：公益社団法人 日本顕微鏡学会】

図 2 スパーシファイ変換を行った係数ベクトル t の次元が 2の場合（t = (t1, t2)T）における L0, L1, L2 ノルムの等値線．座標軸から離れた疎でない点は，L1, L0 ノルムでは原点からの実際の長さよりも等値線が中に入り込んでいるため値が格段に大きく評価され，解として選ばれにくくなる．

図3

TV ノルムの中で使用される水平差分と垂直差分の定義．

下のように定式化される．
まず，これまでは 3 次元画像を表すのが x であったが，以降では説明を複雑にしないため 3 次元物体の一断面のトモグファクトや雑音の低減効果が小さくなる．最適な β を自動決

ラフィー画像再構成を考え，2 次元画像を x で表す．図 3 に定する研究も存在するが，ほとんど全ての場合においてユー示すように，注目する画素 j の周囲で水平方向と垂直方向のザが経験的に β を決定しているのが現状である．筆者らは，差分をとる演算をベクトルの内積で表現したものを，各々投影データ y の値が［0.0, 1.0］の範囲に入るようにスケー

x・ejH, x・ejV で表す（ejH, ejV の定義は図 3 を参照）．このとき，リングした後に画像再構成を行っているが，典型的な β の値

画像 x の TV ノルムは次式で定義される．

として 0.001 ～ 10.0 の範囲で上手く動作する．

J

次に，式（3），（4）の解を求める反復法について述べる．px TV = ∑= 1, K = J で W が直交行列の場合に限って，反復ソフトしき

j=1

い値処理（Iterative Soft-Thresholding）法と呼ばれる単純な構(xie ) + (xie )H 2jV 2j（5）即ち，画像の各画素 j の周囲で水平差分と垂直差分の 2 乗和2,10～12）．

を求めルートをとったものを全ての画素 j について加算した

一方，p = 1, K > J のより一般の場合に解を厳密に求めるにのが TV ノルムである注 3）．これは，次式のように考えてみは，一般化反復ソフトしきい値処理法や Chambolle-Pock ア

ると，連続画像 f(x, y) の濃度勾配の大きさを表す関数 g(x, y)ルゴリズムなど数学的に高度な最新の反復法が必要となの L1 ノルムと解釈することができ，式（3），（4）の評価関造の反復法を用いて，厳密な解を求めることができる

る

14,15）

数を拡張したものである注 4）．

．

CS の理論的側面で興味深い問題は，測定データ y に雑音が含まれない単にデータ数が少ない状況において，不足した測定データから数学的に厳密な画像復元が本当に可能かということである．Donoho や Candes らは CS の最初の論文におい2⎛ ∂f ( x, y) ⎞⎟ ⎛⎜ ∂f ( x, y) ⎞⎟⎟+⎜g ( x, y) = ⎜⎜⎟⎝⎜ ∂x ⎠⎟ ⎜⎝ ∂y ⎠⎟⎟2x TV ≈ ∫∫ g ( x, y) dxdy = L1 norm of g ( x, y)（6）て，システム行列 A が RIP（Restricted Isometry Property）と呼ばれる条件を満足すれば，L1 ノルムを用いた CS により数そして，雑音の有無の各々の場合について，画像再構成を次学的に厳密な画像復元が可能であることを証明した 3,4）．RIPの最適化問題として定式化する．

条件については，興味がある読者は原著論文を参照されたい．
2.2

トータルバリエーション正則化

minimize x TV subject to Ax = y（7）minimize β x TV + ( Ax − y 2 )（8）

CS の中にも多数の異なる手法が存在するが，CT 画像再構成では『トータルバリエーション（Total Variation, TV）正2則化』と呼ばれる手法が最もよく用いられる．発見者の名前にちなんで ROF（Rudin, Osher, Fatemi）モデルとも呼ばれ，
16）

TV 正則化は，画像の濃度変化が平らな部分のみを平滑化し．TV 正則て，濃度勾配が大きなエッジ部を強力に保存する性質をもつ．

化は上述の分析型 CS の特別な場合に相当し，具体的には以この性質は，式（6）において濃度勾配の大きさ |g(x, y)| を

元々は雑音除去を目的に提案された手法である

50

顕微鏡 Vol. 51, No. 1（2016）【著作権者：公益社団法人 日本顕微鏡学会】

2 乗した L2 ノルムを使わずに L1 ノルムを使って評価していることに由来する．エッジ部では |g(x, y)| の値が大きくなるが，これを 2 乗して L2 ノルムで評価すれば（2 乗することで）過度に大きく評価されそれを最小化しているためエッジはぼやけて復元される．一方，2 乗せずに L1 ノルムで評価すれば小さ目に評価され最小化しても上手くエッジが保存される．その効果は実際の画像処理問題において驚くほど大きくCS 最大の長所であり，『エッジ保存平滑化』または『L1 マジック』などと呼ばれる．
式（7），（8）の解を求める反復法の構築については，多くの研究が存在する 14～19）．この問題の困難さは，TV ノルムは一項でもルートの中が零になるような画像 x において微分不可能な関数であることに起因する 注 5）．初期の研究では，
TV ノルムを次式の微分可能な関数で近似して最急降下法など微分可能最適化の手法を適用するのが主流であった．
J ⎛22⎞x TV ≈ ∑ ⎜⎜ (xieHj ) + (xie Vj ) + ε − ε ⎟⎟⎜⎠⎟
j =1 ⎝（9）

図 4 胸部 CT 撮影を模擬した数値ファントムの 16 方向投影データからの再構成画像．

ただし，ε > 0 は小さな定数である．しかし，この手法ではTV 正則化が本来もっている上述のエッジを強力に保存する性質が働きにくい問題があり，2000 年代後半以降は微分不可能最適化の手法で厳密に最小化する研究が主流に方向転換した．この周辺の研究の数は膨大であるが，高度な内容になり割愛する．
図 4 に，医療用胸部 CT 画像を模擬した数値ファントム（モデル画像）を，16 方向投影データから TV 正則化により再構成した画像の例を示す．CT 画像再構成の標準的手法であるフィルタ補正逆投影（FBP）法と ART 法で再構成した画像と比較しているが，TV 正則化の再構成画像はほぼ完全であり驚きの画質改善効果があることが分かる．
なお，以上の CS の説明では測定した投影データが不足し

図 5 2 方向投影データから TV 正則化で正しく再構成されない例（右画像を TV 正則化で再構成すると TV 値がより小さい左画像に誤って再構成される）．ている状況を想定したが，式（4），（8）の定式化を用いれば線量を落として測定した低 SN 比の投影データからの画像再水平と垂直の 2 方向投影データから TV 正則化で画像再構成構成（雑音除去，デノイジング）にも有効である．

を行えば，左画像と右画像の投影データは全く同じでかつ左2.3

CS の問題点と限界

画像の TV 値は右画像の TV 値より小さいため，中央の構造

Donoho や Candes らにより CS の概念が提唱されてから約物が消滅して左画像として誤って再構成される．これは簡単

10 年が経過して，その問題点や限界も次第に明らかになっな例であるが，方向数が増えても同様な例はいくらでも作れてきた．以降では，CS の主要な 2 つの問題点について述べる．

る．筆者の経験でも，測定条件が悪かったり正則化を強くかa）CS を過信しないけ過ぎたりした場合に，CS を適用すると小さな構造物が消CS を用いた CT 画像再構成の研究は 2010 年頃に過熱気味滅する事例によく遭遇した．即ち，RIP 条件が満足される特になり，たった 4 方向の投影データで厳密な画像再構成が可

別な場合を除いて，CS は過信しないで用いることが大切で能と主張する論文まで出現した 17）．これに対して，ART 法あると言える．筆者らは，0.001, 0.005, 0.01, …, 1.0, 5.0, 10.0の考案者である Herman らは，『少数方向投影データからのの刻みで β の値を変化させ画像再構成を行い最も自然な結果

画像再構成は零空間（投影データが全て零になる非零画像集を選び，その画像が十分な画質でない場合は CS の限界と判合が作る空間）の特徴づけが 1970 年代に行われており，厳断する，という考え方をとっている．

密な再構成が可能なはずがない』，という批判を行っていb）滑らかな濃度変化・テクスチャー・エッジの再現性

る 18）．即ち，システム行列 A は厳密な画像再構成のためのCS 画像再構成は，エッジ部を除く同一領域内の濃度値がRIP 条件を満足しない．図 5 に，Herman らの主張を分かり完全に一様である『区分的に一様（piecewise constant）』なやすく説明する簡単な例を示す。この図において，右画像の画像モデルを仮定することに相当する．よって，対象画像に最近の研究と技術

圧縮センシングを用いたトモグラフィー画像再構成

51【著作権者：公益社団法人 日本顕微鏡学会】

滑らかな濃度変化やテクスチャーが存在する場合，滑らかなの TV 正則化を用いた画像再構成法の収束が遅い理由は，同変化が失われ階段状のガクガクした関数として誤って復元さ時反復型の構造をもつためと考えられる．そこで，Row-れる問題点が発生する．それに加えて，測定条件が悪い少数Action 型の反復法で収束を高速化してかつ厳密に TV 正則化

データや低 SN 比の場合には，エッジ座標が誤って推定されの解が求められるように改良した反復法が ISER 高速化の骨エッジ部がギザギザに復元される問題点が発生する．図 6 に，組みである．ISER の具体的なアルゴリズム導出は，文献をTV 正則化によりこれらの問題点が発生した顕著な例を示す．

参照されたい 19）．

左側の 2 種類の数値ファントムを低 SN 比の 32 方向投影デー

図 7，8 に，ISER による電子線トモグラフィー画像再構

タから再構成したのが右画像である．上段の例では，滑らか成の例を示す．試料は Xe イオン照射を施した FePt/Al2O3 薄な濃度変化を無理に階段状関数で近似した再構成画像になっ膜（磁性材料）で，映っている構造物は Fe-Pt ナノ粒子であている．この問題点を解決する手法の開発は，現在におけるる 22）．TEM で測定した明視野の連続傾斜像から画像再構成

20）．一方，図 6 の下段の例を行った．図 7 は，角度範囲が –65 ～ 65 度に制限された 32

では，エッジ座標が誤って推定されエッジ部がギザギザに再方向連続傾斜像から再構成した一断面の画像を表している

CS の大きなテーマになっている構成されている．

3. アルゴリズム ISER（Iterative SEries Reduction）の紹介
最近筆者らは，ISER という名称の TV 正則化に基づく画像再構成アルゴリズムを開発した．以降では，ISER について紹介する．
まず，ISER とは，既存の TV 正則化は反復法の収束が遅く計算量が膨大になる問題点を解決して，約 100 倍の収束の高速化を図った新手法である．具体的には，以下の考え方で高速化を実現している．式（7），（8）の形に定式化した CT
画像再構成に用いる反復法の構造は，全ての投影データを同時に用いて画像更新を行う同時反復型（SIRT 法はこれに属する）と投影データを部分データ（ブロック）の集合に分割してブロック毎に画像更新を行う Row-Action 型またはブロック反復型（ART 法はこれに属する）の 2 種類が存在す
る 2）．Row-Action 型の反復法は，投影データのアクセス順序を工夫すると高速に収束することが知られている 2,21）．既存

図 6 CS で滑らかな濃度変化やエッジが上手く再現できない画質限界を示す例（低 SN 比の 32 方向投影データから TV 正則化で再構成すると，滑らかな濃度変化は階段状に，エッジはギザギザに再構成される）．

52

図 7 Fe-Pt 試料の電子線トモグラフィー実連続傾斜像（測定角度範囲 –65 ～ 65 度，方向数 32）から再構成した断面像の比較（構造物はナノ粒子，傾斜像は白黒反転して表示）．

図 8 Fe-Pt 試料の電子線トモグラフィー実連続傾斜像（測定角度範囲 –65 ～ 65 度，方向数 13）から再構成した 3 次元画像のボリュームレンダリング画像（構造物はナノ粒子）．

顕微鏡 Vol. 51, No. 1（2016）【著作権者：公益社団法人 日本顕微鏡学会】

（傾斜像は白黒反転して表示している）．FBP 法と ART 法の再構成画像と比較しているが，ISER では滑らかな濃度変化やテクスチャーが失われている反面，角度欠損による縦軸方向の形の歪み，粒子境界のぼやけ，ストリークアーティファ
クトはなく粒子形状が鮮明に再構成されている．図 8 は，13方向連続傾斜像から再構成した 3 次元画像をボリュームレンダリングした画像を表している．FBP 法で顕著に現れるストリークアーティファクトが ISER では少なく，全データの131 方向連続傾斜像から再構成した画像に近く粒子形状が再構成できている．
上述の実データの画像再構成（32 方向連続傾斜像）における計算時間は，Core i7-3770 3.4 GHz CPU を搭載したデス
クトップ PC を用いて，FBP 法が 35（min），高速化を行わない同時反復型の TV 正則化である一般化反復ソフトしきい値処理法が 1,959（hour），ISER が 21.67（hour）であった．
FBP 法や SIRT 法と比較すると計算時間が長いものの，図 7，8 で示した絶大な画質改善効果をもつ CS 画像再構成を，ISER の開発によって実用可能な計算時間にまで大幅に短縮することに成功した．
4.

まとめ

本論文では，電子線トモグラフィー分野において今後研究開発や実用化が進むと予想される，CS を用いた CT 画像再構成について解説した．また，最近筆者らが開発した電子線トモグラフィー用途の CS を用いた画像再構成アルゴリズム
ISER を紹介した．読者の方々の役に立てば幸いである．
謝辞

図 7，8 の電子線トモグラフィー連続傾斜像は，九州大学松村晶教授（工学研究院）から提供を受けたものであり，深謝いたします．また，本研究の一部は，JSPS 科研費 15K06103，
JST 先端計測「マテリアル開発系リアルタイム電子線トモグラフィーシステムの開発」の助成を受けたものです．
注
注 1） 電子線トモグラフィーにおいて M × M 画素の L 枚の連続傾斜像から N × N × N 画素の 3 次元画像を再構成する場合，ラスタ走査で N × N × N 個の画素値を一列に並べた J = N × N × N
次 元 ベ ク ト ル が x，M × M × L 個 の 測 定 値 を 一 列 に 並 べ た I =M × M × L 次元ベクトルが y になる．画素数 N は，画像再構成に先立ってユーザがどの程度の空間分解能のトモグラムが必要かを考慮して決定する．
注 2） Lp ノルムは p ≧ 1 の場合には極小解が存在しない最小化が容易な関数となり，このような関数を凸関数と呼ぶ．一方、0 ≦ p < 1の場合には極小解が存在する関数となり，このような関数を非凸関数と呼び厳密な最小化は大変難しくなる．例えば，p = 1 で次元が 2の場合，式（3）は直線 c1 t1 + c2 t2 = y 上で |t1| + |t2| を最小化する問題，式（4）は β(|t1| + |t2|) と (c1 t1 + c2 t2 – y)2 の和を最小化する問題となり，凸関数の最小化問題である．一般の J 次元でも同じであり，p = 1 を用いれば極小解に陥る心配は不要である．
注 3） 3 次元画像再構成の場合には，1）2 次元断面画像の TV ノル

最近の研究と技術

圧縮センシングを用いたトモグラフィー画像再構成

ムを全スライスについて加算する手法，2）スライス間の差分も考慮して式（5）を拡張した 3 次元画像の TV ノルムを使用する手法，の 2 種類がある．筆者らは簡便な 1）の手法を用いている．
注 4） 式（6）において g(x, y) ≧ 0 であるため第 2 式の絶対値記号は不要であるが，L1 ノルムであることを気づかせるため意図的につけた．
注 5） ||x||TV の第 j 項は [(xj + 1 – xj)2 + (xj + N – xj)2]1/2 であり，xj + 1 に
関 す る 偏 導 関 数 は (xj + 1 – xj)/[(xj + 1 – xj)2 + (xj + N – xj)2]1/2 と な り xj =xj + 1 = xj + N の座標（つまり x・ejH = 0, x・ejV = 0 となる画像 x）において不連続（どの方向からその座標に近づくかで微係数が –1 と+1 の場合が存在）となり微分できない．同様に計算すると，同じ座標で xj + N と xj の偏導関数も不連続となり微分できない．
文献

1）Rangayyan, R., Dhawan, A.P. and Gordon, R.: Appl. Opt., 24, 4000–4012 (1985)
2）工藤博幸：第 2 章『投影からの画像再構成』，日本医用画像工学会（編），医用画像工学ハンドブック，日本医用画像工学会，
39–84（2012）
3）Donoho, D.L.: IEEE Trans Inf Theory, 52, 1289–1306 (2006)
4）Candes, E.J., Romberg, J. and Tao, T.: IEEE Trans Inf Theory, 52,
489–509 (2006)5）Leary, R., Saghi, Z., Midgley, P.A. et al.: Ultramicroscopy, 131, 70–91(2013)
6）Saghi, Z., Divitinia, G., Winter, B. et al.: Ultramicroscopy, 160,230–238 (2016)
7）Venkatakrishnan, S.V., Drummy, L.F., Jackson, M.A. et al.: IEEE
Trans Image Process, 22, 4532–4544 (2013)8）Monsegue, N., Jin, X., Echigo, T. et al.: Microscopy and Micro­
analysis, 13, 1362–1367 (2012)9）Batenburg, K.J., Bals, S., Sijbers, J. et al.: Ultramicroscopy, 109,730–740 (2009)10）Elder, Y.C. and Kutyniok, G.: Compressed sensing: theory andapplications, Cambridge University Press (2012)11）Elad, M.: Sparse and redundant representations: from theory toapplications in signal and image processing, Springer (2010)12）Rish, I. and Grabarnik, G.: Sparse modeling: theory, algorithms,and applications, CRC Press (2014)13）Elad, M., Milanfar, P. and Rubinstein, R.: Inverse Problems, 23,947–968 (2007)14）Loris, I. and Verhoeven, C.: Inverse Problems, 27, Article ID 125007(2011)
15）Chambolle, A. and Pock, T.: J. Math. Imaging Vis., 40, 120–145(2011)16）Rudin, L.I., Osher, S. and Fatemi, E.: Physica D, 60, 259–268 (1992)
17）Trzasko, J. and Manduca, A.: IEEE Trans Med Imaging, 28, 106–
121 (2009)18）Herman, G.T. and Davidi, R.: Inverse Problems, 24, Article ID
045011 (2008)
19）山嵜 深，高木景太，根本拓也，他：電子情報通信学会技術研
究報告 MI2014-95（2015）20）Guo, W., Qin, J. and Yin, W.: SIAM J Imaging Sci, 7, 1309–1334(2014)
21）Herman, G.T. and Meyer, L.B.: IEEE Trans Med Imaging, 12,600–609 (1993)22）Shirai, M., Tsumori, K., Kutsuwada, M. et al.: Nucl. Instrument.
Methods in Phys. Res. Section B, 267, 1787–1791 (2009)53【著作権者：公益社団法人 日本顕微鏡学会】