TOP > 国内特許検索 > 画像処理装置、コンピュータプログラム及び画像補完方法 > 明細書

明細書 :画像処理装置、コンピュータプログラム及び画像補完方法

発行国 日本国特許庁(JP)
公報種別 再公表特許(A1)
発行日 令和2年4月30日(2020.4.30)
発明の名称または考案の名称 画像処理装置、コンピュータプログラム及び画像補完方法
国際特許分類 G06T   1/00        (2006.01)
FI G06T 1/00 295
国際予備審査の請求
全頁数 43
出願番号 特願2019-525492 (P2019-525492)
国際出願番号 PCT/JP2018/022618
国際公開番号 WO2018/230615
国際出願日 平成30年6月13日(2018.6.13)
国際公開日 平成30年12月20日(2018.12.20)
優先権出願番号 2017117213
優先日 平成29年6月14日(2017.6.14)
優先権主張国 日本国(JP)
指定国 AP(BW , GH , GM , KE , LR , LS , MW , MZ , NA , RW , SD , SL , ST , SZ , TZ , UG , ZM , ZW) , EA(AM , AZ , BY , KG , KZ , RU , TJ , TM) , EP(AL , AT , BE , BG , CH , CY , CZ , DE , DK , EE , ES , FI , FR , GB , GR , HR , HU , IE , IS , IT , LT , LU , LV , MC , MK , MT , NL , NO , PL , PT , RO , RS , SE , SI , SK , SM , TR) , OA(BF , BJ , CF , CG , CI , CM , GA , GN , GQ , GW , KM , ML , MR , NE , SN , TD , TG) , AE , AG , AL , AM , AO , AT , AU , AZ , BA , BB , BG , BH , BN , BR , BW , BY , BZ , CA , CH , CL , CN , CO , CR , CU , CZ , DE , DJ , DK , DM , DO , DZ , EC , EE , EG , ES , FI , GB , GD , GE , GH , GM , GT , HN , HR , HU , ID , IL , IN , IR , IS , JO , JP , KE , KG , KH , KN , KP , KR , KW , KZ , LA , LC , LK , LR , LS , LU , LY , MA , MD , ME , MG , MK , MN , MW , MX , MY , MZ , NA , NG , NI , NO , NZ , OM , PA , PE , PG , PH , PL , PT , QA , RO , RS , RU , RW , SA , SC , SD , SE , SG , SK , SL , SM , ST , SV , SY , TH , TJ , TM , TN , TR , TT
発明者または考案者 【氏名】曽我部 舞奈
【氏名】大関 真之
【氏名】瀬原 淳子
出願人 【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
個別代理人の代理人 【識別番号】100114557、【弁理士】、【氏名又は名称】河野 英仁
【識別番号】100078868、【弁理士】、【氏名又は名称】河野 登夫
審査請求 未請求
テーマコード 5B057
Fターム 5B057CA08
5B057CA12
5B057CA16
5B057CB08
5B057CB12
5B057CB16
5B057CD06
要約 撮影時間を短縮することができる画像処理装置、コンピュータプログラム及び画像補完方法を提供する。
画像処理装置は、z方向に沿って撮像された複数のxy平面画像の画像データを取得する取得部と、取得部で取得した画像データに基づいて、当該複数のxy平面画像及びz方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定する推定部と、補完後画像の隣接する画素の画素値の差分を算出する第1差分算出部と、補完後画像の各画素の画素値の和を算出する和算出部とを備え、推定部は、第1差分算出部で算出した差分及び和算出部で算出した和に基づいて補完後画像の画像データを推定する。
特許請求の範囲 【請求項1】
実空間画像を補完する画像処理装置であって、
z方向に沿って撮像された複数のxy平面画像の画像データを取得する取得部と、
該取得部で取得した画像データに基づいて、前記複数のxy平面画像及び前記z方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定する推定部と、
前記補完後画像の隣接する画素の画素値の差分を算出する第1差分算出部と、
前記補完後画像の各画素の画素値の和を算出する和算出部と
を備え、
前記推定部は、
前記第1差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する画像処理装置。
【請求項2】
前記取得部で取得した前記複数のxy平面画像の画像データと、前記補完後画像の画像データのうち前記複数のxy平面画像に対応する画像データとの差分を算出する第2差分算出部を備え、
前記推定部は、
前記第2差分算出部で算出した差分に基づいて前記補完後画像の画像データを推定する請求項1に記載の画像処理装置。
【請求項3】
前記第1差分算出部で算出した差分、前記和算出部で算出した和及び前記第2差分算出部で算出した差分の合計を、前記補完後画像の画像データの推定値を変数とする所定関数として算出する関数算出部を備え、
前記推定部は、
前記所定関数が小さくなるように前記推定値を更新して、前記補完後画像の画像データを推定する請求項2に記載の画像処理装置。
【請求項4】
前記第1差分算出部は、
前記補完後画像のz方向に沿って隣接する画素の画素値の差分を算出する請求項1から請求項3のいずれか一項に記載の画像処理装置。
【請求項5】
前記第1差分算出部は、
前記補完後画像の注目画素の重み付けされた画素値と、前記注目画素を間にしてz方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する請求項4に記載の画像処理装置。
【請求項6】
前記第1差分算出部は、
前記補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分を算出する請求項1から請求項5のいずれか一項に記載の画像処理装置。
【請求項7】
前記第1差分算出部は、
前記補完後画像の注目画素の重み付けされた画素値と、前記注目画素を間にしてx方向及びy方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する請求項6に記載の画像処理装置。
【請求項8】
前記第1差分算出部は、
前記補完後画像のz方向に沿って隣接する画素の画素値の差分のL1ノルムを算出し、
前記推定部は、
前記第1差分算出部で算出したL1ノルムが最小になるように前記補完後画像の画像データを推定する請求項4又は請求項5に記載の画像処理装置。
【請求項9】
前記第1差分算出部は、
前記補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分のL2ノルムを算出し、
前記推定部は、
前記第1差分算出部で算出したL2ノルムが最小になるように前記補完後画像の画像データを推定する請求項6又は請求項7に記載の画像処理装置。
【請求項10】
前記推定部が推定した画像データによって得られた前記補完後画像のxz平面画像又はyz平面画像に対してカーブレット変換を施す変換部を備える請求項1から請求項9のいずれか一項に記載の画像処理装置。
【請求項11】
前記和算出部は、
前記補完後画像に対してカーブレット変換を施して得られた各画素の画素値の和を算出する請求項1から請求項9のいずれか一項に記載の画像処理装置。
【請求項12】
複数の高解像度画像と該複数の高解像度画像を疑似的に劣化させた複数の低解像度画像とを対応付けた辞書データを記憶する記憶部と、
該記憶部に記憶した複数の低解像度画像の中から前記推定部が推定した画像データによって得られた前記補完後画像を構成する複数の低解像度画像を特定する特定部と、
該特定部で特定した複数の低解像度画像に対応する複数の高解像度画像に基づいて高解像度補完後画像を生成する生成部と
を備える請求項1から請求項11のいずれか一項に記載の画像処理装置。
【請求項13】
実空間画像を補完する画像処理装置であって、
時系列に撮像された複数の2次元画像の画像データを取得する取得部と、
該取得部で取得した画像データに基づいて、前記複数の2次元画像及び時間方向に沿って位置する補完対象2次元画像を含む補完後画像の画像データを推定する推定部と、
前記補完後画像の隣接する画素の画素値の差分を算出する第1差分算出部と、
前記補完後画像の各画素の画素値の和を算出する和算出部と
を備え、
前記推定部は、
前記第1差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する画像処理装置。
【請求項14】
コンピュータに、実空間画像を補完させるためのコンピュータプログラムであって、
コンピュータに、
z方向に沿って撮像された複数のxy平面画像の画像データを取得する処理と、
取得した画像データに基づいて、前記複数のxy平面画像及び前記z方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定する処理と、
前記補完後画像の隣接する画素の画素値の差分を算出する処理と、
前記補完後画像の各画素の画素値の和を算出する処理と、
算出した差分及び算出した和に基づいて前記補完後画像の画像データを推定する処理と
を実行させるコンピュータプログラム。
【請求項15】
実空間画像を補完する画像補完方法であって、
z方向に沿って撮像された複数のxy平面画像の画像データを取得部が取得し、
取得された画像データに基づいて、前記複数のxy平面画像及び前記z方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定部が推定し、
前記補完後画像の隣接する画素の画素値の差分を第1差分算出部が算出し、
前記補完後画像の各画素の画素値の和を和算出部が算出し、
算出された差分及び算出された和に基づいて前記補完後画像の画像データを前記推定部が推定する画像補完方法。
発明の詳細な説明 【技術分野】
【0001】
本発明は、画像処理装置、コンピュータプログラム及び画像補完方法に関する。
【背景技術】
【0002】
生体イメージングは、医学、生物学の分野で注目されている新しい技術である。生体イメージングでは、生体組織の長時間イメージングによって、生命現象をリアルタイムで捉える。生体組織の立体的画像を得るため、例えば、60μm~100μmの間を数十枚程度の光学切片画像を取得する必要がある。
【0003】
生体組織の立体的画像を得るためには、例えば、多光子励起顕微鏡、共焦点顕微鏡、蛍光顕微鏡などの計測機器が用いられる。非特許文献1には、2光子励起顕微鏡についての改善が開示されている。
【先行技術文献】
【0004】

【非特許文献1】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
【発明の概要】
【発明が解決しようとする課題】
【0005】
生体組織の微細形態を捕らえるためには、多数枚の光学切片画像を取得する必要があり、また高解像度の撮像を必要とする。このため、1回当たりの撮影時間(所要枚数の光学切片画像の撮影時間)が長くなるという問題があった。
【0006】
本開示は斯かる事情に鑑みてなされたものであり、撮影時間を短縮することができる画像処理装置、コンピュータプログラム及び画像補完方法を提供することを目的とする。
【課題を解決するための手段】
【0007】
本開示の実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、z方向に沿って撮像された複数のxy平面画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数のxy平面画像及び前記z方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像の隣接する画素の画素値の差分を算出する第1差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第1差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。
【0008】
本開示の実施の形態に係る画像処理装置は、実空間画像を補完する画像処理装置であって、時系列に撮像された複数の2次元画像の画像データを取得する取得部と、該取得部で取得した画像データに基づいて、前記複数の2次元画像及び時間方向に沿って位置する補完対象2次元画像を含む補完後画像の画像データを推定する推定部と、前記補完後画像の隣接する画素の画素値の差分を算出する第1差分算出部と、前記補完後画像の各画素の画素値の和を算出する和算出部とを備え、前記推定部は、前記第1差分算出部で算出した差分及び前記和算出部で算出した和に基づいて前記補完後画像の画像データを推定する。
【0009】
本開示の実施の形態に係るコンピュータプログラムは、コンピュータに、実空間画像を補完させるためのコンピュータプログラムであって、コンピュータに、z方向に沿って撮像された複数のxy平面画像の画像データを取得する処理と、取得した画像データに基づいて、前記複数のxy平面画像及び前記z方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定する処理と、前記補完後画像の隣接する画素の画素値の差分を算出する処理と、前記補完後画像の各画素の画素値の和を算出する処理と、算出した差分及び算出した和に基づいて前記補完後画像の画像データを推定する処理とを実行させる。
【0010】
本開示の実施の形態に係る画像補完方法は、実空間画像を補完する画像補完方法であって、z方向に沿って撮像された複数のxy平面画像の画像データを取得部が取得し、取得された画像データに基づいて、前記複数のxy平面画像及び前記z方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定部が推定し、前記補完後画像の隣接する画素の画素値の差分を第1差分算出部が算出し、前記補完後画像の各画素の画素値の和を和算出部が算出し、算出された差分及び算出された和に基づいて前記補完後画像の画像データを前記推定部が推定する。
【発明の効果】
【0011】
本開示によれば、撮影時間を短縮することができる。
【図面の簡単な説明】
【0012】
【図1】本実施の形態の画像処理装置の構成の一例を示すブロック図である。
【図2】本実施の形態の画像処理装置による補完後画像の一例を示す模式図である。
【図3】取得した画像データのベクトルの一例を示す模式図である。
【図4】補完後画像の画像データのベクトルの一例を示す模式図である。
【図5】観測行列の一例を示す模式図である。
【図6】x方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。
【図7】y方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。
【図8】z方向に沿って隣接する画素の画素値の差分を算出する差分行列の一例を示す模式図である。
【図9】x方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。
【図10】y方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。
【図11】z方向に沿って隣接する画素の画素値の差分を算出する差分行列の他の例を示す模式図である。
【図12】各行列の大きさの一例を示す説明図である。
【図13】本実施の形態の画像処理装置による高速画像取得の様子の一例を示す模式図である。
【図14】本実施の形態の画像処理装置による変換前後の画像の一例を示す説明図である。
【図15】本実施の形態の画像処理装置によるカーブレット変換の一例を示す模式図である。
【図16】本実施の形態の画像処理装置によるカーブレット変換後の画像の一例を示す説明図である。
【図17】本実施の形態の画像処理装置による補完後画像の第1例を示す説明図である。
【図18】本実施の形態の画像処理装置による補完後画像の第2例を示す説明図である。
【図19】本実施の形態の画像処理装置による補完後画像の精度評価の一例を示す説明図である。
【図20】本実施の形態の画像処理装置による補完後画像の高解像度化の様子の一例を示す模式図である。
【図21】本実施の形態の画像処理装置による撮影時間の一例を示す説明図である。
【図22】本実施の形態の画像処理装置による補完後画像の第3例を示す説明図である。
【図23】本実施の形態の画像処理装置によるパラメータ決定の処理の手順の一例を示すフローチャートである。
【図24】本実施の形態の画像処理装置による画像補完処理の手順の一例を示すフローチャートである。
【図25】本実施の形態の画像処理装置によるパラメータ決定の処理の手順の他の例を示すフローチャートである。
【図26】本実施の形態の画像処理装置による画像補完処理の手順の他の例を示すフローチャートである。
【図27】本実施の形態の画像処理装置の構成の他の例を示す説明図である。
【図28】第2実施形態の画像処理装置による補完後画像の一例を示す模式図である。
【発明を実施するための形態】
【0013】
(第1実施形態)
以下、本発明をその実施の形態を示す図面に基づいて説明する。図1は本実施の形態の画像処理装置100の構成の一例を示すブロック図である。画像処理装置100は、計測機器200に接続されている。計測機器200は、生体組織の立体的画像などの生物画像を撮影することができる機器であり、例えば、多光子励起顕微鏡、共焦点顕微鏡、蛍光顕微鏡、ライトシート型顕微鏡、イオンコンダクタンス顕微鏡、光緩衝断層計などが含まれる。

【0014】
画像処理装置100は、装置全体を制御する制御部10、画像データ取得部11、パラメータ取得部12、第1差分算出部13、第2差分算出部14、和算出部15、推定部16、補完画像生成部17、変換部18、記憶部19、特定部20及び関数算出部21などを備える。

【0015】
画像データ取得部11は、取得部としての機能を有し、生体組織のz方向に沿って撮像された複数のxy平面画像(取得画像、観測画像とも称する)の画像データ(観測値とも称する)を計測機器200から取得する。

【0016】
一つのxy平面画像(光学切片画像とも称する)を、n×m画素の画像とし、z方向に沿って撮影されたxy平面画像の数(スタック数)をzとすると、画像データは、n×m×z行1列の行列(あるいは、n×m×z個の成分を有するベクトル)で表すことができる。なお、例えば、n×m×zは、n*m*zの如く表記することもできる。取得した画像データをベクトルYで表す。画像データは、各画素の画素値を含む、画素値は、例えば、0~255で表される輝度値とすることができる。すなわち、画像データは、例えば、各画素の輝度値のデータであるといえる。

【0017】
推定部16は、画像データ取得部11で取得した画像データに基づいて、当該複数の取得画像(xy平面画像)及びz方向に沿って位置する補完対象画像(補完対象xy平面画像)を含む補完後画像の画像データを推定する。

【0018】
補完対象画像は、撮影された複数の取得画像によって補完される画像であって、実際には撮影されていない画像である。補完対象画像の数(スタック数)をΔzとすると、補完後画像の数(スタック数)z′は、z′=z+Δzとなる。補完後画像の画像データは、n×m×z′行1列の行列(あるいは、n×m×z′個の成分を有するベクトル)で表すことができる。補完後画像の画像データをベクトルXで表す。

【0019】
図2は本実施の形態の画像処理装置100による補完後画像の一例を示す模式図である。図中左側の図は取得画像(観測画像)を示し、右側の図は補完後画像を示す。図2に示すように、取得画像G1、G2、G3は、生体組織のxy平面上の画像であり、z方向に沿って3枚撮影されており、スタック数=3となる。なお、図2では、便宜上、取得画像のスタック数を3としているが、スタック数は3に限定されない。

【0020】
補完後画像は、取得画像G1、G2、G3及び補完対象画像C1、C2、C3を含み、生体組織のxy平面上の画像であり、z方向に沿って6枚あり、スタック数=6となる。補完対象画像C1、C2、C3は、実際に撮影された画像ではなく、本実施の形態の画像処理装置100によって補完される画像である。補完対象画像C1は、取得画像G1とG2との間にあり、補完対象画像C2は、取得画像G2とG3との間にあり、補完対象画像C3は、取得画像G3よりも深い位置にある。なお、画像の枚数、補完対象画像の位置などは図2の例に限定されるものではない。

【0021】
図3は取得した画像データのベクトルYの一例を示す模式図である。便宜上、取得画像が3×3画素で構成されるとする。取得画像G1(z=1)の画素は9個存在し、画素番号(ピクセル番号)を1~9とする。同様に、取得画像G2(z=2)の画素番号(ピクセル番号)を10~18とし、取得画像G3(z=3)の画素番号(ピクセル番号)を19~27とする。画像データのベクトルYは、1~27の27次元ベクトルで表される。ベクトルYの各成分は、それぞれの画素番号における観測値となる。

【0022】
図4は補完後画像の画像データのベクトルXの一例を示す模式図である。補完後画像は、取得画像G1~G3及び補完対象画像C1、C2、C3で構成される。図4に示すように、補完後画像のスタックz=1、z=3、z=5の画像が、取得画像G1(z=1)、取得画像G2(z=2)、取得画像G3(z=3)に対応している。この対応関係を決定するのが、後述の観測行列Aである。補完後画像は、6枚の画像で構成されるので、ベクトルXは、1~54の54次元ベクトルで表される。ベクトルYの各成分は、それぞれの画素番号における推定値となる。

【0023】
図5は観測行列Aの一例を示す模式図である。観測行列Aは、取得した画像データのベクトルYの各成分が、補完後画像の画像データのベクトルXのいずれの成分に対応するかを決定するための行列である。ベクトルYが27次元であり、ベクトルXが54次元である場合、観測行列Aは、27行54列の行列となる。ベクトルYの成分とベクトルXの成分とが対応している場合、観測行列Aの成分は1となり、対応していない場合、観測行列Aの成分は0となる。

【0024】
例えば、図4に例示したように、ベクトルYのy1~y9が、それぞれベクトルXのx1~x9に対応している場合、観測行列Aの(1、1)、(2、2)、…(9、9)成分は1となる。また、ベクトルYのy10~y18が、それぞれベクトルXのx19~x27に対応している場合、観測行列Aの(10、19)、(11、20)、…(18、27)成分は1となる。また、ベクトルYのy19~y27が、それぞれベクトルXのx37~x45に対応している場合、観測行列Aの(19、37)、(29、38)、…(27、45)成分は1となる。

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

【0026】
第1実施として、関数算出部21は、式(1)を用いてコスト関数を算出することができる。argminは、コスト関数が最小となるベクトルXを求めるというような意味である。λ1~λ4は、パラメータである。λ1~λ4は、生物画像がもつ特徴によって適宜設定することができる。

【0027】
【数1】
JP2018230615A1_000003t.gif

【0028】
第2差分算出部14は、式(1)の第1項、||Y-A・X||2 ([Y-A・X]のL2ノルムと称する)を算出する。すなわち、第2差分算出部14は、取得した複数の取得画像(観測画像)の画像データ(ベクトルY)と、補完後画像の画像データのうち当該複数の取得画像に対応する画像データ(A・X)との差分を算出する。

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

【0030】
第1差分算出部13は、式(1)の第2項、第3項及び第4項を算出する。すなわち、第1差分算出部13は、補完後画像の隣接する画素の画素値の差分を算出する。隣接する画素は、例えば、x方向(第2項に対応)、y方向(第3項に対応)及びz方向(第4項に対応)に沿った画素とすることができる。

【0031】
第2項のTvxは、x方向に沿って隣接する画素の画素値の差分を算出するための差分行列であり、第3項のTvyは、y方向に沿って隣接する画素の画素値の差分を算出するための差分行列であり、第4項のTvzは、z方向に沿って隣接する画素の画素値の差分を算出するための差分行列である。

【0032】
生体組織を撮影した生物画像は、隣接する画素の画素値(例えば、輝度値)の差分のスパース性(画素値の差分がゼロとなる特性)を有するという特徴がある。隣接する画素の画素値の差分のスパース性に基づき、第1差分算出部13で算出する、x方向、y方向及びz方向に沿った隣接する画素の画素値の差分を小さくするという手法を用いることができる。

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

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

【0035】
図6はx方向に沿って隣接する画素の画素値の差分を算出する差分行列Tvxの一例を示す模式図である。z=1の補完後画像を考える。x方向の差分は、ピクセル1と2との差、ピクセル2と3との差、ピクセル4と5との差、ピクセル5と6との差、ピクセル7と8との差、ピクセル8と9との差となる。

【0036】
差分行列Tvxの(1、1)成分を1とし、(1、2)成分を-1とすることにより、ピクセル1と2との差を算出することができる。他のピクセルの差も同様である。一番右端のピクセル3、6、9は差をとる相手がいないので、差分行列Tvxの3、6、9行目の成分はない。また、他のスタック(z=1以外)についても同様にすることで、差分行列Tvxを決定することができる。

【0037】
図7はy方向に沿って隣接する画素の画素値の差分を算出する差分行列Tvyの一例を示す模式図である。z=1の補完後画像を考える。y方向の差分は、ピクセル1と4との差、ピクセル2と5との差、ピクセル3と6との差、ピクセル4と7との差、ピクセル5と8との差、ピクセル6と9との差となる。

【0038】
差分行列Tvyの(1、1)成分を1とし、(1、4)成分を-1とすることにより、ピクセル1と4との差を算出することができる。他のピクセルの差も同様である。一番下のピクセル7、8、9は差をとる相手がいないので、差分行列Tvxの7~9行目の成分はない。また、他のスタック(z=1以外)についても同様にすることで、差分行列Tvyを決定することができる。

【0039】
図8はz方向に沿って隣接する画素の画素値の差分を算出する差分行列Tvzの一例を示す模式図である。z=1~6の補完後画像を考える。z方向の差分は、z=1の画像のピクセル1とz=2の画像のピクセル10との差、z=1の画像のピクセル2とz=2の画像のピクセル11との差、…、z=1の画像のピクセル9とz=2の画像のピクセル18との差となる。

【0040】
差分行列Tvzの(1、1)成分を1とし、(1、10)成分を-1とすることにより、z=1の画像のピクセル1とz=2の画像のピクセル10との差を算出することができる。他のピクセルの差も同様である。一番下のz=6の画像のピクセル46~54は差をとる相手がいないので、差分行列Tvxの対応する成分はない。また、他のスタック同士についても同様にすることで、差分行列Tvzを決定することができる。

【0041】
一般的に、ベクトルXを式(2)で表すと、(x1、x2、…xn)は、ベクトルXの成分である。ベクトルXのL1ノルムは、式(3)で表される。また、ベクトルXのL2ノルムは、式(4)で表される。

【0042】
【数2】
JP2018230615A1_000004t.gif

【0043】
第1差分算出部13は、式(1)の第2項及び第3項の如く、補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分のL2ノルムを算出する。また、第1差分算出部13は、式(1)の第4項の如く、補完後画像のz方向に沿って隣接する画素の画素値の差分のL1ノルムを算出する。

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

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

【0046】
和算出部15は、式(1)の第5項、||X||1 ([X]のL1ノルムと称する)を算出する。すなわち、和算出部15は、補完後画像の各画素の画素値の和を算出する。生体組織を撮影した生物画像は、観測対象以外の部分のスパース性(画素値がゼロとなる特性)を有するという特徴がある。観測対象以外の部分のスパース性に基づき、和算出部15で算出した和を小さくするという手法を用いることができる。

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

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

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

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

【0051】
図9はx方向に沿って隣接する画素の画素値の差分を算出する差分行列Tvxの他の例を示す模式図である。なお、便宜上、差分行列を3行×3列の行列として説明する。z=1の補完後画像を考える。図9に示すように、x方向の差分は、ピクセル2の重み付け(図9では、重み付け係数が2)された画素値と、ピクセル1及びピクセル3それぞれの重み付け(図9では、重み付け係数が1)された画素値との差となる。同様に、x方向の差分は、ピクセル5の重み付け(図9では、重み付け係数が2)された画素値と、ピクセル4及びピクセル6それぞれの重み付け(図9では、重み付け係数が1)された画素値との差となる。

【0052】
差分行列Tvxの(2、2)成分を2とし、(2、1)成分及び(2、3)成分をそれぞれ-1とすることにより、補完後画像のピクセル2の重み付けされた画素値と、ピクセル2を間にしてx方向に沿って隣接するピクセル1及びピクセル3の重み付けされた画素値との差分を算出することができる。他のピクセルも同様である。一番右端のピクセル3、6、9、及び一番左端のピクセル1、4、7は差をとる相手がいないので、x方向に隣接する画素の画素値との差分を計算している。また、他のスタック(z=1以外)についても同様にすることで、差分行列Tvxを決定することができる。

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

【0054】
すなわち、補完後画像の注目画素の画素値をP(i+1)とし、重み付け係数をα1とする。注目画素を間にしてx方向に沿って隣接する一方の画素の画素値をP(i)とし、重み付け係数をα0とする。注目画素を間にしてx方向に沿って隣接する他方の画素の画素値をP(i+2)とし、重み付け係数をα2とする。

【0055】
第1差分算出部13は、α1×P(i+1)-{α0×P(i)+α2×P(i+2)}という式を用いて差分を算出する。ここで、α1-(α0+α2)=0とすることができ、例えば、α1=2、α0=α2=1とすると、第1差分算出部13は、2×P(i+1)-{P(i)+P(i+2)}という式により差分を算出することができる。図9の例では、α0=α2=1であり、α1=2である。

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

【0057】
図10はy方向に沿って隣接する画素の画素値の差分を算出する差分行列Tvyの他の例を示す模式図である。z=1の補完後画像を考える。図10に示すように、y方向の差分は、ピクセル4の重み付け(図10では、重み付け係数が2)された画素値と、ピクセル1及びピクセル7それぞれの重み付け(図10では、重み付け係数が1)された画素値との差となる。同様に、y方向の差分は、ピクセル5の重み付け(図10では、重み付け係数が2)された画素値と、ピクセル2及びピクセル8それぞれの重み付け(図10では、重み付け係数が1)された画素値との差となる。

【0058】
差分行列Tvyの(4、4)成分を2とし、(4、1)成分及び(4、7)成分をそれぞれ-1とすることにより、補完後画像のピクセル4の重み付けされた画素値と、ピクセル4を間にしてy方向に沿って隣接するピクセル1及びピクセル7の重み付けされた画素値との差分を算出することができる。他のピクセルも同様である。一番上端のピクセル1、2、3、及び一番下端のピクセル7、8、9は差をとる相手がいないので、y方向に隣接する画素の画素値との差分を計算している。また、他のスタック(z=1以外)についても同様にすることで、差分行列Tvyを決定することができる。

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

【0060】
すなわち、補完後画像の注目画素の画素値をP(i+1)とし、重み付け係数をα1とする。注目画素を間にしてy方向に沿って隣接する一方の画素の画素値をP(i)とし、重み付け係数をα0とする。注目画素を間にしてy方向に沿って隣接する他方の画素の画素値をP(i+2)とし、重み付け係数をα2とする。

【0061】
第1差分算出部13は、α1×P(i+1)-{α0×P(i)+α2×P(i+2)}という式を用いて差分を算出する。ここで、α1-(α0+α2)=0とすることができ、例えば、α1=2、α0=α2=1とすると、第1差分算出部13は、2×P(i+1)-{P(i)+P(i+2)}という式により差分を算出することができる。図10の例では、α0=α2=1であり、α1=2である。

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

【0063】
図11はz方向に沿って隣接する画素の画素値の差分を算出する差分行列Tvzの他の例を示す模式図である。z=1~6の補完後画像を考える。z方向の差分は、ピクセル10の重み付け(図11では、重み付け係数が2)された画素値と、ピクセル1及びピクセル19それぞれの重み付け(図11では、重み付け係数が1)された画素値との差となる。

【0064】
差分行列Tvzの(10、10)成分を2とし、(10、1)成分及び(10、19)成分をそれぞれ-1とすることにより、補完後画像のピクセル10の重み付けされた画素値と、ピクセル10を間にしてz方向に沿って隣接するピクセル1及びピクセル19の重み付けされた画素値との差分を算出することができる。他のピクセルも同様である。一番上のz=1の画像のピクセル1~9、及び一番下のz=6の画像のピクセル46~54は差をとる相手がいないので、z方向に隣接する画素の画素値との差分を計算している。また、他のスタック同士についても同様にすることで、差分行列Tvzを決定することができる。

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

【0066】
すなわち、補完後画像のz方向に沿ったスタックz′を、z′=k、k+1、k+2とする。スタックz′=k+1の補完後画像上の注目画素の画素値をP(k+1)とし、重みけ付け係数をα1とする。同様に、スタックz′=kの補完後画像上の当該注目画素に対応する位置の画素の画素値をP(k)とし、重みけ付け係数をα0とし、スタックz′=k+2の補完後画像上の当該注目画素に対応する位置の画素の画素値をP(k+2)とし、重みけ付け係数をα2とする。

【0067】
第1差分算出部13は、α1×P(k+1)-{α0×P(k)+α2×P(k+2)}という式を用いて差分を算出する。ここで、α1-(α0+α2)=0とすることができ、例えば、α1=2、α0=α2=1とすると、第1差分算出部13は、2×P(k+1)-{P(k)+P(k+2)}という式により差分を算出することができる。図11の例では、α0=α2=1であり、α1=2である。

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

【0069】
図12は各行列の大きさの一例を示す説明図である。取得した画像のスタック数をzとし、補完後画像のスタック数をz′とする。画像の大きさをn×m画素とする。ベクトルY(観測値)は、n*m*z次元、すなわちn*m*z行1列の行列に相当する。観測行列Aは、n*m*z行n*m*z′列の行列となる。ベクトルX(推定値)は、n*m*z′次元、すなわちn*m*z′行1列の行列に相当する。差分行列Tvxは、(n-1)*m*z′行n*m*z′列の行列となり、差分行列Tvyは、n*(m-1)*z′行n*m*z′列の行列となり、差分行列Tvzは、n*m*(z-1)′行n*m*z′列の行列となる。

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

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

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

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

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

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

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

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

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

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

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

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

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

【0083】
xyz方向の補完後画像に対するカーブレット変換は、xy平面画像に対するカーブレット変換、xz平面画像に対するカーブレット変換、yz平面画像に対するカーブレット変換のいずれかを用いることができる。

【0084】
図15は本実施の形態の画像処理装置100によるカーブレット変換の一例を示す模式図である。図中、符号G1~G3、C1~C3は、図2の例と同様である。すなわち、画像G1~G3は取得画像であり、画像C1~C3は補完対象画像である。上段の図は、ベクトルYの画像データを、再構成して、xy平面画像をスタックz=1からz=6まで並べたものである。xy平面画像に対するカーブレット変換は、スタックz=1からz=6までの6枚のxy平面画像それぞれに対して施される。なお、便宜上、画像の枚数を6としているが、実際は、さらに多くの枚数の画像を用いる。

【0085】
下段の図は、ベクトルYの画像データを、再構成して、xz平面画像をスタックy=1からy=6まで並べたものである。xz平面画像に対するカーブレット変換は、スタックy=1からy=6までの6枚のxy平面画像それぞれに対して施される。なお、yz平面画像はxz平面画像と同様であるので、省略している。

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

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

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

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

【0090】
図17は本実施の形態の画像処理装置100による補完後画像の第1例を示す説明図である。図中、上段の図は元画像(実際に撮影した画像)をz方向に沿って3枚並べたものである。中段の図は、比較例として一般的な画像の平滑化を使用して補完した場合を示し、3つの画像のうち、中央の画像が補完した画像である。下段の図は、本実施の形態のL1ノルム及びL2ノルムの組合せを使用して補完した場合を示し、3つの画像のうち、中央の画像が補完した画像である。右端の図は、補完した画像及び元画像の拡大像である。

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

【0092】
図18は本実施の形態の画像処理装置100による補完後画像の第2例を示す説明図である。4つ並んだ各図において、上側の図はz方向の横断面画像を示し、下側の図はxy平面画像の一つを示す。図中、左から2番目の画像は元画像(60スタック)を示し、左端の図は、z-スタック半減画像、すなわち、60スタックに半分である30スタックだけを実際に撮影した画像である。右から2番目の図は30スタックの画像から比較例(一般的な画像の平滑化を使用)による補完画像を示し、右端の図は本実施の形態による補完画像を示す。

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

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

【0095】
図19は本実施の形態の画像処理装置100による補完後画像の精度評価の一例を示す説明図である。NO.1からNO.4の4種類の画像の評価結果を示す。評価指標は、MES(Mean Square Error)及びS/N Ratio(Signal to Noise Ratio)を用いている。MESは、元画像との画素値の違いを表し、0に近いほど元画像に近いことを示す。S/N Ratioは、元画像と比較して、どの程度シグナルが増加し、ノイズが減少したかを表し、大きな値ほど良好な画像であることを示す。

【0096】
図19に示すように、NO.1からNO.4の各画像において、比較例(一般的な画像の平滑化を使用)に比べて、本実施の形態では、MSEが小さく、S/N Ratioが大きいことが分かる。

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

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

【0099】
画像データ取得部11で取得した画像データに係る取得画像の解像度を、例えば、512×512とする(低解像度とする)。推定部16によって推定された推定値に基づく補完後画像は、512×512(低解像度補完画像)となる。

【0100】
特定部20は、記憶部19に記憶した複数の低解像度画像の中から、推定部16が推定した画像データによって得られた補完対象xy平面画像を構成する複数の低解像度画像を特定する。例えば、低解像度補完画像を、所定の大きさの画素領域(例えば、30×30画素で構成される画素領域)毎に走査して、一致又は類似する辞書データ内の低解像度画像を特定する。

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

【0102】
これにより、例えば、512×512で撮影した低解像度補完画像を、1024×1024の高解像度補完画像にすることができる。

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

【0104】
図21は本実施の形態の画像処理装置100による撮影時間の一例を示す説明図である。右側の比較例は、計測機器200だけを用いて3次元の生体組織を撮影した場合のデータを示す。すなわち、1024×1024の解像度で60スタック(60枚)の光学切片画像を撮影したときのデータ量は、127.9MBであり、撮影時間は6分10秒であった。

【0105】
次に、本実施の形態により、1024×1024の解像度で30スタック(30枚)の光学切片画像を撮影したときのデータ量は、65.0MBであり、撮影時間は2分40秒であった。この場合、取得画像は30枚であり、補完対象画像は30枚、補完後画像は、比較例と同様に60枚となる。すなわち、比較例と同じ60枚の画像を得るのに、2分40秒しか要しないことが分かり、撮影時間が大幅に短縮されている。

【0106】
また、512×512の解像度で30スタック(30枚)の光学切片画像を撮影したときのデータ量は、16.3MBであり、撮影時間は1分20秒であった。この場合、取得画像は30枚の低解像度画像であり、高解像度化された補完対象画像は30枚、高解像度化された補完後画像は、比較例と同様に60枚となる。すなわち、比較例と同じ60枚の画像を得るのに、1分20秒しか要しないことが分かり、撮影時間がさらに短縮されている。また、データ量は、比較例が127.9MBであるのに対し、本実施の形態では、16.3MBであり、1撮影当たりのデータ量を削減することができる。

【0107】
また、補完後画像の高解像度化(超解像)の他の例として、以下の手法を用いることもできる。すなわち、図2に示したように、xy平面画像をz方向に対して補完して得られた補完後画像(xyz方向の立体状)を回転し、当該補完画像のyz平面画像をx方向に対して補完して補完後画像を得る。得られた補完後画像をさらに回転し、当該補完画像のxz平面画像をy方向に対して補完することによって、疑似的に3次元(3D)の解像度を上げることができる。なお、補完方向の順番は、z方向、x方向、y方向の順番に限定されない。

【0108】
次に、コスト関数の他の例について説明する。前述の例(実施例1と称する)では、コスト関数として、式(1)を用いた。他の例(実施例2と称する)では、コスト関数として式(5)を用いる。

【0109】
【数3】
JP2018230615A1_000005t.gif

【0110】
実施例2では、第1差分算出部13は、式(5)の第2項の如く、補完後画像のz方向沿って隣接する画素の画素値の差分のL1ノルムを算出する。すなわち、補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分については算出せず、コスト関数に入れていない。実施例2においても、実施例1と同様に、撮影時間を短縮することができ、時間方向の分解能を向上させることができる。なお、実施例1の方が、滑らかな画像を得る可能性が高くなる。

【0111】
また、コスト関数の他の例として、実施例3では、第1差分算出部13は、式(6)の第2項及び第3項の如く、補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分のL1ノルムを算出する。

【0112】
【数4】
JP2018230615A1_000006t.gif

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

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

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

【0116】
しかし、第1実施例の方は、補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分のL2ノルムを算出するので、より生物画像(例えば、細胞画像)らしい滑らかな画像を再構成することができる。

【0117】
図23は本実施の形態の画像処理装置100によるパラメータ決定の処理の手順の一例を示すフローチャートである。パラメータは、コスト関数で用いる、λ1~λ4である。以下では、便宜上、処理の主体を制御部10として説明する。制御部10は、パラメータ決定のための元画像の画像データを取得する(S11)。元画像とは、図2の例では、画像G1~G3、及び補完対象画像C1~C3の全てを含む。すなわち、パラメータを決定するため、まず生体組織の補完対象画像に相当する部分も含めて撮像する。

【0118】
制御部10は、パラメータの初期値(例えば、0、1、0.5など適宜の値とすることができる)を設定し(S12)、取得した画像データのうち、補完対象画像の画像データを除く(S13)。制御部10は、推定値を初期値に設定する(S14)。初期値は、0でも1でもよく、適宜の値を設定することができる。

【0119】
制御部10は、式(1)、式(5)又は式(6)のいずれかで表されるコスト関数が小さくなるように推定値(ベクトルXの各成分)を更新する(S15)。

【0120】
制御部10は、コスト関数が最小であるか否かを判定し(S16)、コスト関数が最小である場合(S16でYES)、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する(S17)。

【0121】
制御部10は、元画像の画像データ及び補完後画像の画像データに基づいて、MSE(Mean Square Error)を算出し(S18)、算出したMSEが閾値以下であるか否かを判定する(S19)。算出したMSEが閾値以下である場合(S19でYES)、制御部10は、パラメータを決定し(S20)、処理を終了する。

【0122】
算出したMSEが閾値以下でない場合(S19でNO)、制御部10は、パラメータを変更し(S21)、ステップS15以降の処理を行う。コスト関数が最小でない場合(S16でNO)、制御部10は、推定値(補完後画像)に対してカーブレット変換を行い(S22)、変換後の推定値をコスト関数に代入して(S23)、ステップS15以降の処理を行う。

【0123】
図24は本実施の形態の画像処理装置100による画像補完処理の手順の一例を示すフローチャートである。制御部10は、画像データを取得し(S31)、パラメータを取得する(S32)。パラメータは、例えば、λ1~λ4、観測行列Aなどである。なお、パラメータλ1~λ4は、図23で示す処理の結果、決定されたパラメータを用いることができる。

【0124】
制御部10は、推定値を初期値に設定する(S33)。すなわち、ベクトルXの各成分を0に設定する。なお、推定値の初期値は、0に限定されず、1など適宜設定することができる。制御部10は、式(1)、式(5)又は式(6)のいずれかで表されるコスト関数が小さくなるように推定値(ベクトルXの各成分)を更新する(S34)。

【0125】
制御部10は、コスト関数が最小であるか否かを判定し(S35)、最小でない場合(S35でNO)、推定値(補完後画像)に対してカーブレット変換を行い(S36)、変換後の推定値をコスト関数に代入して(S37)、ステップS34以降の処理を行う。コスト関数が最小である場合(S35でYES)、制御部10は、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する(S38)。

【0126】
制御部10は、高解像度化が必要であるか否かを判定し(S39)、高解像度化が必要な場合(S39でYES)、辞書データを用いて補完後画像を高解像度化し(S40)、処理を終了する。高解像度化が必要でない場合(S39でNO)、制御部10は、ステップS40の処理を行うことなく処理を終了する。なお、ステップS39、S40の処理は必須ではなく省略してもよい。

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

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

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

【0130】
具体的には、関数算出部21は、式(7)を用いてコスト関数を算出することができる。

【0131】
【数5】
JP2018230615A1_000007t.gif

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

【0133】
図25は本実施の形態の画像処理装置100によるパラメータ決定の処理の手順の他の例を示すフローチャートである。パラメータは、コスト関数で用いる、λ1~λ4である。制御部10は、パラメータ決定のための元画像の画像データを取得する(S51)。元画像とは、図2の例では、画像G1~G3、及び補完対象画像C1~C3の全てを含む。すなわち、パラメータを決定するため、まず生体組織の補完対象画像に相当する部分も含めて撮像する。

【0134】
制御部10は、パラメータの初期値(例えば、0、1、0.5など適宜の値とすることができる)を設定し(S52)、取得した画像データのうち、補完対象画像の画像データを除く(S53)。制御部10は、推定値を初期値に設定する(S54)。初期値は、0でも1でもよく、適宜の値を設定することができる。

【0135】
制御部10は、式(7)で表されるコスト関数が小さくなるように推定値(ベクトルXの各成分)を更新する(S55)。

【0136】
制御部10は、コスト関数が最小であるか否かを判定し(S56)、コスト関数が最小である場合(S56でYES)、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する(S57)。

【0137】
制御部10は、元画像の画像データ及び補完後画像の画像データに基づいて、MSE(Mean Square Error)を算出し(S58)、算出したMSEが閾値以下であるか否かを判定する(S59)。算出したMSEが閾値以下である場合(S59でYES)、制御部10は、パラメータを決定し(S60)、処理を終了する。

【0138】
算出したMSEが閾値以下でない場合(S59でNO)、制御部10は、パラメータを変更し(S61)、ステップS55以降の処理を行う。コスト関数が最小でない場合(S56でNO)、制御部10は、ステップS55以降の処理を行う。

【0139】
図26は本実施の形態の画像処理装置100による画像補完処理の手順の他の例を示すフローチャートである。制御部10は、画像データを取得し(S71)、パラメータを取得する(S72)。パラメータは、例えば、λ1~λ4、観測行列Aなどである。なお、パラメータλ1~λ4は、図25で示す処理の結果、決定されたパラメータを用いることができる。

【0140】
制御部10は、推定値を初期値に設定する(S73)。すなわち、ベクトルXの各成分を0に設定する。なお、推定値の初期値は、0に限定されず、1など適宜設定することができる。制御部10は、式(7)で表されるコスト関数が小さくなるように推定値(ベクトルXの各成分)を更新する(S74)。

【0141】
制御部10は、コスト関数が最小であるか否かを判定し(S75)、最小でない場合(S75でNO)、ステップS74以降の処理を行う。コスト関数が最小である場合(S75でYES)、制御部10は、コスト関数を最小にしたときの推定値に基づいて、補完後画像を生成する(S76)。

【0142】
制御部10は、高解像度化が必要であるか否かを判定し(S77)、高解像度化が必要な場合(S77でYES)、辞書データを用いて補完後画像を高解像度化し(S78)、処理を終了する。高解像度化が必要でない場合(S77でNO)、制御部10は、ステップS78の処理を行うことなく処理を終了する。なお、ステップS77、S78の処理は必須ではなく省略してもよい。

【0143】
本実施の形態の画像処理装置100は、CPU(プロセッサ)、RAM(メモリ)などを備えたコンピュータを用いて実現することもできる。すなわち、図23から図26に示すような、各処理の手順を定めたコンピュータプログラムをコンピュータに備えられたRAM(メモリ)にロードし、コンピュータプログラムをCPU(プロセッサ)で実行することにより、画像処理装置100を実現することができる。

【0144】
図27は本実施の形態の画像処理装置の構成の他の例を示す説明図である。図27において、符号300は、通常のコンピュータである。コンピュータ300は、制御部30、入力部40、出力部50、外部I/F(インタフェース)部60などを備える。制御部30は、CPU31、ROM32、RAM33、I/F(インタフェース)34などを備える。また、コンピュータ300は、辞書データベース70に接続することができる。

【0145】
入力部40は、計測機器200から画像データを取得する。また、入力部40は、パラメータを取得する。出力部50は、画像補完処理結果である画像データを出力する。I/F34は、制御部30と、入力部40、出力部50及び外部I/F部60それぞれとの間のインタフェース機能を有する。

【0146】
外部I/F部60は、コンピュータプログラムを記録した記録媒体M(例えば、DVDなどのメディア)からコンピュータプログラムを読み取ることが可能である。

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

【0148】
上述の実施の形態では、図2に例示したように、取得画像G1~G3それぞれに対応させて一つの画像を補完する構成であったが、これに限定されない。例えば、取得画像G1~G3それぞれに対応させて複数(例えば、二つ)の画像を補完することもできる。例えば、まず、一つの取得画像に対応させて一つの画像を補完し、次に、当該一つの取得画像と補完した画像との間の画像を補完することによって、一つの取得画像に対応させて二つの画像を補完することができる。

【0149】
(第2実施形態)
前述の第1実施形態では、xyz方向の3次元画像において、例えば、z方向に沿って配置される画像の一部を補完するものであったが、第1実施形態の構成は、2次元画像が時系列に並んだ場合において、時間方向に沿って配置される画像の一部を補完する場合にも適用することができる。以下、第2実施形態について説明する。

【0150】
図28は第2実施形態の画像処理装置100による補完後画像の一例を示す模式図である。画像処理装置100の構成は第1実施形態の場合と同様である。図中上側の図は取得画像(観測画像)を示し、下側の図は補完後画像を示す。図28に示すように、取得画像G1、G2、G3は、生体組織のxy平面上の画像であり、時系列(時点t1、t3、t5)に3枚撮影されており、スタック数=3となる。なお、図28では、便宜上、取得画像のスタック数を3としているが、スタック数は3に限定されない。

【0151】
補完後画像は、取得画像G1、G2、G3及び補完対象画像C1、C2、C3を含み、生体組織のxy平面上の画像であり、時間方向に沿って6枚あり、スタック数=6となる。補完対象画像C1、C2、C3は、実際に撮影された画像ではなく、本実施の形態の画像処理装置100によって補完される画像である。補完対象画像C1は、時点t2(時点1と時点t3との間の時点)に対応する画像であり、取得画像G1とG2との間にある。補完対象画像C2は、時点t4(時点3と時点t5との間の時点)に対応する画像であり、取得画像G2とG3との間にある。補完対象画像C3は、時点t6(時点5の後の時点)に対応する画像であり、取得画像G3の後にある。なお、画像の枚数、補完対象画像の時点などは図23の例に限定されるものではない。

【0152】
画像データ取得部11は、時系列に撮像された複数の2次元画像(例えば、xy平面画像)の画像データを取得する。一つの2次元画像を、n×m画素の画像とし、時系列(時間方向)に撮影された2次元画像の数(スタック数)をzとすると、画像データは、n×m×z行1列の行列(あるいは、n×m×z個の成分を有するベクトル)で表すことができる。取得した画像データをベクトルYで表す。

【0153】
推定部16は、画像データ取得部11で取得した画像データに基づいて、当該複数の2次元画像及び時間方向に沿って位置する補完対象2次元画像を含む補完後画像の画像データを推定する。補完対象2次元画像は、撮影された複数の2次元画像によって補完される画像であって、実際には撮影されていない画像である。補完対象2次元画像の数(スタック数)をΔzとすると、補完後画像の数(スタック数)z′は、z′=z+Δzとなる。補完後画像の画像データは、n×m×z′行1列の行列(あるいは、n×m×z′個の成分を有するベクトル)で表すことができる。補完後画像の画像データをベクトルXで表す。

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

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

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

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

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

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

【0160】
取得部は、z方向に沿って撮像された複数のxy平面画像の画像データを取得する。一つのxy平面画像を、n×m画素の画像とし、z方向に沿って撮影されたxy平面画像の数(スタック数)をzとすると、画像データは、n×m×z行1列の行列(あるいは、n×m×z個の成分を有するベクトル)で表すことができる。取得した画像データをベクトルYで表す。

【0161】
推定部は、取得した画像データに基づいて、当該複数のxy平面画像及びz方向に沿って位置する補完対象xy平面画像を含む補完後画像の画像データを推定する。補完対象xy平面画像は、撮影された複数のxy平面画像によって補完される画像であって、実際には撮影されていない画像である。補完対象xy平面画像の数(スタック数)をΔzとすると、補完後画像の数(スタック数)z′は、z′=z+Δzとなる。補完後画像の画像データは、n×m×z′行1列の行列(あるいは、n×m×z′個の成分を有するベクトル)で表すことができる。補完後画像の画像データをベクトルXで表す。

【0162】
第1差分算出部は、補完後画像の隣接する画素の画素値の差分を算出する。隣接する画素は、例えば、x方向、y方向及びz方向に沿った画素とすることができる。

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

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

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

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

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

【0168】
第2差分算出部は、取得した複数のxy平面画像の画像データと、補完後画像の画像データのうち当該複数のxy平面画像に対応する画像データとの差分を算出する。例えば、取得した画像のスタック数z=3とし、補完後画像のスタック数z′=6とする。z=1の取得画像が、z′=1の補完後画像に対応し、z=2の取得画像が、z′=3の補完後画像に対応し、z=3の取得画像が、z′=5の補完後画像に対応する場合、かかる対応付けを決定する行列A(観測行列とも称する)を考える。この場合、行列Aは、3行6列の行列となり、画像データが対応する箇所は1となり、対応しない箇所は0となる。

【0169】
推定部は、第2差分算出部で算出した差分に基づいて補完後画像の画像データを推定する。上述の例では、例えば、||Y-A・X||2 ([Y-A・X]のL2ノルムと称する)を最小二乗法によって小さくするようにして補完後画像の画像データを推定することができる。

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

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

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

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

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

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

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

【0177】
第1差分算出部は、補完後画像の注目画素の重み付けされた画素値と、当該注目画素を間にしてz方向に沿って隣接する各画素の重み付けされた画素値との差分を算出する。例えば、補完後画像のz方向に沿ったスタックz′を、z′=k、k+1、k+2とする。スタックz′=k+1の補完後画像上の注目画素の画素値をP(k+1)とし、重みけ付け係数をα1とする。同様に、スタックz′=kの補完後画像上の当該注目画素に対応する位置の画素の画素値をP(k)とし、重みけ付け係数をα0とし、スタックz′=k+2の補完後画像上の当該注目画素に対応する位置の画素の画素値をP(k+2)とし、重みけ付け係数をα2とする。

【0178】
第1差分算出部は、α1×P(k+1)-{α0×P(k)+α2×P(k+2)}という式を用いて差分を算出する。ここで、α1-(α0+α2)=0とすることができ、例えば、α1=2、α0=α2=1とすると、第1差分算出部は、2×P(k+1)-{P(k)+P(k+2)}という式により差分を算出することができる。これにより、補完後画像の任意の画素の画素値と、z方向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、z方向に沿った偏りを抑制して、3次元構造の生体組織に適した生物画像を補完することができる。

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

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

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

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

【0183】
例えば、補完後画像の注目画素の画素値をP(i+1)とし、重み付け係数をα1とする。注目画素を間にしてx方向に沿って隣接する一方の画素の画素値をP(i)とし、重み付け係数をα0とする。注目画素を間にしてx方向に沿って隣接する他方の画素の画素値をP(i+2)とし、重み付け係数をα2とする。

【0184】
第1差分算出部は、α1×P(i+1)-{α0×P(i)+α2×P(i+2)}という式を用いて差分を算出する。ここで、α1-(α0+α2)=0とすることができ、例えば、α1=2、α0=α2=1とすると、第1差分算出部は、2×P(i+1)-{P(i)+P(i+2)}という式により差分を算出することができる。これにより、補完後画像の任意の画素の画素値と、x向に沿って当該任意の画素の両側に隣接する画素の画素との差分を用いるので、x方向に沿った偏りを抑制して、3次元構造の生体組織に適した生物画像を補完することができる。なお、y方向についても、x方向と同様である。

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

【0186】
第1差分算出部は、補完後画像のz方向に沿って隣接する画素の画素値の差分のL1ノルムを算出する。例えば、ベクトルXの成分を(x1、x2、…xn)とすると、ベクトルXのL1ノルムは、{|x1|+|x2|+…+|xn|}という式で表すことができる。

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

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

【0189】
第1差分算出部は、補完後画像のx方向及びy方向に沿って隣接する画素の画素値の差分のL2ノルムを算出する。例えば、ベクトルXの成分を(x1、x2、…xn)とすると、ベクトルXのL2ノルムは、√{(x1)2 +(x2)2 +…(xn)2 }という式で表すことができる。

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

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

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

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

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

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

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

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

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

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

【0200】
取得部は、時系列に撮像された複数の2次元画像(例えば、xy平面画像)の画像データを取得する。一つの2次元画像を、n×m画素の画像とし、時系列(時間方向)に撮影された2次元画像の数(スタック数)をzとすると、画像データは、n×m×z行1列の行列(あるいは、n×m×z個の成分を有するベクトル)で表すことができる。取得した画像データをベクトルYで表す。

【0201】
推定部は、取得した画像データに基づいて、当該複数の2次元画像及び時間方向に沿って位置する補完対象2次元画像を含む補完後画像の画像データを推定する。補完対象2次元画像は、撮影された複数の2次元画像によって補完される画像であって、実際には撮影されていない画像である。補完対象2次元画像の数(スタック数)をΔzとすると、補完後画像の数(スタック数)z′は、z′=z+Δzとなる。補完後画像の画像データは、n×m×z′行1列の行列(あるいは、n×m×z′個の成分を有するベクトル)で表すことができる。補完後画像の画像データをベクトルXで表す。

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

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

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

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

【0206】
推定値が得られることによって、スタック数zの画像(取得した画像)から、スタック数z′の補完後画像を得ることができる。これにより、1枚の2次元画像を撮影する間に、生体組織の動的な変化がある場合でも、当該撮影する間の生体組織の変化を補完後画像(より具体的には、補完対象2次元画像)によって捉えることができるので、生命現象を一層リアルタイムで捉えることが可能となる。
【符号の説明】
【0207】
10 制御部
11 画像データ取得部
12 パラメータ取得部
13 第1差分算出部
14 第2差分算出部
15 和算出部
16 推定部
17 補完画像生成部
18 変換部
19 記憶部
20 特定部
21 関数算出部
30 制御部
31 CPU
32 ROM
33 RAM
34 I/F
40 入力部
50 出力部
60 外部I/F部
70 辞書データベース
100 画像処理装置
200 計測機器
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16
【図18】
17
【図19】
18
【図20】
19
【図21】
20
【図22】
21
【図23】
22
【図24】
23
【図25】
24
【図26】
25
【図27】
26
【図28】
27