JPWO2018230615A1 - 画像処理装置、コンピュータプログラム及び画像補完方法 - Google Patents

Abstract

撮影時間を短縮することができる画像処理装置、コンピュータプログラム及び画像補完方法を提供する。画像処理装置は、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する取得部と、取得部で取得した画像データに基づいて、当該複数のｘｙ平面画像及びｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する推定部と、補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、補完後画像の各画素の画素値の和を算出する和算出部とを備え、推定部は、第１差分算出部で算出した差分及び和算出部で算出した和に基づいて補完後画像の画像データを推定する。

Description

本発明は、画像処理装置、コンピュータプログラム及び画像補完方法に関する。

生体イメージングは、医学、生物学の分野で注目されている新しい技術である。生体イメージングでは、生体組織の長時間イメージングによって、生命現象をリアルタイムで捉える。生体組織の立体的画像を得るため、例えば、６０μｍ〜１００μｍの間を数十枚程度の光学切片画像を取得する必要がある。

生体組織の立体的画像を得るためには、例えば、多光子励起顕微鏡、共焦点顕微鏡、蛍光顕微鏡などの計測機器が用いられる。非特許文献１には、２光子励起顕微鏡についての改善が開示されている。

Kaspar Podgorski, "Ultra-Bright and -Stable Red and Near-Infrared Squaraine Fluorophores for In Vivo Two-Photon Imaging", PLoS ONE 7(12): e51980, December 14, 2012

生体組織の微細形態を捕らえるためには、多数枚の光学切片画像を取得する必要があり、また高解像度の撮像を必要とする。このため、１回当たりの撮影時間（所要枚数の光学切片画像の撮影時間）が長くなるという問題があった。

本開示は斯かる事情に鑑みてなされたものであり、撮影時間を短縮することができる画像処理装置、コンピュータプログラム及び画像補完方法を提供することを目的とする。

本開示の実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第１差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。

本開示の実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、時系列に撮像された複数の２次元画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数の２次元画像及び時間方向に沿って位置する補完対象２次元画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第１差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。

本開示の実施の形態に係るコンピュータプログラムは、コンピュータに、実空間画像を補完させるためのコンピュータプログラムであって、コンピュータに、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する処理と、取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する処理と、前記補完後画像の隣接する画素の画素値の差分を算出する処理と、前記補完後画像の各画素の画素値の和を算出する処理と、算出した差分及び算出した和に基づいて前記補完後画像の画像データを推定する処理とを実行させる。

本開示の実施の形態に係る画像補完方法は、実空間画像を補完する画像補完方法であって、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得部が取得し、取得された画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定部が推定し、前記補完後画像の隣接する画素の画素値の差分を第１差分算出部が算出し、前記補完後画像の各画素の画素値の和を和算出部が算出し、算出された差分及び算出された和に基づいて前記補完後画像の画像データを前記推定部が推定する。

本開示によれば、撮影時間を短縮することができる。

本実施の形態の画像処理装置の構成の一例を示すブロック図である。 本実施の形態の画像処理装置による補完後画像の一例を示す模式図である。 取得した画像データのベクトルの一例を示す模式図である。 補完後画像の画像データのベクトルの一例を示す模式図である。 観測行列の一例を示す模式図である。 ｘ方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。 ｙ方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。 ｚ方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。 ｘ方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。 ｙ方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。 ｚ方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。 各行列の大きさの一例を示す説明図である。 本実施の形態の画像処理装置による高速画像取得の様子の一例を示す模式図である。 本実施の形態の画像処理装置による変換前後の画像の一例を示す説明図である。 本実施の形態の画像処理装置によるカーブレット変換の一例を示す模式図である。 本実施の形態の画像処理装置によるカーブレット変換後の画像の一例を示す説明図である。 本実施の形態の画像処理装置による補完後画像の第１例を示す説明図である。 本実施の形態の画像処理装置による補完後画像の第２例を示す説明図である。 本実施の形態の画像処理装置による補完後画像の精度評価の一例を示す説明図である。 本実施の形態の画像処理装置による補完後画像の高解像度化の様子の一例を示す模式図である。 本実施の形態の画像処理装置による撮影時間の一例を示す説明図である。 本実施の形態の画像処理装置による補完後画像の第３例を示す説明図である。 本実施の形態の画像処理装置によるパラメータ決定の処理の手順の一例を示すフローチャートである。 本実施の形態の画像処理装置による画像補完処理の手順の一例を示すフローチャートである。 本実施の形態の画像処理装置によるパラメータ決定の処理の手順の他の例を示すフローチャートである。 本実施の形態の画像処理装置による画像補完処理の手順の他の例を示すフローチャートである。 本実施の形態の画像処理装置の構成の他の例を示す説明図である。 第２実施形態の画像処理装置による補完後画像の一例を示す模式図である。

（第１実施形態）
以下、本発明をその実施の形態を示す図面に基づいて説明する。図１は本実施の形態の画像処理装置１００の構成の一例を示すブロック図である。画像処理装置１００は、計測機器２００に接続されている。計測機器２００は、生体組織の立体的画像などの生物画像を撮影することができる機器であり、例えば、多光子励起顕微鏡、共焦点顕微鏡、蛍光顕微鏡、ライトシート型顕微鏡、イオンコンダクタンス顕微鏡、光緩衝断層計などが含まれる。

画像処理装置１００は、装置全体を制御する制御部１０、画像データ取得部１１、パラメータ取得部１２、第１差分算出部１３、第２差分算出部１４、和算出部１５、推定部１６、補完画像生成部１７、変換部１８、記憶部１９、特定部２０及び関数算出部２１などを備える。

画像データ取得部１１は、取得部としての機能を有し、生体組織のｚ方向に沿って撮像された複数のｘｙ平面画像（取得画像、観測画像とも称する）の画像データ（観測値とも称する）を計測機器２００から取得する。

一つのｘｙ平面画像（光学切片画像とも称する）を、ｎ×ｍ画素の画像とし、ｚ方向に沿って撮影されたｘｙ平面画像の数（スタック数）をｚとすると、画像データは、ｎ×ｍ×ｚ行１列の行列（あるいは、ｎ×ｍ×ｚ個の成分を有するベクトル）で表すことができる。なお、例えば、ｎ×ｍ×ｚは、ｎ＊ｍ＊ｚの如く表記することもできる。取得した画像データをベクトルＹで表す。画像データは、各画素の画素値を含む、画素値は、例えば、０〜２５５で表される輝度値とすることができる。すなわち、画像データは、例えば、各画素の輝度値のデータであるといえる。

推定部１６は、画像データ取得部１１で取得した画像データに基づいて、当該複数の取得画像（ｘｙ平面画像）及びｚ方向に沿って位置する補完対象画像（補完対象ｘｙ平面画像）を含む補完後画像の画像データを推定する。

補完対象画像は、撮影された複数の取得画像によって補完される画像であって、実際には撮影されていない画像である。補完対象画像の数（スタック数）をΔｚとすると、補完後画像の数（スタック数）ｚ′は、ｚ′＝ｚ＋Δｚとなる。補完後画像の画像データは、ｎ×ｍ×ｚ′行１列の行列（あるいは、ｎ×ｍ×ｚ′個の成分を有するベクトル）で表すことができる。補完後画像の画像データをベクトルＸで表す。

図２は本実施の形態の画像処理装置１００による補完後画像の一例を示す模式図である。図中左側の図は取得画像（観測画像）を示し、右側の図は補完後画像を示す。図２に示すように、取得画像Ｇ１、Ｇ２、Ｇ３は、生体組織のｘｙ平面上の画像であり、ｚ方向に沿って３枚撮影されており、スタック数＝３となる。なお、図２では、便宜上、取得画像のスタック数を３としているが、スタック数は３に限定されない。

補完後画像は、取得画像Ｇ１、Ｇ２、Ｇ３及び補完対象画像Ｃ１、Ｃ２、Ｃ３を含み、生体組織のｘｙ平面上の画像であり、ｚ方向に沿って６枚あり、スタック数＝６となる。補完対象画像Ｃ１、Ｃ２、Ｃ３は、実際に撮影された画像ではなく、本実施の形態の画像処理装置１００によって補完される画像である。補完対象画像Ｃ１は、取得画像Ｇ１とＧ２との間にあり、補完対象画像Ｃ２は、取得画像Ｇ２とＧ３との間にあり、補完対象画像Ｃ３は、取得画像Ｇ３よりも深い位置にある。なお、画像の枚数、補完対象画像の位置などは図２の例に限定されるものではない。

図３は取得した画像データのベクトルＹの一例を示す模式図である。便宜上、取得画像が３×３画素で構成されるとする。取得画像Ｇ１（ｚ＝１）の画素は９個存在し、画素番号（ピクセル番号）を１〜９とする。同様に、取得画像Ｇ２（ｚ＝２）の画素番号（ピクセル番号）を１０〜１８とし、取得画像Ｇ３（ｚ＝３）の画素番号（ピクセル番号）を１９〜２７とする。画像データのベクトルＹは、１〜２７の２７次元ベクトルで表される。ベクトルＹの各成分は、それぞれの画素番号における観測値となる。

図４は補完後画像の画像データのベクトルＸの一例を示す模式図である。補完後画像は、取得画像Ｇ１〜Ｇ３及び補完対象画像Ｃ１、Ｃ２、Ｃ３で構成される。図４に示すように、補完後画像のスタックｚ＝１、ｚ＝３、ｚ＝５の画像が、取得画像Ｇ１（ｚ＝１）、取得画像Ｇ２（ｚ＝２）、取得画像Ｇ３（ｚ＝３）に対応している。この対応関係を決定するのが、後述の観測行列Ａである。補完後画像は、６枚の画像で構成されるので、ベクトルＸは、１〜５４の５４次元ベクトルで表される。ベクトルＹの各成分は、それぞれの画素番号における推定値となる。

図５は観測行列Ａの一例を示す模式図である。観測行列Ａは、取得した画像データのベクトルＹの各成分が、補完後画像の画像データのベクトルＸのいずれの成分に対応するかを決定するための行列である。ベクトルＹが２７次元であり、ベクトルＸが５４次元である場合、観測行列Ａは、２７行５４列の行列となる。ベクトルＹの成分とベクトルＸの成分とが対応している場合、観測行列Ａの成分は１となり、対応していない場合、観測行列Ａの成分は０となる。

例えば、図４に例示したように、ベクトルＹのｙ１〜ｙ９が、それぞれベクトルＸのｘ１〜ｘ９に対応している場合、観測行列Ａの（１、１）、（２、２）、…（９、９）成分は１となる。また、ベクトルＹのｙ１０〜ｙ１８が、それぞれベクトルＸのｘ１９〜ｘ２７に対応している場合、観測行列Ａの（１０、１９）、（１１、２０）、…（１８、２７）成分は１となる。また、ベクトルＹのｙ１９〜ｙ２７が、それぞれベクトルＸのｘ３７〜ｘ４５に対応している場合、観測行列Ａの（１９、３７）、（２９、３８）、…（２７、４５）成分は１となる。

より具体的には、推定部１６は、関数算出部２１が算出する所定関数としてのコスト関数を用い、ベクトルＸの各成分（推定値）の初期値（例えば、０）を当該コスト関数に入力してコスト関数が最小化となる推定値を更新する。推定部１６は、更新した推定値を使って再度コスト関数を最小化する処理を、コスト関数が最小となるまで繰り返す。これにより推定値を求めることができる。コスト関数が最小になったか否かは、繰り返し算出したコスト関数が、ある値に収束していくかで判定することができ、ある値に収束した場合、最小になったと判定することができる。

第１実施として、関数算出部２１は、式（１）を用いてコスト関数を算出することができる。ａｒｇｍｉｎは、コスト関数が最小となるベクトルＸを求めるというような意味である。λ１〜λ４は、パラメータである。λ１〜λ４は、生物画像がもつ特徴によって適宜設定することができる。

第２差分算出部１４は、式（１）の第１項、||Ｙ-Ａ・Ｘ||₂ （［Ｙ-Ａ・Ｘ］のＬ２ノルムと称する）を算出する。すなわち、第２差分算出部１４は、取得した複数の取得画像（観測画像）の画像データ（ベクトルＹ）と、補完後画像の画像データのうち当該複数の取得画像に対応する画像データ（Ａ・Ｘ）との差分を算出する。

第２差分算出部１４で算出した差分を小さくすることにより、取得した画像データと、補完後画像（具体的には、補完後画像のうちの取得画像）の画像データとの差を小さくしつつ補完後画像の画像データを推定することができる。第１項を最小化することは、取得データ（観測データ）と推定データ（再構成データ）とを一致させることに繋がり、補完後画像の元の画像に対する精度を向上させることができる。

第１差分算出部１３は、式（１）の第２項、第３項及び第４項を算出する。すなわち、第１差分算出部１３は、補完後画像の隣接する画素の画素値の差分を算出する。隣接する画素は、例えば、ｘ方向（第２項に対応）、ｙ方向（第３項に対応）及びｚ方向（第４項に対応）に沿った画素とすることができる。

第２項のＴｖｘは、ｘ方向に沿って隣接する画素の画素値の差分を算出するための差分行列であり、第３項のＴｖｙは、ｙ方向に沿って隣接する画素の画素値の差分を算出するための差分行列であり、第４項のＴｖｚは、ｚ方向に沿って隣接する画素の画素値の差分を算出するための差分行列である。

生体組織を撮影した生物画像は、隣接する画素の画素値（例えば、輝度値）の差分のスパース性（画素値の差分がゼロとなる特性）を有するという特徴がある。隣接する画素の画素値の差分のスパース性に基づき、第１差分算出部１３で算出する、ｘ方向、ｙ方向及びｚ方向に沿った隣接する画素の画素値の差分を小さくするという手法を用いることができる。

ｘ方向及びｙ方向に沿って隣接する画素の画素値の差分を算出し、算出した差分が小さくなるように補完後画像の画像データを推定するので、生体組織に適した生物画像を補完することができる。また、ｚ方向に沿って隣接する画素の画素値の差分を算出し、算出した差分が小さくなるように補完後画像の画像データを推定するので、３次元構造の生体組織に適した生物画像を補完することができる。

次に差分行列について具体的に説明する。

図６はｘ方向に沿って隣接する画素の画素値の差分を算出する差分行列Ｔｖｘの一例を示す模式図である。ｚ＝１の補完後画像を考える。ｘ方向の差分は、ピクセル１と２との差、ピクセル２と３との差、ピクセル４と５との差、ピクセル５と６との差、ピクセル７と８との差、ピクセル８と９との差となる。

差分行列Ｔｖｘの（１、１）成分を１とし、（１、２）成分を-１とすることにより、ピクセル１と２との差を算出することができる。他のピクセルの差も同様である。一番右端のピクセル３、６、９は差をとる相手がいないので、差分行列Ｔｖｘの３、６、９行目の成分はない。また、他のスタック（ｚ＝１以外）についても同様にすることで、差分行列Ｔｖｘを決定することができる。

図７はｙ方向に沿って隣接する画素の画素値の差分を算出する差分行列Ｔｖｙの一例を示す模式図である。ｚ＝１の補完後画像を考える。ｙ方向の差分は、ピクセル１と４との差、ピクセル２と５との差、ピクセル３と６との差、ピクセル４と７との差、ピクセル５と８との差、ピクセル６と９との差となる。

差分行列Ｔｖｙの（１、１）成分を１とし、（１、４）成分を−１とすることにより、ピクセル１と４との差を算出することができる。他のピクセルの差も同様である。一番下のピクセル７、８、９は差をとる相手がいないので、差分行列Ｔｖｘの７〜９行目の成分はない。また、他のスタック（ｚ＝１以外）についても同様にすることで、差分行列Ｔｖｙを決定することができる。

図８はｚ方向に沿って隣接する画素の画素値の差分を算出する差分行列Ｔｖｚの一例を示す模式図である。ｚ＝１〜６の補完後画像を考える。ｚ方向の差分は、ｚ＝１の画像のピクセル１とｚ＝２の画像のピクセル１０との差、ｚ＝１の画像のピクセル２とｚ＝２の画像のピクセル１１との差、…、ｚ＝１の画像のピクセル９とｚ＝２の画像のピクセル１８との差となる。

差分行列Ｔｖｚの（１、１）成分を１とし、（１、１０）成分を−１とすることにより、ｚ＝１の画像のピクセル１とｚ＝２の画像のピクセル１０との差を算出することができる。他のピクセルの差も同様である。一番下のｚ＝６の画像のピクセル４６〜５４は差をとる相手がいないので、差分行列Ｔｖｘの対応する成分はない。また、他のスタック同士についても同様にすることで、差分行列Ｔｖｚを決定することができる。

一般的に、ベクトルＸを式（２）で表すと、（ｘ１、ｘ２、…ｘｎ）は、ベクトルＸの成分である。ベクトルＸのＬ１ノルムは、式（３）で表される。また、ベクトルＸのＬ２ノルムは、式（４）で表される。

第１差分算出部１３は、式（１）の第２項及び第３項の如く、補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分のＬ２ノルムを算出する。また、第１差分算出部１３は、式（１）の第４項の如く、補完後画像のｚ方向に沿って隣接する画素の画素値の差分のＬ１ノルムを算出する。

推定部１６は、第１差分算出部１３で算出したＬ２ノルムが最小になるように補完後画像の画像データを推定する。取得した画像は、ｘｙ平面画像であるので、ｘ方向及びｙ方向に沿って隣接する画素の画素値の情報はもともと存在している。そこで、ｘ方向及びｙ方向については、Ｌ２ノルムを用いることにより、存在する情報を用いて滑らかな生体画像を得ることができる。

また、推定部１６は、第１差分算出部１３で算出したＬ１ノルムが最小になるように補完後画像の画像データを推定する。ｚ方向に沿って隣接する画素の画素値の情報は、補完対象のスタックでは存在しない。そこで、ｚ方向については、Ｌ１ノルムを用いることにより、情報が存在しない場合でも、スパース解（複数の数値の組で構成される解の中にゼロが含まれる）が得られるようにして、より正確な補完を可能にして、滑らかな生体画像を得ることができる。

和算出部１５は、式（１）の第５項、||Ｘ||₁ （［Ｘ］のＬ１ノルムと称する）を算出する。すなわち、和算出部１５は、補完後画像の各画素の画素値の和を算出する。生体組織を撮影した生物画像は、観測対象以外の部分のスパース性（画素値がゼロとなる特性）を有するという特徴がある。観測対象以外の部分のスパース性に基づき、和算出部１５で算出した和を小さくするという手法を用いることができる。

例えば、遺伝学的に蛍光タンパクを特定の細胞で発現させるシステム（計測機器２００を含む）を使用する場合、当該システムは厳密に制御されているため、理論的にはその細胞以外の部位は光ることはない。すなわち、光っている部分以外はノイズ又は自家蛍光であると考えられるので、光っている部分以外の部分は輝度値０にし、真の信号（光っている部分からの信号）はそのままにするという選択をさせることができる。

推定部１６は、第１差分算出部１３で算出した差分及び和算出部１５で算出した和に基づいて補完後画像の画像データを推定する。生体組織を撮影した生物画像には、隣接する画素の画素値（例えば、輝度値）の差分のスパース性（画素値の差分がゼロとなる特性）、及び観測対象以外の部分のスパース性（画素値がゼロとなる特性）を有するという特徴がある。このようなスパース性に注目することによって、取得した少ない画像データ（例えば、スタック数がｚの画像）から多くの情報（例えば、スタック数ｚ′の画像）を抽出することができる。

推定値が得られることによって、スタック数ｚの画像（取得した画像）から、スタック数ｚ′の補完後画像を得ることができる。これにより、ｚ枚（ｚ＜ｚ′）の光学切片画像を取得するだけでｚ′枚の光学切片画像から得られる情報を抽出することができるので、従来よりも短い撮影時間で生体組織の微細形態を捕らえることができる。

差分行列の一例を図６、図７及び図８を用いて説明したが、差分行列は図６から図８の例に限定されない。以下では、差分行列の他の例について説明する。

図９はｘ方向に沿って隣接する画素の画素値の差分を算出する差分行列Ｔｖｘの他の例を示す模式図である。なお、便宜上、差分行列を３行×３列の行列として説明する。ｚ＝１の補完後画像を考える。図９に示すように、ｘ方向の差分は、ピクセル２の重み付け（図９では、重み付け係数が２）された画素値と、ピクセル１及びピクセル３それぞれの重み付け（図９では、重み付け係数が１）された画素値との差となる。同様に、ｘ方向の差分は、ピクセル５の重み付け（図９では、重み付け係数が２）された画素値と、ピクセル４及びピクセル６それぞれの重み付け（図９では、重み付け係数が１）された画素値との差となる。

差分行列Ｔｖｘの（２、２）成分を２とし、（２、１）成分及び（２、３）成分をそれぞれ-１とすることにより、補完後画像のピクセル２の重み付けされた画素値と、ピクセル２を間にしてｘ方向に沿って隣接するピクセル１及びピクセル３の重み付けされた画素値との差分を算出することができる。他のピクセルも同様である。一番右端のピクセル３、６、９、及び一番左端のピクセル１、４、７は差をとる相手がいないので、ｘ方向に隣接する画素の画素値との差分を計算している。また、他のスタック（ｚ＝１以外）についても同様にすることで、差分行列Ｔｖｘを決定することができる。

上述のように、第１差分算出部１３は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてｘ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出することができる。

すなわち、補完後画像の注目画素の画素値をＰ（ｉ＋１）とし、重み付け係数をα１とする。注目画素を間にしてｘ方向に沿って隣接する一方の画素の画素値をＰ（ｉ）とし、重み付け係数をα０とする。注目画素を間にしてｘ方向に沿って隣接する他方の画素の画素値をＰ（ｉ＋２）とし、重み付け係数をα２とする。

第１差分算出部１３は、α１×Ｐ（ｉ＋１）−{α０×Ｐ（ｉ）＋α２×Ｐ（ｉ＋２）}という式を用いて差分を算出する。ここで、α１−（α０＋α２）＝０とすることができ、例えば、α１＝２、α０＝α２＝１とすると、第１差分算出部１３は、２×Ｐ（ｉ＋１）−{Ｐ（ｉ）＋Ｐ（ｉ＋２）}という式により差分を算出することができる。図９の例では、α０＝α２＝１であり、α１＝２である。

これにより、補完後画像の任意の画素の画素値と、ｘ向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、ｘ方向に沿った偏りを抑制して、３次元構造の生体組織に適した生物画像を補完することができる。

図１０はｙ方向に沿って隣接する画素の画素値の差分を算出する差分行列Ｔｖｙの他の例を示す模式図である。ｚ＝１の補完後画像を考える。図１０に示すように、ｙ方向の差分は、ピクセル４の重み付け（図１０では、重み付け係数が２）された画素値と、ピクセル１及びピクセル７それぞれの重み付け（図１０では、重み付け係数が１）された画素値との差となる。同様に、ｙ方向の差分は、ピクセル５の重み付け（図１０では、重み付け係数が２）された画素値と、ピクセル２及びピクセル８それぞれの重み付け（図１０では、重み付け係数が１）された画素値との差となる。

差分行列Ｔｖｙの（４、４）成分を２とし、（４、１）成分及び（４、７）成分をそれぞれ-１とすることにより、補完後画像のピクセル４の重み付けされた画素値と、ピクセル４を間にしてｙ方向に沿って隣接するピクセル１及びピクセル７の重み付けされた画素値との差分を算出することができる。他のピクセルも同様である。一番上端のピクセル１、２、３、及び一番下端のピクセル７、８、９は差をとる相手がいないので、ｙ方向に隣接する画素の画素値との差分を計算している。また、他のスタック（ｚ＝１以外）についても同様にすることで、差分行列Ｔｖｙを決定することができる。

上述のように、第１差分算出部１３は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてｙ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出することができる。

すなわち、補完後画像の注目画素の画素値をＰ（ｉ＋１）とし、重み付け係数をα１とする。注目画素を間にしてｙ方向に沿って隣接する一方の画素の画素値をＰ（ｉ）とし、重み付け係数をα０とする。注目画素を間にしてｙ方向に沿って隣接する他方の画素の画素値をＰ（ｉ＋２）とし、重み付け係数をα２とする。

第１差分算出部１３は、α１×Ｐ（ｉ＋１）−{α０×Ｐ（ｉ）＋α２×Ｐ（ｉ＋２）}という式を用いて差分を算出する。ここで、α１−（α０＋α２）＝０とすることができ、例えば、α１＝２、α０＝α２＝１とすると、第１差分算出部１３は、２×Ｐ（ｉ＋１）−{Ｐ（ｉ）＋Ｐ（ｉ＋２）}という式により差分を算出することができる。図１０の例では、α０＝α２＝１であり、α１＝２である。

これにより、補完後画像の任意の画素の画素値と、ｙ向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、ｙ方向に沿った偏りを抑制して、３次元構造の生体組織に適した生物画像を補完することができる。

図１１はｚ方向に沿って隣接する画素の画素値の差分を算出する差分行列Ｔｖｚの他の例を示す模式図である。ｚ＝１〜６の補完後画像を考える。ｚ方向の差分は、ピクセル１０の重み付け（図１１では、重み付け係数が２）された画素値と、ピクセル１及びピクセル１９それぞれの重み付け（図１１では、重み付け係数が１）された画素値との差となる。

差分行列Ｔｖｚの（１０、１０）成分を２とし、（１０、１）成分及び（１０、１９）成分をそれぞれ-１とすることにより、補完後画像のピクセル１０の重み付けされた画素値と、ピクセル１０を間にしてｚ方向に沿って隣接するピクセル１及びピクセル１９の重み付けされた画素値との差分を算出することができる。他のピクセルも同様である。一番上のｚ＝１の画像のピクセル１〜９、及び一番下のｚ＝６の画像のピクセル４６〜５４は差をとる相手がいないので、ｚ方向に隣接する画素の画素値との差分を計算している。また、他のスタック同士についても同様にすることで、差分行列Ｔｖｚを決定することができる。

上述のように、第１差分算出部１３は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてｚ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出することができる。

すなわち、補完後画像のｚ方向に沿ったスタックｚ′を、ｚ′＝ｋ、ｋ＋１、ｋ＋２とする。スタックｚ′＝ｋ＋１の補完後画像上の注目画素の画素値をＰ（ｋ＋１）とし、重みけ付け係数をα１とする。同様に、スタックｚ′＝ｋの補完後画像上の当該注目画素に対応する位置の画素の画素値をＰ（ｋ）とし、重みけ付け係数をα０とし、スタックｚ′＝ｋ＋２の補完後画像上の当該注目画素に対応する位置の画素の画素値をＰ（ｋ＋２）とし、重みけ付け係数をα２とする。

第１差分算出部１３は、α１×Ｐ（ｋ＋１）−{α０×Ｐ（ｋ）＋α２×Ｐ（ｋ＋２）}という式を用いて差分を算出する。ここで、α１−（α０＋α２）＝０とすることができ、例えば、α１＝２、α０＝α２＝１とすると、第１差分算出部１３は、２×Ｐ（ｋ＋１）−{Ｐ（ｋ）＋Ｐ（ｋ＋２）}という式により差分を算出することができる。図１１の例では、α０＝α２＝１であり、α１＝２である。

これにより、補完後画像の任意の画素の画素値と、ｚ方向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、ｚ方向に沿った偏りを抑制して、３次元構造の生体組織に適した生物画像を補完することができる。

図１２は各行列の大きさの一例を示す説明図である。取得した画像のスタック数をｚとし、補完後画像のスタック数をｚ′とする。画像の大きさをｎ×ｍ画素とする。ベクトルＹ（観測値）は、ｎ＊ｍ＊ｚ次元、すなわちｎ＊ｍ＊ｚ行１列の行列に相当する。観測行列Ａは、ｎ＊ｍ＊ｚ行ｎ＊ｍ＊ｚ′列の行列となる。ベクトルＸ（推定値）は、ｎ＊ｍ＊ｚ′次元、すなわちｎ＊ｍ＊ｚ′行１列の行列に相当する。差分行列Ｔｖｘは、（ｎ-１）＊ｍ＊ｚ′行ｎ＊ｍ＊ｚ′列の行列となり、差分行列Ｔｖｙは、ｎ＊（ｍ-１）＊ｚ′行ｎ＊ｍ＊ｚ′列の行列となり、差分行列Ｔｖｚは、ｎ＊ｍ＊（ｚ-１）′行ｎ＊ｍ＊ｚ′列の行列となる。

上述のように、本実施の形態の画像処理装置１００によれば、撮影時間を短縮することができ、時間方向の分解能を向上させることができるので、高速で画像を取得することができる。

図１３は本実施の形態の画像処理装置１００による高速画像取得の様子の一例を示す模式図である。左側の図は比較例の場合を示し、右側の図は本実施の形態の場合を示す。比較例は、計測機器２００だけを用いて３次元の生体組織を撮影した場合である。図中、曲線で囲まれる領域は、生体組織（例えば、筋細胞など）の撮影対象を示す。また、立方体状の領域は、１回の撮影で撮影される３次元領域を示す。

図１３に示すように、比較例の場合、１回の撮影時間（例えば、１５分など）の間では、当然１箇所しか撮影することができない。一方、本実施の形態では、画像を高速で取得することができるので、当該撮影時間（例えば、１５分）の間に、例えば、４箇所の部分の画像を取得することができる。すなわち、従来であれば、生体組織の一点しか観察することができないところ、本実施の形態によれば、時間方向の分解能が向上するので、生体組織の多点を観察することができる。このように、撮影時間の短縮によって、多点タイムラプス（複数視野の観察）が一層容易になり、例えば、ライブイメージングに利用するマウス匹数の削減が可能となる。

また、撮影時間の短縮を実現し、蛍光強度を維持しながら、必要な情報を過不足なく取得することができる。従来にあっては、顕微鏡のスキャンスピードを増加させて解像度の低い画像を得ることで妥協してきたものが、解像度を犠牲にすることなく撮影時間を短くすることができる。すなわち、生体組織の３Ｄ解析に十分使用できる程度の詳細画像を従来よりも短時間で撮影することができる。

また、撮影時間を短縮することができるので、長時間のレーザ照射が不要となり、計測機器２００で使用するレーザ照射による蛍光の退色、レーザによる細胞及び蛍光タンパク等へのダメージを少なくすることができる。

また、ディープラーニングやＧＡＮ（Generative Adversarial Network）のような情報処理技術に比べて、本実施の形態による処理は負荷が軽く、本実施の形態による処理を実行させるためのコンピュータがハイスペックである必要はなく、本実施の形態の適用（普及）は容易である。

また、３次元構造の生体組織を撮影して得られた生体画像は、特有の性質（人物、風景、建物など撮影した一般画像にはない性質）を有する。特有の性質は、例えば、３次元方向に沿って隣接する画素の画素値の差分のスパース性、観測対象以外の部分のスパース性などである。本実施の形態では、これらの特有の性質（スパース性）を事前情報としたベイズ推定を用いることによって、元画像を再現することができる正確な画像補完が可能となる。

また、本実施の形態では、例えば、波数空間画像などのように、一つの波の情報の中に非常に多くの情報が入っている結果、それらを組み合わせて詳細な画像を表現するシステムとは異なり、実空間画像の一つのピクセルに画素値（例えば、輝度値）という極めて少ない情報だけを頼りに、上述の事前情報と組み合わせることによって、正確な画像（実際に取得していない画像）を補完することができる。

次に、画像を補完する際のノイズ除去について説明する。

変換部１８は、推定部１６が推定した画像データによって得られた補完対象画像に対してカーブレット変換を施す。カーブレット変換は、例えば、生体組織（例えば、細胞など）がもつ「曲線の集合」という画像構造に注目し、曲線の集まりとして捉えられる画像を抽出する手法である。

予め種々の曲線（生体組織を撮影して得られる曲線部分）を辞書データとして収集しておく。補完対象画像を所定の大きさの画素領域（例えば、３０×３０画素で構成される画素領域）毎に走査して、辞書データに含まれる曲線に適合するか否かの判定を行い、曲線に適合した部分だけを抽出する。これにより、補完対象画像に含まれるノイズ（例えば、生体組織に相当しない、格子状など）を除去して補完精度を向上させることができる。

図１４は本実施の形態の画像処理装置１００による変換前後の画像の一例を示す説明図である。左側の図はカーブレット変換前の画像であり、右側の図は変換後の画像である。画像の中央部に存在す明るい部分及び当該部分から伸びている細い足状の部分は、辞書データとして存在する曲線の集合に相当するので、そのまま残されている。一方、変換前の画像の左側に存在する格子点状の部分は、辞書データに存在しないので、ノイズとして除去されている。なお、カーブレット変換は、一例であって、これに限定されない。例えば、曲線に適合した部分だけを抽出して、ノイズを除去することができるような、他の処理を用いることができる。

上述のように、本実施の形態によれば、恣意的でない、ノイズの除去が可能となる。

ｘｙｚ方向の補完後画像に対するカーブレット変換は、ｘｙ平面画像に対するカーブレット変換、ｘｚ平面画像に対するカーブレット変換、ｙｚ平面画像に対するカーブレット変換のいずれかを用いることができる。

図１５は本実施の形態の画像処理装置１００によるカーブレット変換の一例を示す模式図である。図中、符号Ｇ１〜Ｇ３、Ｃ１〜Ｃ３は、図２の例と同様である。すなわち、画像Ｇ１〜Ｇ３は取得画像であり、画像Ｃ１〜Ｃ３は補完対象画像である。上段の図は、ベクトルＹの画像データを、再構成して、ｘｙ平面画像をスタックｚ＝１からｚ＝６まで並べたものである。ｘｙ平面画像に対するカーブレット変換は、スタックｚ＝１からｚ＝６までの６枚のｘｙ平面画像それぞれに対して施される。なお、便宜上、画像の枚数を６としているが、実際は、さらに多くの枚数の画像を用いる。

下段の図は、ベクトルＹの画像データを、再構成して、ｘｚ平面画像をスタックｙ＝１からｙ＝６まで並べたものである。ｘｚ平面画像に対するカーブレット変換は、スタックｙ＝１からｙ＝６までの６枚のｘｙ平面画像それぞれに対して施される。なお、ｙｚ平面画像はｘｚ平面画像と同様であるので、省略している。

図１６は本実施の形態の画像処理装置１００によるカーブレット変換後の画像の一例を示す説明図である。左側の図は、ｘｙ平面画像に対してカーブレット変換を行った後の画像を示し、右側の図は、ｘｚ平面画像に対してカーブレット変換を行った後の画像を示す。なお、ｙｚ平面画像の場合はｘｚ平面画像の場合と同様であるので、省略している。

図中、中央部に左下から右上の方向に向かって白味を帯びて見える領域が、対象の生体組織である。図に示すように、ｘｙ平面画像に対してカーブレット変換を行った後の画像に比べて、ｘｚ平面画像に対してカーブレット変換を行った後の画像の方が、画像の精細度を良いことが分かる。

図１５に示したように、ｘｙ平面画像では、６枚の画像のうち、３枚は取得した画像であるので、当該画像では、元々画像データが滑らかであると言える。これに対して、ｘｚ平面画像では、取得した画像の一部と補完対象画像の一部とが１枚の画像の中に交互に現れるため、画像データはあまり滑らかでなく、カーブレット変換が比較的強く作用すると考えられる。このため、ｘｚ平面画像又はｙｚ平面画像に対してカーブレット変換を行った後の画像の方が、画像の精細度が良くなると考えられる。なお、ｘｙ平面画像に対してカーブレット変換を行ってよい。

次に、本実施の形態の画像処理装置１００による画像の補完結果について説明する。

図１７は本実施の形態の画像処理装置１００による補完後画像の第１例を示す説明図である。図中、上段の図は元画像（実際に撮影した画像）をｚ方向に沿って３枚並べたものである。中段の図は、比較例として一般的な画像の平滑化を使用して補完した場合を示し、３つの画像のうち、中央の画像が補完した画像である。下段の図は、本実施の形態のＬ１ノルム及びＬ２ノルムの組合せを使用して補完した場合を示し、３つの画像のうち、中央の画像が補完した画像である。右端の図は、補完した画像及び元画像の拡大像である。

図１７に示すように、比較例の場合、元の画像に比べて全体に輝度が高くなっている。また、拡大像の中央にある細胞（やや明るい部分）から右下に向かって伸びている細い足状の部分については、元の画像よりも下方に向かって長くなっている。一方、本実施の形態の場合には、足状の部分を正確に補完することができている。

図１８は本実施の形態の画像処理装置１００による補完後画像の第２例を示す説明図である。４つ並んだ各図において、上側の図はｚ方向の横断面画像を示し、下側の図はｘｙ平面画像の一つを示す。図中、左から２番目の画像は元画像（６０スタック）を示し、左端の図は、ｚ−スタック半減画像、すなわち、６０スタックに半分である３０スタックだけを実際に撮影した画像である。右から２番目の図は３０スタックの画像から比較例（一般的な画像の平滑化を使用）による補完画像を示し、右端の図は本実施の形態による補完画像を示す。

ｚ方向の横断面画像を見ると、比較例の場合には、元画像に比べてノイズ（点状の明るい箇所）が多いことが分かる。また、本実施の形態では、元画像を正確に補完していることが分かる。また、スタック数を半減した場合、左端の図が示すように、ｚ方向の横断面画像を見ると、元画像では観測することができる細胞（白い部分）の大部分を観測することができないことが分かる。

このように、本実施の形態によれば、横方向から見ても生体組織の輪郭をはっきり観測することができ、生体組織（例えば、細胞形態）を明瞭に把握することができる。

図１９は本実施の形態の画像処理装置１００による補完後画像の精度評価の一例を示す説明図である。ＮＯ．１からＮＯ．４の４種類の画像の評価結果を示す。評価指標は、ＭＥＳ（Mean Square Error）及びＳ／Ｎ Ratio（Signal to Noise Ratio）を用いている。ＭＥＳは、元画像との画素値の違いを表し、０に近いほど元画像に近いことを示す。Ｓ／Ｎ Ratioは、元画像と比較して、どの程度シグナルが増加し、ノイズが減少したかを表し、大きな値ほど良好な画像であることを示す。

図１９に示すように、ＮＯ．１からＮＯ．４の各画像において、比較例（一般的な画像の平滑化を使用）に比べて、本実施の形態では、ＭＳＥが小さく、Ｓ／Ｎ Ratioが大きいことが分かる。

次に、補完後画像を高解像度化する方法について説明する。

図２０は本実施の形態の画像処理装置１００による補完後画像の高解像度化の様子の一例を示す模式図である。記憶部１９は、複数の高解像度画像と当該複数の高解像度画像を疑似的に劣化させた複数の低解像度画像とを対応付けた辞書データを記憶する。複数の高解像度画像としては、予め種々の生体組織を撮影した画像を用いることができる。高解像度画像は、例えば、６０×６０画素の画像とし、低解像度画像は、３０×３０画素の画像とすることができるが、これに限定されない。なお、本実施の形態では、画像処理装置１００が辞書データを備える構成であるが、画像処理装置１００から外部の辞書データベースに接続する構成であってもよい。

画像データ取得部１１で取得した画像データに係る取得画像の解像度を、例えば、５１２×５１２とする（低解像度とする）。推定部１６によって推定された推定値に基づく補完後画像は、５１２×５１２（低解像度補完画像）となる。

特定部２０は、記憶部１９に記憶した複数の低解像度画像の中から、推定部１６が推定した画像データによって得られた補完対象ｘｙ平面画像を構成する複数の低解像度画像を特定する。例えば、低解像度補完画像を、所定の大きさの画素領域（例えば、３０×３０画素で構成される画素領域）毎に走査して、一致又は類似する辞書データ内の低解像度画像を特定する。

補完画像生成部１７は、生成部としての機能を有し、特定部２０で特定した複数の低解像度画像に対応する複数の高解像度画像に基づいて高解像度補完対象画像を生成する。

これにより、例えば、５１２×５１２で撮影した低解像度補完画像を、１０２４×１０２４の高解像度補完画像にすることができる。

上述の辞書データは、生体組織を撮影した大量の画像を教師データとして、コンピュータ（学習モジュール）に学習させることによって作成することができる。これによって、辞書データには、細胞などの生体組織に特有な形態（ヒゲ状の部分、分岐する部分、細胞の塊部分など）が表現された画像が記録される。なお、生体組織に特有のスパース性によって、蓄積される画像データは、比較的単純化することができ、辞書データの規模は小さくて済む。

図２１は本実施の形態の画像処理装置１００による撮影時間の一例を示す説明図である。右側の比較例は、計測機器２００だけを用いて３次元の生体組織を撮影した場合のデータを示す。すなわち、１０２４×１０２４の解像度で６０スタック（６０枚）の光学切片画像を撮影したときのデータ量は、１２７．９ＭＢであり、撮影時間は６分１０秒であった。

次に、本実施の形態により、１０２４×１０２４の解像度で３０スタック（３０枚）の光学切片画像を撮影したときのデータ量は、６５．０ＭＢであり、撮影時間は２分４０秒であった。この場合、取得画像は３０枚であり、補完対象画像は３０枚、補完後画像は、比較例と同様に６０枚となる。すなわち、比較例と同じ６０枚の画像を得るのに、２分４０秒しか要しないことが分かり、撮影時間が大幅に短縮されている。

また、５１２×５１２の解像度で３０スタック（３０枚）の光学切片画像を撮影したときのデータ量は、１６．３ＭＢであり、撮影時間は１分２０秒であった。この場合、取得画像は３０枚の低解像度画像であり、高解像度化された補完対象画像は３０枚、高解像度化された補完後画像は、比較例と同様に６０枚となる。すなわち、比較例と同じ６０枚の画像を得るのに、１分２０秒しか要しないことが分かり、撮影時間がさらに短縮されている。また、データ量は、比較例が１２７．９ＭＢであるのに対し、本実施の形態では、１６．３ＭＢであり、１撮影当たりのデータ量を削減することができる。

また、補完後画像の高解像度化（超解像）の他の例として、以下の手法を用いることもできる。すなわち、図２に示したように、ｘｙ平面画像をｚ方向に対して補完して得られた補完後画像（ｘｙｚ方向の立体状）を回転し、当該補完画像のｙｚ平面画像をｘ方向に対して補完して補完後画像を得る。得られた補完後画像をさらに回転し、当該補完画像のｘｚ平面画像をｙ方向に対して補完することによって、疑似的に３次元（３Ｄ）の解像度を上げることができる。なお、補完方向の順番は、ｚ方向、ｘ方向、ｙ方向の順番に限定されない。

次に、コスト関数の他の例について説明する。前述の例（実施例１と称する）では、コスト関数として、式（１）を用いた。他の例（実施例２と称する）では、コスト関数として式（５）を用いる。

実施例２では、第１差分算出部１３は、式（５）の第２項の如く、補完後画像のｚ方向沿って隣接する画素の画素値の差分のＬ１ノルムを算出する。すなわち、補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分については算出せず、コスト関数に入れていない。実施例２においても、実施例１と同様に、撮影時間を短縮することができ、時間方向の分解能を向上させることができる。なお、実施例１の方が、滑らかな画像を得る可能性が高くなる。

また、コスト関数の他の例として、実施例３では、第１差分算出部１３は、式（６）の第２項及び第３項の如く、補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分のＬ１ノルムを算出する。

図２２は本実施の形態の画像処理装置１００による補完後画像の第３例を示す説明図である。図２２の例は、血管画像の一例を示す。元画像は、ｚ方向に７枚（スタック数＝７）撮影してある。元画像の各スタックでは、左上から右下に向かって細長い白い部分を観測することができる。これが血管を示している。取得画像は、元画像のうち、左側から右側に向かって奇数番目の画像が取得されているとする。すなわち、元画像のうち、左側から右側に向かって偶数番目の画像が補完対象画像である。

最下段の比較例（一般的な画像の平滑化を使用）では、補完対象画像を含めて、各スタックの画像の左側から下側に亘って格子状の白い点がノイズとして現れている。

一方、実施例１及び実施例３では、血管が元画像の場合と同じように再現されていることが分かる。実施例１を用いるか実施例３を用いるかは、撮影対象（生体組織の部位の違い等）によって、適宜決定することができる。

しかし、第１実施例の方は、補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分のＬ２ノルムを算出するので、より生物画像（例えば、細胞画像）らしい滑らかな画像を再構成することができる。

図２３は本実施の形態の画像処理装置１００によるパラメータ決定の処理の手順の一例を示すフローチャートである。パラメータは、コスト関数で用いる、λ１〜λ４である。以下では、便宜上、処理の主体を制御部１０として説明する。制御部１０は、パラメータ決定のための元画像の画像データを取得する（Ｓ１１）。元画像とは、図２の例では、画像Ｇ１〜Ｇ３、及び補完対象画像Ｃ１〜Ｃ３の全てを含む。すなわち、パラメータを決定するため、まず生体組織の補完対象画像に相当する部分も含めて撮像する。

制御部１０は、パラメータの初期値（例えば、０、１、０．５など適宜の値とすることができる）を設定し（Ｓ１２）、取得した画像データのうち、補完対象画像の画像データを除く（Ｓ１３）。制御部１０は、推定値を初期値に設定する（Ｓ１４）。初期値は、０でも１でもよく、適宜の値を設定することができる。

制御部１０は、式（１）、式（５）又は式（６）のいずれかで表されるコスト関数が小さくなるように推定値（ベクトルＸの各成分）を更新する（Ｓ１５）。

制御部１０は、コスト関数が最小であるか否かを判定し（Ｓ１６）、コスト関数が最小である場合（Ｓ１６でＹＥＳ）、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する（Ｓ１７）。

制御部１０は、元画像の画像データ及び補完後画像の画像データに基づいて、ＭＳＥ（Mean Square Error）を算出し（Ｓ１８）、算出したＭＳＥが閾値以下であるか否かを判定する（Ｓ１９）。算出したＭＳＥが閾値以下である場合（Ｓ１９でＹＥＳ）、制御部１０は、パラメータを決定し（Ｓ２０）、処理を終了する。

算出したＭＳＥが閾値以下でない場合（Ｓ１９でＮＯ）、制御部１０は、パラメータを変更し（Ｓ２１）、ステップＳ１５以降の処理を行う。コスト関数が最小でない場合（Ｓ１６でＮＯ）、制御部１０は、推定値（補完後画像）に対してカーブレット変換を行い（Ｓ２２）、変換後の推定値をコスト関数に代入して（Ｓ２３）、ステップＳ１５以降の処理を行う。

図２４は本実施の形態の画像処理装置１００による画像補完処理の手順の一例を示すフローチャートである。制御部１０は、画像データを取得し（Ｓ３１）、パラメータを取得する（Ｓ３２）。パラメータは、例えば、λ１〜λ４、観測行列Ａなどである。なお、パラメータλ１〜λ４は、図２３で示す処理の結果、決定されたパラメータを用いることができる。

制御部１０は、推定値を初期値に設定する（Ｓ３３）。すなわち、ベクトルＸの各成分を０に設定する。なお、推定値の初期値は、０に限定されず、１など適宜設定することができる。制御部１０は、式（１）、式（５）又は式（６）のいずれかで表されるコスト関数が小さくなるように推定値（ベクトルＸの各成分）を更新する（Ｓ３４）。

制御部１０は、コスト関数が最小であるか否かを判定し（Ｓ３５）、最小でない場合（Ｓ３５でＮＯ）、推定値（補完後画像）に対してカーブレット変換を行い（Ｓ３６）、変換後の推定値をコスト関数に代入して（Ｓ３７）、ステップＳ３４以降の処理を行う。コスト関数が最小である場合（Ｓ３５でＹＥＳ）、制御部１０は、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する（Ｓ３８）。

制御部１０は、高解像度化が必要であるか否かを判定し（Ｓ３９）、高解像度化が必要な場合（Ｓ３９でＹＥＳ）、辞書データを用いて補完後画像を高解像度化し（Ｓ４０）、処理を終了する。高解像度化が必要でない場合（Ｓ３９でＮＯ）、制御部１０は、ステップＳ４０の処理を行うことなく処理を終了する。なお、ステップＳ３９、Ｓ４０の処理は必須ではなく省略してもよい。

上述の処理では、推定値の更新を繰り返す都度、カーブレット変換を実行するが、これに限定されるものではない。例えば、コスト関数が最小になった後で、得られた推定値に対してカーブレット変換を行ってもよい。

上述の例では、推定部１６が推定した画像データによって得られた補完対象画像に対してカーブレット変換を施す場合に、式（１）全体に対してカーブレット変換を施す場合について説明したが、これに限定されない。

例えば、推定部１６が推定した画像データによって得られた補完対象画像に対してカーブレット変換を施す場合に、和算出部１５が、補完後画像に対してカーブレット変換を施して得られた各画素の画素値の和を算出してもよい。

具体的には、関数算出部２１は、式（７）を用いてコスト関数を算出することができる。

この場合、和算出部１５は、補完後画像に対してカーブレット変換を施して得られた各画素の画素値の和を算出する。補完後画像の画像データをベクトルＸで表す。ベクトルＸに対するカーブレット（Curvelet）変換をＣ（Ｘ）で表す。なお、カーブレット係数は、予め準備されたマザーカーブレット（例えば、円弧を拡大又は縮小したような曲線）を回転処理又は平行移動処理を行うことによって生成することができる。和算出部１５は、例えば、||Ｃ（Ｘ）||₁ という式（Ｃ（Ｘ）のＬ１ノルム）で和を算出することができる。カーブレット変換は、エッジのある画像をスパース（疎）に表現できるので、例えば、算出した和を小さくすることにより、本来的に生体組織の画像が有さないエッジ等を効果的に除去することができる。以下、式（７）を用いる場合の処理について説明する。

図２５は本実施の形態の画像処理装置１００によるパラメータ決定の処理の手順の他の例を示すフローチャートである。パラメータは、コスト関数で用いる、λ１〜λ４である。制御部１０は、パラメータ決定のための元画像の画像データを取得する（Ｓ５１）。元画像とは、図２の例では、画像Ｇ１〜Ｇ３、及び補完対象画像Ｃ１〜Ｃ３の全てを含む。すなわち、パラメータを決定するため、まず生体組織の補完対象画像に相当する部分も含めて撮像する。

制御部１０は、パラメータの初期値（例えば、０、１、０．５など適宜の値とすることができる）を設定し（Ｓ５２）、取得した画像データのうち、補完対象画像の画像データを除く（Ｓ５３）。制御部１０は、推定値を初期値に設定する（Ｓ５４）。初期値は、０でも１でもよく、適宜の値を設定することができる。

制御部１０は、式（７）で表されるコスト関数が小さくなるように推定値（ベクトルＸの各成分）を更新する（Ｓ５５）。

制御部１０は、コスト関数が最小であるか否かを判定し（Ｓ５６）、コスト関数が最小である場合（Ｓ５６でＹＥＳ）、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する（Ｓ５７）。

制御部１０は、元画像の画像データ及び補完後画像の画像データに基づいて、ＭＳＥ（Mean Square Error）を算出し（Ｓ５８）、算出したＭＳＥが閾値以下であるか否かを判定する（Ｓ５９）。算出したＭＳＥが閾値以下である場合（Ｓ５９でＹＥＳ）、制御部１０は、パラメータを決定し（Ｓ６０）、処理を終了する。

算出したＭＳＥが閾値以下でない場合（Ｓ５９でＮＯ）、制御部１０は、パラメータを変更し（Ｓ６１）、ステップＳ５５以降の処理を行う。コスト関数が最小でない場合（Ｓ５６でＮＯ）、制御部１０は、ステップＳ５５以降の処理を行う。

図２６は本実施の形態の画像処理装置１００による画像補完処理の手順の他の例を示すフローチャートである。制御部１０は、画像データを取得し（Ｓ７１）、パラメータを取得する（Ｓ７２）。パラメータは、例えば、λ１〜λ４、観測行列Ａなどである。なお、パラメータλ１〜λ４は、図２５で示す処理の結果、決定されたパラメータを用いることができる。

制御部１０は、推定値を初期値に設定する（Ｓ７３）。すなわち、ベクトルＸの各成分を０に設定する。なお、推定値の初期値は、０に限定されず、１など適宜設定することができる。制御部１０は、式（７）で表されるコスト関数が小さくなるように推定値（ベクトルＸの各成分）を更新する（Ｓ７４）。

制御部１０は、コスト関数が最小であるか否かを判定し（Ｓ７５）、最小でない場合（Ｓ７５でＮＯ）、ステップＳ７４以降の処理を行う。コスト関数が最小である場合（Ｓ７５でＹＥＳ）、制御部１０は、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する（Ｓ７６）。

制御部１０は、高解像度化が必要であるか否かを判定し（Ｓ７７）、高解像度化が必要な場合（Ｓ７７でＹＥＳ）、辞書データを用いて補完後画像を高解像度化し（Ｓ７８）、処理を終了する。高解像度化が必要でない場合（Ｓ７７でＮＯ）、制御部１０は、ステップＳ７８の処理を行うことなく処理を終了する。なお、ステップＳ７７、Ｓ７８の処理は必須ではなく省略してもよい。

本実施の形態の画像処理装置１００は、ＣＰＵ（プロセッサ）、ＲＡＭ（メモリ）などを備えたコンピュータを用いて実現することもできる。すなわち、図２３から図２６に示すような、各処理の手順を定めたコンピュータプログラムをコンピュータに備えられたＲＡＭ（メモリ）にロードし、コンピュータプログラムをＣＰＵ（プロセッサ）で実行することにより、画像処理装置１００を実現することができる。

図２７は本実施の形態の画像処理装置の構成の他の例を示す説明図である。図２７において、符号３００は、通常のコンピュータである。コンピュータ３００は、制御部３０、入力部４０、出力部５０、外部Ｉ／Ｆ（インタフェース）部６０などを備える。制御部３０は、ＣＰＵ３１、ＲＯＭ３２、ＲＡＭ３３、Ｉ／Ｆ（インタフェース）３４などを備える。また、コンピュータ３００は、辞書データベース７０に接続することができる。

入力部４０は、計測機器２００から画像データを取得する。また、入力部４０は、パラメータを取得する。出力部５０は、画像補完処理結果である画像データを出力する。Ｉ／Ｆ３４は、制御部３０と、入力部４０、出力部５０及び外部Ｉ／Ｆ部６０それぞれとの間のインタフェース機能を有する。

外部Ｉ／Ｆ部６０は、コンピュータプログラムを記録した記録媒体Ｍ（例えば、ＤＶＤなどのメディア）からコンピュータプログラムを読み取ることが可能である。

なお、図示していないが、記録媒体Ｍに記録されたコンピュータプログラムは、持ち運びが自由なメディアに記録されたものに限定されるものではなく、インターネット又は他の通信回線を通じて伝送されるコンピュータプログラムも含めることができる。また、コンピュータには、複数のプロセッサを搭載した１台のコンピュータ、あるいは、通信ネットワークを介して接続された複数台のコンピュータで構成されるコンピュータシステムも含まれる。

上述の実施の形態では、図２に例示したように、取得画像Ｇ１〜Ｇ３それぞれに対応させて一つの画像を補完する構成であったが、これに限定されない。例えば、取得画像Ｇ１〜Ｇ３それぞれに対応させて複数（例えば、二つ）の画像を補完することもできる。例えば、まず、一つの取得画像に対応させて一つの画像を補完し、次に、当該一つの取得画像と補完した画像との間の画像を補完することによって、一つの取得画像に対応させて二つの画像を補完することができる。

（第２実施形態）
前述の第１実施形態では、ｘｙｚ方向の３次元画像において、例えば、ｚ方向に沿って配置される画像の一部を補完するものであったが、第１実施形態の構成は、２次元画像が時系列に並んだ場合において、時間方向に沿って配置される画像の一部を補完する場合にも適用することができる。以下、第２実施形態について説明する。

図２８は第２実施形態の画像処理装置１００による補完後画像の一例を示す模式図である。画像処理装置１００の構成は第１実施形態の場合と同様である。図中上側の図は取得画像（観測画像）を示し、下側の図は補完後画像を示す。図２８に示すように、取得画像Ｇ１、Ｇ２、Ｇ３は、生体組織のｘｙ平面上の画像であり、時系列（時点ｔ１、ｔ３、ｔ５）に３枚撮影されており、スタック数＝３となる。なお、図２８では、便宜上、取得画像のスタック数を３としているが、スタック数は３に限定されない。

補完後画像は、取得画像Ｇ１、Ｇ２、Ｇ３及び補完対象画像Ｃ１、Ｃ２、Ｃ３を含み、生体組織のｘｙ平面上の画像であり、時間方向に沿って６枚あり、スタック数＝６となる。補完対象画像Ｃ１、Ｃ２、Ｃ３は、実際に撮影された画像ではなく、本実施の形態の画像処理装置１００によって補完される画像である。補完対象画像Ｃ１は、時点ｔ２（時点１と時点ｔ３との間の時点）に対応する画像であり、取得画像Ｇ１とＧ２との間にある。補完対象画像Ｃ２は、時点ｔ４（時点３と時点ｔ５との間の時点）に対応する画像であり、取得画像Ｇ２とＧ３との間にある。補完対象画像Ｃ３は、時点ｔ６（時点５の後の時点）に対応する画像であり、取得画像Ｇ３の後にある。なお、画像の枚数、補完対象画像の時点などは図２３の例に限定されるものではない。

画像データ取得部１１は、時系列に撮像された複数の２次元画像（例えば、ｘｙ平面画像）の画像データを取得する。一つの２次元画像を、ｎ×ｍ画素の画像とし、時系列（時間方向）に撮影された２次元画像の数（スタック数）をｚとすると、画像データは、ｎ×ｍ×ｚ行１列の行列（あるいは、ｎ×ｍ×ｚ個の成分を有するベクトル）で表すことができる。取得した画像データをベクトルＹで表す。

推定部１６は、画像データ取得部１１で取得した画像データに基づいて、当該複数の２次元画像及び時間方向に沿って位置する補完対象２次元画像を含む補完後画像の画像データを推定する。補完対象２次元画像は、撮影された複数の２次元画像によって補完される画像であって、実際には撮影されていない画像である。補完対象２次元画像の数（スタック数）をΔｚとすると、補完後画像の数（スタック数）ｚ′は、ｚ′＝ｚ＋Δｚとなる。補完後画像の画像データは、ｎ×ｍ×ｚ′行１列の行列（あるいは、ｎ×ｍ×ｚ′個の成分を有するベクトル）で表すことができる。補完後画像の画像データをベクトルＸで表す。

より具体的には、推定部１６は、第１差分算出部１３で算出した差分及び和算出部１５で算出した和に基づいて補完後画像の画像データを推定する。生体組織を撮影した生物画像には、隣接する画素の画素値（例えば、輝度値）の差分のスパース性（画素値の差分がゼロとなる特性）、及び観測対象以外の部分のスパース性（画素値がゼロとなる特性）を有するという特徴がある。このようなスパース性に注目することによって、取得した少ない画像データ（例えば、スタック数がｚの画像）から多くの情報（例えば、スタック数ｚ′の画像）を抽出することができる。

隣接する画素の画素値の差分のスパース性に基づき、第１差分算出部１３で算出した差分が小さくなる。また、観測対象以外の部分のスパース性に基づき和算出部１５で算出した和が小さくなる。第１差分算出部１３で算出した差分及び和算出部１５で算出した和をコスト関数とし、まず、ベクトルＸの各成分（推定値）の初期値（例えば、０）をコスト関数に入力してコスト関数が最小化となる推定値を更新し、更新した推定値を使って再度コスト関数を最小化する処理を繰り返すことによって、推定値を求めることができる。なお、コスト関数を算出する際に、第１実施形態と同様に、第２差分算出部１４で算出した差分を含めることができる。

推定値が得られることによって、スタック数ｚの画像（取得した画像）から、スタック数ｚ′の補完後画像を得ることができる。これにより、１枚の２次元画像を撮影する間に、生体組織の動的な変化がある場合でも、当該撮影する間の生体組織の変化を補完後画像（より具体的には、補完対象２次元画像）によって捉えることができるので、生命現象を一層リアルタイムで捉えることが可能となる。

本実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第１差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。

本実施の形態に係るコンピュータプログラムは、コンピュータに、実空間画像を補完させるためのコンピュータプログラムであって、コンピュータに、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する処理と、取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する処理と、前記補完後画像の隣接する画素の画素値の差分を算出する処理と、前記補完後画像の各画素の画素値の和を算出する処理と、算出した差分及び算出した和に基づいて前記補完後画像の画像データを推定する処理とを実行させる。

本実施の形態に係る画像補完方法は、実空間画像を補完する画像補完方法であって、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得部が取得し、取得された画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定部が推定し、前記補完後画像の隣接する画素の画素値の差分を第１差分算出部が算出し、前記補完後画像の各画素の画素値の和を和算出部が算出し、算出された差分及び算出された和に基づいて前記補完後画像の画像データを前記推定部が推定する。

取得部は、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する。一つのｘｙ平面画像を、ｎ×ｍ画素の画像とし、ｚ方向に沿って撮影されたｘｙ平面画像の数（スタック数）をｚとすると、画像データは、ｎ×ｍ×ｚ行１列の行列（あるいは、ｎ×ｍ×ｚ個の成分を有するベクトル）で表すことができる。取得した画像データをベクトルＹで表す。

推定部は、取得した画像データに基づいて、当該複数のｘｙ平面画像及びｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する。補完対象ｘｙ平面画像は、撮影された複数のｘｙ平面画像によって補完される画像であって、実際には撮影されていない画像である。補完対象ｘｙ平面画像の数（スタック数）をΔｚとすると、補完後画像の数（スタック数）ｚ′は、ｚ′＝ｚ＋Δｚとなる。補完後画像の画像データは、ｎ×ｍ×ｚ′行１列の行列（あるいは、ｎ×ｍ×ｚ′個の成分を有するベクトル）で表すことができる。補完後画像の画像データをベクトルＸで表す。

第１差分算出部は、補完後画像の隣接する画素の画素値の差分を算出する。隣接する画素は、例えば、ｘ方向、ｙ方向及びｚ方向に沿った画素とすることができる。

和算出部は、補完後画像の各画素の画素値の和を算出する。

推定部は、第１差分算出部で算出した差分及び和算出部で算出した和に基づいて補完後画像の画像データを推定する。生体組織を撮影した生物画像には、隣接する画素の画素値（例えば、輝度値）の差分のスパース性（画素値の差分がゼロとなる特性）、及び観測対象以外の部分のスパース性（画素値がゼロとなる特性）を有するという特徴がある。このようなスパース性に注目することによって、取得した少ない画像データ（例えば、スタック数がｚの画像）から多くの情報（例えば、スタック数ｚ′の画像）を抽出することができる。

隣接する画素の画素値の差分のスパース性に基づき、第１差分算出部で算出した差分が小さくなる。また、観測対象以外の部分のスパース性に基づき和算出部で算出した和が小さくなる。第１差分算出部で算出した差分及び和算出部で算出した和をコスト関数とし、まず、ベクトルＸの各成分（推定値）の初期値（例えば、０）をコスト関数に入力してコスト関数が最小化となる推定値を更新し、更新した推定値を使って再度コスト関数を最小化する処理を繰り返すことによって、推定値を求めることができる。

推定値が得られることによって、スタック数ｚの画像（取得した画像）から、スタック数ｚ′の補完後画像を得ることができる。これにより、ｚ枚（ｚ＜ｚ′）の光学切片画像を取得するだけでｚ′枚の光学切片画像から得られる情報を抽出することができるので、従来よりも短い撮影時間で生体組織の微細形態を捕らえることができる。

本実施の形態に係る画像処理装置は、前記取得部で取得した前記複数のｘｙ平面画像の画像データと、前記補完後画像の画像データのうち前記複数のｘｙ平面画像に対応する画像データとの差分を算出する第２差分算出部を備え、前記推定部は、前記第２差分算出部で算出した差分に基づいて前記補完後画像の画像データを推定する。

第２差分算出部は、取得した複数のｘｙ平面画像の画像データと、補完後画像の画像データのうち当該複数のｘｙ平面画像に対応する画像データとの差分を算出する。例えば、取得した画像のスタック数ｚ＝３とし、補完後画像のスタック数ｚ′＝６とする。ｚ＝１の取得画像が、ｚ′＝１の補完後画像に対応し、ｚ＝２の取得画像が、ｚ′＝３の補完後画像に対応し、ｚ＝３の取得画像が、ｚ′＝５の補完後画像に対応する場合、かかる対応付けを決定する行列Ａ（観測行列とも称する）を考える。この場合、行列Ａは、３行６列の行列となり、画像データが対応する箇所は１となり、対応しない箇所は０となる。

推定部は、第２差分算出部で算出した差分に基づいて補完後画像の画像データを推定する。上述の例では、例えば、||Ｙ-Ａ・Ｘ||₂ （［Ｙ-Ａ・Ｘ］のＬ２ノルムと称する）を最小二乗法によって小さくするようにして補完後画像の画像データを推定することができる。

これにより、取得した画像データと、補完後画像（具体的には、補完後画像のうちの取得した画像）の画像データとの差に基づいて補完後画像の画像データを推定することができる。これにより、補完後画像の元の画像に対する精度を向上させることができる。

本実施の形態に係る画像処理装置は、前記第１差分算出部で算出した差分、前記和算出部で算出した和及び前記第２差分算出部で算出した差分の合計を、前記補完後画像の画像データの推定値を変数とする所定関数として算出する関数算出部を備え、前記推定部は、前記所定関数が小さくなるように前記推定値を更新して、前記補完後画像の画像データを推定する。

関数算出部は、第１差分算出部で算出した差分、和算出部で算出した和及び第２差分算出部で算出した差分の合計を、補完後画像の画像データの推定値を変数とする所定関数として算出する。所定関数は、コスト関数とも称する。

推定部は、所定関数が小さくなるように推定値を更新して、補完後画像の画像データを推定する。例えば、推定値の初期値を所定関数に代入して、所定関数が小さくなる推定値を求め、求めた推定値をさらに所定関数に代入して所定関数が小さくなる推定値を求めることによって、推定値を更新する。所定関数を小さくする処理を繰り返すことによって、所定関数が最小になったときの推定値を補完後画像の画像データとすることができる。

本実施の形態に係る画像処理装置において、前記第１差分算出部は、前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分を算出する。

第１差分算出部は、補完後画像のｚ方向に沿って隣接する画素の画素値の差分を算出する。ｚ方向に沿って隣接する画素の画素値の差分を算出し、算出した差分が小さくなるように補完後画像の画像データを推定するので、３次元構造の生体組織に適した生物画像を補完することができる。

本実施の形態に係る画像処理装置において、第１差分算出部は、前記補完後画像の注目画素の重み付けされた画素値と、前記注目画素を間にしてｚ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する。

第１差分算出部は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてｚ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する。例えば、補完後画像のｚ方向に沿ったスタックｚ′を、ｚ′＝ｋ、ｋ＋１、ｋ＋２とする。スタックｚ′＝ｋ＋１の補完後画像上の注目画素の画素値をＰ（ｋ＋１）とし、重みけ付け係数をα１とする。同様に、スタックｚ′＝ｋの補完後画像上の当該注目画素に対応する位置の画素の画素値をＰ（ｋ）とし、重みけ付け係数をα０とし、スタックｚ′＝ｋ＋２の補完後画像上の当該注目画素に対応する位置の画素の画素値をＰ（ｋ＋２）とし、重みけ付け係数をα２とする。

第１差分算出部は、α１×Ｐ（ｋ＋１）−{α０×Ｐ（ｋ）＋α２×Ｐ（ｋ＋２）}という式を用いて差分を算出する。ここで、α１−（α０＋α２）＝０とすることができ、例えば、α１＝２、α０＝α２＝１とすると、第１差分算出部は、２×Ｐ（ｋ＋１）−{Ｐ（ｋ）＋Ｐ（ｋ＋２）}という式により差分を算出することができる。これにより、補完後画像の任意の画素の画素値と、ｚ方向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、ｚ方向に沿った偏りを抑制して、３次元構造の生体組織に適した生物画像を補完することができる。

本実施の形態に係る画像処理装置において、前記第１差分算出部は、前記補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分を算出する。

第１差分算出部は、補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分を算出する。ｘ方向及びｙ方向に沿って隣接する画素の画素値の差分を算出し、算出した差分が小さくなるように補完後画像の画像データを推定するので、生体組織に適した生物画像を補完することができる。

本実施の形態に係る画像処理装置において、前記第１差分算出部は、前記補完後画像の注目画素の重み付けされた画素値と、前記注目画素を間にしてｘ方向及びｙ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する。

第１差分算出部は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてｘ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する。また、第１差分算出部は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてｙ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する。

例えば、補完後画像の注目画素の画素値をＰ（ｉ＋１）とし、重み付け係数をα１とする。注目画素を間にしてｘ方向に沿って隣接する一方の画素の画素値をＰ（ｉ）とし、重み付け係数をα０とする。注目画素を間にしてｘ方向に沿って隣接する他方の画素の画素値をＰ（ｉ＋２）とし、重み付け係数をα２とする。

第１差分算出部は、α１×Ｐ（ｉ＋１）−{α０×Ｐ（ｉ）＋α２×Ｐ（ｉ＋２）}という式を用いて差分を算出する。ここで、α１−（α０＋α２）＝０とすることができ、例えば、α１＝２、α０＝α２＝１とすると、第１差分算出部は、２×Ｐ（ｉ＋１）−{Ｐ（ｉ）＋Ｐ（ｉ＋２）}という式により差分を算出することができる。これにより、補完後画像の任意の画素の画素値と、ｘ向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、ｘ方向に沿った偏りを抑制して、３次元構造の生体組織に適した生物画像を補完することができる。なお、ｙ方向についても、ｘ方向と同様である。

本実施の形態に係る画像処理装置において、前記第１差分算出部は、前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分のＬ１ノルムを算出し、前記推定部は、前記第１差分算出部で算出したＬ１ノルムが最小になるように前記補完後画像の画像データを推定する。

第１差分算出部は、補完後画像のｚ方向に沿って隣接する画素の画素値の差分のＬ１ノルムを算出する。例えば、ベクトルＸの成分を（ｘ１、ｘ２、…ｘｎ）とすると、ベクトルＸのＬ１ノルムは、{|ｘ１|＋|ｘ２|＋…＋|ｘｎ|}という式で表すことができる。

推定部は、第１差分算出部で算出したＬ１ノルムが最小になるように補完後画像の画像データを推定する。ｚ方向に沿って隣接する画素の画素値の情報は、補完対象のスタックでは存在しない。そこで、Ｌ１ノルムを用いることにより、情報が存在しない場合でも、スパース解（解の中にゼロが含まれる）が得られるようにして、より正確な補完を可能にして、滑らかな生体画像を得ることができる。

本実施の形態に係る画像処理装置において、前記第１差分算出部は、前記補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分のＬ２ノルムを算出し、前記推定部は、前記第１差分算出部で算出したＬ２ノルムが最小になるように前記補完後画像の画像データを推定する。

第１差分算出部は、補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分のＬ２ノルムを算出する。例えば、ベクトルＸの成分を（ｘ１、ｘ２、…ｘｎ）とすると、ベクトルＸのＬ２ノルムは、√{（ｘ１）² ＋（ｘ２）² ＋…（ｘｎ）² }という式で表すことができる。

推定部は、第１差分算出部で算出したＬ２ノルムが最小になるように補完後画像の画像データを推定する。取得した画像は、ｘｙ平面画像であるので、ｘ方向及びｙ方向に沿って隣接する画素の画素値の情報はもともと存在している。そこで、Ｌ２ノルムを用いることにより、存在する情報を用いて滑らかな生体画像を得ることができる。

本実施の形態に係る画像処理装置は、前記推定部が推定した画像データによって得られた前記補完後画像のｘｚ平面画像又はｙｚ平面画像に対してカーブレット変換を施す変換部を備える。

変換部は、推定部が推定した画像データによって得られた補完後画像のｘｚ平面画像又はｙｚ平面画像に対してカーブレット変換を施す。カーブレット変換は、曲線の集まりとして捉えられる画像を抽出する手法である。予め種々の曲線（生体組織を撮影して得られる曲線部分）を辞書データとして収集しておく。補完後画像のｘｚ平面画像又はｙｚ平面画像を所定の大きさの画素領域（例えば、３０×３０画素で構成される画素領域）毎に走査して、辞書データに含まれる曲線に適合するか否かの判定を行い、曲線に適合した部分だけを抽出する。これにより、補完後画像のｘｚ平面画像又はｙｚ平面画像に含まれるノイズ（例えば、生体組織に相当しない、点など）を除去して補完精度を向上させることができる。なお、補完後画像のｘｚ平面画像に対してカーブレット変換を施すこともできる。

本実施の形態に係る画像処理装置において、前記和算出部は、前記補完後画像に対してカーブレット変換を施して得られた各画素の画素値の和を算出する。

和算出部は、補完後画像に対してカーブレット変換を施して得られた各画素の画素値の和を算出する。補完後画像の画像データをベクトルＸで表す。ベクトルＸに対するカーブレット（Curvelet）変換をＣ（Ｘ）で表す。和算出部は、例えば、||Ｃ（Ｘ）||₁ という式（Ｃ（Ｘ）のＬ１ノルム）で和を算出することができる。カーブレット変換は、エッジのある画像をスパース（疎）に表現できるので、例えば、算出した和を小さくすることにより、本来的に生体組織の画像が有さないエッジ等を効果的に除去することができる。

本実施の形態に係る画像処理装置は、複数の高解像度画像と該複数の高解像度画像を疑似的に劣化させた複数の低解像度画像とを対応付けた辞書データを記憶する記憶部と、該記憶部に記憶した複数の低解像度画像の中から前記推定部が推定した画像データによって得られた前記補完後画像を構成する複数の低解像度画像を特定する特定部と、該特定部で特定した複数の低解像度画像に対応する複数の高解像度画像に基づいて高解像度補完後画像を生成する生成部とを備える。

記憶部は、複数の高解像度画像と当該複数の高解像度画像を疑似的に劣化させた複数の低解像度画像とを対応付けた辞書データを記憶する。複数の高解像度画像としては、予め種々の生体組織を撮影した画像を用いることができる。高解像度画像は、例えば、６０×６０画素の画像とし、低解像度画像は、３０×３０画素の画像とすることができるが、これに限定されない。

特定部は、記憶部に記憶した複数の低解像度画像の中から推定部が推定した画像データによって得られた補完後画像を構成する複数の低解像度画像を特定する。例えば、補完後画像を、所定の大きさの画素領域（例えば、３０×３０画素で構成される画素領域）毎に走査して、一致又は類似する低解像度画像を特定する。

生成部は、特定部で特定した複数の低解像度画像に対応する複数の高解像度画像に基づいて高解像度補完画像を生成する。これにより、例えば、５１２×５１２で撮影した補完後画像を、１０２４×１０２４の高解像度画像にすることができる。

本実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、時系列に撮像された複数の２次元画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数の２次元画像及び時間方向に沿って位置する補完対象２次元画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第１差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。

取得部は、時系列に撮像された複数の２次元画像（例えば、ｘｙ平面画像）の画像データを取得する。一つの２次元画像を、ｎ×ｍ画素の画像とし、時系列（時間方向）に撮影された２次元画像の数（スタック数）をｚとすると、画像データは、ｎ×ｍ×ｚ行１列の行列（あるいは、ｎ×ｍ×ｚ個の成分を有するベクトル）で表すことができる。取得した画像データをベクトルＹで表す。

推定部は、取得した画像データに基づいて、当該複数の２次元画像及び時間方向に沿って位置する補完対象２次元画像を含む補完後画像の画像データを推定する。補完対象２次元画像は、撮影された複数の２次元画像によって補完される画像であって、実際には撮影されていない画像である。補完対象２次元画像の数（スタック数）をΔｚとすると、補完後画像の数（スタック数）ｚ′は、ｚ′＝ｚ＋Δｚとなる。補完後画像の画像データは、ｎ×ｍ×ｚ′行１列の行列（あるいは、ｎ×ｍ×ｚ′個の成分を有するベクトル）で表すことができる。補完後画像の画像データをベクトルＸで表す。

第１差分算出部は、補完後画像の隣接する画素の画素値の差分を算出する。隣接する画素は、例えば、ｘ方向、ｙ方向及び時間方向に沿った画素とすることができる。

和算出部は、補完後画像の各画素の画素値の和を算出する。

推定部は、第１差分算出部で算出した差分及び和算出部で算出した和に基づいて補完後画像の画像データを推定する。生体組織を撮影した生物画像には、隣接する画素の画素値（例えば、輝度値）の差分のスパース性（画素値の差分がゼロとなる特性）、及び観測対象以外の部分のスパース性（画素値がゼロとなる特性）を有するという特徴がある。このようなスパース性に注目することによって、取得した少ない画像データ（例えば、スタック数がｚの画像）から多くの情報（例えば、スタック数ｚ′の画像）を抽出することができる。

隣接する画素の画素値の差分のスパース性に基づき、第１差分算出部で算出した差分が小さくなる。また、観測対象以外の部分のスパース性に基づき和算出部で算出した和が小さくなる。第１差分算出部で算出した差分及び和算出部で算出した和をコスト関数とし、まず、ベクトルＸの各成分（推定値）の初期値（例えば、０）をコスト関数に入力してコスト関数が最小化となる推定値を更新し、更新した推定値を使って再度コスト関数を最小化する処理を繰り返すことによって、推定値を求めることができる。

推定値が得られることによって、スタック数ｚの画像（取得した画像）から、スタック数ｚ′の補完後画像を得ることができる。これにより、１枚の２次元画像を撮影する間に、生体組織の動的な変化がある場合でも、当該撮影する間の生体組織の変化を補完後画像（より具体的には、補完対象２次元画像）によって捉えることができるので、生命現象を一層リアルタイムで捉えることが可能となる。

１０ 制御部
１１ 画像データ取得部
１２ パラメータ取得部
１３ 第１差分算出部
１４ 第２差分算出部
１５ 和算出部
１６ 推定部
１７ 補完画像生成部
１８ 変換部
１９ 記憶部
２０ 特定部
２１ 関数算出部
３０ 制御部
３１ ＣＰＵ
３２ ＲＯＭ
３３ ＲＡＭ
３４ Ｉ／Ｆ
４０ 入力部
５０ 出力部
６０ 外部Ｉ／Ｆ部
７０ 辞書データベース
１００ 画像処理装置
２００ 計測機器

【０００２】
とができる画像処理装置、コンピュータプログラム及び画像補完方法を提供することを目的とする。
課題を解決するための手段
［０００７］
本開示の実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像のｘ方向に沿って隣接する画素の画素値の差分のノルム、前記補完後画像のｙ方向に沿って隣接する画素の画素値のノルム、及び前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分のノルムの合計を算出する第１差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第１差分算出部で算出した合計及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。
［０００８］
本開示の実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、時系列に撮像された複数の２次元画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数の２次元画像及び時間方向に沿って位置する補完対象２次元画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像のｘ方向に沿って隣接する画素の画素値の差分のノルム、前記補完後画像のｙ方向に沿って隣接する画素の画素値のノルム、及び前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分のノルムの合計を算出する第１差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第１差分算出部で算出した合計及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。
［０００９］
本開示の実施の形態に係るコンピュータプログラムは、コンピュータに、実空間画像を補完させるためのコンピュータプログラムであって、コンピュータに、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する処理と、取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する処理と、前記補完後画像のｘ方向に沿って隣接する画素の画素値の差分のノルム、前記補完後画像のｙ方向に沿って隣接する画素の画素値のノルム、及び前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分のノルムの合計を算出する処理と、前記補完後画像の各画素の画素値の和を算出する処理と、

【０００３】
算出した合計及び算出した和に基づいて前記補完後画像の画像データを推定する処理とを実行させる。
［００１０］
本開示の実施の形態に係る画像補完方法は、実空間画像を補完する画像補完方法であって、ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得部が取得し、取得された画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定部が推定し、前記補完後画像のｘ方向に沿って隣接する画素の画素値の差分のノルム、前記補完後画像のｙ方向に沿って隣接する画素の画素値のノルム、及び前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分のノルムの合計を第１差分算出部が算出し、前記補完後画像の各画素の画素値の和を和算出部が算出し、算出された合計及び算出された和に基づいて前記補完後画像の画像データを前記推定部が推定する。
発明の効果
［００１１］
本開示によれば、撮影時間を短縮することができる。
図面の簡単な説明
［００１２］
［図１］本実施の形態の画像処理装置の構成の一例を示すブロック図である。
［図２］本実施の形態の画像処理装置による補完後画像の一例を示す模式図である。
［図３］取得した画像データのベクトルの一例を示す模式図である。
［図４］補完後画像の画像データのベクトルの一例を示す模式図である。
［図５］観測行列の一例を示す模式図である。
［図６］ｘ方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。
［図７］ｙ方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。
［図８］ｚ方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。
［図９］ｘ方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。
［図１０］ｙ方向に沿って隣接する画素の画素値の差分を算出する差分行列の他

Claims (15)

1. 実空間画像を補完する画像処理装置であって、
   ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する取得部と、
   該取得部で取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する推定部と、
   前記補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、
   前記補完後画像の各画素の画素値の和を算出する和算出部と
   を備え、
   前記推定部は、
   前記第１差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する画像処理装置。
2. 前記取得部で取得した前記複数のｘｙ平面画像の画像データと、前記補完後画像の画像データのうち前記複数のｘｙ平面画像に対応する画像データとの差分を算出する第２差分算出部を備え、
   前記推定部は、
   前記第２差分算出部で算出した差分に基づいて前記補完後画像の画像データを推定する請求項１に記載の画像処理装置。
3. 前記第１差分算出部で算出した差分、前記和算出部で算出した和及び前記第２差分算出部で算出した差分の合計を、前記補完後画像の画像データの推定値を変数とする所定関数として算出する関数算出部を備え、
   前記推定部は、
   前記所定関数が小さくなるように前記推定値を更新して、前記補完後画像の画像データを推定する請求項２に記載の画像処理装置。
4. 前記第１差分算出部は、
   前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分を算出する請求項１から請求項３のいずれか一項に記載の画像処理装置。
5. 前記第１差分算出部は、
   前記補完後画像の注目画素の重み付けされた画素値と、前記注目画素を間にしてｚ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する請求項４に記載の画像処理装置。
6. 前記第１差分算出部は、
   前記補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分を算出する請求項１から請求項５のいずれか一項に記載の画像処理装置。
7. 前記第１差分算出部は、
   前記補完後画像の注目画素の重み付けされた画素値と、前記注目画素を間にしてｘ方向及びｙ方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する請求項６に記載の画像処理装置。
8. 前記第１差分算出部は、
   前記補完後画像のｚ方向に沿って隣接する画素の画素値の差分のＬ１ノルムを算出し、
   前記推定部は、
   前記第１差分算出部で算出したＬ１ノルムが最小になるように前記補完後画像の画像データを推定する請求項４又は請求項５に記載の画像処理装置。
9. 前記第１差分算出部は、
   前記補完後画像のｘ方向及びｙ方向に沿って隣接する画素の画素値の差分のＬ２ノルムを算出し、
   前記推定部は、
   前記第１差分算出部で算出したＬ２ノルムが最小になるように前記補完後画像の画像データを推定する請求項６又は請求項７に記載の画像処理装置。
10. 前記推定部が推定した画像データによって得られた前記補完後画像のｘｚ平面画像又はｙｚ平面画像に対してカーブレット変換を施す変換部を備える請求項１から請求項９のいずれか一項に記載の画像処理装置。
11. 前記和算出部は、
    前記補完後画像に対してカーブレット変換を施して得られた各画素の画素値の和を算出する請求項１から請求項９のいずれか一項に記載の画像処理装置。
12. 複数の高解像度画像と該複数の高解像度画像を疑似的に劣化させた複数の低解像度画像とを対応付けた辞書データを記憶する記憶部と、
    該記憶部に記憶した複数の低解像度画像の中から前記推定部が推定した画像データによって得られた前記補完後画像を構成する複数の低解像度画像を特定する特定部と、
    該特定部で特定した複数の低解像度画像に対応する複数の高解像度画像に基づいて高解像度補完後画像を生成する生成部と
    を備える請求項１から請求項１１のいずれか一項に記載の画像処理装置。
13. 実空間画像を補完する画像処理装置であって、
    時系列に撮像された複数の２次元画像の画像データを取得する取得部と、
    該取得部で取得した画像データに基づいて、前記複数の２次元画像及び時間方向に沿って位置する補完対象２次元画像を含む補完後画像の画像データを推定する推定部と、
    前記補完後画像の隣接する画素の画素値の差分を算出する第１差分算出部と、
    前記補完後画像の各画素の画素値の和を算出する和算出部と
    を備え、
    前記推定部は、
    前記第１差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する画像処理装置。
14. コンピュータに、実空間画像を補完させるためのコンピュータプログラムであって、
    コンピュータに、
    ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得する処理と、
    取得した画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定する処理と、
    前記補完後画像の隣接する画素の画素値の差分を算出する処理と、
    前記補完後画像の各画素の画素値の和を算出する処理と、
    算出した差分及び算出した和に基づいて前記補完後画像の画像データを推定する処理と
    を実行させるコンピュータプログラム。
15. 実空間画像を補完する画像補完方法であって、
    ｚ方向に沿って撮像された複数のｘｙ平面画像の画像データを取得部が取得し、
    取得された画像データに基づいて、前記複数のｘｙ平面画像及び前記ｚ方向に沿って位置する補完対象ｘｙ平面画像を含む補完後画像の画像データを推定部が推定し、
    前記補完後画像の隣接する画素の画素値の差分を第１差分算出部が算出し、
    前記補完後画像の各画素の画素値の和を和算出部が算出し、
    算出された差分及び算出された和に基づいて前記補完後画像の画像データを前記推定部が推定する画像補完方法。